一种消除自由飞配平载荷下飞机变形中刚性位移的方法，基本特征在于，给配平载荷作用的结构有限元模型施加任意合理的静定约束，利用工程商用软件计算得到结构位移数据，借助一致紧支径向基插值函数在选定的区域上重构其位移场，构造变形梯度矩阵，对变形梯度矩阵实施极分解运算得到对称矩阵和正交矩阵的乘积，对称矩阵描述了结构的伸长，正交矩阵反映了结构的旋转。根据极分解获得的刚性旋转矩阵，构造结构整体的刚性位移，继而从总位移中消除结构的整体刚性位移，从而得到自由结构的真实弹性位移。本发明方法不改变结构总刚矩阵的基本特性，运算量小、计算简单，可高效地消除在计算静变形过程中由于静定约束设置引起的结构整体刚性位移。
1.一种消除自由飞配平载荷下飞机变形中刚性位移的方法，其特征在于，计算流程如下：S1：根据实际结构的性质在有限元分析软件中建立结构有限元模型；
S2：对步骤S1中建立的有限元模型施加自平衡载荷及静定约束，然后利用NASTRAN软件求解线性静变形；
S3：根据求解得到的结点位移，利用一致紧支径向基插值函数在选定的区域上重构其位移场，并由此求解位移梯度；
S5：根据结构刚体旋转性质，利用刚性旋转矩阵构造结构的整体刚性位移；
S6：从总位移中消除结构的整体刚性位移获得其真实弹性变形。
2.根据权利要求1所述的一种消除自由飞配平载荷下飞机变形中刚性位移的方法，其特征在于：从总位移中消除结构刚性位移获得真实弹性变形，用公式表示为：ud(X)＝u(X)-ur(X) (1)式中，ud为消除刚性位移后自由结构的弹性变形；u为静定约束下结构总位移；ur为结构整体的刚性位移；X表示结构中任意一点的总体坐标，坐标原点位于结构质心处，坐标轴方向符合工程应用习惯；
ur(X)＝R(X-Xc) (2)式中，R为刚性旋转的正交阵，即RTR＝I；Xc为约束点坐标；
该刚性旋转矩阵可通过对反映结构整体刚性位移的变形梯度矩阵进行极分解得到，式中，F为结构的变形梯度，x为结构体变形后的任一点坐标；R为正交转动张量，U为对称张量，R和U可以通过下式获得：对于二维平面情况，若结构整体绕着垂直纸面向外的z轴旋转了θ角，则其旋转矩阵可以写为
3.根据权利要求2所述的一种消除自由飞配平载荷下飞机变形中刚性位移的方法，其特征在于：求解结构变形梯度需要知道结构中任意一点的位移，本申请采用一致紧支径向基函数插值来计算任意一点的位移，一致紧支径向基插值函数的形式如下：式中，s(X)是结构节点变形的插值函数；N为选定的插值节点个数；αi为插值系数；
为所采用的基函数形式；||·||为2-范数；p(X)为X的线性多项式，其表达式为：p(X)＝γ0+γ1X+γ2Y+γ3Z (10)式中，X点处的坐标为(X，Y，Z)；γi(i＝0，...，3)为插值系数；
式(9)及式(10)中插值系数需满足如下条件：式中，qj(X)为次数不超过p(X)的多项式，可以取为1、X、Y、Z；通过选择插值基函数，根据结构节点坐标和已知的位移分量，即可求解式(9)中的未知系数；将已知结构节点坐标及位移分量代入式(9)，由式(10)和式(11)可以得到：us＝Css·αs (12)式中，us为节点单方向位移分量构成的列阵，αs为待求系数列阵，具体形式可以写为：式中，si(i＝1,…,N)为各节点单向位移分量，工程上用ui，vi，wi表示节点i沿总体坐标系的各位移分量；
式中，Css为可逆矩阵， 为以节点间距离为自变量的一致紧支径向基函数；
将Css代入式(12)即可确定结构节点和位移形成的插值曲面；将任意一点坐标代入式(9)，得到该点的位移，即：式中，s为该点位移分量u、v或w；us列阵根据所求分量确定；Aas的具体形式如下：
4.根据权利要求3所述的一种消除自由飞配平载荷下飞机变形中刚性位移的方法，其特征在于：式(3)中结构任意一点的变形梯度具体形式可以写为：式中，x为结构中任意一点X变形后的坐标，x的分量为x、y、z；X的分量为X、Y、Z；
式(17)中的变形梯度表达式不易显式表达，因此写为关于位移的表达；若u为任意一点的位移，其分量为u、v、w，由于当前构型由初始构型变形所得，故而有：x＝X+u (18)将式(18)代入式(17)，可以得到变形梯度关于位移的表达式：式中，I为3×3的单位矩阵；
根据式(9)可知，其位移导数可以通过下式求得：式中，DXas、DYas和DZas为式(16)对坐标X、Y、Z求导所得，即：式中， 和 为对X、Y和Z的偏导数；
将式(20)代入式(19)中，便可以得到结构的变形梯度具体表达式，其反应了结构整体的弹性与刚性变形。
5.根据权利要求4所述的一种消除自由飞配平载荷下飞机变形中刚性位移的方法，其特征在于：将结构的变形梯度按照式(4)进行极分解得到反映结构整体刚性运动的旋转矩阵R；根据式(2)计算各个节点的刚性位移：uri＝R(Xi-Xc) i＝1,…,N (22)式中，uri为第i个节点的刚性位移；Xi为第i各节点的坐标；N为节点个数；
将静定约束下结构各个节点的总位移减去其整体刚性位移，便可以得到消除整体刚性位移后的静弹性变形位移，即自由结构受自平衡载荷的静变形：udi＝ui-uri i＝1,…,N (23)式中，udi为结构消除整体刚性位移后各节点的位移，即无约束状态下结构的位移；
ui为结构在静定约束下由有限元分析软件计算得到的总位移；uri为由式(22)计算得到的各个节点的刚性位移。
一种消除自由飞配平载荷下飞机变形中刚性位移的方法
技术领域
[0001] 本发明适用于飞机自由飞状态下受配平载荷作用的真实弹性位移精确高效解算法，更具体地，针对受配平载荷作用的飞机，结合有限元分析软件，精确快速地消除因静定约束引起的结构整体刚性位移，得到其自由飞状态下的真实弹性变形。
背景技术
[0002] 工程结构的静力学线弹性数值分析是工程设计的基础技术工作。在实际工程结构分析中，经常涉及自由结构，即约束不足或完全无约束的情况，如自由飞状态下的飞机。对自由飞状态下的飞机，其受到的气动载荷与质量惯性载荷形成配平载荷，即两种载荷互相平衡，此时飞机以匀速运动方式飞行；实际上，飞机升力面上的气动载荷与飞行姿态以及弹性变形等相关。飞机的地面静力试验或飞机在平衡载荷作用下的数值计算均需在飞机上施加静定约束，虽该静定约束不产生支反力，但却会引入飞机整体的刚性位移，其静态变形并不能真实反映自由飞状态下的飞机弹性位移。因此，消除因约束引入的整体刚性位移，对于获其真实弹性位移具有重要的工程意义。
[0003] 一般地，基于位移法的有限元数值计算必须建立完整适定的约束，或利用刚性位移模态消除结构总刚度矩阵的奇异性。处于无约束自由飞状态下的飞机，受到配平载荷作用，因缺少求解所需的静定约束，现有商用软件一般采用矩阵缩聚或模态技术求解其静变形。文献1“MSC/NASTRAN Linear Static Analysis:User’s Guide.Version 68,The MacNeal-Schwendler Corporation.1994”利用刚性惯性释放技术获得平衡外部气动载荷的质量惯性力，使结构处于一种无约束静力平衡状态；其次引入总位移自由度的质量加权平均作为自由结构平衡方程的额外约束进行求解。但该方法得到的位移是弹性与刚性位移的一种混合解，并且破坏了自由结构总刚矩阵的稀疏性。文献2“ZAERO User’s Manual.Version 8.2,ZONA Technology Inc.March 2008”利用有限阶弹性模态的线性组合模拟结构弹性变形，并将气动载荷、惯性载荷及弹性内力构成的平衡方程变换到模态坐标系下进行求解，从而获得自由结构的弹性变形解。该方法利用有限模态子空间代替整个构型空间，因而存在模型的截断误差。文献3“A localized FETI method for structural parallel analysis[J].Structure&Environment Engineering,2007.8.Vol.34,No.4”提出由无约束有限元模型建立刚性模态阵，再用投影矩阵构造方法消除漂浮子域的奇异性，从而获得无约束结构的位移解。该方法由于投影矩阵是对称满阵，破坏了原结构总刚矩阵的稀疏性，给大规模求解带来不便。
发明内容
[0004] 本发明针对自由飞状态配平载荷作用下飞机静变形弹性位移计算的上述不足，提出了一种消除结构整体刚性位移的精确高效方法。基本特征在于，给配平载荷作用的结构有限元模型施加任意合理的静定约束，利用工程商用软件计算得到结构位移数据，借助一致紧支径向基插值函数在选定的区域上重构其位移场，进而计算其变形梯度矩阵，对变形梯度矩阵实施极分解运算得到对称矩阵和正交矩阵的乘积，对称矩阵描述了结构的伸长，正交矩阵反映了结构的旋转。根据极分解获得的刚性旋转矩阵，构造结构整体的刚性位移，继而从总位移中消除结构的整体刚性位移，从而得到自由结构的真实弹性位移。
[0005] 计算流程如下：
[0006] S1：根据实际结构的性质在有限元分析软件中建立结构有限元模型；
[0007] S2：对步骤S1中建立的有限元模型施加自平衡载荷及静定约束，然后利用NASTRAN软件求解线性静变形；
[0008] S3：根据求解得到的结点位移，利用一致紧支径向基插值函数在选定的区域上重构其位移场，并由此求解位移梯度；
[0009] S4：对位移梯度进行极分解获得刚性旋转矩阵；
[0010] S5：根据结构刚体旋转性质，利用刚性旋转矩阵构造结构的整体刚性位移；
[0011] S6：从总位移中消除结构的整体刚性位移获得其真实弹性变形。
[0012] 本发明技术对于自由飞状态下飞机配平载荷作用的真实弹性位移计算具有重要工程应用价值及意义。具体方法包括以下步骤：
[0013] 步骤1：原始数据及理论计算关系。针对飞机结构，利用商用有限元分析软件建立结构静力学分析的有限元模型，并设置恰当的静定约束，在自平衡载荷作用下，该静定约束不产生支反力。调用结构线性静态分析求解器，即可获得结构在静定约束下的静变形。静定约束下结构的静变形包含了整体刚性位移，因而从总位移中消除结构刚性位移便可得到自由结构真实的弹性变形：
[0014] ud(X)＝u(X)-ur(X) (1)
[0015] 式中，ud为消除刚性位移后自由结构的弹性变形；u为静定约束下结构总位移；ur为结构整体的刚性位移；X表示结构中任意一点的总体坐标，坐标原点位于结构质心处，坐标轴方向符合工程应用习惯。
[0016] 由于静定约束必须只能施加在飞机壳体结构上的点，这仅限制了结构的整体刚性平动位移，由于约束点不在质心处，则必然引起结构体的刚性转动。根据结构体刚性转动定义，式(1)中结构整体刚性位移ur可以写作：
[0017] ur(X)＝R(X-Xc) (2)
[0018] 式中，R为刚性旋转的正交阵，即RTR＝I；Xc为约束点坐标。该刚性旋转矩阵可通过对反映结构整体刚性位移的变形梯度矩阵进行极分解得到。
[0019] 对于任意一个非奇异二阶张量，即该二阶张量行列式值不为零，可以通过极分解得到一个二阶对称张量和一个二阶正交张量。变形梯度表示了当前构型与初始构型的关系，对任意一个由初始构型变形得到的当前构型，按照逆变换可以得到初始构型，因而变形梯度可逆，也就是说其行列式值不为零。因此，可以对变形梯度进行极分解，其形式可以写为：
[0020]
[0021] 式中，F为结构的变形梯度，x为结构体变形后的任一点坐标；R为正交转动张量，反映了结构的刚性转动；U为对称张量，反映了结构的变形伸长。R和U可以通过下式获得：
[0022]
[0023] 对于二维平面情况，若结构整体绕着垂直纸面向外的z轴旋转了θ角，则其旋转矩阵可以写为
[0024]
[0025] 步骤2：结构中任意一点位移。求解结构变形梯度需要知道结构中任意一点的位移。利用商用有限元分析软件求解得到的是离散后各个单元节点的位移，为了获得结构中任意一点的位移，需要利用节点位移插值得到。在有限元法中，任意一点位移可以通过该点所在单元节点插值得到。
[0026] 对四节点平面等参单元，假设母单元为边长1×1的正方形，母单元坐标原点为其中心点，实际结构单元内任意一点位移通过下式得到：
[0027]
[0028] 式中，ui,vi(i＝1,…,4)为单元节点位移；Ni(i＝1,…,4)为等参元形函数，可以表示为：
[0029]
[0030] 式中，ξ,η为母单元中任意一点坐标，-1≤ξ≤1，-1≤η≤1，ξi,ηi为节点i的坐标，实际结构单元中任意一点与母单元中一点坐标相对应，其具体关系可以表示为：
[0031]
[0032] 通过以上单个单元内插值可得到单元内任一点的位移，但需要知道的单元类型，对不同类型的单元，插值函数并不一致，并且通常商用有限元分析软件中单元类型所对应的插值形函数一般无法得到具体形式，因而式(6)无法显式表达。同时，对单个单元插值并不能反映整体的刚性位移，从而无法正确消除刚性位移。
[0033] 步骤3：一致紧支径向基函数插值。为了克服使用单个单元进行插值的缺点，本发明采用一致紧支径向基函数插值来计算任意一点的位移。一致紧支径向基插值函数的形式如下：
[0034]
[0035] 式中，s(X)是结构节点变形的插值函数；N为选定的插值节点个数；αi为插值系数；为所采用的基函数形式；║·║为2-范数。p(X)为X的线性多项式，其表达式为：
[0036] p(X)＝γ0+γ1X+γ2Y+γ3Z (10)
[0037] 式中，X点处的坐标为(X,Y,Z)；γi(i＝0,…,3)为插值系数。
[0038] 式(9)及式(10)中插值系数需满足如下条件：
[0039]
[0040] 式中，qj(X)为次数不超过p(X)的多项式，可以取为1、X、Y、Z。
[0041] 通过选择插值基函数，根据结构节点坐标和已知的位移分量，即可求解式(9)中的未知系数。将已知结构节点坐标及位移分量代入式(9)，由式(10)和式(11)可以得到：
[0042] us＝Css·αs (12)
[0043] 式中，us为节点单方向位移分量构成的列阵，αs为待求系数列阵，具体形式可以写为：
[0044]
[0045] 式中，si(i＝1,…,N)为各节点单向位移分量，工程上用ui，vi，wi表示节点i沿总体坐标系的各位移分量。
[0046] 式(12)中矩阵Css具体形式如下：
[0047]
[0048] 式中，Css为可逆矩阵， 为以节点间距离为自变量的一致紧支径向基函数。
[0049] 将Css代入式(12)即可确定结构节点和位移形成的插值曲面。将任意一点坐标代入式(9)，得到该点的位移，即：
[0050]
[0051] 式中，s为该点位移分量u、v或w；us列阵根据所求分量确定；Aas的具体形式如下：
[0052]
[0053] 步骤4：构造变形梯度矩阵。式(3)中结构任意一点的变形梯度具体形式可以写为：
[0054]
[0055] 式中，x为结构中任意一点X变形后的坐标，x的分量为x、y、z；X的分量为X、Y、Z。
[0056] 式(17)中的变形梯度表达式不易显式表达，因此写为关于位移的表达。若u为任意一点的位移，其分量为u、v、w，由于当前构型由初始构型变形所得，故而有：
[0057] x＝X+u (18)
[0058] 将式(18)代入式(17)，可以得到变形梯度关于位移的表达式：
[0059]
[0060] 式中，I为3×3的单位矩阵。
[0061] 由步骤3可以得到任意一点的位移，根据式(9)可知，其位移导数可以通过下式求得：
[0062]
[0063] 式中，DXas、DYas和DZas为式(16)对坐标X、Y、Z求导所得，即：
[0064]
[0065] 式中， 和 为对X、Y和Z的偏导数。
[0066] 将式(20)代入式(19)中，便可以得到结构的变形梯度具体表达式，其反应了结构整体的弹性与刚性变形。
[0067] 步骤5：结构整体刚性位移的消除。根据步骤4中构造的变形梯度按式(4)进行极分解得到反映结构整体刚性运动的旋转矩阵R。根据式(2)计算各个节点的刚性位移：
[0068] uri＝R(Xi-Xc)i＝1,…,N (22)
[0069] 式中，uri为第i个节点的刚性位移；Xi为第i各节点的坐标；N为节点个数。
[0070] 将步骤1中静定约束下结构各个节点的总位移减去其整体刚性位移，便可以得到消除整体刚性位移后的静弹性变形位移，即自由结构受自平衡载荷的静变形：
[0071] udi＝ui-uri i＝1,…,N (23)
[0072] 式中，udi为结构消除整体刚性位移后各节点的位移，即无约束状态下结构的位移；ui为结构在静定约束下由有限元分析软件计算得到的总位移；uri为由式(22)计算得到的各个节点的刚性位移。
附图说明
[0073] 图1是本发明一种消除自由飞配平载荷下飞机变形中刚性位移方法的流程图；
[0074] 图2是一平面梁几何模型及有限元网格图；
[0075] 图3是静定约束下平面梁静变形图；
[0076] 图4是平面梁刚性位移图；
[0077] 图5是无约束平面梁自平衡载荷静变形图；
[0078] 图6是一无人机机翼几何及骨架结构图，其中(a)为机翼平面几何形状图，(b)为机翼结构内部骨架结构图；
[0079] 图7是机翼结构有限元网格及离散质量模型图，其中(a)和(b)为机翼结构有限元网格图，(c)为机翼结构离散质量分布图；
[0080] 图8是静定约束下机翼静变形图；
[0081] 图9是机翼刚性位移图；
[0082] 图10是无约束机翼自平衡载荷静变形图。
具体实施方式
[0083] 关于本发明的特点及流程，下面将针对具体结构，结合附图和实施例对本发明进一步说明。
[0084] 实施例1
[0085] 本实施例是对一受自平衡载荷作用的二维平面梁静变形消除刚性位移的过程。具体过程包括以下步骤：
[0086] 步骤S11，在PATRAN中建立二维平面梁模型，在长度方向划分10个网格，高度方向划分1个网格，其模型及有限元网格如图2所示，平面梁模型几何尺寸及材料参数如表1所示。
[0087] 表1平面梁几何及材料参数
[0088]
[0089] 步骤S12，在节点3和节点9上施加竖直向上载荷，其大小为1000N，节点17上施加竖直向下的载荷，其大小为2000N，这些施加的载荷对平面梁形成自平衡载荷系。在平面梁左端施加静定约束，同时约束节点1的u、v方向及节点12的u方向位移。调用NASTRAN软件对步骤1中建立的平面梁有限元模型进行线性静力求解，其静变形如图3所示，表2列出了各个节点的位移，在该自平衡载荷系作用下平面梁自由端最大挠度达到了梁长度的27％。
[0090] 表2平面梁各个节点位移
[0091]
[0092] 步骤S13，选取C2函数作为紧支径向基函数，根据步骤1中得到梁模型各个节点坐标，代入式(14)和式(21)得到矩阵Css、DXas和DYas；利用NASTRAN求得的各个节点位移，代入式(13)，分别形成关于u和v的位移列阵。利用式(19)和式(20)计算得到平面梁质心的变形梯度矩阵F：
[0093]
[0094] 步骤S14，根据式(4)求解计算旋转矩阵R：
[0095]
[0096] 其代表了平面梁的刚性转动。对二维平面梁模型，根据式(5)可知，平面梁绕垂直纸面向外的z轴旋转角为16.1°。
[0097] 得到刚性旋转矩阵R后，根据式(22)可以得到各个节点的刚性位移，刚性位移图如图4所示。
[0098] 步骤S15，由步骤S12所得的各个节点总位移与步骤S14所得的刚性位移，按照式(23)消除由于静定约束导致的刚性位移，得到无约束平面梁在自平衡载荷作用下静变形，其变形图如图5所示。
[0099] 从步骤S13到步骤S15，可完成结构变形梯度矩阵以及旋转矩阵计算，并以此从总体位移中消除刚性位移。对比图3和图5可以看出，本发明算法通过消除静定约束引起的刚性位移，获得了自由无约束平面梁在自平衡载荷下的真实静变形；由于载荷为对称载荷，其静变形应为对称变形，图5显示消除刚性位移后静变形为对称变形；无约束平面梁静变形与未消除刚性位移的静变形有较大差别，无约束平面梁消除刚性位移后挠度会减小。图4显示了由于引入静定约束而导致结构的刚性位移，该刚性位移过程中平面梁保持了面积不变，符合结构刚性位移要求。
[0100] 表3不同挠度消除刚性位移后静变形对称误差
[0101]
[0102] 表3给出了一端约束平面梁在不同自平衡载荷大小下两端对称性误差，结果显示平面梁在该自平衡载荷下真实静变形具有很高的对称性；随着载荷逐渐增加，梁自由端挠度逐渐增加，直至挠度达到长度的27％，虽然对称误差有所增加，但总体符合精度要求，显示了本发明算法的精确性。
[0103] 实施例2
[0104] 上述实施例1是对简单平面梁的分析，为验证本发明算法对复杂结构的适用性，下面对一薄壁结构的无人机机翼进行操作说明。具体过程包括以下步骤：
[0105] 步骤S21，利用CATIA及Hypermesh软件对无人机机翼的CAD模型进行建模并进行网格划分。无人机机翼几何外形及内部骨架结构如图6所示，可划分为中央翼、左右外翼，分割点为A3L及A3R，图中标注了关键点的位置。该模型关于中面A0-B0对称，表4给出了右半机翼尺寸参数，其半展长为8900mm。
[0106] 表4机翼平面尺寸参数
[0107]
[0108] 无人机机翼有限元网格如图7所示。本算例采用了较小尺度的板壳和梁两类单元的有限元网格剖分方法。结构质量采用集中质量方法建模，即将每个结构单元的分布质量均分到有限元模型的节点上。根据表5所示的结构材料及质量密度输入有限元模型参数。
[0109] 表5模型材料性能参数
[0110]
[0111] 步骤S22，在机翼上施加竖直向上的均布载荷，其大小为300N，利用Nastran软件惯性释放技术获得平衡该均布载荷的质量惯性力，结构此时处于自平衡状态，在节点712上施加静定约束后进行线性静力求解，其静变形如图8所示，翼稍一端最大挠度为4930mm，达到了展长的27.7％。
[0112] 步骤S23，选取C2函数作为紧支径向基函数，根据步骤1中得到无人机机翼模型各个节点坐标，代入式(14)和式(21)得到矩阵Css、DXas和DYas；利用步骤S22中NASTRAN求得的各个节点位移，代入式(13)，分别形成关于u,v和w的位移列阵。利用式(19)和式(20)计算得到机翼质心的变形梯度矩阵F：
[0113]
[0114] 步骤S24，根据式(4)求解计算旋转矩阵R：
[0115]
[0116] 其代表了无人机机翼的刚性转动。
[0117] 得到刚性旋转矩阵R后，根据式(22)可以得到各个节点的刚性位移，刚性位移图如图9所示。
[0118] 步骤S25，由步骤S22所得的各个节点总位移与步骤S24所得的刚性位移，按照式(23)消除由于静定约束导致的刚性位移，得到无约束机翼在自平衡载荷作用下的静变形，其变形图如图10所示，得到的机翼静变形关于中面对称。
[0119] 从图中可以得到三点结论：
[0120] 1.由于静定约束的设置，对称的自平衡载荷作用下结构静变形不再对称，该静变形与真实弹性变形具有较大区别，因而已不能真实地反映自由结构的弹性位移。
[0121] 2.将自平衡载荷下结构静变形中的刚性位移消除后可得到自由结构在该配平载荷下的真实弹性位移，即使自由端具有较大的挠度，依然可以保证消除真实的刚性变形，误差符合要求，显示了本算法合理性、精确性。
[0122] 3.对比现有商用有限元分析软件可知，本发明算法过程简单，易于操作，具有高效的效率。
[0123] 由此可见，本发明算法整体高效便捷，具有较高的数值精度，且算法易于在现代常规个人电脑上执行。可为受配平载荷作用的现代飞机在自由飞状态下的真实静变形提供精度可靠的有效工程计算方法。
