非均匀纤维增强复合材料分布优化设计方法

摘要

本发明公开非均匀纤维增强复合材料分布优化设计方法，包括步骤：第一步、设计问题初始化；第二步、形成有限元刚度方程；第三步、求解有限单元刚度方程及应力、应变；第四步：优化材料的方向和密度；第五步：迭代收敛性检查；第六步：优化结果的后处理。本发明具有较高的优化计算效率，适合于任意不均匀分布的纤维增强复合材料优化设计方法。

说明书

技术领域

本发明属于非均匀各向异性材料优化设计技术领域，特别一种涉及非均匀纤维增强复合材料分布优化设计方法

背景技术

在航天航空、材料、机械、土木以及船舶和水利等领域经常需要各种复杂的纤维增强复合材料设计。纤维增强复合材料(Fiber Reinforced Polymer,或Fiber ReinforcedPlastic,简称FRP)是由增强纤维材料,如玻璃纤维,碳纤维,芳纶纤维等与基体材料经过缠绕、模压或拉挤等成型工艺而形成的复合材料，这些材料及其结构不仅需要满足使用功能的要求，而且希望使用的材料尽可能少。如何优化设计各种不同性能的纤维增强复合材料是上述各领域需要解决的关键问题，然而，解决这些问题并不是一件容易的事情，由于材料科学和制造技术的发展，使得我们能够制造材料范围不断扩展，但是如何科学设计的问题确没有能够合理解决。这属于材料与结构的拓扑优化设计问题。

目前的各种连续体结构拓扑优化方法主要是集中在设计均匀等厚带孔板或匀质带空洞的空间连续体结构，如均匀化方法(Homogenization Method)，进化结构优化方法(Evolutionary Structural Optimization Method，ESO)，水平集方法(Level SetMethod)，独立连续映射(Independent Continuous Mapping,ICM)优化方法，各向同性固体罚中间密度方法(Solid Isotropic Material with Penalization,SIMP)。这些方法都不能用于非均匀各向异性纤维增强复合材料优化设计。

自由材料优化设计方法(Free Material Optimization)用于一般非均匀各向异性材料优化设计，但是，其假设的各向异性材料优化设计变量是弹性张量的所有分量。这种材料更多情况下只是一种理论上的抽象材料，一般难于找到对应的真实材料，甚至其质量的含义都很模糊，难以定义，所以目前只是一种理论探讨，不能解决实际工程优化设计问题。

常规的符合材料优化设计的方法只能优化设计几层的铺设角度，在每一层内，而纤维的方向只能是同一个角度，每一层的厚度也是确定的，故常规方法不能对任意位置的材料分布优化。

针对上诉问题，本发明人提出一种非均匀纤维增强复合材料分布优化设计方法。

发明内容

本发明的目的在于提供非均匀纤维增强复合材料分布优化设计方法，以解决非均匀各向异性纤维增强复合材料拓扑优化设计问题。

为了实现上述目的，本发明采用如下技术方案：

非均匀纤维增强复合材料分布优化设计方法，包括以下步骤：

第一步、设计问题初始化

根据实际使用要求设定一个设计域，采用有限元方法将设计域划分有限单元，建立有限元模型，再根据实际工作状况对有限单元模型施加位移约束条件和荷载；

初始化纤维增强复合材料在结点位置的方向和密度；

第二步、形成有限元刚度方程

根据当前纤维增强复合材料的分布形成结点位置的弹性矩阵，由结点位置的弹性矩阵，通过形函数插值得到有限单元内任意位置的弹性矩阵，再由弹性矩阵计算出有限单元刚度矩阵，接着由有限单元刚度矩阵按照结点对应位置累加得到结构刚度矩阵，最后由结构刚度矩阵得到有限单元刚度方程；

第三步、求解有限单元刚度方程及应力、应变

求解有限单元刚度方程，得到结构在结点位置的位移，并据此计算出结构单元内的应力和应变，再根据结构单元内的应变计算出结点位置的应变；由应力分量计算出主应力的方向以及主应力的方向的应变值；

第四步：优化材料的方向和密度

将纤维增强复合材料的两个主方向调整到主应力方向；根据满应力准则的应力比公式调整纤维增强复合材料在结点位置的密度值；

第五步：迭代收敛性检查

比较纤维增强复合材料在结点位置的方向和密度在两次连续迭代的相对变化量，当相对变化量小于一个预设的数值时，优化迭代结束；否则返回第二步；

第六步：优化结果的后处理

采用有限单元的形函数，根据纤维增强复合材料在结点位置的方向和密度插值得到各单元内部任意位置的纤维增强复合材料分布。

在所述第一步中，当初始化的纤维增强复合材料采用的是正交纤维增强复合材料模型时，初始角度为零，初始密度为任意值。

所述第一步中，采用有限元方法将设计域划分有限单元具体是：对于规则区域选用矩形规则单元进行划分；对于不规则区域，则采用三角形单元、四结点四边形单元或八结点四边形单元进行划分。

所述第二步至第五步具体是：

第二步：

设初始化的纤维增强复合材料在结点位置的初始角度为零，初始密度为任意值，设迭代指标为i＝0，则

式(1)中上角标0表示初始迭代指标，下角标j表示结点编号，J为结点总数，b表示纤维增强复合材料的两个主轴编号；

根据当前结构纤维增强复合材料的分布形成结点j位置的弹性矩阵，该弹性矩阵如下：

式(2)中E是弹性模量，上角标i为迭代指标，和是材料当前迭代下在结点位置的方向和密度，s_br是常数矩阵的分量，且

g_r(α_j)是材料方向的函数矩阵的分量，且

g(α)＝[cos2α sin2α 1](4)

A_r是常数矩阵，且

由结点位置的弹性矩阵式，通过形函数插值得到有限单元内任意点位置的弹性矩阵如下：

式(6)中ξ,η为单元内部任意位置的局部坐标，N_j(ξ,η)是形函数；

由弹性矩阵式可以计算出单元刚度矩阵

式(8)中H_ejr是与设计变量无关的常数矩阵；

由单元刚度矩阵式按照结点对应位置累加得到结构刚度矩阵，

最终形成刚度方程，

KⁱUⁱ＝F(11)

式中Uⁱ和F分别是当前结构有限元方法的结点位移列向量和结点力列向量；

第三步：

依照下式求解结构刚度方程得到结构在结点位置的位移

Uⁱ＝(Kⁱ)^-1F(12)

并据此计算出结构单元内的应力和应变

式(13)中εⁱ是应变向量，和分别是纤维增强复合材料沿x和y坐标轴方向的正应变和剪应变分量，B和是几何矩阵和单元结点位移向量，根据结构单元内的应变可以计算出结点位置的应变，

式中表示结点位置的应变，εⁱ表示围绕结点j的各单元在结点位置的应变，n_j表示围绕结点j的单元数量，S_j为围绕结点j的单元集合，由线弹性假设的胡克定律可以计算出应力

和分别是纤维增强复合材料在结点j位置沿x和y坐标轴方向的正应力和剪应力分量，由应力分量可以计算出主应力的方向

以及沿主应力方向的应变值，

第四步：

将纤维增强复合材料的两个主方向调整到主应力方向，

根据满应力准则的应力比公式

调整纤维增强复合材料在结点位置的密度值，式(19)中是结点j位置在第i次迭代时的应变，ε_p是材料的允许最大应变；

第五步：

通过式(22)比较纤维增强复合材料在结点位置的方向和密度在两次连续迭代的相对变化量，当变化量小于一个事先给定的数值时，优化迭代结束，否则迭代指标增加1，i＝i+1，并返回第二步重新计算

式(22)中δ为事先指定的一个很小的值；

第六步：

采用有限元的形函数，根据纤维增强复合材料在结点位置的方向和密度插值得到各单元内部任意位置的纤维增强复合材料分布，如下：

得到的结构材料体积为

式(24)中

对于规则矩形单元，

式中V_e是单元体积。

所述第二步中，所述当有限单元采用四边形单元时，形函数为

N_j(ξ,η)＝(1+ξ_jξ)(1+η_jη)/4，j＝1,2,3,4(7)

式(7)中ξ_j和η_j是结点在局部坐标系下的坐标值；

当有限单元采用三角形单元时，则使用面积坐标作为形函数。

若划分所述有限单元采用规则单元网格，则H_ejr与有限单元无关，且通过公式(9)计算；

所述第四步中，在调整纤维增强复合材料在结点位置的密度值时，通过式(20)来限制材料在结点位置的密度不要过低

式中是本次迭代得到材料所有结点位置的密度最大值

采用上述方案后，本发明有益效果是：本发明基于纤维增强复合材料模型，在满足应力约束的条件下实现结构的材料体积最小化。本方法首先假设在设计域内均匀连续分布一种纤维增强复合材料，采用有限元方法进行结构分析，得到应力，应变分布场。根据结点位置的主应力方向以及主应力方向的应变大小，采用满应力准则的应力比公式优化纤维增强复合材料在结点位置的分布方向和密度。由形函数插值得到任意位置材料分布。经过反复迭代，直到收敛。最后形成的纤维增强复合材料优化分布场用材料在结点位置的分布方向和密度表示，不需要分层。

本发明具有较高的优化计算效率，而且，相比常规的复合增强材料优化设计方法仅能够设计几种分层结构的优化分布，每层的内部完全一致，不能优化，本发明则并没有采用分层设计，不受分层限制，材料允许任意分布，具有更大的优化设计空间。

下面结合附图对本发明做进一步的说明。

附图说明

附图1是本发明非均匀纤维增强复合材料分布优化方法的流程图；

附图2是本发明一个矩形梁结构的初始设计域；

附图3是采用本发明方法设计获得的最佳材料分布。

具体实施方式

如图1所示本实施例揭示的非均匀纤维增强复合材料分布优化设计方法，具体包括以下步骤：

第一步、设计问题初始化

根据实际使用要求设定一个设计域，采用有限元方法将设计域划分有限单元，形成有限单元模型，采用有限元方法将设计域划分有限单元模型具体是：对于规则区域选用矩形规则单元进行划分；对于不规则区域，则采用三角形单元、四结点四边形单元或八结点四边形单元进行划分。再根据实际工作状况对有限单元模型施加位移约束条件和荷载；

初始化纤维增强复合材料在结点位置的方向和密度；当初始化的纤维增强复合材料采用的是正交纤维增强复合材料模型时，初始角度为零，初始密度为任意值。

第二步、形成有限元刚度方程

设初始化的纤维增强复合材料在结点位置的初始角度为零，初始密度为任意值，设迭代指标为i＝0，则

式(1)中上角标0表示初始迭代指标，下角标j表示结点编号，J为结点总数，b表示纤维增强复合材料的两个主轴编号；

根据当前结构纤维增强复合材料的分布形成结点j位置的弹性矩阵，该弹性矩阵如下：

式(2)中E是弹性模量，上角标i为迭代指标，和是材料当前迭代下在结点位置的方向和密度，s_br是常数矩阵的分量，且

g_r(α_j)是材料方向的函数矩阵的分量，且

g(α)＝[cos2α sin2α 1](4)

A_r是常数矩阵，且

由结点位置的弹性矩阵式，通过形函数插值得到有限单元内任意位置的弹性矩阵如下：

式(6)中ξ,η为单元内部任意位置的局部坐标，N_j(ξ,η)是形函数；当有限单元采用四边形单元时，形函数为

N_j(ξ,η)＝(1+ξ_jξ)(1+η_jη)/4，j＝1,2,3,4(7)

式(7)中ξ_j和η_j是结点在局部坐标系下的坐标值；

当有限单元采用三角形单元时，则使用面积坐标作为形函数。

由弹性矩阵式可以计算出单元刚度矩阵

式(8)中H_ejr是与设计变量无关的常数矩阵；若划分所述有限单元采用规则单元网格，则H_ejr与有限单元无关，且通过公式(9)计算；

由单元刚度矩阵式按照结点对应位置累加得到结构刚度矩阵，

最终形成刚度方程，

KⁱUⁱ＝F(11)

式中Uⁱ和F分别是当前结构有限元方法的结点位移列向量和结点力列向量；

第三步、求解有限单元刚度方程及应力、应变

依照下式求解结构刚度方程得到结构在结点位置的位移

Uⁱ＝(Kⁱ)^-1F(12)

并据此计算出结构单元内的应力和应变

式(13)中εⁱ是应变向量，和分别是纤维增强复合材料沿x和y坐标轴方向的正应变和剪应变分量，B和是几何矩阵和单元结点位移向量，根据结构单元内的应变可以计算出结点位置的应变，

式中表示结点位置的应变，εⁱ表示围绕结点j的各单元在结点位置的应变，n_j表示围绕结点j的单元数量，S_j为围绕结点j的单元集合，由线弹性假设的胡克定律可以计算出应力

和分别是纤维增强复合材料在结点j位置沿x和y坐标轴方向的正应力和剪应力分量，由应力分量可以计算出主应力的方向

以及沿主应力方向的应变值，

第四步：优化材料的方向和密度

将纤维增强复合材料的两个主方向调整到主应力方向，

根据满应力准则的应力比公式

调整纤维增强复合材料在结点位置的密度值，式(19)中是结点j位置在第i次迭代时的应变，ε_p是材料的允许最大应变；

在调整纤维增强复合材料在结点位置的密度值时，通过式(20)来限制材料在结点位置的密度不要过低

式中是本次迭代得到材料所有结点位置的密度最大值

第五步：迭代收敛性检查

通过式(22)比较纤维增强复合材料在结点位置的方向和密度在两次连续迭代的相对变化量，当变化量小于一个事先给定的数值时，优化迭代结束，否则迭代指标增加1，i＝i+1，并返回第二步重新计算

式(22)中δ为事先指定的一个很小的值；

第六步：优化结果的后处理

采用有限元的形函数，根据纤维增强复合材料在结点位置的方向和密度插值得到各单元内部任意位置的纤维增强复合材料分布，如下：

得到的结构材料体积为

式(24)中

对于规则矩形单元，

式中V_e是单元体积。

因此，本发明的非均匀纤维增强复合材料分布优化设计方法解决了以下问题：

(a)纤维增强复合材料模型的构造方法以及数学表达式；

(b)有限单元计算的刚度矩阵及其密度数学表达式；

(c)非均匀各向异性纤维增强复合连续分布场的描述方法；

(d)纤维增强复合材料优化方法；

(e)纤维增强复合材料优化分布的可视化方法。

以下是本发明用于悬臂矩形梁的拓扑优化设计的应用实例：

一个长1.6米，高1米，厚0.01米的矩形设计域如图2所示。左边固定，右边中点有一个100kN竖直向下的集中力作用。弹性模量E＝210GPa，允许应力σ_p＝160MPa。优化设计纤维增强复合材料分布。最后的材料优化分布可视化图形为图3所示。

优化设计步骤如下：

1.将设计域划分为2032×640个四结点矩形单元，左边所有结点施加位移约束，右边中间结点施加向下的结点力；

2.形成结构刚度矩阵，并进行有限元分析，得到结点位移列向量，据此，逐个单元计算结点位置应变分量并累加后取平均值得到结点位置应变；

3.根据材料在结点位置的分布计算弹性矩阵，由结点位置的应变和弹性矩阵计算得到结点位置的主应力方向和主应力方向的应变；

4.根据满应力准则调整材料在结点位置的方向和角度；

5.验证收敛条件，如果未收敛则返回2，否则进入下一步；

6.计算结构的材料体积，并插值得到材料在任意位置的方向和密度；

7.优化结果可视化输出。

上述说明示出并描述了本发明的优选实施例，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。