摩擦摆隔震系统

摩擦摆隔震系统1. 模型2. 运动方程3. 状态转换4. 运动方程求解5. LS-DYNA 求解验证5.1位移对比5.2 加速度对比6. 简单的多层结构FPS计算程序

最近接触到一种隔震方式叫摩擦摆隔震系统（Friction Pendulum System）。

1. 模型

其最简单的模型如下图所示：

fps model

$\mu$ $m_1$ $m_2$ $m_7$ $c_i、k_i$ 为各层刚度与阻尼系数。

2. 运动方程

首先，滑块水平、竖直方向坐标与转角之间的关系：

x_1=rsin(\theta), \dot{x}_1=rcos(\theta)\dot{\theta}, \ddot{x}_1=-rsin(\theta)\dot{\theta}^{2}+rcos(\theta)\ddot{\theta}\\ y_1=r-rcos(\theta), \dot{y}_1=rsin(\theta)\dot{\theta}, \ddot{y}_1=rcos(\theta)\dot{\theta}^2+rsin(\theta)\ddot{\theta}

model

$F$ $W$ ，其值分别为：

m_1\ddot{x_1}+c_2(\dot{x_1}-\dot{x_2})+k_2(x_1-x_2)=F\\ W=\sum_{i=1}^{n} m_ig+\sum_{i=1}^{n} m_i\ddot{y_1}

$T$ $P$ ，其值分别为：

T=Fcos(\theta)+Wsin(\theta)\\ P=Wcos(\theta)-Fsin(\theta)

或者

F=Tcos(\theta)-Psin(\theta)\\ W=Tsin(\theta)+Pcos(\theta)

滑动状态下，根据库伦摩擦公式：

T=-P\mu \cdot sgn(\dot\theta)

则有：

W=P\mu sin(\theta)+Pcos(\theta)=P(\mu sin(\theta)+cos(\theta))\\ F=P\mu cos(\theta)-Psin(\theta)=P(\mu cos(\theta)-sin(\theta))\\ =\frac {W}{\mu sin(\theta)+cos(\theta)}(\mu cos(\theta)-sin(\theta))

把方程1和方程6带入方程2，即可得到滑块的运动微分方程：

\ddot \theta+\frac{A_2}{A_1}\dot{\theta}^2+\frac{A_3}{A_1}-\frac{1}{m_1rA_1}[k_2(rsin\theta-x_2)+c_2(r\dot{\theta}cos\theta-\dot{x_2})](\mu sgn(\dot{\theta})sin\theta-cos\theta)\\ =\frac{\ddot{a_g}}{rA_1}(\mu sgn(\dot \theta)sin\theta-cos\theta)

$sgn(\dot\theta)$ 为符号函数，用于描述摩擦力与运动方向的关系。而

A_1=1+\frac{\sum_{i=2}^n m_i}{m_1}({\mu}sgn(\dot\theta)cos(\theta)+sin\theta)sin\theta

A_2=\mu sgn(\dot\theta)+\frac{\sum_{i=2}^n m_i}{m_1}(\mu sgn(\dot\theta)cos\theta+sin\theta)cos\theta

A_3=\frac{g}{r}\frac{\sum_{i=1}^n m_i}{m_1}(\mu sgn(\dot\theta)cos\theta+sin\theta)

中间层的运动微分方程为：

m_i\ddot{x_i}+k_i(x_i-x_{i-1})+k_{i+1}(x_i-x_{i+1})+c_i(\dot {x_i}-\dot x_{i-1})+c_{i+1}(\dot x_i-\dot x_{i+1})=-m_i\ddot{a_g}

$m_n$ 的运动微分方程为：

m_n\ddot x_n+k_n(x_n-x_{n-1})+c_n(\dot y_n-\dot y_{n-1})=-m_n\ddot a_g

由方程7、11、12组成的方程组即为该摩擦摆隔震系统的微分方程组。

3. 状态转换

$T<P\mu$ 时，且滑块速度为0时，滑块则不会滑动，相当于固定基础，此时方程7将不存在。

4. 运动方程求解

求解该方程组的Matlab代码如下:


113
1
function [ rt,rd,rv,ra,rrd] = solve_fps(floordata,r,u,dt,ag)
2
% *SOLVE_FPS 摩擦摆隔振系统求解函数*
3
4
% floordata：楼层参数,其格式如下，
5
%     {[m1],[k1],[c1]
6
%      [m1],[k1],[c1]
7
%      [m2],[k2],[c2]
8
%      ...
9
%      [mn],[kn],[cn]}
10
%
11
%     r：摩擦摆半径；
12
%     u：摩擦系数；
13
%     dt：地震时程采样时间间隔；
14
%     ag：地震时程数据，单位为g；
15
% 返回值：
16
%   rt：时刻数组；
17
%   rd：响应位移时程；
18
%   rv：响应速度；
19
%   ra：响应加速度；
20
%   rrd：层间位移；
21
22
nfloor=size(floordata,1);   %楼层数
23
ag=ag*9.8;                  %乘以g
24
tg=0:dt:dt*length(ag)-dt;   %时间向量
25
m=floordata(:,1);           %m=各层质量
26
m=cell2mat(m);
27
k=floordata(2:end,2);       %k=各层刚度，
28
k=cell2mat(k);
29
k=[1;k];
30
k=k*1000;                   %原刚度单位为kN，固乘以1000；
31
c=floordata(2:end,3);       %c=各层阻尼系数
32
c=cell2mat(c);
33
c=[1;c];
34
35
m1=m(1);                    %滑块质量
36
M=sum(m(2:end));            %上层总质量
37
38
init=zeros(nfloor*2,1);     %初始值0
39
[rt,ry]=ode45(@wffc_fps,tg,init);   %龙格库塔法求解
40
41
rd=ry(:,1:2:2*nfloor);      %各层位移
42
rv=ry(:,2:2:2*nfloor);      %各层速度
43
44
%根据计算得到的位移、速度，带入方程，计算相对加速度响应
45
ra=zeros(length(rt),nfloor);
46
for tn=1:length(rt)
47
    dy=wffc_fps(rt(tn),ry(tn,:));
48
    ra(tn,:)=dy(2:2:2*nfloor);
49
end
50
51
ra(:,1)=r*cos(rd(:,1)).*ra(:,1)-r*sin(rd(:,1)).*rv(:,1).^2;
52
rv(:,1)=r*cos(rd(:,1)).*rv(:,1);
53
rd(:,1)=r*sin(rd(:,1));
54
55
56
% %相对加速度加上地面加速度，得到绝对加速度
57
for i=1:nfloor
58
    ra(:,i)=ra(:,i)+ag;
59
end
60
61
rrd=rd(:,2:end)-rd(:,1:end-1);
62
rrd=[rd(:,1),rrd];
63
64
    function dy=wffc_fps(t,y)
65
        % n 自由度系统微分方程组
66
        dy=zeros(nfloor*2,1);
67
        g=9.8;
68
        uy=interp1(tg,ag,t,'linear'); %根据线性插值，获取对应时间的地震加速度值
69
70
        % 中间参数
71
        sgn=sign(y(2));
72
        siny=sin(y(1));
73
        cosy=cos(y(1));
74
        A1=1+(M/m1)*(u*sgn*cosy+siny)*siny;
75
        A2=u*sgn+(M/m1)*(u*sgn*cosy+siny)*cosy;
76
        A3=(g/r)*((M+m1)/m1)*(u*sgn*cosy+siny);
77
78
        %微分方程组
79
        dy(1)=y(2);
80
        dy(2)=uy*(u*sgn*siny-cosy)/r/A1+(k(2)*(r*siny-y(3))+c(2)*(r*y(2)*cosy-y(4)))*(u*sgn*siny-cosy)/m1/r/A1-A3/A1-A2/A1*y(2)^2;
81
82
        %计算切向力与法向力，用于判断滑块是否滑动
83
        G=(m1+M)*g+r*M*(cosy*y(2)^2+siny*dy(2));
84
        P=-G*siny-(k(2)*(r*siny-y(3))+c(2)*(r*y(2)*cosy-y(4)))*cosy-m1*uy*cosy;
85
        N=G*cosy-(m1*uy+k(2)*(r*siny-y(3))+c(2)*(r*y(2)*cosy-y(4)))*siny+m1*r*y(2)^2;
86
        F=u*N;
87
        %如果滑块没有滑动，则直接设置其速度与加速度为0
88
        if abs(y(2))<1e-5&&abs(P)<=abs(F)
89
            dy(1)=0;
90
            dy(2)=0;
91
        end
92
93
        %如果仅有两层，则第2层就是顶层。
94
        if nfloor==2
95
            dy(3)=y(4);
96
            dy(4)=(-m(2)*uy-c(2)*(y(4)-r*y(2)*cosy)-k(2)*(y(3)-r*siny))/m(2);
97
        end
98
        
99
        if nfloor>2
100
            dy(3)=y(4);
101
            dy(4)=(-m(2)*uy-c(3)*(y(4)-y(6))-c(2)*(y(4)-r*y(2)*cosy)-k(3)*(y(3)-y(5))-k(2)*(y(3)-r*siny))/m(2);
102
103
            for fn=3:nfloor
104
                dy(fn*2-1)=y(fn*2);
105
                if fn==nfloor
106
                    dy(fn*2)=(-m(fn)*uy-c(fn)*(y(2*fn)-y(2*fn-2))-k(fn)*(y(2*fn-1)-y(2*fn-3)))/m(fn);
107
                else
108
                    dy(fn*2)=(-m(fn)*uy-c(fn+1)*(y(2*fn)-y(2*fn+2))-c(fn)*(y(2*fn)-y(2*fn-2))-k(fn+1)*(y(2*fn-1)-y(2*fn+1))-k(fn)*(y(2*fn-1)-y(2*fn-3)))/m(fn);
109
                end
110
            end
111
        end
112
    end
113
end

5. LS-DYNA 求解验证

为了验证上述求解的正确性，用DYNA做了一个小例子，进行时程分析。

DYNA模型

如图所示，基座半径为1m，滑块质量186kg，摩擦系数为0.1。基础之上共有5层，每层质量都与基座相同。每层都用质量点表示，各层之间通过离散梁单元连接，梁单元使用材料模型为LINEAR_ELASTIC_DISCRETE_BEAM，设置其沿水平方向的刚度即横向刚度为186000N/m，阻尼系数为450，沿重力方向刚度取大三个数量级，阻尼系数也大两个数量级。

通过下面的对比可以看出，DYNA的计算结果几乎与MATLAB算出来的一致。

有几点可能的原因导致结算结果的轻微偏差：

MATLAB计算时，需要对进行滑块的状态转换，判断滑块不滑动的充要条件是，切向力小于滑动摩擦力且速度为0。为提高迭代速度，代码中取滑块速度小于1e-5时则取其速度为0，这里会导致轻微偏差；
在DYNA中模拟摩擦摆最关键的是滑块与基座之间的接触设置。由于有限元软件中接触不可避免的产生接触面间的相互渗透，当滑块与基座的节点在计算过程中渗透量越大时，其摩擦系数并不是理想的设置的摩擦系数，导致较大偏差。为使仿真过程中的等效摩擦系数是期望的值，需要尽可能提高接触界面间的刚度，使其渗透减小。

5.1位移对比

第一层位移对比

第二层位移对比

第6层位移对比

5.2 加速度对比

第一层加速度对比

第六层加速度对比

6. 简单的多层结构FPS计算程序

根据上述内容，可编写如下程序，加入了没有隔振系统的计算结果，输出任意层的位移、加速度、层间位移等，可以对比隔振效果。