一种气动载荷转换为结构载荷的方法与流程

[0001]
本发明涉及计算机仿真技术领域，尤其涉及一种气动载荷转换为结构载荷的方法。


背景技术：

[0002]
飞行器在实际使用过程中承受着各种各样的载荷，大体可以分为两大类：一类为空气动力载荷，另一类为惯性力载荷。高速运行飞行器受到气动载荷的影响比较大，在飞机结构有限元分析中，空气动力载荷最终在结构设计中往往需要转化为结构上的受力，也即通常需要将气动载荷转换为结构有限元的节点载荷，基于该种需求，提供一种将气动载荷转换为结构节点载荷的方法具有非常重要的工程意义。


技术实现要素：

[0003]
针对上述现有技术的缺点，本发明的目的旨在提供一种高效、快速、准确的气动载荷转换为结构载荷的方法。
[0004]
本发明实施例提供的一种气动载荷转换为结构载荷的方法，该方法包括：
[0005]
气动面载荷预处理，形成气动面cp数据；结构数据预处理，形成结构节点空间位置数据；
[0006]
根据形成的气动面cp数据和结构节点空间位置数据，插值计算得到结构节点cp数据；
[0007]
获得结构节点cp数据后，在单元面上积分形成单元载荷数据，对单元载荷数据在结构节点上进行分载计算，得到结构节点力载荷。
[0008]
优选地，上述气动载荷转换为结构载荷的方法中，所述气动面类型包括扁平面和曲面；
[0009]
所述气动面为扁平面时，气动面载荷预处理流程包括：
[0010]
读取输入的气动面轮廓数据，形成分载网格面；
[0011]
读取气动分布载荷cp，根据双线性插值方法，形成分载网格面cp数据；
[0012]
所述气动面为曲面时，气动面载荷预处理流程包括：
[0013]
读取预存的气动面节点cp数据。
[0014]
优选地，上述气动载荷转换为结构载荷的方法中，所述根据双线性插值方法，形成分载网格面cp数据的流程如下：
[0015]
假设分载网格面的四个角点的坐标和cp数据分别为：(x1，y1，cp1)、(x2，y2，cp2)、(x3，y3，cp3)和(x4，y4，cp4)；
[0016]
由四个角点插值生成双线性曲面，则该面内任一点的参数坐标和cp数据为：
[0017]
x(u，w)＝x1+u(x2-x1)+[x3+u(x4-x3)-x1-u(x2-x1)]w
[0018]
y(u，w)＝y1+u(y2-y1)+[y3+u(y4-y3)-y1-u(y2-y1)]w
[0019]
cp(u，w)＝cp1(1-u)(1-w)+cp3(1-u)w+cp2u(1-w)+cp4uw
[0020]
式中：u∈[0，1]；w∈[0，1]，u、w均为参数坐标。
[0021]
优选地，上述气动载荷转换为结构载荷的方法中，所述插值计算得到结构节点cp数据的具体流程包括：
[0022]
4.1)根据结构节点的空间位置数据，确定在设定范围内与所有结构节点关联的气动cp节点；
[0023]
4.2)根据步骤4.1)确定的气动cp节点的空间位置数据和cp数据，通过样条基函数法或径向基函数法，计算得到局部气动cp数据分布场的插值系数；
[0024]
4.3)根据结构节点的空间位置数据和步骤4.2)中得到的局部气动cp数据分布场的插值系数，计算得到结构节点cp数据。
[0025]
优选地，上述气动载荷转换为结构载荷的方法中，所述对单元载荷数据在结构节点上进行分载计算的具体流程包括：
[0026]
5.1)获取结构单元的节点编号及其空间坐标；
[0027]
5.2)根据结构单元的类型，计算结构单元法向量；
[0028]
5.3)设置结构单元积分点数目；
[0029]
5.4)根据积分点参数坐标和结构单元节点的物理坐标，计算得到结构单元高斯点的物理坐标；
[0030]
5.5)根据高斯点物理坐标和插值系数计算得到高斯点的cp数据；
[0031]
5.6)根据结构单元法向量计算任一空间平面旋转至xy平面的旋转矩阵，并将空间单元旋转到xy平面；
[0032]
5.7)通过高斯积分计算得到空间单元旋转后的压心和压力；
[0033]
5.8)在旋转后的单元上进行分载；
[0034]
5.9)将分载得到的基于xy平面单元的节点力旋转回原空间；
[0035]
5.10)对属于不同结构单元的同一节点的节点载荷进行叠加，得到结构节点力载荷；
[0036]
5.11)判断所有结构单元是否均完成处理，如没完成处理，则继续重复步骤5.1)-5.10)。
[0037]
优选地，上述气动载荷转换为结构载荷的方法中，所述结构单元的类型包括三角形和四边形，所述结构单元的类型为三角形时，通过平面单元两个相连的单元边向量叉乘获得垂直于平面的结构单元法向量；所述结构单元的类型为四边形时，通过对角线向量叉乘获得结构单元法向量。
[0038]
优选地，上述气动载荷转换为结构载荷的方法中，所述在旋转后的单元上进行分载计算具体采用最小优化法，或非线性规划法，或多点排直接分配法进行计算。
[0039]
优选地，上述气动载荷转换为结构载荷的方法中，所述通过高斯积分计算得到空间单元旋转后的压心和压力的公式如下：
[0040][0041]
其中，f为计算得到的压力，x
f
，y
f
，z
f
是压力作用点，即压心的坐标。
[0042]
优选地，上述气动载荷转换为结构载荷的方法中，还包括在得到结构节点力载荷后，判断载荷转换前后总压和压心是否一致，如一致，则输出节点载荷文件，如不一致，则对节点力载荷进行修正。
[0043]
优选地，上述气动载荷转换为结构载荷的方法中，所述对节点力载荷进行修正的具体流程包括：
[0044]
10.1)读取气动面提供的总压和压心数据；
[0045]
10.2)基于步骤10.1)的总压和压心数据，计算各个方向上的合力矩信息；
[0046]
10.3)针对单个方向建立线性方程组，并求解系数；
[0047]
10.4)根据步骤10.3)获得的系数，结合结构节点坐标计算得到修正后的单方向分量载荷；
[0048]
10.5)重复步骤10.1)-10.4)，直到x、y、z三个方向分量载荷均处理完毕；
[0049]
10.6)输出修正后的节点力载荷。
[0050]
与现有技术相比，本发明气动载荷转换为结构载荷的方法包括：气动面载荷预处理，形成气动面cp数据；结构数据预处理，形成结构节点空间位置数据；根据形成的气动面cp数据和结构节点空间位置数据，插值计算得到结构节点cp数据；获得结构节点cp数据后，在单元面上积分形成单元载荷数据，对单元载荷数据在结构节点上进行分载计算，得到结构节点力载荷本发明通过读取气动载荷数据，继而形成压力场，通过读取结构有限元数据，并在形成的压力场内进行插值，形成单元面压力载荷，并最终转化为节点力载荷。本发明通过读取气动载荷数据，继而形成压力场，通过读取结构有限元数据，并在形成的压力场内进行插值，形成单元面压力载荷，并最终转化为节点力载荷。本发明能够快速有效地转化气动载荷为结构载荷，并保证转换前后的总载荷大小和压心相同。
附图说明
[0051]
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简要介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0052]
图1为本发明提供的一种气动载荷转换为结构载荷的方法流程图；
[0053]
图2为本发明提供的插值计算得到结构节点cp数据的方法流程图；
[0054]
图3为本发明提供的对单元载荷数据在结构节点上进行分载计算的流程图；
[0055]
图4为本发明提供的物理坐标系下的四边形单元结构示意图；
[0056]
图5为图4等参坐标系下的四边形单元结构示意图；
[0057]
图6为本发明提供的对节点力载荷进行修正的方法流程图。
具体实施方式
[0058]
为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步地详细描述，显然，所描述的实施例仅仅是本发明一部份实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
[0059]
下面结合说明书附图对本发明实施例作进一步详细描述。
[0060]
如图1所示，本发明实施例公开了一种气动载荷转换为结构载荷的方法，该方法包括：
[0061]
气动面载荷预处理，形成气动面cp数据；结构数据预处理，形成结构节点空间位置数据；
[0062]
根据形成的气动面cp数据和结构节点空间位置数据，插值计算得到结构节点cp数据；
[0063]
获得结构节点cp数据后，在单元面上积分形成单元载荷数据，对单元载荷数据在结构节点上进行分载计算，得到结构节点力载荷。
[0064]
优选地，本发明实施例结构数据的预处理可通过节点文件和单元文件读入结构数据。
[0065]
本发明采用python语言进行开发，提供方便的文本输入和文本内容显示窗口，通过读取气动载荷数据，继而形成压力场，通过读取结构有限元数据，并在形成的压力场内进行插值，形成单元面压力载荷，并最终转化为节点力载荷。本发明能够快速有效地将气动载荷转化为结构载荷，并保证转换前后的总载荷大小和压心相同。
[0066]
优选地，本发明实施例气动面类型包括扁平面和曲面。本发明实施中，气动部门会根据气动面的几何特点给出气动载荷，若气动面接近平面，即气动面为扁平面时，会给出气动面轮廓数据和分布载荷cp，通过对气动数据的处理，会得到加载面上的cp值；若气动面是大曲率结构，即气动面为曲面时，会给出详细的气动面网格，可以直接读取预存的数据。
[0067]
具体地，所述气动面为扁平面时，气动面载荷预处理流程包括：读取输入的气动面轮廓数据，形成分载网格面；读取气动分布载荷cp，根据双线性插值方法，形成分载网格面cp数据。所述气动面为曲面时，气动面载荷预处理流程包括：读取预存的气动面节点cp数据。
[0068]
本发明分载网格面cp数据是一种特殊的气动面cp数据，其数据量相对气动面cp数据少，一般的气动面cp数据由cfd计算得出，数据点较多，扁平面cp数据一般由工程方法得出，数据点较少。
[0069]
优选地，本发明所述根据双线性插值方法，形成分载网格面cp数据的流程如下：
[0070]
假设分载网格面的四个角点的坐标和cp数据分别为：(x1，y1，cp1)、(x2，y2，cp2)、(x3，y3，cp3)和(x4，y4，cp4)；
[0071]
由四个角点插值生成双线性曲面，则该面内任一点的参数坐标和cp数据为：
[0072]
x(u，w)＝x1+u(x2-x1)+[x3+u(x4-x3)-x1-u(x2-x1)]w
[0073]
y(u，w)＝y1+u(y2-y1)+[y3+u(y4-y3)-y1-u(y2-y1)]w
[0074]
cp(u，w)＝cp1(1-u)(1-w)+cp3(1-u)w+cp2u(1-w)+cp4uw
[0075]
式中：u∈[0，1]；w∈[0，1]，u、w均为参数坐标。
[0076]
优选地，如图2所示，本发明实施例在预处理完成后，就可以开始对结构节点的cp数据进行插值计算，具体地，本发明所述插值计算得到结构节点cp数据的具体流程包括：
[0077]
4.1)根据结构节点的空间位置数据，确定在设定范围内与所有结构节点关联的气动cp节点；
[0078]
4.2)根据步骤4.1)确定的气动cp节点的空间位置数据和cp数据，通过样条基函数法或径向基函数法，计算得到局部气动cp数据分布场的插值系数；
[0079]
4.3)根据结构节点的空间位置数据和步骤4.2)中得到的局部气动cp数据分布场的插值系数，计算得到结构节点cp数据。
[0080]
本发明需要根据结构节点的坐标位置，找到在一定范围内与之关联的气动cp节点。具体实施中，重复步骤4.1)直到所有结构节点均找到与之关联的气动cp节点。
[0081]
本发明实施例通过径向奇函数法计算得到局部气动cp数据分布场的插值系数的方式如下：
[0082]
在定义域x∈r
d
上，所有形如φ(x-c)＝φ(||x-c||)及其线性组合而成的函数空间称为由函数φ导出的径向基空间。当x
j
两两不同时，φ(x-x
j
)线性无关，故可选其作为某函数子空间的一组基，当x
j
几乎充满r时，由φ(x-x
j
)生成的函数空间几乎可以逼近所有函数。常用的基函数有：gauss函数markoff函数multi-quadric函数φ(r)＝(c2+r2)-β
等等。
[0083]
本发明实施例采用multi-quadric函数插值，即为对于给定数据点集{(x
j
,f
j
)}
n
，寻找函数使其满足插值条件k＝1,l,n。针对multi-quadric函数插值的唯一性定理：对于mq基函数插值问题解存在且唯一，同时增加一个多项式基函数，因此取插值函数为：
[0084][0085]
考虑到求解的复杂性，取p
j
(x)为常数项，同时已知参数c取0，将线性方程组简化，因此插值函数为：
[0086][0087]
根据无穷远边界条件进行约束：以及已知条件：s(x
j
)＝f
j
，即可确定权重a
j
和b。
[0088]
本发明实施例通过样条基函数法计算得到局部气动cp数据分布场的插值系数的方式如下：
[0089]
在n维欧式空间内域d上有n个向量(每个向量代表一个节点坐标)及其对应的函数
值，即已知x
i
＝[x1,x2,l,x
n
]
t
,(i＝1,2,l,n)。设向量x的函数w＝w(x)(称为真实函数)是单值连续的，根据已知条件构成一个逼近函数，公式如下：
[0090][0091]
式中，c1,c2,l,c
n+1+n
为待定系数，n是x的维度；ε为给定常量，对于一般平坦曲面ε＝10-2
:1，对于有奇性的曲面可取ε＝10-5
:10-6
；x
pi
为第i个已知节点的第p维坐标；x
pk
为待求节点x
k
的第p维坐标。
[0092]
根据无穷远处边界条件可得：根据n个已知节点函数值可确定待定系数。
[0093]
优选地，如图3所示，本发明实施例在获得结构节点的cp数据后，就可以开始对结构节点进行分载计算，具体地，本发明所述对单元载荷数据在结构节点上进行分载计算的具体流程包括：
[0094]
5.1)获取结构单元的节点编号及其空间坐标；
[0095]
5.2)根据结构单元的类型，计算结构单元法向量；
[0096]
5.3)设置结构单元积分点数目；
[0097]
5.4)根据积分点参数坐标和结构单元节点的物理坐标，计算得到结构单元高斯点的物理坐标；
[0098]
5.5)根据高斯点物理坐标和插值系数计算得到高斯点的cp数据；
[0099]
5.6)根据结构单元法向量计算任一空间平面旋转至xy平面的旋转矩阵，并将空间单元旋转到xy平面；
[0100]
5.7)通过高斯积分计算得到空间单元旋转后的压心和压力；
[0101]
5.8)在旋转后的单元上进行分载；
[0102]
5.9)将分载得到的基于xy平面单元的节点力旋转回原空间；
[0103]
5.10)对属于不同结构单元的同一节点的节点载荷进行叠加，得到结构节点力载荷；
[0104]
5.11)判断所有结构单元是否均完成处理，如没完成处理，则继续重复步骤5.1)-5.10)。
[0105]
本发明所述结构单元的类型包括三角形和四边形，所述结构单元的类型为三角形时，通过平面单元两个相连的单元边向量叉乘获得垂直于平面的结构单元法向量；所述结构单元的类型为四边形时，通过对角线向量叉乘获得结构单元法向量。
[0106]
如图4所示为本发明实施例物理坐标系(x,y)下的四边形单元，图5为转换到等参坐标系(ξ,η)下的四边形单元。可见，在等参坐标下，单元形状非常规则，有利于数值积分运算，因此在有限元计算中，常常采用等参元的高斯积分方法对单元进行积分计算。即将一个函数的定积分转化为若干个点的函数值的加权组合。
[0107]
高斯点根据分布函数多项式的次数选择高斯点的个数，本发明实施例在程序中内置了最多3x 3个高斯点的坐标和权系数。
[0108]
假定在等参单元上具有分布函数f(ξ,η),对单元进行积分计算即需要计算如下所示定积分：
[0109][0110]
通过高斯积分，可以将上式写为：
[0111][0112]
其中，h
i
,h
j
是对应的高斯点权系数，可以通过查表获得，f(ξ
i
,η
i
)为积分点(ξ
i
,η
i
)的函数值，(ξ
i
,η
i
)为高斯点的参数坐标，也可以通过查表获得。
[0113]
优选地，本发明所述在旋转后的单元上进行分载计算可以选择使用最小优化法，或非线性规划法，或多点排直接分配法等三种方法进行计算。
[0114]
本发明实施例采用多点排直接分配法进行计算的方式如下：假设有限元节点与气动点之间有一个无形的梁，它是以气动点一端为固支的悬臂梁，其自由端上的有限元节点分配到载荷p
j
时产生的变形能为：
[0115][0116]
于是全系统的变形能为：
[0117][0118]
载荷的转换要对系统影响最小，即在满足静力等效的条件下，系统的变形能最小，即：
[0119][0120]
其中，n是单元中的节点数。
[0121]
采用拉格朗日(lagrange)乘子法建立拉格朗日函数：
[0122][0123]
式中，λ、λ
x
、λ
y
为拉格朗日乘子。为使得f(λλ
x
λ
y
)最小，有：
[0124]
[0125]
因此有：
[0126][0127]
将上式代入问题方程组，可以得到如下线性代数方程组：
[0128][0129]
对上式求解，即可得到λ、λ
x
、λ
y
，代入即可求得分载后的载荷p
j
。
[0130]
本发明实施例采用该方法满足了静力等效原则和传力路线不变，静力等效原则保证了总载荷和总压心不变，而传力路线不变主要体现在保证载荷真实传递。
[0131]
本发明实施例采用最小优化法和非线性规划法进行计算的方式如下：
[0132]
需要求解的问题可以描述为：对于一个处于平面上的单元，找到一组作用在节点上的力载荷满足与作用在压心上的合力和合力矩等效，并使得和或者最小。
[0133][0134]
具体实施中，对于上述优化问题，可以直接使用最小优化方法和非线性规划两种方法进行求解。
[0135]
实施中，本发明通过高斯积分计算得到空间单元旋转后的压心和压力的公式如下：
[0136][0137]
其中，f为计算得到的压力，x
f
，y
f
，z
f
是压力作用点，即压心的坐标。
[0138]
优选地，本发明气动载荷转换为结构载荷的方法中，还包括在得到结构节点力载荷后，判断载荷转换前后总压和压心是否一致，如一致，则输出节点载荷文件，如不一致，则对节点力载荷进行修正。
[0139]
本发明实施中，在获得结构节点的节点力载荷后，需要依据总压和压心与未转换前的气动载荷进行比较，如果两者一致，则不用修正，如果差距较大则需要对节点力载荷进行修正。
[0140]
如图1、6所示，本发明对节点力载荷进行修正的具体流程包括：
[0141]
10.1)读取气动面提供的总压和压心数据；
[0142]
10.2)基于步骤10.1)的总压和压心数据，计算各个方向上的合力矩信息；
[0143]
10.3)针对单个方向建立线性方程组，并求解系数；
[0144]
10.4)根据步骤10.3)获得的系数，结合结构节点坐标计算得到修正后的单方向分量载荷；
[0145]
10.5)重复步骤10.1)-10.4)，直到x、y、z三个方向分量载荷均处理完毕；
[0146]
10.6)输出修正后的节点力载荷。
[0147]
本发明实施中在求解系数时，假设节点载荷在该方向的分量的修正量与节点其它两个方向的坐标满足线性关系，从而建立线性方程组，求解得到系数。
[0148]
具体地，本发明实施例根据求解系数，结合结构节点坐标计算得到修正后的分量载荷的修正方式如下：对于节点载荷的修正，假设每个节点载荷修正量与节点坐标为线性关系，通过待定系数法确定。
[0149]
假设已知总载为(f
x
,f
y
,f
z
)，压心为(x,y,z)，可以在x、y、z三个方向对节点载荷的三个分量进行修正。
[0150]
对于x方向，合力为f
x
，合力矩为m
xz
＝f
x
gz、m
xy
＝f
x
gy，因此处理后的节点载荷在x方向的分量p
ix
(i＝1,2,l,n)在x方向合力及合力矩需与其总合力及力矩平衡。即：
[0151][0152]
其中，y
i
、z
i
为节点yz坐标，假设节点载荷在x方向的分量的修正量与节点坐标满足线性关系，即：δp
ix
＝a0+a1y
i
+a2z
i
，则可以通过解线性方程组求出系数a0、a1、a2。
[0153]
同理，可按照上述方法求出节点载荷在y方向、z方向上分量的修正量。
[0154]
综上，本发明提供的气动载荷转换为结构载荷的方法，提高了气动载荷处理的效率和准确性。
[0155]
尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例做出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。
[0156]
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围
之内，则本发明也意图包含这些改动和变型在内。