固体力学中的数值方法_Newmark算法

Newmark 算法程序计算报告

1.Newmark 算法理论基础

在~t t t +∆的时间区域内，Newmark 积分方法采用下列的假设，即

()2112t t t t t t t t

t t t t t t t t

δδαα+∆+∆+∆+∆=+-+∆⎡⎤⎣⎦⎡⎤⎛⎫

=+∆+-+∆ ⎪⎢⎥⎝⎭⎣⎦

a a a a a a a a a （1.1） 其中α和δ是按积分精度和稳定性要求决定的参数。另一方面，α和δ取不同的值则代表了不同的数值积分方案。当1/6,1/2αδ==时，（1.1）式相应于线性加速度法，因为这时他们可以由下式得到

()/ (0)t t t t t t t t ττ+∆+∆=+-∆≤≤∆a a a a （1.2）

当1/4,1/2αδ==时，Newmark 算法相应于常平均加速度法这样一种无条件稳定的积分方案。此时，t ∆内的加速度为

()1

2

t t t t t +∆+∆=

+a a a （1.3） 因此，将（1.1）式可以得到

()2

11112t t t t t t t t t ααα+∆+∆⎛⎫

=

---- ⎪∆∆⎝⎭

a a a a a （1.4） 代入到动力学平衡方程中可以得到

2

2111112 112t t t t t t t t t t t t t t t t δαααααδδδααα+∆+∆⎡⎤

⎛⎫⎛⎫++=+++-+ ⎪ ⎪⎢⎥∆∆∆∆⎝⎭⎝⎭⎣⎦⎡⎤⎛⎫⎛⎫+-+-∆ ⎪ ⎪⎢⎥

∆⎝⎭⎝⎭⎣⎦

K M C a Q M a a a C a a a （1.5）

2.Newmark 算法计算步骤

（1）形成刚度矩阵 K 、质量矩阵M 、阻尼矩阵 C 。 （2）给定0 a ，0a ，和0a

（3）选择时间步长t ∆及参数α和δ，并计算积分常数。 这里要求：()2

0.5,0.250.5δαδ≥≥+

()01232

4567111,,,121,2,1

,2c c c c t t t t c c c t c t δαααα

δδδδαα=

===-∆∆∆∆⎛⎫=-=-=∆-=∆ ⎪⎝⎭

（4）形成有效刚度矩阵01ˆˆc c + K :K =K +M C （5）三角分解ˆˆT

K

:K =LDL 对于每一个时间步长

（1）计算时间t t +∆的有效载荷

()()023145ˆt t t t t t t t t t c c c C c c c +∆+∆=++++++Q Q M a a a a a a

（2）求解时间t t +∆的位移

ˆT t t t t

+∆+∆=LDL a Q （3）计算时间t t +∆的加速度和速度

()02367t t t t t t t t t t t t t

c c c c c +∆+∆+∆+∆=---=++a a a a a a a a a

3.程序设计思路

（1） 读入质量矩阵、刚度矩阵、阻尼矩阵、载荷列向量；读入初始状态值，如初始位移、初始速度、初始加速度；读入控制变量，如时间步长、时间步； （2） 计算8个积分常数；

（3） 计算有效刚度矩阵ˆK

并对有效刚度矩阵ˆK 进行LU 分解； （4） 求解时间t t +∆的有效载荷向量ˆt t

+∆Q ； （5） 求解时间t t +∆的位移t t +∆a ；

（6） 求解时间t t +∆的的加速度t t +∆a 和速度t t +∆a （7）计算下一时刻的有效载荷向量并循环。

4.问题描述与计算结果

4.1 问题描述

一个三自由度系统。它的运动方程是

100210003014200010226-⎡⎤⎡⎤⎛⎫ ⎪⎢⎥⎢⎥-- ⎪⎢⎥⎢⎥

⎪⎢⎥⎢⎥-⎣⎦⎣⎦⎝⎭

a +a =

初始条件：当0t =时，000,0==a a 。

已知此系统的固有频率是

：

123ωωω=

==。相应的振动周期是：12310.89, 4.44, 3.628T T T ===.

4.2 计算结果

取时间步长为3/100.363t T ∆==计算得到 时间

t ∆

2t ∆

3t ∆

4t ∆

5t ∆

1a 0.0002 0.0027 0.0155 0.0589 0.1701 2a 0.0078 0.0594 0.2275 0.5997 1.2402 3a

0.3714

1.3973

2.8413

4.3935

5.7660

时间

6t ∆

7t ∆

8t ∆

9t ∆

10t ∆

1a 0.4012 0.8066 1.4227 2.2454 3.2137 2a 2.1579 3.2934 4.5282 5.7147 6.7165 3a

6.7784

7.4035

7.7604

8.0568

8.5045

附件：

DIMENSION M(3,3),C(3,3),K(3,3),CC(8) !定义质量，阻尼比，刚度，积分参数数组

DIMENSION A0(3),A1(3),A2(3),A0T(3),A1T(3),A2T(3),Q(3),P(3) !定义位移，速度，加速度,力向量DIMENSION EQK(3,3),EQQ(3) !定义有效刚度矩阵和有效载荷矩阵

DIMENSION EQC1(3),EQC2(3) !定义参数向量

PARAMETER(AA=0.25,DERTA=0.50)

DO I=1,3

A0(I)=0

A1(I)=0

A2(I)=0

END DO

A2(3)=6 !设置初始位移，初始速度，初始加速度

OPEN(5,FILE='INPUT.dat',STATUS='OLD')

OPEN(6,FILE='OUTPUT.DAT')

READ(5,*) DT !读入时间步长

READ(5,*) NUM !读入时间步数

WRITE(6,100)

100FORMAT(/2X,'dt=')

WRITE(6,110)(DT)

110FORMAT(4X,F7.4)

READ(5,*)(M(I,1),M(I,2),M(I,3),I=1,3) !读入质量矩阵

WRITE(6,120)

120FORMAT(/2X,'MASS=')

WRITE(6,130)((M(I,J),J=1,3),I=1,3)

130FORMAT(4X,3I4)

READ(5,*)(C(I,1),C(I,2),C(I,3),I=1,3) !读入阻尼矩阵

WRITE(6,140)

140FORMAT(/2X,'C=')

WRITE(6,130)((C(I,J),J=1,3),I=1,3)

READ(5,*)(K(I,1),K(I,2),K(I,3),I=1,3) !读入刚度矩阵

WRITE(6,150)

150FORMAT(/2X,'K=')

WRITE(6,130)((K(I,J),J=1,3),I=1,3)

CC(1)=1/(AA*DT*DT)

CC(2)=DERTA/(AA*DT)

CC(3)=1/(AA*DT)

CC(4)=1/(2*AA)-1

CC(5)=DERTA/AA-1

CC(6)=DT/2*(DERTA/AA-2)

CC(7)=DT*(1-DERTA)

CC(8)=DERTA*DT !计算积分常数