最近在上的一门课,用到了Newmark-beta逐步积分法,其他文章都描述的很模糊,于是自己搞懂了之后特此记录下来,仅是个人的一些理解,可能有错误,欢迎交流。
有这么一道题
其中的参数: 梁长 L m 25
杨氏弹性模量 E N/m^2 2.89e9
泊松比 σ 0.2
转动惯量 I m^4 2.9
每延米质量 m kg/m 2303
悬挂刚度 k_v kN/m 1595
悬吊质量 M_v kg 5750
行车速度 V km/h 100
要求画出梁的跨中位移和加速度曲线
课上讲了一个例题,只不过是移动力是恒定的,而此题的力随时间是变化的。
用到的解法是振型叠加法(/分解法)
经过一番推导会得到一个方程(过程不必赘述,想要资料请留言邮箱)
移动力作用下简支梁第n阶广义动力平衡方程
m L 2 q ¨ n ( t ) + m L 2 ω n 2 q n ( t ) = p ( t ) s i n n π x t L (a) \frac{mL}{2}\ddot q_n (t)+\frac{mL}{2} \omega_n^2 q_n (t)=p(t)sin\frac{n\pi x_t}{L} \tag{a} 2mLq¨n(t)+2mLωn2qn(t)=p(t)sinLnπxt(a) 其 中 , m 为 每 延 米 质 量 , k g / m ; L 为 梁 长 , m ; ω n = n 2 π 2 E I m L 4 ; x t = v t ; p ( t ) 为 移 动 集 中 力 。 其中,\\[2ex]m为每延米质量,kg/m;\\[2ex] L为梁长,m;\\[2ex] \omega_n=n^2 \pi ^2 \sqrt{\frac{EI}{mL^4}};\\[2ex] x_t =vt;\\[2ex] p(t)为移动集中力。 其中,m为每延米质量,kg/m;L为梁长,m;ωn=n2π2mL4EI;xt=vt;p(t)为移动集中力。
移动集中力p(t)可以由下式计算出
p
(
t
)
=
M
v
g
+
k
v
[
z
(
t
)
−
∑
i
=
1
∞
q
i
(
t
)
s
i
n
i
π
v
t
L
]
(b)
p(t)=M_v g+k_v [z(t)-\sum_{i=1}^{\infty}q_i(t)sin\frac{i\pi vt}{L}] \tag{b}
p(t)=Mvg+kv[z(t)−i=1∑∞qi(t)sinLiπvt](b)
其
中
,
M
v
为
悬
吊
质
量
,
k
g
;
k
v
为
悬
挂
刚
度
,
k
N
/
m
;
z
(
t
)
为
质
量
块
位
移
,
m
。
其中,\\[2ex]M_v为悬吊质量,kg;\\[2ex] k_v为悬挂刚度,kN/m;\\[2ex] z(t)为质量块位移,m。
其中,Mv为悬吊质量,kg;kv为悬挂刚度,kN/m;z(t)为质量块位移,m。
对质量块进行受力分析会有
M
v
z
¨
(
t
)
+
k
v
[
z
(
t
)
−
∑
i
=
1
∞
q
i
(
t
)
s
i
n
i
π
v
t
L
]
=
0
(c)
M_v \ddot z(t)+k_v[z(t)-\sum_{i=1}^\infty q_i(t)sin\frac{i\pi vt}{L}]=0 \tag{c}
Mvz¨(t)+kv[z(t)−i=1∑∞qi(t)sinLiπvt]=0(c)
将上述(a)、(b)、(c)三个式子结合,可以推出如下方程
{
M
v
z
¨
(
t
)
+
k
v
z
(
t
)
−
k
v
∑
i
=
1
∞
q
i
(
t
)
s
i
n
i
π
v
t
L
=
0
q
¨
n
(
t
)
+
[
ω
n
2
q
n
(
t
)
+
2
m
L
∑
i
=
1
∞
k
v
q
i
(
t
)
s
i
n
i
π
v
t
L
s
i
n
n
π
v
t
L
]
−
2
m
L
k
v
z
(
t
)
s
i
n
n
π
v
t
L
=
2
m
L
M
v
g
s
i
n
n
π
v
t
L
\begin{cases} M_v \ddot z(t)+k_vz(t)-k_v \sum_{i=1}^\infty q_i(t)sin\frac{i\pi vt}{L}=0 \\[4ex] \ddot q_n(t)+[\omega_n^2 q_n(t)+\frac{2}{mL}\sum_{i=1}^\infty k_v q_i(t)sin\frac{i\pi vt}{L}sin\frac{n\pi vt}{L}]-\frac{2}{mL}k_v z(t)sin\frac{n\pi vt}{L}\\[2ex]=\frac{2}{mL}M_vgsin\frac{n\pi vt}{L} \end{cases}
⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧Mvz¨(t)+kvz(t)−kv∑i=1∞qi(t)sinLiπvt=0q¨n(t)+[ωn2qn(t)+mL2∑i=1∞kvqi(t)sinLiπvtsinLnπvt]−mL2kvz(t)sinLnπvt=mL2MvgsinLnπvt
整理成矩阵形式就可以求解了。
OK,式子都有了,关键问题就是如何求解。这是一个微分方程问题,老师讲授此部分时推荐的newmark-beta逐步积分法。当时对此方法听了个大概,感觉跟数值分析课程上讲的差不多,没怎么注意。但是当自己真正要用其进行matlab编程的时候犯了难,网上此法求解过程描述的也有,但是不够清楚。。也可能是本人太笨,最后找到一篇学位论文,大致明白了原理,编程也就不是问题了。
论文链接放出来,感兴趣的朋友可以自己看看,在3.1.1节
基于Duhamel积分的地震作用下车桥耦合振动分析
在Newmark-beta中用到的字母个人认为应该是 γ 和 β \gamma和\beta γ和β,学位论文中可能是为了降重用了 α 和 β \alpha和\beta α和β
个人建议在使用字母时一定不要擅自乱用,一般的字母或形式都有其含义,乱用很容易造成麻烦,也显得不是很专业。
下面简单说明一下该方法的推导过程
首先会设系统振动的微分方程为
M
z
¨
+
C
z
˙
+
K
z
=
F
(1)
M\ddot z+C\dot z+Kz=F \tag{1}
Mz¨+Cz˙+Kz=F(1)
newmark-beta法会假设:
z
˙
t
+
Δ
t
=
z
˙
t
+
[
(
1
−
β
)
z
¨
t
+
β
z
˙
t
+
Δ
t
]
Δ
t
(2)
\dot z_{t+\Delta t}=\dot z_t +[(1-\beta)\ddot z_t +\beta \dot z_{t+\Delta t}]\Delta t \tag{2}
z˙t+Δt=z˙t+[(1−β)z¨t+βz˙t+Δt]Δt(2)
z
t
+
Δ
t
=
z
t
+
z
˙
t
Δ
t
+
[
(
1
/
2
−
γ
)
+
γ
z
¨
t
+
Δ
t
]
Δ
t
2
(3)
z_{t+\Delta t}=z_t + \dot z_t \Delta t+[(1/2-\gamma)+\gamma \ddot z_{t+\Delta t}]\Delta t^2 \tag{3}
zt+Δt=zt+z˙tΔt+[(1/2−γ)+γz¨t+Δt]Δt2(3)
其中
γ
和
β
\gamma和\beta
γ和β的取值对该算法有较大的影响,当满足以下条件时,算法无条件稳定
β
>
=
1
/
2
\beta >= 1/2
β>=1/2
γ
>
=
(
1
/
2
+
β
)
2
/
4
\gamma>=(1/2+\beta)^2 /4
γ>=(1/2+β)2/4一般取值
γ
=
0.25
,
β
=
0.5
\gamma=0.25,\beta=0.5
γ=0.25,β=0.5也称作平均加速度法。
根据式(2)和(3),可以推导出用于迭代的计算式
z
¨
t
+
Δ
t
=
1
γ
Δ
t
2
(
z
t
+
Δ
t
+
z
t
)
−
1
γ
Δ
t
z
˙
t
−
(
1
2
γ
−
1
)
z
¨
t
(4)
\ddot z_{t+\Delta t}=\frac{1}{\gamma \Delta t^2}(z_{t+\Delta t}+z_t)-\frac{1}{\gamma \Delta t}\dot z_t-(\frac{1}{2\gamma}-1)\ddot z_t \tag{4}
z¨t+Δt=γΔt21(zt+Δt+zt)−γΔt1z˙t−(2γ1−1)z¨t(4)
z
˙
t
+
Δ
t
=
β
γ
Δ
t
(
z
t
+
Δ
t
−
z
t
)
+
(
1
−
β
γ
z
˙
t
)
+
(
1
−
β
2
γ
Δ
t
z
¨
t
)
(5)
\dot z_{t+\Delta t}=\frac{\beta}{\gamma \Delta t}(z_{t+\Delta t}-z_t)+(1-\frac{\beta}{\gamma}\dot z_t)+(1-\frac{\beta}{2\gamma}\Delta t \ddot z_t) \tag{5}
z˙t+Δt=γΔtβ(zt+Δt−zt)+(1−γβz˙t)+(1−2γβΔtz¨t)(5)然后
t
+
Δ
t
t+\Delta t
t+Δt时刻系统的振动微分方程如下
M
z
¨
t
+
Δ
t
+
C
z
˙
t
+
Δ
t
+
K
z
t
+
Δ
t
=
F
t
+
Δ
t
(6)
M\ddot z_{t+\Delta t}+C\dot z_{t+\Delta t}+Kz_{t+\Delta t}=F_{t+\Delta t} \tag{6}
Mz¨t+Δt+Cz˙t+Δt+Kzt+Δt=Ft+Δt(6)把式(4)、(5)代入(6)中就会得到
K
^
z
t
+
Δ
t
=
F
^
t
+
Δ
t
(7)
\hat K z_{t+\Delta t}=\hat F_{t+\Delta t}\tag{7}
K^zt+Δt=F^t+Δt(7)
其中,
K
^
=
K
+
1
γ
Δ
t
2
M
+
β
γ
Δ
t
C
(8)
\hat K=K+\frac{1}{\gamma \Delta t^2}M+\frac{\beta}{\gamma \Delta t}C\tag{8}
K^=K+γΔt21M+γΔtβC(8)
F
^
t
+
Δ
t
=
F
t
+
Δ
t
+
M
[
1
γ
Δ
t
2
z
t
+
1
γ
Δ
t
z
˙
t
+
(
1
2
γ
−
1
)
z
¨
t
]
+
C
[
β
γ
Δ
t
z
t
+
(
β
γ
−
1
)
z
˙
t
+
(
β
γ
−
1
)
Δ
t
z
¨
t
]
(9)
\hat F_{t+\Delta t}=F_{t+\Delta t}+M[\frac{1}{\gamma \Delta t^2}z_t+\frac{1}{\gamma \Delta t}\dot z_t+(\frac{1}{2\gamma}-1)\ddot z_t]\\[2ex]+C[\frac{\beta}{\gamma \Delta t}z_t+(\frac{\beta}{\gamma}-1)\dot z_t+(\frac{\beta}{\gamma-1})\Delta t\ddot z_t] \tag{9}
F^t+Δt=Ft+Δt+M[γΔt21zt+γΔt1z˙t+(2γ1−1)z¨t]+C[γΔtβzt+(γβ−1)z˙t+(γ−1β)Δtz¨t](9)
由式(9)直接可以算出
z
t
+
Δ
t
z_{t+\Delta t}
zt+Δt,即下一时刻的位移,然后通过式(4)和式(2)就可以算出加速度和速度,注意一定是先算(4),因为(2)中要用到(4)的结果,为什么不用(5)算速度呢,个人想法是提高精度,类似赛德尔迭代一样。在编程的时候会把各矩阵或向量的系数用
α
0
至
α
7
\alpha_0至\alpha_7
α0至α7代替,好像一般也有固定的顺序,建议还是遵守规则,按照参考中的来。
下面就附上matlab程序,
设定了两个函数,在另外两个文件中
function w = w_(n)
L=25;
E=2.89e9;
I=2.9;
m=2303;
w=(n^2)*(pi^2)*((E*I)/(m*L^4))^(1/2);
end
function phy = phy_(n,t)
v=100/3.6;
L=25;
phy=sin(n*pi*v*t/L);
end
开始计算
clear;
%初始参数
L = 25;
E = 2.89e9;
sigma = 0.2;
I = 2.9;
m = 2303;
k_v = 1595;
M_v = 5750;
v = 100/3.6;
%生成参数
rho_k = 2 * k_v / (m*L);
rho_F = 2 * M_v * 9.81 / (m*L);
M = diag([1,1,1,1,1,M_v]);%质量矩阵
T = L / v;%通过时长
nt = 1000;%步数
dt = T / nt;%步长
gama = 0.25;beta = 0.5;
a0 = 1 / (gama * dt^2);
a1 = beta / (gama * dt);
a2 = 1 / (gama * dt);
a3 = 1 / (2 * gama) - 1;
a4 = beta / gama -1 ;
a5 = (beta / (2 * gama) - 1) * dt;
a6 = dt * (1 - beta);
a7 = dt * beta;
d = zeros(6,nt);%初始化广义位移,6自由度,q1-q5,z_t
v = zeros(6,nt);%速度
a = zeros(6,nt);%加速度
%newmark-beta法
for i = 2:nt %第1步已赋初值,从第2步开始到1000步
t = (i - 1) * dt; %响应时间
K_0 = [w_(1)^2 + rho_k * phy_(1,t) * phy_(1,t),rho_k * phy_(1,t) * phy_(2,t),...
rho_k * phy_(1,t) * phy_(3,t),rho_k * phy_(1,t) * phy_(4,t),...
rho_k * phy_(1,t) * phy_(5,t),-rho_k * phy_(1,t);%第一行
rho_k * phy_(2,t) * phy_(1,t),w_(2)^2 + rho_k * phy_(2,t) * phy_(2,t),...
rho_k * phy_(2,t) * phy_(3,t),rho_k * phy_(2,t) * phy_(4,t),...
rho_k * phy_(2,t) * phy_(5,t),-rho_k * phy_(2,t);%第二行
rho_k * phy_(3,t) * phy_(1,t),rho_k * phy_(3,t) * phy_(2,t),...
w_(3)^2 + rho_k * phy_(3,t) * phy_(3,t),rho_k * phy_(3,t) * phy_(4,t),...
rho_k * phy_(3,t) * phy_(5,t),-rho_k * phy_(3,t);%第三行
rho_k * phy_(4,t) * phy_(1,t),rho_k * phy_(4,t) * phy_(2,t),...
rho_k * phy_(4,t) * phy_(3,t),w_(4)^2 + rho_k * phy_(4,t) * phy_(4,t),...
rho_k * phy_(4,t) * phy_(5,t),-rho_k * phy_(4,t);%第四行
rho_k * phy_(5,t) * phy_(1,t),rho_k * phy_(5,t) * phy_(2,t),...
rho_k * phy_(5,t) * phy_(3,t),rho_k * phy_(5,t) * phy_(4,t),...
w_(5)^2 + rho_k * phy_(5,t) * phy_(5,t),-rho_k * phy_(5,t);%第五行
-k_v * phy_(1,t),-k_v * phy_(2,t),-k_v * phy_(3,t),-k_v * phy_(4,t),-k_v * phy_(5,t),k_v;%第六行
];%初始广义刚度矩阵6*6
K = K_0 + a0 * M;%用于计算的刚度矩阵
F_0 = [rho_F * phy_(1,t);rho_F * phy_(2,t);rho_F * phy_(3,t)
rho_F * phy_(4,t);rho_F * phy_(5,t);0];%广义载荷矩阵
F = F_0 + M * (a0 * d(:,i-1) + a2 * v(:,i-1) + a3 * a(:,i-1));%用于计算的载荷矩阵
d(:,i) = K \ F;%计算位移
a(:,i) = a0 * (d(:,i) - d(:,i-1)) - a2 * v(:,i-1) - a3 * a(:,i-1);%加速度
v(:,i) = v(:,i-1) + a6 * a(:,i-1) + a7 * a(:,i);%速度
end
t0=linspace(0,T,nt);
figure
plot(t0,d(1,:) * sin(1 * pi / 2) + d(2,:) * sin(2 * pi / 2) + d(3,:) * sin(3 * pi / 2) + d(4,:) * sin(4 * pi / 2) + d(5,:) * sin(5 * pi / 2))
ylabel('位移/m');
grid on;
xlabel('时间/s');
title('跨中位移响应时程');
figure
plot(t0,a(1,:) * sin(1 * pi / 2) + a(2,:) * sin(2 * pi / 2) + a(3,:) * sin(3 * pi / 2) + a(4,:) * sin(4 * pi / 2) + a(5,:) * sin(5 * pi / 2));
ylabel('加速度/m/s^2');
grid on;
xlabel('时间/s');
title('跨中加速度响应时程');
上述程序中
d就是位移矩阵,前面几行是梁的N阶振型坐标,最后一行是质量块的坐标
v是速度矩阵,a是加速度矩阵。
最后放一下结果图
关于阶数,个人认为五阶已经是足够了很保守,三阶就ok,二阶应该也可以。
看一下各阶数的比较,一阶还是主要部分,二阶稍微能察觉到一点,后面几阶数量级太小,几乎为零。