CN107357974A - 非均匀纤维增强复合材料分布优化设计方法 - Google Patents
非均匀纤维增强复合材料分布优化设计方法 Download PDFInfo
- Publication number
- CN107357974A CN107357974A CN201710494613.9A CN201710494613A CN107357974A CN 107357974 A CN107357974 A CN 107357974A CN 201710494613 A CN201710494613 A CN 201710494613A CN 107357974 A CN107357974 A CN 107357974A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msubsup
- msub
- mtd
- formula
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 239000002131 composite material Substances 0.000 title claims abstract description 86
- 239000000835 fiber Substances 0.000 title claims abstract description 85
- 238000013461 design Methods 0.000 title claims abstract description 58
- 238000000034 method Methods 0.000 title claims abstract description 45
- 238000005457 optimization Methods 0.000 title claims abstract description 45
- 238000009826 distribution Methods 0.000 title claims abstract description 41
- 239000000463 material Substances 0.000 claims abstract description 52
- 230000002708 enhancing effect Effects 0.000 claims abstract description 7
- 238000007689 inspection Methods 0.000 claims abstract description 4
- 238000012805 post-processing Methods 0.000 claims abstract description 4
- 239000011159 matrix material Substances 0.000 claims description 49
- 238000006073 displacement reaction Methods 0.000 claims description 16
- 238000010276 construction Methods 0.000 claims description 10
- 230000001186 cumulative effect Effects 0.000 claims description 5
- 101100234408 Danio rerio kif7 gene Proteins 0.000 claims description 3
- 101100221620 Drosophila melanogaster cos gene Proteins 0.000 claims description 3
- 230000005483 Hooke's law Effects 0.000 claims description 3
- 101100398237 Xenopus tropicalis kif11 gene Proteins 0.000 claims description 3
- 230000001788 irregular Effects 0.000 claims description 3
- 238000010008 shearing Methods 0.000 claims description 3
- 238000005516 engineering process Methods 0.000 description 3
- 229920002430 Fibre-reinforced plastic Polymers 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 2
- 239000003733 fiber-reinforced composite Substances 0.000 description 2
- 239000011151 fibre-reinforced plastic Substances 0.000 description 2
- 238000000265 homogenisation Methods 0.000 description 2
- 230000001965 increasing effect Effects 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000000465 moulding Methods 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 229920000049 Carbon (fiber) Polymers 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 229920006231 aramid fiber Polymers 0.000 description 1
- 239000004917 carbon fiber Substances 0.000 description 1
- 239000012141 concentrate Substances 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007717 exclusion Effects 0.000 description 1
- 239000003365 glass fiber Substances 0.000 description 1
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000012783 reinforcing fiber Substances 0.000 description 1
- 239000012779 reinforcing material Substances 0.000 description 1
- 238000012916 structural analysis Methods 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/04—Constraint-based CAD
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/06—Power analysis or power optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Nitrogen And Oxygen Or Sulfur-Condensed Heterocyclic Ring Systems (AREA)
- Moulding By Coating Moulds (AREA)
Abstract
本发明公开非均匀纤维增强复合材料分布优化设计方法，包括步骤：第一步、设计问题初始化；第二步、形成有限元刚度方程；第三步、求解有限单元刚度方程及应力、应变；第四步：优化材料的方向和密度；第五步：迭代收敛性检查；第六步：优化结果的后处理。本发明具有较高的优化计算效率，适合于任意不均匀分布的纤维增强复合材料优化设计方法。
Description
技术领域
本发明属于非均匀各向异性材料优化设计技术领域，特别一种涉及非均匀纤维增强复合材料分布优化设计方法
背景技术
在航天航空、材料、机械、土木以及船舶和水利等领域经常需要各种复杂的纤维增强复合材料设计。纤维增强复合材料(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为迭代指标，和是材料当前迭代下在结点位置的方向和密度，sbr是常数矩阵的分量，且
gr(αj)是材料方向的函数矩阵的分量，且
g(α)＝[cos2α sin2α 1] (4)
Ar是常数矩阵，且
由结点位置的弹性矩阵式，通过形函数插值得到有限单元内任意点位置的弹性矩阵如下：
式(6)中ξ,η为单元内部任意位置的局部坐标，Nj(ξ,η)是形函数；
由弹性矩阵式可以计算出单元刚度矩阵
式(8)中Hejr是与设计变量无关的常数矩阵；
由单元刚度矩阵式按照结点对应位置累加得到结构刚度矩阵，
最终形成刚度方程，
KiUi＝F (11)
式中Ui和F分别是当前结构有限元方法的结点位移列向量和结点力列向量；
第三步：
依照下式求解结构刚度方程得到结构在结点位置的位移
Ui＝(Ki)-1F (12)
并据此计算出结构单元内的应力和应变
式(13)中εi是应变向量，和分别是纤维增强复合材料沿x和y坐标轴方向的正应变和剪应变分量，B和是几何矩阵和单元结点位移向量，根据结构单元内的应变可以计算出结点位置的应变，
式中表示结点位置的应变，εi表示围绕结点j的各单元在结点位置的应变，nj表示围绕结点j的单元数量，Sj为围绕结点j的单元集合，由线弹性假设的胡克定律可以计算出应力
和分别是纤维增强复合材料在结点j位置沿x和y坐标轴方向的正应力和剪应力分量，由应力分量可以计算出主应力的方向
以及沿主应力方向的应变值，
第四步：
将纤维增强复合材料的两个主方向调整到主应力方向，
根据满应力准则的应力比公式
调整纤维增强复合材料在结点位置的密度值，式(19)中是结点j位置在第i次迭代时的应变，εp是材料的允许最大应变；
第五步：
通过式(22)比较纤维增强复合材料在结点位置的方向和密度在两次连续迭代的相对变化量，当变化量小于一个事先给定的数值时，优化迭代结束，否则迭代指标增加1，i＝i+1，并返回第二步重新计算
式(22)中δ为事先指定的一个很小的值；
第六步：
采用有限元的形函数，根据纤维增强复合材料在结点位置的方向和密度插值得到各单元内部任意位置的纤维增强复合材料分布，如下：
得到的结构材料体积为
式(24)中
对于规则矩形单元，
式中Ve是单元体积。
所述第二步中，所述当有限单元采用四边形单元时，形函数为
Nj(ξ,η)＝(1+ξjξ)(1+ηjη)/4，j＝1,2,3,4 (7)
式(7)中ξj和ηj是结点在局部坐标系下的坐标值；
当有限单元采用三角形单元时，则使用面积坐标作为形函数。
若划分所述有限单元采用规则单元网格，则Hejr与有限单元无关，且通过公式(9)计算；
所述第四步中，在调整纤维增强复合材料在结点位置的密度值时，通过式(20)来限制材料在结点位置的密度不要过低
式中是本次迭代得到材料所有结点位置的密度最大值
采用上述方案后，本发明有益效果是：本发明基于纤维增强复合材料模型，在满足应力约束的条件下实现结构的材料体积最小化。本方法首先假设在设计域内均匀连续分布一种纤维增强复合材料，采用有限元方法进行结构分析，得到应力，应变分布场。根据结点位置的主应力方向以及主应力方向的应变大小，采用满应力准则的应力比公式优化纤维增强复合材料在结点位置的分布方向和密度。由形函数插值得到任意位置材料分布。经过反复迭代，直到收敛。最后形成的纤维增强复合材料优化分布场用材料在结点位置的分布方向和密度表示，不需要分层。
本发明具有较高的优化计算效率，而且，相比常规的复合增强材料优化设计方法仅能够设计几种分层结构的优化分布，每层的内部完全一致，不能优化，本发明则并没有采用分层设计，不受分层限制，材料允许任意分布，具有更大的优化设计空间。
下面结合附图对本发明做进一步的说明。
附图说明
附图1是本发明非均匀纤维增强复合材料分布优化方法的流程图；
附图2是本发明一个矩形梁结构的初始设计域；
附图3是采用本发明方法设计获得的最佳材料分布。
具体实施方式
如图1所示本实施例揭示的非均匀纤维增强复合材料分布优化设计方法，具体包括以下步骤：
第一步、设计问题初始化
根据实际使用要求设定一个设计域，采用有限元方法将设计域划分有限单元，形成有限单元模型，采用有限元方法将设计域划分有限单元模型具体是：对于规则区域选用矩形规则单元进行划分；对于不规则区域，则采用三角形单元、四结点四边形单元或八结点四边形单元进行划分。再根据实际工作状况对有限单元模型施加位移约束条件和荷载；
初始化纤维增强复合材料在结点位置的方向和密度；当初始化的纤维增强复合材料采用的是正交纤维增强复合材料模型时，初始角度为零，初始密度为任意值。
第二步、形成有限元刚度方程
设初始化的纤维增强复合材料在结点位置的初始角度为零，初始密度为任意值，设迭代指标为i＝0，则
式(1)中上角标0表示初始迭代指标，下角标j表示结点编号，J为结点总数，b表示纤维增强复合材料的两个主轴编号；
根据当前结构纤维增强复合材料的分布形成结点j位置的弹性矩阵，该弹性矩阵如下：
式(2)中E是弹性模量，上角标i为迭代指标，和是材料当前迭代下在结点位置的方向和密度，sbr是常数矩阵的分量，且
gr(αj)是材料方向的函数矩阵的分量，且
g(α)＝[cos2α sin2α 1] (4)
Ar是常数矩阵，且
由结点位置的弹性矩阵式，通过形函数插值得到有限单元内任意位置的弹性矩阵如下：
式(6)中ξ,η为单元内部任意位置的局部坐标，Nj(ξ,η)是形函数；当有限单元采用四边形单元时，形函数为
Nj(ξ,η)＝(1+ξjξ)(1+ηjη)/4，j＝1,2,3,4 (7)
式(7)中ξj和ηj是结点在局部坐标系下的坐标值；
当有限单元采用三角形单元时，则使用面积坐标作为形函数。
由弹性矩阵式可以计算出单元刚度矩阵
式(8)中Hejr是与设计变量无关的常数矩阵；若划分所述有限单元采用规则单元网格，则Hejr与有限单元无关，且通过公式(9)计算；
由单元刚度矩阵式按照结点对应位置累加得到结构刚度矩阵，
最终形成刚度方程，
KiUi＝F (11)
式中Ui和F分别是当前结构有限元方法的结点位移列向量和结点力列向量；
第三步、求解有限单元刚度方程及应力、应变
依照下式求解结构刚度方程得到结构在结点位置的位移
Ui＝(Ki)-1F (12)
并据此计算出结构单元内的应力和应变
式(13)中εi是应变向量，和分别是纤维增强复合材料沿x和y坐标轴方向的正应变和剪应变分量，B和是几何矩阵和单元结点位移向量，根据结构单元内的应变可以计算出结点位置的应变，
式中表示结点位置的应变，εi表示围绕结点j的各单元在结点位置的应变，nj表示围绕结点j的单元数量，Sj为围绕结点j的单元集合，由线弹性假设的胡克定律可以计算出应力
和分别是纤维增强复合材料在结点j位置沿x和y坐标轴方向的正应力和剪应力分量，由应力分量可以计算出主应力的方向
以及沿主应力方向的应变值，
第四步：优化材料的方向和密度
将纤维增强复合材料的两个主方向调整到主应力方向，
根据满应力准则的应力比公式
调整纤维增强复合材料在结点位置的密度值，式(19)中是结点j位置在第i次迭代时的应变，εp是材料的允许最大应变；
在调整纤维增强复合材料在结点位置的密度值时，通过式(20)来限制材料在结点位置的密度不要过低
式中是本次迭代得到材料所有结点位置的密度最大值
第五步：迭代收敛性检查
通过式(22)比较纤维增强复合材料在结点位置的方向和密度在两次连续迭代的相对变化量，当变化量小于一个事先给定的数值时，优化迭代结束，否则迭代指标增加1，i＝i+1，并返回第二步重新计算
式(22)中δ为事先指定的一个很小的值；
第六步：优化结果的后处理
采用有限元的形函数，根据纤维增强复合材料在结点位置的方向和密度插值得到各单元内部任意位置的纤维增强复合材料分布，如下：
得到的结构材料体积为
式(24)中
对于规则矩形单元，
式中Ve是单元体积。
因此，本发明的非均匀纤维增强复合材料分布优化设计方法解决了以下问题：
(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.优化结果可视化输出。
上述说明示出并描述了本发明的优选实施例，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。
Claims (7)
1.非均匀纤维增强复合材料分布优化设计方法，其特征在于，包括以下步骤：
第一步、设计问题初始化
根据实际使用要求设定一个设计域，采用有限元方法将设计域划分有限单元，建立有限元模型，再根据实际工作状况对有限单元模型施加位移约束条件和荷载；
初始化纤维增强复合材料在结点位置的方向和密度；
第二步、形成有限元刚度方程
根据当前纤维增强复合材料的分布形成结点位置的弹性矩阵，由结点位置的弹性矩阵，通过形函数插值得到有限单元内任意位置的弹性矩阵，再由弹性矩阵计算出有限单元刚度矩阵，接着由有限单元刚度矩阵按照结点对应位置累加得到结构刚度矩阵，最后由结构刚度矩阵得到有限单元刚度方程；
第三步、求解有限单元刚度方程及应力、应变
求解有限单元刚度方程，得到结构在结点位置的位移，并据此计算出结构单元内的应力和应变，再根据结构单元内的应变计算出结点位置的应变；由应力分量计算出主应力的方向以及主应力的方向的应变值；
第四步：优化材料的方向和密度
将纤维增强复合材料的两个主方向调整到主应力方向；根据满应力准则的应力比公式调整纤维增强复合材料在结点位置的密度值；
第五步：迭代收敛性检查
比较纤维增强复合材料在结点位置的方向和密度在两次连续迭代的相对变化量，当相对变化量小于一个预设的数值时，优化迭代结束；否则返回第二步；
第六步：优化结果的后处理
采用有限单元的形函数，根据纤维增强复合材料在结点位置的方向和密度插值得到各单元内部任意位置的纤维增强复合材料分布。
2.如权利要求1所述的非均匀纤维增强复合材料分布优化设计方法，其特征在于：所述第一步中，采用有限元方法将设计域划分有限单元具体是：对于规则区域选用矩形规则单元进行划分；对于不规则区域，则采用三角形单元、四结点四边形单元或八结点四边形单元进行划分。
3.如权利要求1所述的非均匀纤维增强复合材料分布优化设计方法，其特征在于：在所述第一步中，当初始化的纤维增强复合材料采用的是正交纤维增强复合材料模型时，初始角度为零，初始密度为任意值。
4.如权利要求1所述的非均匀纤维增强复合材料分布优化设计方法，其特征在于：所述第二步至第五步具体是：
第二步：
设初始化的纤维增强复合材料在结点位置的初始角度为零，初始密度为任意值，设迭代指标为i＝0，则
<mrow>
<msubsup>
<mi>a</mi>
<mi>j</mi>
<mn>0</mn>
</msubsup>
<mo>=</mo>
<mn>0</mn>
<mo>,</mo>
<msubsup>
<mi>t</mi>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
<mn>0</mn>
</msubsup>
<mo>=</mo>
<mn>0</mn>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mo>...</mo>
<mo>,</mo>
<mi>J</mi>
<mo>;</mo>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
式(1)中上角标0表示初始迭代指标，下角标j表示结点编号，J为结点总数，b表示纤维增强复合材料的两个主轴编号；
根据当前结构中纤维增强复合材料的分布形成结点j位置的弹性矩阵，该弹性矩阵如下：
<mrow>
<msubsup>
<mi>D</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>=</mo>
<mi>D</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>t</mi>
<mrow>
<mn>1</mn>
<mi>j</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>,</mo>
<msubsup>
<mi>t</mi>
<mrow>
<mn>2</mn>
<mi>j</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>,</mo>
<msubsup>
<mi>&alpha;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>E</mi>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mn>2</mn>
</munderover>
<msubsup>
<mi>t</mi>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
<mi>i</mi>
</msubsup>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>r</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mn>3</mn>
</munderover>
<msub>
<mi>s</mi>
<mrow>
<mi>b</mi>
<mi>r</mi>
</mrow>
</msub>
<msub>
<mi>g</mi>
<mi>r</mi>
</msub>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&alpha;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>)</mo>
</mrow>
<msub>
<mi>A</mi>
<mi>r</mi>
</msub>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>2</mn>
<mo>)</mo>
</mrow>
</mrow>
式(2)中E是弹性模量，上角标i为迭代指标，和是材料当前迭代下在结点位置的方向和密度，sbr是常数矩阵的分量，且
<mrow>
<mi>s</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>3</mn>
<mo>)</mo>
</mrow>
</mrow>
gr(αj)是材料方向的函数矩阵的分量，且
g(α)＝[cos2α sin2α 1] (4)
Ar是常数矩阵，且
<mrow>
<msub>
<mi>A</mi>
<mn>1</mn>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<mi>d</mi>
<mi>i</mi>
<mi>a</mi>
<mi>g</mi>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
<msub>
<mi>A</mi>
<mn>2</mn>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mn>4</mn>
</mfrac>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
<msub>
<mi>A</mi>
<mn>3</mn>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<mi>d</mi>
<mi>i</mi>
<mi>a</mi>
<mi>g</mi>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>5</mn>
<mo>)</mo>
</mrow>
</mrow>
由结点位置的弹性矩阵式，通过形函数插值得到有限单元内任意点位置的弹性矩阵如下：
<mrow>
<msub>
<mi>D</mi>
<mi>e</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>&xi;</mi>
<mo>,</mo>
<mi>&eta;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>&Element;</mo>
<msub>
<mi>S</mi>
<mi>e</mi>
</msub>
</mrow>
</munder>
<msub>
<mi>N</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>&xi;</mi>
<mo>,</mo>
<mi>&eta;</mi>
<mo>)</mo>
</mrow>
<msub>
<mi>D</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>6</mn>
<mo>)</mo>
</mrow>
</mrow>
式(6)中ξ,η为单元内部任意位置的局部坐标，Nj(ξ,η)是形函数；
由弹性矩阵式可以计算出单元刚度矩阵
<mrow>
<msubsup>
<mi>k</mi>
<mi>e</mi>
<mi>i</mi>
</msubsup>
<mo>=</mo>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>&Element;</mo>
<msub>
<mi>S</mi>
<mi>e</mi>
</msub>
</mrow>
</munder>
<munder>
<mo>&Sigma;</mo>
<mi>b</mi>
</munder>
<msubsup>
<mi>t</mi>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
<mi>i</mi>
</msubsup>
<munder>
<mo>&Sigma;</mo>
<mi>r</mi>
</munder>
<msub>
<mi>s</mi>
<mrow>
<mi>b</mi>
<mi>r</mi>
</mrow>
</msub>
<msub>
<mi>g</mi>
<mi>r</mi>
</msub>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&alpha;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>)</mo>
</mrow>
<msub>
<mi>H</mi>
<mrow>
<mi>e</mi>
<mi>j</mi>
<mi>r</mi>
</mrow>
</msub>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>8</mn>
<mo>)</mo>
</mrow>
</mrow>
式(8)中Hejr是与设计变量无关的常数矩阵；
由单元刚度矩阵式按照结点对应位置累加得到结构刚度矩阵，
<mrow>
<msup>
<mi>K</mi>
<mi>i</mi>
</msup>
<mo>=</mo>
<munder>
<mo>&Sigma;</mo>
<mi>e</mi>
</munder>
<msubsup>
<mi>k</mi>
<mi>e</mi>
<mi>i</mi>
</msubsup>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>10</mn>
<mo>)</mo>
</mrow>
</mrow>
最终形成刚度方程，
KiUi＝F (11)
式中Ui和F分别是当前结构有限元方法的结点位移列向量和结点力列向量；
第三步：
依照下式求解结构刚度方程得到结构在结点位置的位移
Ui＝(Ki)-1F (12)
并据此计算出结构单元内的应力和应变
<mrow>
<msup>
<mi>&epsiv;</mi>
<mi>i</mi>
</msup>
<mo>=</mo>
<msup>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msubsup>
<mi>&epsiv;</mi>
<mi>x</mi>
<mi>i</mi>
</msubsup>
</mtd>
<mtd>
<msubsup>
<mi>&epsiv;</mi>
<mi>y</mi>
<mi>i</mi>
</msubsup>
</mtd>
<mtd>
<msubsup>
<mi>&gamma;</mi>
<mrow>
<mi>x</mi>
<mi>y</mi>
</mrow>
<mi>i</mi>
</msubsup>
</mtd>
</mtr>
</mtable>
</mfenced>
<mi>T</mi>
</msup>
<mo>=</mo>
<msubsup>
<mi>BU</mi>
<mi>e</mi>
<mi>i</mi>
</msubsup>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>13</mn>
<mo>)</mo>
</mrow>
</mrow>
式(13)中εi是应变向量，和分别是纤维增强复合材料沿x和y坐标轴方向的正应变和剪应变分量，B和是几何矩阵和单元结点位移向量，根据结构单元内的应变可以计算出结点位置的应变，
<mrow>
<msubsup>
<mi>&epsiv;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>=</mo>
<msup>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msubsup>
<mi>&epsiv;</mi>
<mrow>
<mi>j</mi>
<mi>x</mi>
</mrow>
<mi>i</mi>
</msubsup>
</mtd>
<mtd>
<msubsup>
<mi>&epsiv;</mi>
<mrow>
<mi>j</mi>
<mi>y</mi>
</mrow>
<mi>i</mi>
</msubsup>
</mtd>
<mtd>
<msubsup>
<mi>&gamma;</mi>
<mrow>
<mi>j</mi>
<mi>x</mi>
<mi>y</mi>
</mrow>
<mi>i</mi>
</msubsup>
</mtd>
</mtr>
</mtable>
</mfenced>
<mi>T</mi>
</msup>
<mo>=</mo>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mi>e</mi>
<mo>&Element;</mo>
<msub>
<mi>S</mi>
<mi>j</mi>
</msub>
</mrow>
</munder>
<msup>
<mi>&epsiv;</mi>
<mi>i</mi>
</msup>
<mo>/</mo>
<msub>
<mi>n</mi>
<mi>j</mi>
</msub>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mi>J</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>14</mn>
<mo>)</mo>
</mrow>
</mrow>
式中表示结点位置的应变，εi表示围绕结点j的各单元在结点位置的应变，nj表示围绕结点j的单元数量，Sj为围绕结点j的单元集合，由线弹性假设的胡克定律可以计算出应力
<mrow>
<msubsup>
<mi>&sigma;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>=</mo>
<msup>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msubsup>
<mi>&sigma;</mi>
<mrow>
<mi>j</mi>
<mi>x</mi>
</mrow>
<mi>i</mi>
</msubsup>
</mtd>
<mtd>
<msubsup>
<mi>&sigma;</mi>
<mrow>
<mi>j</mi>
<mi>y</mi>
</mrow>
<mi>i</mi>
</msubsup>
</mtd>
<mtd>
<msubsup>
<mi>&tau;</mi>
<mrow>
<mi>j</mi>
<mi>x</mi>
<mi>y</mi>
</mrow>
<mi>i</mi>
</msubsup>
</mtd>
</mtr>
</mtable>
</mfenced>
<mi>T</mi>
</msup>
<mo>=</mo>
<mi>D</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>t</mi>
<mrow>
<mn>1</mn>
<mi>j</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>,</mo>
<msubsup>
<mi>t</mi>
<mrow>
<mn>2</mn>
<mi>j</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>,</mo>
<msubsup>
<mi>&alpha;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>)</mo>
</mrow>
<msubsup>
<mi>&epsiv;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mi>J</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>15</mn>
<mo>)</mo>
</mrow>
</mrow>
和分别是纤维增强复合材料在结点j位置沿x和y坐标轴方向的正应力和剪应力分量，由应力分量可以计算出主应力的方向
<mrow>
<msubsup>
<mi>&theta;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<msup>
<mi>tan</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mfrac>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&tau;</mi>
<mrow>
<mi>j</mi>
<mi>x</mi>
<mi>y</mi>
</mrow>
<mi>i</mi>
</msubsup>
</mrow>
<mrow>
<msubsup>
<mi>&sigma;</mi>
<mrow>
<mi>j</mi>
<mi>x</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&sigma;</mi>
<mrow>
<mi>j</mi>
<mi>y</mi>
</mrow>
<mi>i</mi>
</msubsup>
</mrow>
</mfrac>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mi>J</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>16</mn>
<mo>)</mo>
</mrow>
</mrow>
以及沿主应力方向的应变值，
<mrow>
<msubsup>
<mi>&epsiv;</mi>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&epsiv;</mi>
<mrow>
<mi>j</mi>
<mi>x</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>&epsiv;</mi>
<mrow>
<mi>j</mi>
<mi>y</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mo>&PlusMinus;</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<mo>&lsqb;</mo>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&epsiv;</mi>
<mrow>
<mi>j</mi>
<mi>x</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&epsiv;</mi>
<mrow>
<mi>j</mi>
<mi>y</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mn>2</mn>
<msubsup>
<mi>&theta;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>&tau;</mi>
<mrow>
<mi>j</mi>
<mi>x</mi>
<mi>y</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mi>s</mi>
<mi>i</mi>
<mi>n</mi>
<mn>2</mn>
<msubsup>
<mi>&theta;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>&rsqb;</mo>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mi>J</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>17</mn>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
第四步：
将纤维增强复合材料的两个主方向调整到主应力方向，
<mrow>
<msubsup>
<mi>&alpha;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>&theta;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mi>J</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>18</mn>
<mo>)</mo>
</mrow>
</mrow>
根据满应力准则的应力比公式
<mrow>
<msubsup>
<mover>
<mi>t</mi>
<mo>&OverBar;</mo>
</mover>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mo>=</mo>
<mo>|</mo>
<msubsup>
<mi>&epsiv;</mi>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>|</mo>
<msubsup>
<mi>t</mi>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
<mi>i</mi>
</msubsup>
<mo>/</mo>
<msub>
<mi>&epsiv;</mi>
<mi>p</mi>
</msub>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mi>J</mi>
<mo>;</mo>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>19</mn>
<mo>)</mo>
</mrow>
</mrow>
调整纤维增强复合材料在结点位置的密度值，式(19)中是结点j位置在第i次迭代时的应变，εp是材料的允许最大应变；
第五步：
通过式(22)比较纤维增强复合材料在结点位置的方向和密度在两次连续迭代的相对变化量，当变化量小于一个事先给定的数值时，优化迭代结束，否则迭代指标增加1，i＝i+1，并返回第二步重新计算
<mrow>
<mo>|</mo>
<mn>1</mn>
<mo>-</mo>
<msubsup>
<mi>t</mi>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mo>/</mo>
<msubsup>
<mi>t</mi>
<mi>m</mi>
<mi>i</mi>
</msubsup>
<mo>|</mo>
<mo>&le;</mo>
<mi>&delta;</mi>
<mo>,</mo>
<mo>|</mo>
<msubsup>
<mi>&alpha;</mi>
<mi>j</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&alpha;</mi>
<mi>j</mi>
<mi>i</mi>
</msubsup>
<mo>|</mo>
<mo>&le;</mo>
<mi>&delta;</mi>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mi>J</mi>
<mo>;</mo>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>22</mn>
<mo>)</mo>
</mrow>
</mrow>
式(22)中δ为事先指定的一个很小的值；
第六步：
采用有限元的形函数，根据纤维增强复合材料在结点位置的方向和密度插值得到各单元内部任意位置的纤维增强复合材料分布，如下：
<mrow>
<mi>&alpha;</mi>
<mrow>
<mo>(</mo>
<mi>&xi;</mi>
<mo>,</mo>
<mi>&eta;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>&Element;</mo>
<msub>
<mi>S</mi>
<mi>e</mi>
</msub>
</mrow>
</munder>
<msub>
<mi>N</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>&xi;</mi>
<mo>,</mo>
<mi>&eta;</mi>
<mo>)</mo>
</mrow>
<msub>
<mi>&alpha;</mi>
<mi>j</mi>
</msub>
</mrow>
<mrow>
<msub>
<mi>t</mi>
<mi>b</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>&xi;</mi>
<mo>,</mo>
<mi>&eta;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>&Element;</mo>
<msub>
<mi>S</mi>
<mi>e</mi>
</msub>
</mrow>
</munder>
<msub>
<mi>N</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>&xi;</mi>
<mo>,</mo>
<mi>&eta;</mi>
<mo>)</mo>
</mrow>
<msub>
<mi>t</mi>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
</msub>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mi>J</mi>
<mo>;</mo>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>23</mn>
<mo>)</mo>
</mrow>
</mrow>
得到的结构材料体积为
<mrow>
<mi>V</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>V</mi>
<mi>e</mi>
</msub>
<munder>
<mo>&Sigma;</mo>
<mi>j</mi>
</munder>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<munder>
<mo>&Sigma;</mo>
<mi>b</mi>
</munder>
<msub>
<mi>t</mi>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
</msub>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>24</mn>
<mo>)</mo>
</mrow>
</mrow>
3
式(24)中
<mrow>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<msub>
<mi>V</mi>
<mi>e</mi>
</msub>
</mfrac>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mi>e</mi>
<mo>&Element;</mo>
<msub>
<mi>S</mi>
<mi>j</mi>
</msub>
</mrow>
</munder>
<msub>
<mo>&Integral;</mo>
<msub>
<mi>V</mi>
<mi>e</mi>
</msub>
</msub>
<msub>
<mi>N</mi>
<mi>j</mi>
</msub>
<mi>d</mi>
<mi>V</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>25</mn>
<mo>)</mo>
</mrow>
</mrow>
对于规则矩形单元，
式中Ve是单元体积。
5.如权利要求4所述的非均匀纤维增强复合材料分布优化设计方法，其特征在于：所述第二步中，所述当有限单元采用四边形单元时，形函数为
Nj(ξ,η)＝(1+ξjξ)(1+ηjη)/4，j＝1,2,3,4 (7)
式(7)中ξj和ηj是结点在局部坐标系下的坐标值；
当有限单元采用三角形单元时，则使用面积坐标作为形函数。
6.如权利要求4所述的非均匀纤维增强复合材料分布优化设计方法，其特征在于：所述第二步中，若划分所述有限单元采用规则单元网格，则Hejr与有限单元无关，且通过公式(9)计算；
<mrow>
<msub>
<mi>H</mi>
<mrow>
<mi>e</mi>
<mi>j</mi>
<mi>r</mi>
</mrow>
</msub>
<mo>=</mo>
<mi>E</mi>
<msub>
<mo>&Integral;</mo>
<msub>
<mi>V</mi>
<mi>e</mi>
</msub>
</msub>
<msub>
<mi>N</mi>
<mi>j</mi>
</msub>
<msup>
<mi>B</mi>
<mi>T</mi>
</msup>
<msub>
<mi>A</mi>
<mi>r</mi>
</msub>
<mi>B</mi>
<mi>d</mi>
<mi>V</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>9</mn>
<mo>)</mo>
</mrow>
<mo>.</mo>
</mrow>
7.如权利要求4所述的非均匀纤维增强复合材料分布优化设计方法，其特征在于：所述第四步中，在调整纤维增强复合材料在结点位置的密度值时，通过式(20)来限制材料在结点位置的密度不要过低
<mrow>
<msubsup>
<mi>t</mi>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mo>=</mo>
<mi>m</mi>
<mi>a</mi>
<mi>x</mi>
<mrow>
<mo>(</mo>
<msup>
<mn>10</mn>
<mrow>
<mo>-</mo>
<mn>5</mn>
</mrow>
</msup>
<mo>&times;</mo>
<msubsup>
<mi>t</mi>
<mi>m</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mo>,</mo>
<msubsup>
<mover>
<mi>t</mi>
<mo>&OverBar;</mo>
</mover>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mo>)</mo>
</mrow>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mi>J</mi>
<mo>;</mo>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>20</mn>
<mo>)</mo>
</mrow>
</mrow>
式中是本次迭代得到材料所有结点位置的密度最大值
<mrow>
<msubsup>
<mi>t</mi>
<mi>m</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mo>=</mo>
<munder>
<mi>max</mi>
<mrow>
<mi>b</mi>
<mo>,</mo>
<mi>j</mi>
</mrow>
</munder>
<mrow>
<mo>(</mo>
<msubsup>
<mover>
<mi>t</mi>
<mo>&OverBar;</mo>
</mover>
<mrow>
<mi>b</mi>
<mi>j</mi>
</mrow>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>21</mn>
<mo>)</mo>
</mrow>
<mo>.</mo>
</mrow>
4
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710207878 | 2017-03-31 | ||
| CN2017102078786 | 2017-03-31 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107357974A true CN107357974A (zh) | 2017-11-17 |
| CN107357974B CN107357974B (zh) | 2020-07-31 |
Family
ID=60272485
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710494613.9A Expired - Fee Related CN107357974B (zh) | 2017-03-31 | 2017-06-26 | 非均匀纤维增强复合材料分布优化设计方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107357974B (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110096756A (zh) * | 2019-04-08 | 2019-08-06 | 河海大学 | 一种考虑荷载不确定性的自由曲面结构形态创建方法 |
| CN110955941A (zh) * | 2019-11-29 | 2020-04-03 | 华中科技大学 | 基于向量场的复合材料结构优化设计方法及设备 |
| WO2020211012A1 (zh) * | 2019-04-17 | 2020-10-22 | 大连理工大学 | 一种面向混杂纤维复合材料板壳结构的快速协同优化方法 |
| CN111832211A (zh) * | 2020-07-27 | 2020-10-27 | 内蒙古工业大学 | 一种复合纤维风力机叶片的刚度优化方法 |
| US11455438B2 (en) | 2018-02-01 | 2022-09-27 | Toyota Motor Engineering & Manufacturing North America, Inc. | Methods for topology optimization using a membership variable |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2005038211A (ja) * | 2003-07-15 | 2005-02-10 | Toyota Motor Corp | 有限要素を用いた数値解析における曲面モデル形状の分析方法 |
| US20060158449A1 (en) * | 2002-05-31 | 2006-07-20 | Ugs Plm Solutions, Inc., A Delaware Corporation | Computerized deformation analyzer |
| CN102043883A (zh) * | 2010-12-29 | 2011-05-04 | 长沙理工大学 | 基于材料破损约束的连续体结构拓扑设计建模及优化设计方法 |
| CN103279591A (zh) * | 2013-04-24 | 2013-09-04 | 西北工业大学 | 基于附加单元的固体重量压力载荷下的拓扑优化设计方法 |
| CN104239624A (zh) * | 2014-09-05 | 2014-12-24 | 西安交通大学 | 一种机床床身内部结构优化设计方法 |
| CN105512385A (zh) * | 2015-12-04 | 2016-04-20 | 吉林大学 | 采用序列幂函数插值方法实现结构多材料拓扑优化的方法 |
-
2017
- 2017-06-26 CN CN201710494613.9A patent/CN107357974B/zh not_active Expired - Fee Related
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20060158449A1 (en) * | 2002-05-31 | 2006-07-20 | Ugs Plm Solutions, Inc., A Delaware Corporation | Computerized deformation analyzer |
| JP2005038211A (ja) * | 2003-07-15 | 2005-02-10 | Toyota Motor Corp | 有限要素を用いた数値解析における曲面モデル形状の分析方法 |
| CN102043883A (zh) * | 2010-12-29 | 2011-05-04 | 长沙理工大学 | 基于材料破损约束的连续体结构拓扑设计建模及优化设计方法 |
| CN103279591A (zh) * | 2013-04-24 | 2013-09-04 | 西北工业大学 | 基于附加单元的固体重量压力载荷下的拓扑优化设计方法 |
| CN104239624A (zh) * | 2014-09-05 | 2014-12-24 | 西安交通大学 | 一种机床床身内部结构优化设计方法 |
| CN105512385A (zh) * | 2015-12-04 | 2016-04-20 | 吉林大学 | 采用序列幂函数插值方法实现结构多材料拓扑优化的方法 |
Non-Patent Citations (5)
| Title |
|---|
| GUEST J ET AL: "Achieving minimum length scale in topology optimization using nodal design variables and projection functions", 《INTERNATIONAL JOURNAL FOR NUMERICAL》 * |
| PRAGER W ET AL: "A note on discretized Michell structures", 《COMPUTER METHODS IN APPLIED》 * |
| 李霞等: "采用类桁架连续体的桁架结构拓扑优化方法", 《华侨大学学报(自然科学版)》 * |
| 辛星等: "基于类桁架连续体框架结构位移约束的拓扑优化", 《福建建筑》 * |
| 郑伟伟等: "应用类桁架模型的连续体拓扑优化方法", 《华侨大学学报（自然科学版）》 * |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11455438B2 (en) | 2018-02-01 | 2022-09-27 | Toyota Motor Engineering & Manufacturing North America, Inc. | Methods for topology optimization using a membership variable |
| CN110096756A (zh) * | 2019-04-08 | 2019-08-06 | 河海大学 | 一种考虑荷载不确定性的自由曲面结构形态创建方法 |
| CN110096756B (zh) * | 2019-04-08 | 2020-01-17 | 河海大学 | 一种考虑荷载不确定性的自由曲面结构形态创建方法 |
| WO2020211012A1 (zh) * | 2019-04-17 | 2020-10-22 | 大连理工大学 | 一种面向混杂纤维复合材料板壳结构的快速协同优化方法 |
| CN110955941A (zh) * | 2019-11-29 | 2020-04-03 | 华中科技大学 | 基于向量场的复合材料结构优化设计方法及设备 |
| CN110955941B (zh) * | 2019-11-29 | 2022-03-18 | 华中科技大学 | 基于向量场的复合材料结构优化设计方法及设备 |
| CN111832211A (zh) * | 2020-07-27 | 2020-10-27 | 内蒙古工业大学 | 一种复合纤维风力机叶片的刚度优化方法 |
| CN111832211B (zh) * | 2020-07-27 | 2023-07-07 | 内蒙古工业大学 | 一种复合纤维风力机叶片的刚度优化方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107357974B (zh) | 2020-07-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107357974A (zh) | 非均匀纤维增强复合材料分布优化设计方法 | |
| Di Sciuva et al. | Bending, free vibration and buckling of functionally graded carbon nanotube-reinforced sandwich plates, using the extended Refined Zigzag Theory | |
| Daniel et al. | Engineering mechanics of composite materials | |
| WO2020211012A1 (zh) | 一种面向混杂纤维复合材料板壳结构的快速协同优化方法 | |
| CN102720295B (zh) | 一种基于索穹顶张拉和承载全过程分析的预应力确定方法 | |
| Mishnaevsky Jr et al. | Hybrid and hierarchical nanoreinforced polymer composites: Computational modelling of structure–properties relationships | |
| CN110083900B (zh) | 一种面向混杂纤维复合材料板壳结构的快速协同优化方法 | |
| Muc et al. | Design of plates with curved fibre format | |
| CN107451309A (zh) | 一种多尺度计算复杂复合材料结构等效热膨胀系数的方法 | |
| CN104699900A (zh) | 一种等厚度复合材料层合板的设计方法 | |
| Ton-That et al. | A combined strain element in static, frequency and buckling analyses of laminated composite plates and shells | |
| Qi-Hua et al. | Crashworthiness design of gradient bionic Al/CFRP hybrid tubes under multiple loading conditions | |
| Van Campen et al. | Design of fiber-steered variable-stiffness laminates based on a given lamination parameters distribution | |
| Rajad et al. | Fiber orientation effect on the behavior of the composite materials of the horizontal axis wind turbine blade (HAWTB) | |
| Wen et al. | Mesh-free micromechanical model for woven fabric composite elastic moduli | |
| Liu et al. | Path-driven shell lattices designed for continuous fiber composite 3D printing | |
| Weaver et al. | Optimisation of variable stiffness plates | |
| McKeown | Optimal composite structures by deflection-variable programming | |
| Fanxiao et al. | Study on structural characteristics of composite smart grille based on principal component analysis | |
| Li et al. | Two-scale and three-scale computational continua models of composite curved beams | |
| Honda et al. | Design for the maximum natural frequency of laminated composite plates by optimally distributed short fibers | |
| Marchal et al. | Fast layer fiber orientation optimization method for continuous fiber-reinforced material extrusion process | |
| Maksimović | Optimum design of composite structures | |
| Kwon | Multiscale and multilevel modeling of composites | |
| Rao et al. | Modal analysis of thin FRP skew symmetric angle-ply laminate with circular cut-out |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200731 |