0 引言

应变片是全机静力试验中一种重要的测量元件,是对全机承载能力最直接的观察手段。在全机静力试验过程中,需要对10000个左右的应变片进行监控,全机应变片分布如图1所示。在试验前,需要计算每个应变片监测点的理论应变值,该值就是应变监测值。试验时,通过应变片测量结果与监测值的对比,确保试验的顺利进行,验证结构强度分析方法和全机有限元模型的正确性。

应变监测值可以由自然网格分析结果结合工程算法计算得到,或者由精细有限元模型分析结果直接对应变片所在单元做结果变换得到。自然网格模型的网格尺寸由典型结构件站位决定,网格尺寸大,结构细节特征在模型简化过程中丢失,宏观上能够准确表达全机刚度,应变监测值需要经过工程算法进行推算,但局部细节处应变计算精度不足。相比自然网格模型,全机精细有限元模型^[1-2]正逐渐被广泛使用,其单元尺寸可达厘米级别,网格数量达到千万级别,结构细节特征更加丰富,应变监测值无需复杂计算便可以得到,且计算精度更高。

然而,全机精细有限元模型规模和计算结果数据量都十分庞大,手工获取应变监测值繁琐且耗时。本文基于全机精细有限元模型,提出了一种获取应变监测值的解决方案,降低了应变监测值数据的准备难度,大大提高了工作效率。

1 基本原理

1.1 应变坐标变换

在全机精细有限元模型求解分析时,采用MSC NASTRAN求解器,该求解器对应变符号作如下规定:正应变以伸长为正,剪应变以使直角变小为正,且剪应变为工程剪应变,如图2所示的ε_y和γ_yz均为正值。

对于二维应变状态,坐标系oxy下应变分量变换至坐标系ox′y′下,变换公式^[3]如下:

式中,θ表示从x轴转向x′轴的角度,以逆时针为正,如图3所示。

全机精细有限元模型主要采用壳单元建模,单元尺寸为10mm左右,应变片尺寸为4mm×8mm,网格尺寸略大于应变片尺寸。假设应变片对应单元的中心点处的应变为应变片的应变。若已知应变片对应单元以及该单元中心点在单元坐标系下的应变状态,单元坐标系为oxy,应变片测量方向沿着坐标系ox′y′的x′轴,由公式(1)可以计算出应变片测量的应变,即

更为精确的方法是根据单元节点应变状态插值得到应变片中心点位置的应变状态,再根据公式(2)计算应变片方向的应变,应变插值方法参见1.2节。

1.2 等参单元应变插值

在二维平面oxy内,对于一阶等参四边形单元,如图4所示,其几何形状的插值公式^[4]为:

式中,$N_1=\frac{1}{4}(1-\zeta )(1-\eta ),N_2=\frac{1}{4}(1+\zeta )(1-\eta ),N_3=\frac{1}{4}(1+\zeta )(1+\eta ),N_4=\frac{1}{4}(1-\zeta )(1+\eta )$。

对于单元所在平面内任意一点P(x₀, y₀),该点的应变状态可以由单元节点应变状态通过形函数Ni(i=1, 2, 3, 4)插值得到,即

式中,ε_I为某个应变状态分量,ε_iI为单元节点i对应的某个应变状态分量,N_i(ξ₀,η₀)为第i个形函数在P点对应的自然坐标值(ξ₀,η₀)下的值,该值代表ε_iI的插值权重。对于单元中心点,自然坐标值为(0.0, 0.0),节点插值权系数均为0.25,所以中心点的应变状态为四个节点应变状态的平均。

根据公式(3),利用牛顿法^[5]求解非线性方程组可以得到P点的自然坐标值(ξ₀,η₀)。将P点坐标带入公式(3)后,并改写为

式中,

$\begin{array}{c}\left\{\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\end{array}\right\}=\frac{1}{4}\left[\begin{array}{cccc}1& 1& 1& 1\\ -1& 1& 1& -1\\ -1& -1& 1& 1\\ 1& -1& 1& -1\end{array}\right]\left\{\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\end{array}\right\},\\ \left\{\begin{array}{c}{b}_0\\ b_1\\ b_2\\ b_3\end{array}\right\}=\frac{1}{4}\left[\begin{array}{cccc}1& 1& 1& 1\\ -1& 1& 1& -1\\ -1& -1& 1& 1\\ 1& -1& 1& -1\end{array}\right]{\left\{\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right\}}_0\end{array}$

则牛顿迭代公式为

在迭代的过程中,当公式(6)计算的残差在容差范围内时,便可得到P点的自然坐标值(ξ₀,η₀),即$∥X_{k+1}-X_k∥<\delta$,其中$\parallel \cdot \parallel$为L2范数,δ为容差。

2 应变监测值获取方法

2.1 获取流程

根据应变坐标变换和应变插值,本文提出了一套基于S模型的应变监测值获取流程,如图5所示。该流程分前处理和后处理两个部分,前处理主要完成:(1) 建立应变片几何模型,称该模型为S模型;(2) 确定S模型对应的最近单元;(3)计算所需的数据,相关数据形成S模型配置文件作为前处理的最终结果。后处理主要完成:(1) 获取最近单元的应变状态;(2) 利用S模型配置文件中相关数据进行应变坐标变换和应变插值。

相比后处理,前处理的工作量更大,但一旦生成一个工况下的S模型配置文件,便可以冻结该文件,在后续多个工况计算结果的后处理中重复使用该配置文件。

2.2 应变片S模型

在模型试验贴片工作之前,需要在三维建模软件中完成应变片几何模型的创建,该数模文件一般称之为S模型,如图6所示。

单片是由4个几何点组成的矩形,花片是由10个几何点组成的二维封闭多边形,如图7所示。对每个应变片定义应变片坐标系,且均为右手笛卡尔正交坐标系,确保应变片坐标系正z轴指向贴片者,应变片坐标系定义规则如下:

(1) 对于单片,应变片坐标系原点为应变片的中心点,正x轴为应变片长边方向,正y轴为应变片短边方向,正z轴为应变片法向;

(2) 对于花片,应变片坐标系原点为应变片的角点,正x轴为0°片测量方向,正y轴为90°片测量方向,坐标系第一象限与正x轴夹角为45°的为45°片测量方向。

2.3 最近单元搜索方法及相关物理量的计算

在应变监测值获取过程中,S模型配置文件中必须包含应变片最近单元编号、应变片在最近单元的单元坐标系下的夹角、应变插值权重系数、应变片距离最近单元参考平面的距离。

在前处理过程中,在应变片坐标系原点沿着±z轴方向寻找距离应变片最近的单元,将此单元作为应变片所粘贴的单元,应变片的测量结果由此单元的应变状态经坐标变换计算得到。

搜索距离应变片最近的单元时,需要将与应变片不平行的单元从搜索集中去除。对单元搜索集内的每个单元,计算应变片z轴与单元所在平面的交点,找到交点在单元内的单元。有限元模型中使用的二维单元包括三角形单元和四边形单元,几何上均为凸多边形。判断点是否在凸多边形内的方法有弧长法^[6]、面积法^[7]、角度法^[8]、射线法^[9]等,其中面积法适用范围广且算法简单,如图8所示。若交点与单元节点依次构成的三角形面积之和等于单元面积,则交点在单元内。

若应变片z轴和单元交点在单元内则计算交点至应变片坐标系原点的距离,按此规则找到单元搜索集内所有距离应变片原点最近的单元,该单元就是应变片所粘贴的单元。整个搜索过程中,考虑单元的偏置,从而真实地反映单元所在的空间位置。

成功找到应变片所粘贴的单元后,将应变片的测量方向投影到该单元坐标系下计算测量方向与单元坐标系x轴的夹角,如图9所示。

若需要通过插值计算应变片测量方向的应变,还应该计算插值点在应变片最近单元内的权重系数。插值点定义为应变片原点沿应变片z轴与最近单元的交点,并按照1.2节计算单元节点插值权重系数N_i(ξ₀,η₀)。若应变片最近单元是三角形单元,则单元节点插值权重系数直接取(1/3, 1/3, 1/3)。

此外,若单元应变输出结果为单元参考面应变和参考面曲率,还需要根据公式(7)进一步计算单元表面的应变:

这里,$\mathrm{}{\epsilon }_x^0$、${\epsilon }_y^0$和${\gamma }_y^0$为单元参考面应变,χ_x、χ_y和χ_xy为单元参考面曲率,z为计算点至参考面的距离。所以,还应该计算应变片所在的单元粘贴面距离单元参考平面的距离,并根据公式(7)计算应变片位置的应变。

3 应变监测值获取工具

根据本文所述方法和流程,在HyperWroks前后处理软件中使用tcl/tk语言^[10]开发出一套应变监测值获取工具,共包含前处理和后处理两个工具,如图10所示。

前处理工具根据应变片S模型生成应变片配置文件。应变片配置文件是在应变片几何搜索和插值计算的基础上,将后处理工具所需相关信息汇总得到的一个列表。为了在后处理工具中计算应变片的测量值,一个完整的应变片配置文件包括应变片号、应变片粘贴的单元、单元具体粘贴面、应变片测量方向与粘贴单元坐标系x轴的夹角、单元节点插值权重系数、单元表面至单元参考面的距离,如图11所示。

应变片配置文件生成之后,使用后处理工具计算应变片测量方向的应变。若应变输出设置为,则程序直接根据单元坐标系下Z1面或Z2面应变,并按照公式(2)计算应变片测量方向应变。若应变输出设置为,则程序首先根据单元坐标系下参考面的应变和曲率,结合公式(7)计算Z1面或Z2面的应变,最后按照公式(2)计算应变片测量方向的应变。当勾选“interpolate with corner data”时,程序使用单元节点应变进行单元面内应变插值。

输出的应变结果如图12所示,图13给出了机翼下壁板某处剖面应变对比结果,试验值和计算值吻合度较高。

此外,采用插值算法提高应变片计算精度是本文的一个重要研究内容。图14给出了机翼下壁板油箱检查口孔壁上一些应变片结果对比,油箱检查口附近有较大的应变梯度,可以看出采用插值方法后应变片的计算精度有明显提高。

4 结论

本文详细介绍了基于精细有限元模型的全机静强度试验应变监测值获取方法和流程,开发了相应的应变监测值获取工具。该方法和工具具有一般通用性,适用于所有结构强度试验。本文所述方法克服了传统方法使用工程方法带来的计算繁琐、效率低、局部计算精度不足的缺点,通过面内和面外插值提高了应变片监测值的计算精度,大大降低了应变片监测值数据的准备难度,有效地辅助了试验顺利进行。

参考文献