参赛队号 # 1753
ai? bixj yjxm ym1 yj1 ym?xjym?xmyj ?yj?ym ?i, j, m?
??ci?1 xj1 xm?xm?xj
上式?i, j, m?表示下标轮换，如i?j,j?m,m?i。以下同此。
将求得的函数?1?t?~?6?t?代入（5-1）式可将位移函数表示成结点位移的函数，即
??u?t??Niui?t??Njuj?t??Nmum?t? ? （5-4）
??v?t??Nivi?t??Njvj?t??Nmvm?t?其中
1Ni?? ?ai?bix?ciy? ?i, j, m2ANi，Nj，Nm称为单元的插值函数或形函数。ai，bi，ci，…，cm是常数，取
决于单元的3个节点坐标。
（5-4）式的矩阵形式是
?u?t???Niu??????v?t???00NiNj00NjNm0?ui?t????vt???i?0???uj?t????? ?Nm??vj?t???um?t??????vm?t??????INiINje?? （5-5） ?NaINm?a a a??ijm?T确定单元位移以后，可以很方便的利用几何方程求得单元的应变。单元的应变
??x???ε???y??LNae?L??Ni????xy?
NjeNm?a????LNiLNjeLNm?a?
7
参赛队号 # 1753
???BiBjee （5-6） Bm?a?Ba?B称为应变矩阵，L是微分算子。
?????x?Bi?LNi??0?????y??0?????Ni0?y???????x????Ni???x0????0?Ni????Ni??y??0???Ni? ??y??Ni??x??同理可得Bj,Bm，所以三结点单元的应变矩阵B
B=??BiBj?bi1?Bm???2A?0?ci?0bjci0bicj0cjbjbm0cm0??cm? （5-7） bm??在应变梯度较大的部位，单元的划分应该适当密集，否则不能反应应变的真实变化
而导致较大的误差。 5.2.1 控制方程的建立
有弹性力学的知识得到三维弹性动力学的基本方程是：
平衡方程：?ij,j?fi??ui,tt??ui,t?0 （在V域内） （5-8）
1几何方程：?ij?(ui,j?uj,i) （在V域内） （5-9）
2物理方程：?ij?Dijkl?kl （在V域内） （5-10）
边界条件：ui?ui （在Su边界上） （5-11）
?ijnj?Ti （在Sσ边界上） （5-12）
初始条件：ui(x,y,z,0)?ui(x,y,z) （5-13） ui,t(x,y,z,0)?ui,t(x,y,z) （5-14） （5-8）式中，ρ是质量密度，μ是阻尼系数，μi,tt和μi,t分别是μi对t的二次导数和一次导数，即分别表示i方向的加速度和速度；—ρμi,tt和—μμi,t分别代表惯性力和阻尼力。它们作为体积力的一部分出现在平衡方程中，是弹性动力学和静力学相区别的基本特点之一。由于在动力学中载荷是时间的函数，以上各式中的各个符号位移、应力、应变是时间的函数。因此，在动力学的定解问题中加入了初始条件（5-13）式。
以上各式联立可以解决动力学问题，但是偏微分方程不能直接解出，必须引用有限元的方法。
在动力学分析中，因为引入了时间坐标，处理的是四维（x,y,z,t）问题。在有限元分析中一般采用部分离散的方法，即只对空间域进行离散。如图5.2所示。
平衡方程（5-8）式及力的边界条件（5-11）、（5-12）式的等效积分形式的伽辽金提法可表示如下：
8
??u(?对上式的第1项??u?ViViij,j?fi??ui,tt??ui,t)dV??S?参赛队号 # 1753
?ui(?ijnj?Ti)ds?0 （5-15）
ij,jdV进行分部积分，并带入物理方程，则从上式可以得到
V?(??ijDijkl?kl??ui?ui,tt??ui?ui,t)dV
???uifidV???uiTidV （5-16）
VS?将空间离散后的位移表达式（5-12）（现在情况下，u1=u，u2=v）代入上式，并注意到结点位移变化δa的任意性，最终得到系统的控制方程（在动力学问题中，又称运动方程）如下：
Ma(t)?Ca(t)+Ka(t)=Q(t) （5-17）
其中a(t)和a(t)分别是系统的结点加速度向量和结点速度向量，M，C，K和Q(t)分别是系统的质量矩阵、阻尼矩阵、刚度矩阵和结点载荷向量，并分别由各自的单元矩阵和向量集成，即
??????M??MeeC??Cee
（5-18）
K??KeeQ??Qee其中
M???NNdVVeeTKe??BDBdVVeeTVeS?TC???NNdVVeeT （5-19）
TQ??NfdV??eNTdVMe，Ce，Ke和Qe分别是单元的质量矩阵、阻尼矩阵、刚度矩阵和载荷向量。
5.2.3 求解方法
本计算过程中采用Newmark方法求解控制方程。
在t～t+Δt的时间区域内，Newmark积分方法采用下列的假设：
?t??t?a?t?[(1??)??t????t??t]?t （5-20） aaa?t?t?[(1??)??t????t??t]?t2 （5-21） at??t?at?aaa其中α和δ是按积分精度和稳定行要求决定的参数。另一方面，α和δ取不同数值
则代表了不同的数值积分方案。当α=1/6和δ=1/2时，（5-20）和（5-21）式相应于线性加速度法，因为这是他们可以由下式，即时间间隔Δt内线性假设的加速度表达式的积分得到。
??t?????t?(??t??t???t)?/?taaaa(0????t) （5-22）
当α=1/4和δ=1/2时，Newmark方法相当于常平均加速度法这样一种无条件稳定的
积分方案。此时，Δt内加速度为：
??at???1(??at??t???at) （5-23） 29
Newmark方法中相同时间t+Δt的位移解答at+Δt是通过满足时间t+Δt的运动方程得
参赛队号 # 1753
到的。即由：
???Mat??t?Cat??t+Kat??t=Qt??t （5-24）
而得到的。为此首先从（5-21）式解得
??at??t?111?(a?a)?a?(?1)??at （5-25）
??t2t??tt??tt2??t、??t计算at??t的将上式代入（5-20）式，然后再一并代入（5-24）式，则得到从at、aa两步递推公式
(K?1??t2M????tC)at??t?Qt??t?M[11??t2at???ta?1t?(2??1)?a?t]? C[???ta?t?(??1)a??t?(2??1)?t?a?t]至此，可将利用Newmark方法逐步求解控制方程的算法步骤归结如下： 1 初始计算
(1) 形成刚度矩阵K、指令矩阵M和阻尼矩阵C。
(2) 给定a0、a?0和?a?0（?a?0由?a?0?M?1(Q0?Ca?0?Ka0)确定） (3) 选择时间步长Δt及参数α和δ，并计算积分常数。
这里要求：??0.50,??0.25(0.5??)2
c?1
??t2,c1??110??t,c2???t,c3?2??1c??t?
4???1,c5?2(??2),c6??t(1??),c7???t(4) 形成有效刚度矩阵K?: K??K?c0M?c1C (5) 三角分解K?:K??LDLT 2 对于每一时间步长
(1) 计算时间t+Δt的有效载荷
Q??Qt??t?M(c0at?c2a?t?c3?a?t)?C(c1at?c4a?t?c5?a?t)
(2) 求解时间t+Δt的位移
LDLTat??t?Q?t??t (3) 计算时间t+Δt的加速度和速度
?a?t??t?c0(at??t?at)?c2a?t?c3?a?t a?t??t?a?t?c6?a?t?c7?a?t??t 10
(5-26)
参赛队号 # 1753
6、模型的求解
6.1 客机模型
在客机迫降过程中，主要是机身下部和水面接触，因此忽略了不与水体接触的部件。在ABAQUS中按照1：1建立了客机主要部件的有限元模型，共分为五部分：机头、机舱、机尾、机翼和舱门。飞机模型的拉格朗日有限单元由二维壳单元和一位梁单元组成。如图6.1所示。
图6.1 客机模型网格图
发动机、设备等在有限元模型中没有办法真实建模,而由壳单元、梁单元和杆单元建立的模型的总体质量与真实模型存在差距.根据质量报告,将缺少的部分采用集中质量单元添加上去。 6.2计算工况
研究客机迫降姿态的影响，为此，把固定客机入水时的速度，本模型中，入水速度均为:X轴方向15 m/s,Z轴方向0.3 m/s。
客机迫降时，客机还可以控制，因此，入水角度不可能较大，应在很小的范围内变化。本文选取了四种工况进行计算，旨在说明迫降姿态对客机影响的趋势。四种姿态的攻角分别为：5°、10°、12°、15°。
11