2011数学中国数学建模网络挑战赛A题特等奖论文要点 下载本文

参赛队号 # 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