CN110162826A - 薄壁结构热气动弹性动响应分析方法 - Google Patents

Abstract

本发明提出一种薄壁结构热气动弹性动响应分析方法，根据薄壁结构和边界条件，将薄壁结构离散为七个自由度壳单元，将加肋离散为七个自由度梁单元，建立薄壁的有限元模型，得到结构有限元网格；将薄壁结构温度场施加到结构有限元网格上，对有限元模型进行线性颤振分析，得到每个减缩频率k对应的非定常气动力矩阵Q_k，拟合出非定常气动力矩阵的时域表达式；对得到的非定常气动力矩阵的时域表达式进行迭代计算，进行薄壁的非线性瞬态响应分析，得到薄壁结构非线性热气动弹性动响应。本发明的方法可以解决复杂薄壁结构的非线性热气动弹性动响应的计算问题。

Description

薄壁结构热气动弹性动响应分析方法

技术领域

本发明涉及热气动弹性响应分析，属于飞行器气动弹性技术技术领域。

背景技术

在超声速、高超声速飞行器上广泛使用大尺寸薄壁结构，如热防护结构、发动机保护罩和大尺寸翼面蒙皮等，这类薄壁结构容易发生弹性力、惯性力、气动力和热应力四者之间相互耦合引发的热气动弹性问题，这种耦合往往导致飞行器的性能下降甚至结构破坏。另外，由于薄壁结构横向挠度在板中面引起面内应力，面内应力在振动中产生一种硬弹性的效应，这是一种典型的结构几何非线性问题。

对于壁板颤振的经典方法如Galerkin方法、Rayleigh-Ritz方法、谐波平衡方法等仅均限与形状规则的矩形薄板，且边界条件为简支边界或者固支边界条件，对于形状不规则、边界条件复杂的情况，经典方法不适用。

采用有限元方法不仅适用于形状不规则、边界条件复杂的薄壁结构，而且能够考虑气动加热、面内内力及不同来流方向等复杂问题。航空航天领域常用的有限元分析软件MSC.Nastran提供了结构线性气动弹性响应的计算功能，但无法进行结构非线性的气动弹性动响应分析。

发明内容

本发明的目的在于克服现有技术不足，提供一种薄壁结构非线性气动弹性动响应分析方法，可以获得薄壁结构的响应和非定常气动力随时间变化历程。

本发明的技术解决方案：

一种薄壁结构热气动弹性动响应分析方法，包括以下步骤：

步骤一，根据薄壁结构和边界条件，将薄壁结构离散为空间壳单元，将加肋离散为梁单元，建立薄壁的有限元模型，所述的壳单元每个节点有七个自由度Δ_C＝[u v w w_,x w_,yw_,xy θ_z]^T，其中，u、v表示壳单元两个方向的面内位移，w、w_,x、w_,y、w_,xy分别表示壳单元四个弯曲位移，θ_z表示壳单元扭转角；

所述的梁单元每个节点有七个自由度，Δ_L＝[u′ v′ w′ θ_x θ_y θ′_z θ_x,x]

其中u′、v′表示梁单元两个方向的面内位移，w'表示梁单元弯曲位移，θ_x，θ_y，θ′_z分别代表三个方向的角位移，θ_x,x表示梁单元的扭转角；

步骤二，根据薄壁结构建立气动面网格，并建立气动面网格和结构有限元模型之间的插值关系，得到结构有限元网格；

步骤三，将薄壁结构温度场施加到步骤二中得到的结构有限元网格上，对有限元模型进行线性颤振分析，得到每个减缩频率k对应的非定常气动力矩阵Q_k并拟合出Q_k的时域表达式；

步骤四，通过对步骤三得到的非定常气动力矩阵的时域表达式进行迭代计算，进行薄壁的非线性瞬态响应分析，得到薄壁结构非线性热气动弹性动响应；

所述的步骤三中非定常气动力矩阵的时域表达式为：其中Q₀-Q_2+j是实待定矩阵，s是拉氏复变量，r_j是由经验给定的正实数,r₁＝0.1～0.25，r₂＝0.3～0.45，r₃＝0.5～0.65，r₂＝0.7～0.9，b为薄壁结构沿来流方向的长度，V为飞行器的速度。

所述的薄壁结构非线性热气动弹性动响应包括薄壁结构非线性热气动弹性气动力和薄壁结构的位移。

所述的薄壁结构的位移求解表达式为：

w＝H_wT_bw_b，u＝H_uT_mw_m，v＝H_vT_mw_m

其中

其中，w_m＝{u,v}^T，w_b＝{w,w_,x,w_,y,w_,xy,θ_z}^T，为薄壁结构沿来流方向的宽度，为薄壁结构沿来流方向的长度；

H_w＝[1 x y x² xy y² x³ x²y xy² y³ x³y x²y² xy³ x³y² x²y³ x³y³]，H_u＝[1 x yxy 0 0 0 0]，H_v＝[0 0 0 0 1 x y xy]，x,y,z分别表示三个坐标方向。

所述的薄壁结构非线性热气动弹性气动力表达式为：

式中，q_a为气流动压，M_∞为飞行器马赫数，γ为比热比，w_,t为弯曲运动速度，H_w,x为矩阵H_w对x的导数。

本发明与现有技术相比的有益效果：

(1)本发明通过构建加肋壁板结构的壳单元和空间梁单元，并在壳单元和空间梁单元上增加了表示扭转自由度的扭转角，可以更真实模拟非线性结构的运动效应，通过扭转自由度的加入，将扭转角引入到了单位刚度矩阵中，应用在有限元模型的建立中，由此方法建立的有限元模型对复杂结构形式和混合边界条件的薄壁具有很好适应性，适于对飞行器薄壁结构的气动弹性动响应分析；

(2)本发明提出了非线性壁板气动弹性动响应瞬态过程的非定常气动力求解方法，可直接用于评估壁板颤振对局部气动特性的影响；

(3)本发明将时域气动力矩阵代入到非线性瞬态响应分析程序中，兼顾了非线性瞬态响应，可以实现薄壁结构几何非线性特性的模拟、同时引入了非定常气动力，实现了薄壁结构的非线性热气动弹性动响应分析。

附图说明

图1为本发明非线性热气动弹性动响应分析方法流程图；

图2为本发明薄壁结构示意图；

图3为本发明薄壁结构壳单元；

图4为本发明薄壁结构空间梁单元；

图5为本发明薄壁结构随时间变化的位移响应图；

图6为本发明薄壁结构气动弹性动响应相位图。

具体实施方式

下面结合具体实例及附图对本发明进行详细说明。

根据图1所示，一种薄壁结构热气动弹性动响应分析方法，包括以下步骤：

步骤一，根据薄壁结构和边界条件建立薄壁的有限元模型，薄壁结构如图2所示，薄壁结构离散为空间壳单元，如图3所示，对于加筋薄壁结构，还需要将加筋离散为空间梁单元，如图4所示。

对于薄壁结构，每个节点有七个自由度Δ_C＝[u v w w_,x w_,y w_,xy θ_z]^T，包括两个面内位移w_m＝{u,v}^T和四个弯曲位移w_b＝{w,w_,x,w_,y,w_,xy,θ_z}^T，单元如图5所示；空间梁单元每个节点有七个自由度Δ_L＝[u′ v′ w′ θ_x θ_y θ′_z θ_x,x]，如图4所示其中θ_z和θ_x,x均代表扭转角。

对于壳单元，位移函数可以表示为：

w＝α₁+α₂x+α₃y+α₄x²+α₅xy+α₆y²+α₇x³+α₈x²y+α₉xy²+

α₁₀y³+α₁₁x³y+α₁₂x²y²+α₁₃xy³+α₁₄x³y²+α₁₅x²y³+α₁₆x³y³

＝H_wα

u＝ζ₁+ζ₂x+ζ₃y+ζ₄xy＝H_uζ

v＝ζ₅+ζ₆x+ζ₇y+ζ₈xy＝H_vζ

其中H_w＝[1 x y x² xy y² x³ x²y xy² y³ x³y x²y² xy³ x³y² x²y³ x³y³]，α＝[α₁α₂ … α₁₆]^T，H_u＝[1 x y xy 0 0 0 0]，H_v＝[0 0 0 0 1 x y xy]，ζ＝[ζ₁ ζ₂ … ζ₈]^T。

为了确定α和ζ，可将图3中结点1,2,3,4的坐标代入下面的表达式

其中，w_i、w_i,x、w_i,y、w_i,xy表示结点i的弯曲位移，H_wi、H_wi,x、H_wi,y、H_wi,xy表示结点i的弯曲位移的计算系数，u_i、v_i表示结点i的面内位移，H_ui、H_vi表示结点i的面内位移的计算系数。

对于空间梁单元，采用线性插值计算，公式如下：

u'＝u₁(1-ξ)+u₂ξ

v'＝v₁F₁+θ_z1lF₂+v₂F₃+θ_z2lF₄

w'＝w₁F₁+θ_y1lF₂+w₂F₃+θ_y2lF₄

θ_x＝θ_x1F₁+(θ_x,x)₁lF₂+θ_x2F₃+(θ_x,x)₂lF₄

其中l为梁单元长度，ξ＝x/l,F₁＝1-3ξ²+2ξ³,F₂＝ξ-2ξ²+ξ³,F₃＝3ξ²-2ξ³,F₄＝ξ³-ξ²，u₁、u₂、v₁、v₂为x和y方向的位移、F₁、F₂、F₃、F₄为系数、θ_z′₁、θ_z′₂、θ_y1、θ_y2、θ_x1、θ_x2为两个节点在三方向的角位移、(θ_x,x)₁、(θ_x,x)₂为两个节点的扭转角，下标1和2分别代表空间梁单元的两个节点。

设F＝(F₁ F₂ F₃ F₄)，则梁单元的弯曲刚度矩阵为：

其中，EI为弯曲刚度。

则梁单元的扭转刚度矩阵：

其中，GJ为扭转刚度。

因为拉压、扭转和双向弯曲互不耦联，从而可由拉压、扭转和两个方向的弯曲单元根据位移编号顺序通过组合来形成单元刚度矩阵，单元刚度矩阵为对称矩阵，表示如下：

其中a_L＝EA/l,b_L＝12EI_z/l³,c_L＝6EI_z/l²,d_L＝12EI_y/l³,e_L＝6EI_y/l²,f_L＝6GJ/5l,g_L＝GJ/10,h_L＝2EI_y/l,i_L＝2EI_z/l,j_L＝2GJl/15,k_L＝-GJl/30，EA为单元轴向拉压刚度，EI_z为xoy平面内的弯曲刚度，EI_y为xoz平面内的弯曲刚度，I_y是筋条对面板中面的惯性矩，设A₁为筋条在yoz平面的截面面积，GJ为扭转刚度，J＝I_y+I_z，通过单元刚度矩阵可以得到结构的整体刚度矩阵K_s，在有限元模型的建立时使用；

步骤二，根据薄壁结构建立气动面网格，并建立气动面网格和结构有限元模型之间的插值关系，插值矩阵为G_ka，在MSC.Nastran的前后处理软件MSC.Patran将壁板结构温度场施加到结构有限元网格上，采用MSC.Nastran求解序列SOL103对有限元模型进行线性颤振分析，得到不同减缩频率k下的非定常气动力矩阵Q_k，利用Roger有理函数拟合的方法拟合出非定常气动力的时域表达式，非定常气动力的有理函数表达式为：

其中Q₀-Q_2+j是实待定矩阵，s是拉氏复变量，r_j是由经验给定的正实数，r₁＝0.2,r₂＝0.4,r₃＝0.6,r₄＝0.8；

在一个实施例中，减缩频率k取为颤振减缩频率附近的若干k值，每个k值对应一个非定常气动力矩阵，矩阵包含实部和虚部，计算得到的一组k下的非定常气动力矩阵Q_k，将简谐条件带入式，并分解为实部、虚部，其中是拉氏复变量、e^iωk是简谐条件的表达方式。在一个实施例中，取7个为一组，则Q_k中每个元素q的待定系数为7个，其中实部为虚部为

采用最小二乘法可求得系数q₀-q₆，由此完成了气动力从频域到时域的转化，得到时域非定常气动力。

壁板气动弹性动响应系统质量矩阵为M＝M_s-M_a，阻尼矩阵为C＝C_s(t)-C_a，刚度矩阵为K＝K_s(t)-K_a，其中M_s、C_s(t)、K_s(t)为结构刚度矩阵、阻尼矩阵、刚度矩阵，为气动力产生的附加质量矩阵，为气动阻尼矩阵，为气动刚度矩阵；

步骤三，在MSC.Nastran的前后处理软件MSC.Patran将壁板结构温度场施加到结构有限元网格上，并进行非线性瞬态响应分析的设置，生成求解序列为SOL129的计算文件.bdf；利用MSC.Nastran的二次开发语言DMAP将壁板系统的质量矩阵、阻尼矩阵和刚度矩阵写入非线性瞬态响应分析的计算文件.bdf中，替换有限元模型自身的质量矩阵、阻尼矩阵和刚度矩阵，作为薄壁气动弹性动响应分析的初始刚度矩阵，其中阻尼矩阵和刚度矩阵均为响应的函数，因此在每一个非线性迭代步进行更新；

步骤四，采用MSC.Nastran求解序列SOL129进行壁板非线性瞬态响应分析，可完成薄壁结构非线性热气动弹性动响应分析。

通过求解可以得到结构任意点的位移：

w＝H_wT_bw_b，u＝H_uT_mw_m，v＝H_vT_mw_m

其中

非定常气动力表示为：

结合结构响应的表达式，可得到耦合非线性结构响应的非定常气动力求解方法为：

式中，q_a＝ρ_aV²/2为气流动压，M_∞为飞行器马赫数，γ为比热比，w_,t为弯曲运动速度，H_w,x为矩阵H_w对x的导数。

本发明未详细说明部分为本领域技术人员公知技术。

Claims (4)

1.一种薄壁结构热气动弹性动响应分析方法，其特征在于：包括以下步骤，

步骤一，根据薄壁结构和边界条件，将薄壁结构离散为空间壳单元，将加肋也离散为壳单元，建立薄壁的有限元模型，所述的壳单元每个节点有七个自由度Δ_C＝[u v w w_,x w_,yw_,xy θ_z]^T，其中，u、v表示壳单元两个方向的面内位移，w、w_,x、w_,y、w_,xy分别表示壳单元四个弯曲位移，θ_z表示壳单元扭转角；

所述的梁单元每个节点有七个自由度，Δ_L＝[u′ v′ w′ θ_x θ_y θ_z′ θ_x,x]

其中u′、v′表示梁单元两个方向的面内位移，w'表示梁单元弯曲位移，θ_x，θ_y，θ_z′分别代表在x,y,z方向的角位移，θ_x,x表示梁单元的扭转角；

步骤二，根据薄壁结构建立气动面网格，建立气动面网格和结构有限元模型之间的插值关系，得到结构有限元网格；

步骤三，将薄壁结构温度场施加到步骤二中得到的结构有限元网格上，对有限元模型进行线性颤振分析，得到每个减缩频率k对应的非定常气动力矩阵Q_k，拟合出非定常气动力矩阵的时域表达式；

步骤四，通过对步骤三得到的非定常气动力矩阵的时域表达式进行迭代计算，得到薄壁结构非线性热气动弹性动响应。

2.根据权利要求1所述的薄壁结构热气动弹性动响应分析方法，其特征在于：非定常气动力矩阵的时域表达式为：

其中Q₀-Q_2+j是实待定矩阵，s是拉氏复变量，r_j是由经验给定的正实数，r₁＝0.1～0.25，r₂＝0.3～0.45，r₃＝0.5～0.65，r₂＝0.7～0.9，b为薄壁结构沿来流方向的长度，V为飞行器的速度。

3.根据权利要求1所述的薄壁结构热气动弹性动响应分析方法，其特征在于：所述的薄壁结构的位移求解表达式为：

w＝H_wT_bw_b，u＝H_uT_mw_m，v＝H_vT_mw_m

其中

其中，w_m＝{u,v}^T，w_b＝{w,w_,x,w_,y,w_,xy,θ_z}^T，为薄壁结构沿来流方向的宽度，为薄壁结构沿来流方向的长度，

H_w＝[1 x y x² xy y² x³ x²y xy² y³ x³y x²y² xy³ x³y² x²y³ x³y³]，H_u＝[1 x y xy 00 0 0]，H_v＝[0 0 0 0 1 x y xy]，x,y,z分别表示三个坐标。

4.根据权利要求1所述的薄壁结构热气动弹性动响应分析方法，其特征在于：所述的薄壁结构非线性热气动弹性气动力表达式为：

式中，q_a为气流动压，M_∞为马赫数，γ为比热比，w_,t为弯曲运动速度，H_w,x为矩阵H_w对x的导数。