有限元法基础讲稿第14讲_经济法基础讲稿
有限元法基础讲稿第14讲由刀豆文库小编整理,希望给你工作、学习、生活带来方便,猜你可能喜欢“经济法基础讲稿”。
青岛大学讲稿
讲
授
内
容
备
注 第10讲(第17周)4.1 结构动力学问题有限元方法
动力学问题在国民经济和科学技术的发展中有着广泛的应用领域。最经常遇到的是结构动力学问题,它有两类研究对象:一类是在运动状态下工作的机械或结构,例如高速旋转的电机、汽轮机、离心压缩机,往复运动的内燃机、冲压机床,以及高速运行的车辆、飞行器等,它们承受着本身惯性及与周围介质或结构相互作用的动力载荷。如何保证它们运行的平稳性及结构的安全性,是极为重要的研究课题。另一类是承受动力载荷作用的工程结构,例如建于地面的高层建筑和厂房,石化厂的反应塔和管道,核电站的安全壳和热交换器,近海工程的海洋石油平台等,它们可能承受强风、水流、地震以及波浪等各种动力载荷的作用。这些结构的破裂、倾覆和垮塌等破坏事故的发生,将给人民的生命财产造成巨大的损失。正确分析和设计这类结构,在理论和实际上也都是具有意义的课题。
动力学研究的另一重要领域是波在介质中的传播问题。它是研究短暂作用于介质边界或内部的载荷所引起的位移和速度的变化,如何在介质中向周围传播,以及在界面上如何反射、折射等的规律。它的研究在结构的抗震设计、人工地震勘探、无损探伤等领域都有广泛的应用背景,因此也是近20多年一直受到工程和科技界密切关注的课题。
现在应用有限单元法和高速电子计算机,已经可以比较正确地进行各种复杂结构的动力计算,本章阐明如何应用有限单元法进行动力分析。
4.1.1 运动方程
结构离散化以后,在运动状态中各节点的动力平衡方程如下
Fi +Fd +P(t)=Fe
(2-2-1)
式中:Fi、Fd、P(t)分别为惯性力、阻尼力和动力荷载,均为向量;Fe为弹性力。
弹性力向量可用节点位移δ和刚度矩阵K表示如下
Fe =K δ
式中:刚度矩阵K的元素Kij为节点j的单位位移在节点i引起的弹性力。
根据达朗贝尔原理,可利用质量矩阵M和节点加速度下
FiMδt22δt22表示惯性力如
式中:质量矩阵的元素Mij为节点j的单位加速度在节点i引起的惯性力。
设结构具有粘滞阻尼,可用阻尼矩阵C和节点速度
FdCδt2δt表示阻尼力如下
式中:阻尼矩阵的元素Cij为节点j的单位速度在节点i引起的阻尼力。
将各力代入式(2-2-1),得到运动方程如下
Mδt22CδtKδP(t)
(2-2-2)
记
δ,δδ δ2tt2则运动方程可写成C δK δP(t)
(2-2-3)M δ,结构各节点的在地震时,设地面加速度为a,结构相对于地面的加速度为δ,在计算惯性力时须用它代替式(2-2-3)中的δ。至于弹性力实际加速度等于a+δ,与地和阻尼力,则分别取决于结构的应变和应变速率,即取决于位移δ和速度δ面加速度无关。
2.2.2 质量矩阵
下面用m表示单元质量矩阵,M表示整体质量矩阵。求出单元质量矩阵后,进行适当的组合即可得到整体质量矩阵。组合方法与由单元刚度矩阵求整体刚度矩阵时相似。
在动力计算中可采用两种质量矩阵,即协调质量矩阵和集中质量矩阵。1.协调质量矩阵
从运动的结构中取出一个微小部分,根据达朗贝尔原理,在它的单位体积上作用的惯性力为
pirt22
式中:ρ为材料的密度。
在对结构进行离散化以后,取出一个单元,并采用如下形式的位移函数
rN δ
e则
piNδt22e
再利用荷载移置的一般公式求得作用于单元节点上的惯性力为
FieNTpidVNNdVTδt22e
即
ee Fimδ可见,单元质量矩阵为
mNNdV
(2-2-4)
T如此计算单元质量矩阵,单元的动能和位能是互相协调的,因此叫做协调质量矩阵。
2.集中质量矩阵
假定单元的质量集中在它的节点上,质量的平移和转动可同样处理。这样得到的质量矩阵是对角线矩阵。
单元集中质量矩阵定义如下:
mdV
(2-2-5)
T式中,为函数i的矩阵,i在分配给节点i的区域内取l,在域外取0。
由于分配给各节点的区域不能交错,所以由上式计算的质量矩阵是对角线的。
3.平面等应变三角形单元集中质量矩阵与协调质量矩阵
设单元重量为W,将它3等分,分配给每一节点,得到单元集中质量矩阵如下
10W0m3g00***0100000010000 001
(2-2-6)
单元协调质量矩阵为
1201W4m3g0140******140120014014012
(2-2-7)
在单元数目相同的条件下,两种质量矩阵给出的计算精度是相差不多的。集中质量矩阵不但本身易于计算,而且由于它是对角线矩阵,可使动力计算简化很 3 多。对于某些问题,如梁、板、壳等。由于可省去转动惯性项,运动方程的自由度数量可显著减少。当采用高次单元时,推导集中质量矩阵是困难的。另外,只要离散化时保持了单元之间的连续性,由协调质量矩阵算得的频率代表结构真实自振频率的上限。
2.2.3 阻尼矩阵
如前所述,结构的质量矩阵[M]和刚度矩阵[K]是由单元质量矩阵[m]和单元刚度矩阵[M]e经过集合而建立起来的。相对来说,阻尼问题比较复杂,结构的阻尼矩阵[C]不是由单元阻尼矩阵经过集合而得到的,而是根据已有的实测资料,由振动过程中结构整体的能量消耗来决定阻尼矩阵的近似值。
1.单自由度体系的阻尼
单自由度体系的自由振动方程为
ck0 m式中:m为质量;c为阻尼系数;k为刚度系数;δ为变位。
上式两边除以m后得到
20 2其中,频率)。
ζ称为阻尼比,ω为体系的自振频率(角k/m,c/2m,设初始条件为:当t=0时,δ=δ0,=v0,符合这些初始条件的解为
expt0cosdtv00d2sindt
(2-2-8)
d1
体系的自振频率为ωd,其振幅随着时间而逐渐衰减。
根据实测资料,大多数结构的阻尼比都是很小的数,较多为ζ=0.01~0.10,一般都小于0.20。可见,阻尼对自振频率的影响是很小的,通常可取ωd=ω。
2.多自由度体系的阻尼
如果假定阻尼力正比于质点运动速度,从运动的结构中取出一微小部分,在它的单位体积上作用的阻尼力为
pdte rNδ式中:α为比例常数;ρ为材料密度;N为形函数。
利用荷载移置的一般公式求得作用于单元e的节点上的阻尼力如下
FdeNTpddVe NNdVδT即
ee FdC δ 4 而
CNNdVm
(2-2-9)
T可见,此时单元阻尼矩阵正比于单元质量矩阵。如果假定阻尼力正比于应变速度,则阻尼应力可表为
ζdDεte DBδ所以作用于单元e的节点上的阻尼力为
FedBTδddVBTeC δe DBdVδ其中
CBTeKe
(2-2-10)DBdVδ可见,此时单元阻尼矩阵正比于单元刚度矩阵K e。
前面已经说过,通常是根据实测资料,由振动过程中结构整体的能量消耗来决定阻尼的近似值,因此不是计算单元阻尼矩阵,而直接计算结构的整体阻尼矩阵C。一般采用如下的线性关系,并称为瑞利(Rayleigh)阻尼,即
CMK
(2-2-11)
其中的系数α和β根据实测资料决定。
现在说明如何计算α和β。设φi和φj为两个振型。对式(2-2-11)的两边先后乘以φi,再前乘以φT得到 jφjCφiφjMφiφjKφi
(2-2-12)
TTT根振型正交性再由式(2-2-12)得到
φjCφi0φjCφiTT2jm pjij ij其中
mpjφjMφj
T令
则
2j2jj
(2-2-13)
φjCφj2jjmpj
T由式(2-2-13)得到
j2j2j
(2-2-14)
实测两个阻尼比即可求解α与β。结构动力学方程主要采用振型叠加法和直接积分法。前者用到振型正交条件,但不同的振型之间不能解耦时(在结构与地基的相互作用问题中,地基的阻尼往
往大于结构本身的阻尼,对于结构和地基应分别给以不同的α与β值),应采用直接积分法求解。
2.2.4 结构自振频率与振型
在式(2-2-3)中,令P(t)=0,得到自由振动方程。在实际工程中,阻尼对结构自振频率和振型的影响不大,因此可进一步忽略阻尼力,得到无阻尼自由振动的运动方程
0
(2-2-15)K δM δ设结构作下述简谐运动
δφcost
把上式代人式(2-2-15),可得到齐次方程
(KM)φ0
(2-2-16)
2在自由振动时,结构中各节点的振幅{Ф}不全为零,所以结构自振频率方程为
KM0
(2-2-17)
结构的刚度矩阵[K]和质量矩阵[M]都是n阶方阵,其中n是节点自由度的数目,所以上式是关于ω2的n次代数方程,由此可求出结构的自振频率
ω1≤ω2≤ω3≤…≤ωn
对于每个自振频率,由式(2-2-16)可确定一组各节点的振幅值 i=[ i1,
Ti2,…, in],它们互相之间应保持固定的比值,但绝对值可任意变化,它们构成一个向量,称为特征向量,在工程上通常称为结构的振型。
因为在每个振型中,各节点的振幅是相对的,其绝对值可取任意数值。在实际工作中,常用以下两种方法之一来决定振型的具体数值:
(1)规准化振型:取 i的某一项,例如取第n项为1,即 in=1,于是
i=[ i1, i2,…,1]T
(2-2-18)
这样的振型称为规准化振型。
(2)正则化振型:选取 ij的数值,使
φiMφi
1(2-2-19)
T2这样的振型称为正则化振型。
设已求得一振型φii1,i2,,in,如令
Tjiij/in
(2-2-20)
则得到的φii1,i2,,in为规准化振型。如令
Tjiij/c
(2-2-21)
cφiMφiTT1/2
则得到的φii1,i2,,in为正则化振型。
令
mpiφiTMφi
(2-2-22)
当M为集中质量矩阵时,则
m10in00m200i10i2 mninmpii12i2ms1s2is
当φi为正则化振型时,有
mpi=1 令
kpiφiTKφiφiTi2Mφii2mpi
(2-2-23)
式中,mpi和kpi分别称为第i阶振型相应的广义质量和广义刚度。由式(2-2-23)得
ikpi/mpi
(2-2-24)
[例2-3]求解K =ω2M的振型,其中
2K10141021,M12014101 2求解说明
频率方程为
20.5KM2214120120.52100
求得三个自振频率为
22212,24,36
将122代入式(2-2-16)中,得到第1振型必须满足的方程组如下
11-12+0=0,-11+212-13=0,11-12+13=0 联立前两个方程解出
11=13,12=13
取13=1,得到规准化的第一振型为
1=[1 1 1]T
用同样方法得到第2、3振型为
2=[-1 0 1]T 3=[1-1 1]T
由式(2-2-21)得到正则化振型如下
1=[1/
21/2
1/
2=[-1 0 1]T 3=[1/2
-1/2
1/
2]T ]T
22.2.5 振型叠加法求解结构的受迫振动
目前,常用的求解结构受迫振动的方法有两种,即振型叠加法和直接积分法。用振型i的线性叠加来表示处于运动状态中的结构位移向量
δφ11tφ22tφnntφt
(2-2-25)
iii1n用φTM前乘上式的两边,由于振型正交性,等式右边的n项中只剩下i=jj这一项,即
φjMδTjtφTMφjjmpjjt
由此得到
itφiMδmpiT
(2-2-26)
i的初始值可表示为 i和φiMδ(0)mpiTi0
(2-2-27)
i0T(0)φiMδmpi
(2-2-28)
现在考虑下列运动方程的求解:
C δK δP(t)M δ把式(2-2-25)代入上式,得到
nniinMφi1iKCφii1φii1iP(t)
对上式两边前乘以φT,并令C=αM+βK,得到 jnφi1TjiMφiφMTji1niKφinφi1TjTKφiiφjP(t)由于振型正交性,得到
iimpiiimpiiφjP(t)mpi22T由于i22ii,上式进一步化为
i2iiiii21mpiφjP(t)i1,2,3,,n
(2-2-29)
T这是二阶常微分方程,这样的方程共有n个,它们是互相独立的。式(2-2-29)在形式上与单自由度体系的运动方程相同。其解答可用数值积分方法计算,也可用Duhamel积分计算如下:
it1dimpiiitt0P* eiitsinditd
(2-2-30)
ei0iii0sinditi0cosditdi其中
2di1i
P*tφiP(t)
T把ηi(t)代人式(2-2-25),即得到所需解答。在用有限元方法进行结构动力分析时,自由度数目n可以达到几百甚至几千,但由于高阶振型对结构动力反应的影响一般都很小,通常只要计算一部分低阶振型就够了。例如,对于地震荷载,一般只要计算前面5~20个振型。对于爆炸和冲击荷载,就需要取更多的振型,有时需取出多达2n/3个振型进行计算,而对于振动激发的动力反应,有时只有一部分中间的振型起作用。
运动方程(2-2-3)是二阶常微分方程组,可用数值积分方法直接求解。应用于动力问题的直接积分方法很多,有线性加速度方法、Wilson方法、Newmark方法等,此不赘述。