CN103366070B - 一种可用于直升机和固定翼飞行器的复合材料梁设计方法 - Google Patents

Abstract

本发明公开了一种可用于直升机和固定翼飞行器的复合材料梁设计方法，属于飞行器结构设计领域。本发明首先根据直升机或者固定翼飞行器的总体设计性能建立复合材料梁优化设计模型，设计变量为每一种材料在每一个单元处的体积分数；然后采用有限元软件建立复合材料梁剖面有限元模型；再次采用复合材料梁优化设计模型对复合材料剖面梁有限元模型进行优化设计，可以采用体积分数过滤方法或者敏度过滤方法，基于序列线性规划算法更新体积分数，得到梁剖面的材料分布规律；最后根据材料分布规律得到梁剖面的拓扑结构。本发明提出的复合材料梁剖面设计方法可以应用于各类型直升机旋翼和固定翼飞行器机翼的剖面设计。

Description

一种可用于直升机和固定翼飞行器的复合材料梁设计方法

技术领域

本发明属于飞行器结构设计技术领域，具体涉及一种可用于直升机和固定翼飞行器的复合材料梁设计方法，可以应用于概念设计阶段直升机旋翼和固定翼机翼的剖面结构设计。

背景技术

直升机旋翼和固定翼飞行器机翼是飞行器结构的重要组成部份，直接承受气动力，并为直升机和固定翼飞行器提供升力和操纵力。它们具有一个共同的特点：均为细长结构，因此在分析中往往被简化为梁结构。梁结构剖面特性对直升机和固定翼飞行器总体性能有重要影响。其剖面刚度特性是进行气动弹性分析和动力学设计的基础，合理的结构刚度能够大大提高动稳定性。剪心位置同样十分重要，作用在剪心的力将不会产生附加的扭矩，所以合理设置剪心和气动中心的相对位置对于减小剖面扭矩具有重要作用。在进行直升机或者固定翼飞行器设计的过程中，往往首先以梁剖面刚度特性、质量特性和剪心位置等作为设计变量对直升机或者固定翼飞行器的总体性能进行优化设计。然而为了达到指定的刚度特性和剪心位置等剖面性能，需要采用合适的材料和结构型式，但是如何选择并合理布置结构才能使得复合材料梁剖面达到指定的特性对飞行器结构设计非常重要。

由于复合材料优良的力学特性和较高的可设计性，在现代直升机和固定翼飞行器设计中已经取代传统的木质和金属材料成为了旋翼和机翼设计的主要材料之一。虽然复合材料的应用为直升机和固定翼飞行器设计提供了更多的选择，但是同样也使得梁结构设计变得更加复杂。传统的复合材料梁设计方法往往在拓扑结构不变的情况下采用复合材料的铺层角度、厚度和材料等作为设计变量，但是通过优化得到的旋翼和机翼剖面结构的拓扑型式没有改变，因此存在设计空间较小，可能难以满足指定的刚体特性和剪心位置等要求。

发明内容

针对传统的复合材料梁设计方法的不足，本发明提出将拓扑优化设计方法应用于复合材料梁设计，通过对材料类型及其分布进行重新设计，从而改变梁剖面结构的材料连接方式和相互关系，能够改变其拓扑结构型式，并得到指定刚度特性和剪心位置的复合材料梁结构。

本发明的基本技术方案是首先根据直升机和固定翼飞行器总体设计要求确定复合材料梁的气动外形、候选的材料类型及其材料常数和剖面特性，由此得到复合材料梁优化设计模型；其次建立复合材料梁剖面有限元模型；然后求解优化复合材料梁剖面有限元模型，得到最优的材料分布规律；最后根据材料分布得到梁剖面的拓扑结构。具体通过如下步骤实现：

（1）根据直升机或固定翼飞行器总体设计要求确定复合材料梁的气动外形、候选的材料类型及其材料常数和剖面特性，由此建立复合材料梁优化设计模型；

（2）建立复合材料梁剖面有限元模型，并且每一种材料在每一个有限元单元均对应一个设计变量，即体积分数；

（3）采用建立的复合材料梁优化设计模型对复合材料梁剖面有限元模型进行优化设计，得到最优的材料分布规律；

（4）由材料分布规律得到最佳的梁剖面的拓扑结构。

根据本发明提供的设计方法得到的复合材料梁考虑了多种材料类型，并且能够得到最佳的材料分布规律，为进一步进行气动弹性分析和动力学设计提供更多的选择，进而最大限度的改善直升机或固定翼飞行器的性能。

本发明提出的方法具有如下优点：

(1)复合材料梁可以考虑任意种类的材料，并为剖面设计提供更多的选择，最大限度提高直升机和固定翼飞行器的性能。

(2)采用刚度特性或者剪心位置等作为优化设计目标，并可以改变梁剖面的拓扑结构，从而实现指定的设计目标。

附图说明

图1是直升机和固定翼飞行器设计流程图；

图2是本发明中复合材料梁优化设计流程图；

图3是复合材料梁上梁剖面内坐标示意图；

图4是纤维平面示意图；

图5是纤维示意图；

图6是方形剖面的纤维方向示意图；

图7是方形剖面的纤维平面方向示意图；

图8是NACA0018剖面纤维方向示意图；

图9是NACA0018剖面纤维平面方向示意图。

图中:

1、复合材料梁；2、坐标系统；3、梁剖面；4、纤维平面；

5、纤维平面角度；6、纤维角度；7、纤维。

具体实施方式

下面结合附图对本发明提出的一种可用于直升机和固定翼飞行器的复合材料梁设计方法进行详细说明。

图1所示为直升机和固定翼飞行器设计的一般流程，首先根据总体性能要求确定气动弹性和动力学特性；然后确定直升机旋翼桨叶或者固定翼飞行器机翼的刚度和质量分布；最后由形状或者尺寸优化得到剖面结构型式和尺寸。在进行剖面设计过程中，往往采用固定剖面拓扑结构，以复合材料铺层角度、厚度或者尺寸作为设计变量的方法，难以满足指定的刚体特性和剪心位置等要求。图2是本发明提出的复合材料梁优化设计方法的流程图，本发明对现有技术中的复合材料梁优化设计方法进行改进，将拓扑优化设计方法引入到剖面刚度和剪心位置优化设计中，具体步骤如下：

（1）根据直升机或固定翼飞行器的总体性能的设计要求确定旋翼或机翼梁剖面的气动外形和剖面特性，并根据设计要求确定材料类型及其材料常数，建立复合材料梁优化设计模型，具体如下：

采用复合材料梁理论进行求解梁剖面特性。气动外形、刚度特性或者剪心位置要求等由直升机或固定翼飞行器的总体性能设计要求确定，根据气动外形对梁剖面进行网格划分，而刚度特性和剪心位置等则作为设计目标函数。根据设计要求确定所需的材料类型及其材料常数。

剪心位置和刚度是复合材料梁的重要参数：减小剪心位置和气动中心位置二者的间距可以减小由气动力产生的附加扭矩，从而改善梁剖面的受力状况；气动弹性和动力学特性则主要由剖面刚度决定。由此综合得到复合材料梁优化模型的目标函数由两部分组成，第一部分由各个刚度与刚度指定值之差加权平方和组成；第二部为剪心位置与气动中心位置的差值的平方。为了突出各个不同刚度和剪心位置的贡献，采用权值的办法，可以根据需要选择不同的权值。梁剖面可以采用m种材料模型，而梁剖面有n个不同的有限元单元组成，每一种材料在每一个有限元单元处均对应了一个设计变量－体积分数ρ_ij，每一个有限元单元的所有材料体积分数ρ_ij之和为1，因此共有nm个设计变量。约束条件中可以考虑质量约束条件，扭转刚度、弯曲刚度和刚度矩阵中的耦合项均可以作为约束条件出现，由此可以建立复合材料梁优化模型如下：

目标函数: $f=\underset{q=1}{\overset{6}{\mathrm{\Σ}}}\underset{p=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{pq}}{(K_{\mathrm{pq}}-K_{\mathrm{pq}0})}^2+{\mathrm{\ω}}_0({(x_{\mathrm{SC}}-x_{\mathrm{AC}})}^2+{(y_{\mathrm{SC}}-y_{\mathrm{AC}})}^2)$

约束条件:W≤W₀；

$\underset{q=1}{\overset{6}{\mathrm{\Σ}}}\underset{p=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\ω}}_{\mathrm{pq}}+{\mathrm{\ω}}_0=1;$

$\begin{array}{cc}0\≤{\mathrm{\ω}}_{\mathrm{pq}}\≤\mathrm{1,0}\≤{\mathrm{\ω}}_0\≤1& \∀1\≤p\≤\mathrm{6,1}\≤q\≤6;\end{array}$

$\begin{array}{cc}{\mathrm{\ρ}}_{\min }\≤{\mathrm{\ρ}}_{\mathrm{ij}}\≤1& \∀1\≤i\≤n,1\≤j\≤m,\end{array}$

$\begin{array}{cc}\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{ij}}=1& \∀1\≤i\≤n;\end{array}$

$K_{22}^L\≤K_{22}\≤K_{22}^U,$ $K_{33}^L\≤K_{33}\≤K_{33}^U;$

$\begin{array}{cc}{K}_{\mathrm{pq}}\≤K_{\mathrm{pq}}^U,& p,q=1,2...6;p\≠q;\end{array}$

其中SC、AC和K分别为剪心、气动中心和刚度矩阵。

W,W₀,n和m分别是设计过程中的复合材料梁质量、设计要求的复合材料梁质量、有限元单元数量和材料种类数量。第j种材料在第i个有限元单元的体积分数值为ρ_ij，所有有限元单元中每一种材料的体积分数ρ_ij的最小值均设置为ρ_min＝0.001。K₂₂，K₃₃为剖面扭转、弯曲刚度，而和则分别为扭转刚度的下界和上界、弯曲刚度的下界和上界；K_pq和分别为刚度矩阵的耦合项及其上界，K_pq0为指定的刚度目标值。x_AC和y_AC分别为气动中心的横坐标值和纵坐标值，x_SC和y_SC分别为剪心的横坐标值和纵坐标值。ω_pq(p,q＝1,2...6)为刚度目标函数的权重系数，ω₀为剪心目标函数的权重系数，这两类权重系数可以根据不同设计要求指定。

在复合材料梁的剖面刚度求解过程中有两个重要问题，首先是单一材料的梁剖面刚度求解过程，其次是考虑多材料影响后的刚度特性。

首先对于单一材料的梁剖面3刚度求解，复合材料梁1结构及其坐标系统2如图3所示，x轴和y轴位于复合材料梁剖面3内，z轴则沿着梁展向。复合材料梁1考虑了剪切效应，内力包括拉力、剪力、弯矩和扭矩等，并且忽略了锥形、预扭和曲度等因素的影响。梁剖面3平面内任意一点的位移由翘曲位移和刚体位移组成，由虚功原理建立复合材料梁1运动方程，并由此得到刚度矩阵，在求解过程中并没有考虑端部效应的影响。

采用本发明的方法首先得到复合材料梁剖面3的柔度矩阵如式（1）所示：

$F=\left[\begin{array}{ccc}{X}^T& X^{\′T}& Y\end{array}\right]\left[\begin{array}{ccc}E& C& R\\ C^T& M& L\\ R^T& L^T& V\end{array}\right]\left[\begin{array}{c}X\\ Y^{\′}\\ Y\end{array}\right]$

其中各个矩阵定义如式（2）所示：

$\begin{array}{cc}V={\∫}_A{Z}^T{S}^T\mathrm{QSZdA}& R={\∫}_A{N}^T{B}^T\mathrm{QSZdA}\\ E={\∫}_A{N}^T{B}^T\mathrm{QBNdA}& C={\∫}_A{N}^T{B}^T\mathrm{QSNdA}\\ L={\∫}_A{Z}^T{S}^T\mathrm{QSNdA}& M={\∫}_A{N}^T{B}^T\mathrm{QSNdA}\end{array}---\left(2\right)$

其中Q为有限元单元的柔度矩阵，A为梁剖面的剖面面积，N为有限元单元的型函数方程，B为应变位移矩阵，S＝[0₃I₃]，其中0₃和I₃分别为3×3的零矩阵和单位矩阵。Z＝[I₃nn^T]，I₃为3×3的单位矩阵，nn如式(3)所示：

$\mathrm{nn}=\left[\begin{array}{ccc}0& -z& y\\ z& 0& -x\\ -y& x& 0\end{array}\right]---\left(3\right)$

式中x,y和z均为复合材料梁坐标系xyz中的坐标值。

S＝[0₃I₃]，其中0₃和I₃为3×3的零矩阵和单位矩阵。

其中X、X'和Y、Y′可以由式（5-6）得到，并且X、X'和Y、Y′满足关系式X'＝dX/dz和Y′＝dY/dz。。得到柔度矩阵F后，梁剖面3刚度矩阵K可以有式（4）得到。

K=F^-1(4)

$\left[\begin{array}{ccc}E& R& D\\ R^T& V& 0\\ D^T& 0& 0\end{array}\right]\left[\begin{array}{c}X\\ Y\\ {\mathrm{\λ}}_1\end{array}\right]=\left[\begin{array}{ccc}{C}^T-C& L& 0\\ L^T& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}X\\ Y\\ {\mathrm{\λ}}_1\end{array}\right]+\left[\begin{array}{c}0\\ I_{6\×6}\\ 0\end{array}\right]---\left(5\right)$

$\left[\begin{array}{ccc}E& R& D\\ R^T& V& 0\\ D^T& 0& 0\end{array}\right]\left[\begin{array}{c}{X}^{\′}\\ Y^{\′}\\ {\mathrm{\λ}}_2\end{array}\right]=\left[\begin{array}{c}0\\ T_r\\ 0\end{array}\right]---\left(6\right)$

其中 $D=\left[\begin{array}{ccc}{I}_3& ...& I_3\\ {\mathrm{nn}}_1& ...& {\mathrm{nn}}_{\mathrm{node}}\end{array}\right],$ 此矩阵的引入可以消除翘曲位移对刚体位移的可能的贡献,nn_node由式（3）可得，具体可通过將节点node处的单位节点坐标值(x,y,z)代入式（3）后得到，node为复合材料梁剖面的节点数目。λ₁和λ₂为拉格朗日乘子，I_6×6为6×6的单位矩阵，T_r如式(7)所示：

$T_r=\left[\begin{array}{cc}{0}_3& t_r\\ 0_3& 0_3\end{array}\right],$ 其中 $t_r=\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right]---\left(7\right)$

对于多材料的梁剖面3刚度求解，主要利用材料柔度矩阵将多材料问题转变为单一材料问题。具体分析如下：

由于梁剖面3每一个有限元单元由m种材料组成，因此需要将各种不同材料的柔度矩阵组装成新的单元柔度矩阵。本发明中不但可以采用各向同性的材料，也可以采用各向异性的材料，需要采用式(8)将各种不同材料的柔度矩阵组装为单元的柔度矩阵Q_i。

$Q_i={{\mathrm{\ρ}}_{i1}}^{\mathrm{\η}}{\stackrel{\‾}{Q}}_1+{{\mathrm{\ρ}}_{i2}}^{\mathrm{\η}}{\stackrel{\‾}{Q}}_2+...+{{\mathrm{\ρ}}_{\mathrm{im}}}^{\mathrm{\η}}{\stackrel{\‾}{Q}}_m---\left(8\right)$

其中Q_i(i＝1,2,3...n)为第i个有限元单元的柔度矩阵，ρ_im为第m种材料在第i个有限元单元处的体积分数，为第m种材料的柔度矩阵，处罚因子η可以设置为3。

各向异性复合材料的纤维平面4和纤维7如图4和图5所示，一般的各向异性材料的柔度矩阵如式（9）所示。

${\stackrel{\‾}{Q}}_o=\left[\begin{array}{cccccc}\frac{1}{E_x}& -\frac{v_{\mathrm{xy}}}{E_x}& 0& 0& 0& -\frac{v_{\mathrm{xz}}}{E_z}\\ -\frac{v_{\mathrm{xy}}}{E_x}& \frac{1}{E_y}& 0& 0& 0& -\frac{v_{\mathrm{yz}}}{E_z}\\ 0& 0& \frac{1}{G_{\mathrm{xy}}}& 0& 0& 0\\ 0& 0& 0& \frac{1}{G_{\mathrm{xz}}}& 0& 0\\ 0& 0& 0& 0& \frac{1}{G_{\mathrm{yz}}}& 0\\ -\frac{v_{\mathrm{xz}}}{E_z}& -\frac{v_{\mathrm{yz}}}{E_z}& 0& 0& 0& \frac{1}{E_z}\end{array}\right]---\left(9\right)$

其中E_x、E_y和E_z为三个方形的弹性模量，v_xy、v_yz和v_zx分别为泊松比，G_xy、G_yz和G_zx分别为xy、yz和zx平面内的剪切刚度。对于不同角度的各向异性材料的柔度矩阵可由式（10）得到。

${\stackrel{\‾}{Q}}_i=\mathrm{\Ω}{\stackrel{\‾}{Q}}_o{\mathrm{\Ω}}^T(i=\mathrm{1,2}...m)---\left(10\right)$

其中Ω为对应于纤维角度6和纤维平面角度5的转换矩阵，和分别为转换后的材料柔度矩阵和原始刚度矩阵。

将各个转换后的单元材料柔度矩阵代入式(8)即可得到多材料的梁剖面3有限元单元的柔度矩阵。得到有限元单元的柔度矩阵之后，即可由公式（1-7）得到复合材料梁剖面的柔度矩阵和刚度矩阵。

由此得到了梁剖面3的柔度矩阵，并且剪心位置(x_SC,y_SC)可以有式（11-12）得到：

$x_{\mathrm{SC}}=-\frac{F_{62}+F_{64}(L-z)}{F_{66}}---\left(11\right)$

$y_{\mathrm{AC}}=\frac{F_{61}+F_{65}(L-z)}{F_{66}}---\left(12\right)$

其中F_ij(i,j＝1,2...6)为复合材料梁剖面的柔度矩阵中(i,j)分量，L和z分别为复合材料梁1长度和梁剖面3的位置。

（2）采用有限元软件建立梁剖面3的有限元模型，导出网格、节点等信息，每一种材料在每一个有限元单元均对应一个设计变量，即体积分数，每一个有限元单元处的所有材料的体积分数之和为1；单元、节点和材料信息等均作为优化设计的基础。

（3）采用建立的复合材料梁优化模型对梁剖面3有限元模型进行优化设计，得到最优的材料分布规律；

为了进行优化设计，敏度分析和优化算法是必不可少的。首先对于敏度分析，敏度是更新单元体积分数的重要依据，本发明中的敏度求解由两个部分组成，首先对式（5-6）设计变量求导，可以得到和然后对梁剖面柔度矩阵（1）设计变量求导，得到的敏度信息如式（13）所示。

$\frac{\∂F}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}={\mathrm{\η\ρ}}_{\mathrm{ij}}^{\mathrm{\η}-1}\left[\begin{array}{ccc}{X}^T& X^{\′T}& Y^T\end{array}\right]\left[\begin{array}{ccc}{E}_0& C_0& R_0\\ {C_0}^T& M_0& L_0\\ {R_0}^T& {L_0}^T& V_0\end{array}\right]\left[\begin{array}{c}X\\ X^{\′}\\ Y\end{array}\right]+\left[\begin{array}{ccc}\frac{{\∂X}^T}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}& \frac{{\∂X}^{\′T}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}& \frac{{\∂Y}^T}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}\end{array}\right]\left[\begin{array}{ccc}E& C& R\\ C^T& M& L\\ R^T& L^T& V\end{array}\right]\left[\begin{array}{c}X\\ X^{\′}\\ Y\end{array}\right]+\left[\begin{array}{ccc}{X}^T& X^{\′T}& Y^T\end{array}\right]\left[\begin{array}{ccc}E& C& R\\ C^T& M& L\\ R^T& L^T& V\end{array}\right]\left[\begin{array}{c}\frac{\∂X}{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}\\ \frac{\∂X^{\′}}{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}\\ \frac{\∂Y}{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}\end{array}\right]$

(13)

上式中，下标为0的矩阵为单元体积分数为1时的矩阵。由此得到了柔度矩阵关于各个有限元单元设计变量的敏度信息，并作为更新有限元单元体积分数的依据。

容易得到梁剖面3刚度和剪心关于设计变量ρ_ij的敏度可由式(15-17)得到：

$\frac{\∂K}{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}=\frac{\∂F^{-1}}{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}=-K\frac{\∂F}{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}K---\left(15\right)$

$\frac{{\∂x}_{\mathrm{SC}}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}=-\frac{\frac{{\∂F}_{62}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}+\frac{{\∂F}_{64}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}(L-z)}{{\left(\frac{{\∂F}_{66}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}\right)}^2}{F}_{66}+\frac{F_{62}+F_{64}(L-z)}{{\left(\frac{{\∂F}_{66}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}\right)}^2}\frac{{\∂F}_{66}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}---\left(16\right)$

$\frac{{\∂y}_{\mathrm{SC}}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}=-\frac{\frac{{\∂F}_{61}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}+\frac{{\∂F}_{65}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}(L-z)}{{\left(\frac{{\∂F}_{66}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}\right)}^2}{F}_{66}+\frac{F_{61}+F_{65}(L-z)}{{\left(\frac{{\∂F}_{66}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}\right)}^2}\frac{{\∂F}_{66}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}---\left(17\right)$

目标函数的敏度求解公式如式(18)所示，将式(15-17)代入即可得到目标函数的敏度信息。

$\frac{\∂f}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}=\underset{p=1}{\overset{6}{\mathrm{\Σ}}}\underset{q=1}{\overset{6}{\mathrm{\Σ}}}2{\mathrm{\ω}}_{\mathrm{pq}}(K_{\mathrm{pq}}-K_{\mathrm{pq}0})\frac{{\∂K}_{\mathrm{pq}}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}+2{\mathrm{\ω}}_0((x_{\mathrm{SC}}-x_{\mathrm{AC}})\frac{{\∂x}_{\mathrm{SC}}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}}+(y_{\mathrm{SC}}-y_{\mathrm{AC}})\frac{{\∂y}_{\mathrm{SC}}}{{\∂\mathrm{\ρ}}_{\mathrm{ij}}})---\left(18\right)$

为了避免棋盘格问题，得到目标函数关于各个有限元单元设计变量的敏度信息之后，可以采用体积分数过滤算法或者敏度过滤算法对体积分数进行处理，并采用序列线性规划算法更新每一个有限元单元的设计变量。

（4）由材料分布规律得到最佳的拓扑结构，为进一步进行详细设计奠定基础。

求得了各种材料在复合材料梁中的分布情况，即可得到梁剖面3的拓扑形状，为进一步进行详细设计奠定基础。

实施例

为了验证本发明提出的复合材料梁设计方法，分别对方形剖面和NACA0018两种剖面形状的固定翼飞行器机翼剖面拓扑结构进行设计。方形剖面由25×25个网格组成，每一个网格均具有9种类型的材料。其中包括了一种各向同性材料和一种各向异性材料，并且各向异性材料考虑了两种纤维平面方向和四种纤维方向，因此共有八种候选的各向异性材料类型。两种材料的材料属性如表1所示:

表1材料属性

方形剖面梁长度为80mm，端部剪切力载荷，最终得到的剖面拓扑结构中的纤维和纤维平面方向分别如图6和图7所示。单元空白表示该单元主要由各向同性材料组成，而单元中的直线长短及其角度则代表对应的各向异性材料的体积分数及其种类。目标函数为x方向的弯曲刚度K₄₄和y方向的剪切刚度K₂₂，目标刚度均设置为0，约束函数考虑了质量约束条件。最终得到的材料分布及其拓扑结构如图6和图7所示，由图可知材料主要分布于剖面两端，主要是由弯曲刚度的要求所致，而目标函数中的剪切刚度则使得纤维角度主要为0°方向。纤维平面角度在靠近中部中性轴的位置主要有45°和-45°。纤维平面在左上和右下方向为0°，其余部位均为90°。

在第二个算例中采用了NACA0018翼型，梁剖面由1266个单元组成，复合材料梁长度为200mm，采用了与前面算例中相同的材料类型和数量，并且选择了和前一算例一样的目标函数和约束函数。最终得到的纤维平面和纤维角度如图8和图9所示。结论与方形剖面类似，但稍有不同的是梁剖面中出现多个不同的传递剪切变形的路径。

算例表明，本发明提出的复合材料梁是可靠的，可以应用于直升机旋翼或者固定翼飞行器机翼设计。

Claims (3)

1.一种可用于直升机和固定翼飞行器的复合材料梁设计方法，其特征在于，包括如下步骤：

第一步，根据直升机或者固定翼飞行器总体性能要求确定复合材料梁的气动外形和梁剖面特性、指定候选的材料类型及其属性，将每一种材料在每一个单元的体积分数度均作为一个设计变量，以梁剖面特性和剪心位置作为设计目标函数，约束条件为质量、体积分数和刚度要求，建立复合材料梁优化模型如下：

目标函数: $f=\underset{q=1}{\overset{6}{\Σ}}\underset{p=1}{\overset{6}{\Σ}}{\mathrm{\ω}}_{pq}{(K_{pq}-K_{pq0})}^2+{\mathrm{\ω}}_0({\left(x_{SC}-x_{AC}\right)}^2+{\left(y_{SC}-y_{AC}\right)}^2)$

约束条件:W≤W₀

$\underset{q=1}{\overset{6}{\Σ}}\underset{p=1}{\overset{6}{\Σ}}{\mathrm{\ω}}_{pq}+{\mathrm{\ω}}_0=1;$

$0\≤{\mathrm{\ω}}_{pq}\≤1,0\≤{\mathrm{\ω}}_0\≤1;\∀1\≤p\≤6,1\≤q\≤6$

${\mathrm{\ρ}}_{min}\≤{\mathrm{\ρ}}_{ij}\≤1,\∀1\≤i\≤n,1\≤j\≤m$

$\underset{j=1}{\overset{m}{\Σ}}{\mathrm{\ρ}}_{ij}=1,\∀1\≤i\≤n;$

$K_{22}^L\≤K_{22}\≤K_{22}^U,K_{33}^L\≤K_{33}\≤K_{33}^U;$

K_pq≤K^U _p，p＝1,2,…,6；q＝1,2,…,6；p≠q；

W,W₀,n和m分别是设计过程中的复合材料梁质量、设计要求的复合材料梁质量、有限元单元数量和材料种类数量；第j种材料在第i个有限元单元的体积分数值为ρ_ij，所有有限元单元中每一种材料的体积分数ρ_ij的最小值均设置为ρ_min＝0.001；K₂₂，K₃₃为剖面扭转、弯曲刚度，而和和则分别为扭转的下界和上界、弯曲刚度的下界和上界；K_pq和分别为刚度矩阵的耦合项及其上界，K_pq0为指定的刚度目标值；x_AC，y_AC，x_SC和y_SC分别为气动中心横坐标值和纵坐标值，剪心的横坐标值和纵坐标值；ω_pq为刚度目标函数的权重系数，ω₀为剪心目标函数的权重系数，p＝1,2,…,6；q＝1,2,…,6；

第二步，采用有限元软件建立梁剖面的有限元模型，导出网格、节点信息；

第三步，采用复合材料梁优化模型对指定的梁剖面有限元模型进行优化设计，得到材料在梁剖面的分布规律；

第四步，根据得到的材料在梁剖面的分布规律得到最终的复合材料梁的拓扑结构；

所述的梁剖面采用任意数量和种类的材料，每一种材料在每一个有限元单元处均有一个设计变量，即体积分数，且每一个有限元单元处的所有材料的体积分数之和为1；梁剖面的刚度求解，首先是单一材料的梁剖面刚度求解过程，其次是考虑多材料影响后的刚度特性，具体如下：

对于单一材料的梁剖面刚度求解，复合材料梁结构及其坐标系统，x轴和y轴位于复合材料梁剖面内，z轴则沿着梁展向；复合材料梁考虑了剪切效应，内力包括拉力、剪力、弯矩和扭矩，并且忽略了锥形、预扭和曲度因素的影响；梁剖面平面内任意一点的位移由翘曲位移和刚体位移组成，由虚功原理建立复合材料梁运动方程，并由此得到刚度矩阵，在求解过程中并没有考虑端部效应的影响；首先得到复合材料梁剖面的柔度矩阵如式(1)所示：

$F=\left[\begin{array}{ccc}{X}^T& X^{\′T}& Y\end{array}\right]\left[\begin{array}{ccc}E& C& R\\ C^T& M& L\\ R^T& L^T& V\end{array}\right]\left[\begin{array}{c}X\\ Y^{\′}\\ Y\end{array}\right]---\left(1\right)$

其中各个矩阵定义如式(2)所示：

V＝∫_AZ^TS^TQSZdAR＝∫_AN^TB^TQSZdA

E＝∫_AN^TB^TQBNdAC＝∫_AN^TB^TQSNdA(2)

L＝∫_AZ^TS^TQSNdAM＝∫_AN^TB^TQSNdA

其中Q为有限元单元的柔度矩阵，A为梁剖面的剖面面积，N为有限元单元的型函数方程，B为应变位移矩阵，S＝[0₃I₃]，其中0₃和I₃分别为3×3的零矩阵和单位矩阵，Z＝[I₃nn^T]，I₃为3×3的单位矩阵，nn如式(3)所示：

$nn=\left[\begin{array}{ccc}0& -z& y\\ z& 0& -x\\ -y& x& 0\end{array}\right]---\left(3\right)$

式中x,y和z均为复合材料梁坐标系xyz中的坐标值；

S＝[0₃I₃]，其中0₃和I₃为3×3的零矩阵和单位矩阵；

其中X、X'和Y、Y′可以由式(5-6)得到，并且X、X'和Y、Y′满足关系式X'＝dX/dz和Y′＝dY/dz；得到柔度矩阵F后，梁剖面刚度矩阵K由式(4)得到：

K＝F^-1(4)

$\left[\begin{array}{ccc}E& R& D\\ R^T& V& 0\\ D^T& 0& 0\end{array}\right]\left[\begin{array}{c}X\\ Y\\ {\mathrm{\λ}}_1\end{array}\right]=\left[\begin{array}{ccc}{C}^T-C& L& 0\\ L^T& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}X\\ Y\\ {\mathrm{\λ}}_1\end{array}\right]+\left[\begin{array}{c}0\\ I_{6\×6}\\ 0\end{array}\right]---\left(5\right)$

$\left[\begin{array}{ccc}E& R& D\\ R^T& V& 0\\ D^T& 0& 0\end{array}\right]\left[\begin{array}{c}{X}^{\′}\\ Y^{\′}\\ {\mathrm{\λ}}_2\end{array}\right]=\left[\begin{array}{c}0\\ T_r\\ 0\end{array}\right]---\left(6\right)$

其中 $D=\left[\begin{array}{ccc}{I}_3& ...& I_3\\ {\mathrm{nn}}_1& ...& {\mathrm{nn}}_{node}\end{array}\right],$ nn_node由式(3)得，具体通过将节点node处的单位节点坐标值(x,y,z)代入式(3)后得到，node为复合材料梁剖面的节点数目，λ₁和λ₂为拉格朗日乘子，I_6×6为6×6的单位矩阵，T_r如式(7)所示：

$T_r=\left[\begin{array}{cc}{0}_3& t_r\\ 0_3& 0_3\end{array}\right],$ 其中 $t_r=\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right]---\left(7\right)$

对于多材料的梁剖面刚度求解，由于梁剖面每一个有限元单元由m种材料组成，因此需要采用式(8)将各种不同材料的柔度矩阵组装为单元的柔度矩阵Q_i：

$Q_i={{\mathrm{\ρ}}_{i1}}^{\mathrm{\η}}{\stackrel{\‾}{Q}}_1+{{\mathrm{\ρ}}_{i2}}^{\mathrm{\η}}{\stackrel{\‾}{Q}}_2+...+{{\mathrm{\ρ}}_{im}}^{\mathrm{\η}}{\stackrel{\‾}{Q}}_m---\left(8\right)$

其中Q_i为第i个有限元单元的柔度矩阵，ρ_im为第m种材料在第i个有限元单元处的体积分数，为第j种材料的柔度矩阵，η为处罚因子，j＝1,2...m；i＝1,2,3...n；

各向异性材料的柔度矩阵如式(9)所示：

${\stackrel{\‾}{Q}}_O=\left[\begin{array}{cccccc}\frac{1}{E_x}& -\frac{v_{xy}}{E_x}& 0& 0& 0& -\frac{v_{xz}}{E_z}\\ -\frac{v_{xy}}{E_x}& \frac{1}{E_y}& 0& 0& 0& -\frac{v_{yz}}{E_z}\\ 0& 0& \frac{1}{G_{xy}}& 0& 0& 0\\ 0& 0& 0& \frac{1}{G_{xz}}& 0& 0\\ 0& 0& 0& 0& \frac{1}{G_{yz}}& 0\\ -\frac{v_{xz}}{E_z}& -\frac{v_{yz}}{E_z}& 0& 0& 0& \frac{1}{E_z}\end{array}\right]---\left(9\right)$

其中E_x、E_y和E_z为三个方形的弹性模量，v_xy、v_yz和v_zx分别为泊松比，G_xy、G_yz和G_zx分别为xy、yz和zx平面内的剪切刚度；对于不同角度的各向异性材料的柔度矩阵由式(10)得到：

${\stackrel{\‾}{Q}}_j=\mathrm{\Ω}{\stackrel{\‾}{Q}}_O{\mathrm{\Ω}}^T---\left(10\right)$

其中Ω为对应于纤维角度和纤维平面角度的转换矩阵，和分别为转换后的材料柔度矩阵和原始刚度矩阵；j＝1,2...m；

将各个转换后的单元材料柔度矩阵代入式(8)即得到多材料的梁剖面有限元单元的柔度矩阵；得到有限元单元的柔度矩阵之后，即由公式(1-7)得到复合材料梁剖面的柔度矩阵和刚度矩阵；

由此得到了梁剖面的柔度矩阵，并且剪心位置(x_SC,y_SC)有式(11-12)得到：

$x_{SC}=-\frac{F_{62}+F_{64}(L-z)}{F_{66}}---\left(11\right)$

$y_{SC}=\frac{F_{61}+F_{65}(L-z)}{F_{66}}---\left(12\right)$

其中F_ij为复合材料梁剖面的柔度矩阵中(i,j)分量，i＝1,2...6；j＝1,2...6；L和z分别为复合材料梁长度和梁剖面的位置。

2.根据权利要求1所述的一种可用于直升机和固定翼飞行器的复合材料梁设计方法，其特征在于：第三步中，为了进行优化设计，得到的敏度信息如式(13)所示：

$\frac{\∂F}{\∂{\mathrm{\ρ}}_{ij}}={\mathrm{\η\ρ}}_{ij}^{\mathrm{\η}-1}\left[\begin{array}{ccc}{X}^T& X^{\′T}& Y^T\end{array}\right]\left[\begin{array}{ccc}{E}_0& C_0& R_0\\ {C_0}^T& M_0& L_0\\ {R_0}^T& {L_0}^T& V_0\end{array}\right]\left[\begin{array}{c}X\\ Y^{\′}\\ Y\end{array}\right]+\left[\begin{array}{ccc}\frac{\∂X^T}{\∂{\mathrm{\ρ}}_{ij}}& \frac{\∂X^{\′T}}{\∂{\mathrm{\ρ}}_{ij}}& \frac{\∂Y^T}{\∂{\mathrm{\ρ}}_{ij}}\end{array}\right]\left[\begin{array}{ccc}E& C& R\\ C^T& M& L\\ R^T& L^T& V\end{array}\right]\left[\begin{array}{c}X\\ Y^{\′}\\ Y\end{array}\right]+\left[\begin{array}{ccc}{X}^T& X^{\′T}& Y^T\end{array}\right]\left[\begin{array}{ccc}E& C& R\\ C^T& M& L\\ R^T& L^T& V\end{array}\right]\left[\begin{array}{c}\frac{\∂X}{\∂{\mathrm{\ρ}}_{ij}}\\ \frac{\∂X^{\′}}{\∂{\mathrm{\ρ}}_{ij}}\\ \frac{\∂Y}{\∂{\mathrm{\ρ}}_{ij}}\end{array}\right]---\left(13\right)$

上式中，下标为0的矩阵为单元体积分数为1时的矩阵；由此得到了柔度矩阵关于各个有限元单元设计变量的敏度信息，并作为更新有限元单元体积分数的依据；

得到梁剖面刚度和剪心关于设计变量ρ_ij的敏度由式(15-17)得到：

$\frac{\∂K}{\∂{\mathrm{\ρ}}_{ij}}=\frac{\∂F^{-1}}{\∂{\mathrm{\ρ}}_{ij}}=-K\frac{\∂F}{\∂{\mathrm{\ρ}}_{ij}}K---\left(15\right)$

$\frac{\∂x_{SC}}{\∂{\mathrm{\ρ}}_{ij}}=-\frac{\frac{\∂F_{62}}{\∂{\mathrm{\ρ}}_{ij}}+\frac{\∂F_{64}}{\∂{\mathrm{\ρ}}_{ij}}(L-z)}{{\left(\frac{\∂F_{66}}{\∂{\mathrm{\ρ}}_{ij}}\right)}^2}{F}_{66}+\frac{F_{62}+F_{64}(L-z)}{{\left(\frac{\∂F_{66}}{\∂{\mathrm{\ρ}}_{ij}}\right)}^2}\frac{\∂F_{66}}{\∂{\mathrm{\ρ}}_{ij}}---\left(16\right)$

$\frac{\∂y_{SC}}{\∂{\mathrm{\ρ}}_{ij}}=-\frac{\frac{\∂F_{61}}{\∂{\mathrm{\ρ}}_{ij}}+\frac{\∂F_{65}}{\∂{\mathrm{\ρ}}_{ij}}(L-z)}{{\left(\frac{\∂F_{66}}{\∂{\mathrm{\ρ}}_{ij}}\right)}^2}{F}_{66}+\frac{F_{61}+F_{65}(L-z)}{{\left(\frac{\∂F_{66}}{\∂{\mathrm{\ρ}}_{ij}}\right)}^2}\frac{\∂F_{66}}{\∂{\mathrm{\ρ}}_{ij}}---\left(17\right)$

目标函数的敏度求解公式如式(18)所示，将式(15-17)代入即得到目标函数的敏度信息；

$\frac{\∂f}{\∂{\mathrm{\ρ}}_{ij}}=\underset{p=1}{\overset{6}{\Σ}}\underset{q=1}{\overset{6}{\Σ}}2{\mathrm{\ω}}_{pq}(K_{pq}-K_{pq0})\frac{\∂K_{pq}}{\∂{\mathrm{\ρ}}_{ij}}+2{\mathrm{\ω}}_0\left(\right(x_{SC}-x_{AC})\frac{\∂x_{SC}}{\∂{\mathrm{\ρ}}_{ij}}+(y_{SC}-y_{AC}\left)\frac{\∂y_{SC}}{\∂{\mathrm{\ρ}}_{ij}}\right)---\left(18\right).$

3.根据权利要求2所述的一种可用于直升机和固定翼飞行器的复合材料梁设计方法，其特征在于：为了避免棋盘格问题，得到目标函数关于各个有限元单元设计变量的敏度信息之后，采用密度过滤算法或者灵敏度过滤算法对体积分数进行处理，并采用序列线性规划算法更新每一个有限元单元的设计变量。