基于网格参数化与B样条参数化结合的复杂曲面加筋优化方法

基于网格参数化与b样条参数化结合的复杂曲面加筋优化方法
技术领域
1.本发明涉及一种复杂曲面加筋的优化方法，特别是一种基于网格参数化与b样条参数化结合的加筋优化方法。


背景技术：

2.薄壁加筋是航空、航天、航海等领域广泛中应用的结构。其原因在于通过加筋可以提高薄壁结构的力学性能，以保证薄壁结构在复杂的工况下能保持结构外形。而在航空航天等工程领域，这类薄壁结构通常是形状极为复杂的曲面，如飞机机头、太空舱外壳、潜水艇外壁等。这种复杂的曲面结构通常难以用参数公式表达，大大增加了曲面加筋布局优化设计的难度。因此针对复杂曲面的加筋布局优化设计成为了国内外学者尤其航空航天领域的热门研究方向。
3.文献“topology optimization in b-spline space[j].qian x.computer methods in applied mechanics&engineering,2013,265(oct.1):15-35.”提出了一种基于b样条的拓扑优化方法。该方法用在伪密度法的基础上引入了b样条参数场，将传统的以单元为基础的伪密度设计变量替换为b样条控制点，以此得到可以表示整个结构拓扑构型的b样条参数场。该方法不需要额外的过滤函数，避免了传统密度法中棋盘格及锯齿状缺陷，所得到的结果光滑且具有高阶连续性，在加筋方面的应用有着突出的优势。然而该方法对设计目标有着严苛的要求，即设计域要严格限定在参数域平面上，这一要求严重限制了该方法在复杂曲面的加筋优化的应用。
[0004]
在已提交的专利“基于b样条参数化的薄壁结构的加筋建模与优化方法，(申请号：cn202011342095.7)”中，实现了将b样条参数化方法应用于薄壁结构加筋优化。然而该方法对曲面的描述采用的是参数映射的方法，即需要得到被优化的曲面的表达公式，这一缺陷限制了该方法的应用范围。对于航空航天等领域广泛存在的复杂曲面，通常难以得到一个有效的表达公式来描述。因此该方法难以应用于这类复杂曲面的加筋优化。
[0005]
文献“a local/global approach to mesh parameterization,liu l,zhang l,xu y,gotsman c,gortler s.computer graphics forum,2008,27(5):1495-1504.”提出了一种高效的网格参数化方法。网格参数化方法是计算机图形学中，将用三角网格表示的复杂曲面展开到平面参数域的一种一一映射方法。这类方法不需要三维曲面的参数表达函数，通过三维曲面上的三角网格与平面参数域的三角网格的顶点一一对应的方法，实现复杂曲面到平面之间的映射。对于不可展曲面、封闭曲面等复杂曲面，可以通过割缝法对三维曲面进行分割，然后进行展开，因此该方法对于任意曲面均可以适用。该方法可以解决b样条参数化方法的参数域和复杂曲面的物理域之间的映射问题，而b样条参数化方法的高阶连续性又可以弥补采用割缝法的网格参数化中引入的不连续问题，两种方法完美契合。


技术实现要素：

[0006]
要解决的技术问题
[0007]
传统的拓扑优化方法的设计变量依托于网格，所得到的薄壁加筋结构存在锯齿状、棋盘格问题等缺陷，为了克服这类问题通常需要额外的过滤函数。本发明提出一种基于网格参数化与b样条参数化结合的复杂曲面加筋优化方法。
[0008]
技术方案
[0009]
一种基于网格参数化与b样条参数化结合的复杂曲面加筋优化方法，其特征在于步骤如下：
[0010]
步骤1：对三维曲面进行三角网格划分，得到三维网格顶点坐标x及坐标编号t；首先判断该三维曲面是否为封闭曲面，如果是封闭曲面，则需要对三维曲面网格分割之后再进行网格参数化；
[0011]
步骤2：计算网格参数化结果的扭曲度；根据已知的三角网格的雅克比矩阵j
t
，计算其相应的最大特征值γ
t
和最小特征值γ
t
，三角网格的扭曲度可以表示为
[0012][0013]
为了避免参数域等比放大缩小对扭曲度的影响，对d
t2
做归一化处理
[0014][0015]
其中，ai和a
i*
分别表示第i个三角网格在三维曲面和参数平面的面积；d
t
＝1表示网格没有被拉伸压缩，d
t
》1表示网格被拉伸，d
t
《1表示网格被压缩；
[0016]
步骤3：对网格参数化结构进行评估，当所有d
t
均处于设定区间时，认为网格参数化结果是扭曲度较小，满足b样条参数化方法要求的；当有d
t
存在不处于设定区间时，认为局部网格变形过大，会导致b样条参数化方法精度降低；因而需要对三维网格进行切割，对切割后的网格进行展开，并再次评估扭曲度；
[0017]
步骤4：在参数域设置控制点作为设计变量设计变量的取值范围为加筋的高度许可范围[0,h
max
]，选取合适的b样条阶数p并建立b样条参数场：
[0018][0019]
其中，n
i,p
(ξ)，n
j,p
(η)表示b样条的基函数；
[0020]
步骤5：在三维曲面上，对三角网格进行法向拉伸，生成n三棱柱实体网格用于有限元分析；三维曲面上三角网格位置的加筋高度，通过将该三角网格在平面参数域对应的坐标带入b样条参数场得到hc＝h(u
t
,v
t
)；三角网格对应的一列三棱柱单元的伪密度根据加筋高度hc与单元上下边界位置h
k-1
、hk之间的关系得到，利用heaviside函数近似为
[0021]
[0022]
其中，为heaviside参数；将所得的伪密度ρk乘到单元刚度矩阵并进行组装，即可得到模型整体的刚度矩阵；在所得的有限元模型上施加载荷及边界条件，得到有限元分析模型；
[0023]
步骤6：根据步骤5所得的有限元模型，建立以结构柔顺度最小为优化目标，加筋体积为约束的优化模型：
[0024][0025]
其中，c为整体加筋结构的柔顺度；k、u、f分别为有限元模型的刚度矩阵、位移向量、载荷向量；表示体积约束；
[0026]
步骤7：对步骤6建立的优化模型进行优化求解，得到优化结果。
[0027]
本发明进一步的技术方案：步骤1利用arap方法，对三角网格进行网格参数化：
[0028][0029]
其中，φ为关于参数坐标u与三角网格相似矩阵l
t
相关的能量函数，而相似矩阵可以表示为求解过程为，首先通过简单的固定边界网格参数化方法得到初始的u0，根据x计算出相应的雅克比矩阵j
t0
；对雅克比矩阵j
t0
进行svd分解得到对应的l
t0
；根据能量方程导数为0，即可以得到公式
[0030][0031]
其中，θ
ij
为三角形ij边所对的角度，j∈g(i)表示与i点相邻的所以点的编号；通过求解方程可以得到新的参数坐标u1；重复上述步骤直到前后两步之间能量变化在一个极小范围内，则认为参数化算法收敛，得到的网格参数化坐标为最优解。
[0032]
本发明进一步的技术方案：步骤3的设定区域为[0.5,2]，具体为：对网格参数化结构进行评估，当所有dt均处于[0.5,2]区间时，认为网格参数化结果是扭曲度较小，满足b样条参数化方法要求的；当存在dt小于0.5或大于2时，认为局部网格变形过大，会导致b样条参数化方法精度降低；因而需要对三维网格进行切割，对切割后的网格进行展开，并再次评估扭曲度。
[0033]
本发明进一步的技术方案：步骤3中采用采用切割法对三维网格进行预处理，对切
割后的网格进行参数化，并再次评估参数化结果的扭曲度；若扭曲度符合许可标准，则输出网格参数化结果，若扭曲度不符合许可标准则再次使用切割法继续对网格进行处理，直至扭曲度符合标准。
[0034]
本发明进一步的技术方案：切割法分为两种形式，一种是切割展平法，一种是分段参数化法；二者的区别在于切割展平法并不将曲面完全分割，而分段参数化会将曲面分割成两段乃至更多；分段参数化方法能够得到扭曲度更小的网格参数化结果，但是同时会引入更多的切缝；因而使用中优先考虑切割展平法，当切割展平法得到的参数化结果扭曲度不满足时进而采用分段参数化方法。
[0035]
本发明进一步的技术方案：步骤7中采用利用全局收敛的移动渐近算法gcmma进行优化求解。
[0036]
有益效果
[0037]
本发明提出的一种基于网格参数化与b样条参数化结合的复杂曲面加筋优化方法，相比现有技术的有益效果是：
[0038]
1、本发明采用三角网格近似的方法描述三维曲面，对于工程应用中最常用的cad模型均可适用。无需对cad模型进行二次建模，实现工程应用中cad模型能够拿来即，减小了模型近似过程的误差。
[0039]
2、本发明采用网格参数化方法实现三维曲面到二维参数面的一一映射，避免了求解曲面的参数表达式的复杂过程。该方法允许采用切割展开的技术，对于封闭曲面、不扩展曲面，可以采用分而治之的策略将其变形最小的展平到参数平面并用于分析优化。对于工程应用中常见的拼接曲面等难以得到参数表达式的复杂模型，该方法能很好的满足使用要求，成功将b样条参数化方法的应用拓展到任意曲面。
[0040]
3、该方法采用b样条参数化方法优化加筋布局，所得到的加筋光滑且连续性好，能够避免锯齿状问题。b样条参数化方法自身具有过滤效果，能够避免传统拓扑优化中常见的棋盘格问题。所得到的优化结果无需二次设计就可以应用到加筋布局设计的指导工作。
[0041]
4、本方法通过网格参数化与b样条参数化两种方法的有机结合，保留了二者的优势，将b样条参数化方法拓展到任意三维曲面，其应用前景不仅仅局限于加筋结构优化，对于其他领域如曲面结构的孔洞形状优化均有实用意义。本方法所得的结果是基于b样条参数场，可以导出为cad的文件模型。因而本方法可以实现从cad到cae再到cad的分析优化流程的无缝衔接。
附图说明
[0042]
附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。
[0043]
图1是本发明方法采用两种割缝法解决封闭曲面问题的切缝位置及展开结果图。
[0044]
图2是本发明方法解决受扭圆柱筒加筋问题的设计域、边界条件和载荷工况。
[0045]
图3是本发明方法处理圆柱筒曲面的arap方法网格参数化结果。
[0046]
图4是本发明方法采用割缝法处理圆柱筒曲面时的割缝位置及网格参数化结果。
[0047]
图5是本发明方法处理圆柱筒问题的优化迭代曲线。
[0048]
图6是本发明方法关于圆柱筒问题优化后的控制点分布及加筋结果。
[0049]
图7是本发明方法解决飞机机头加筋问题的设计域、边界条件和工况。。
[0050]
图8是本发明方法采用arap参数化方法对飞机机头曲面网格的展平结果。
[0051]
图9是本发明方法关于飞机机头问题的优化迭代曲线。
[0052]
图10是本发明方法对飞机机头优化后的控制点分布及加筋效果图。
具体实施方式
[0053]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0054]
本发明提出了利用网格参数化对复杂曲面与b样条参数场进行一一对应，以b样条参数化对加筋高度进行布局及高度同时优化的方法。首先，本方法利用三角网格对复杂曲面结构进行离散化，得到相应的三角网格的三维节点坐标与节点编号顺序。利用网格参数化方法，将三角网格一一映射到平面参数域，得到相应的参数域下的uv坐标。网格参数化过程要避免三角网格出现翻转，网格参数化结果的好坏通过扭曲度参数进行度量。对于展开后扭曲度较大的曲面及封闭曲面，通常会采用切割缝法，将曲面剖分成容易展开的曲面进行网格参数化。参数域下的三角网格与三维曲面的三角网格通过节点编号顺序实现一一对应，这样就得到了复杂三维曲面在二维参数平面的映射结果。然后，在三维曲面的网格基础上，通过网格拉伸的方法，得到三棱柱实体单元，实体单元的总高度为加筋的为最大许可高度h
max
。在平面参数域上，设置控制点作为设计变量，建立b样条参数场，并保证b样条参数场包含网格展开后的二维参数模型。其中，b样条参数场中设计变量的取值范围为0-h
max
，b样条的参数场即表示加筋的高度。在b样条参数化方法中，合理的选取控制点数、网格数及b样条阶数，可以控制过滤半径的大小，避免锯齿状及棋盘格问题。三棱柱实体单元的伪密度由b样条参数场所代表的加筋高度与单元的位置关系决定。三维曲面上三角网格所对应的加筋高度与可以通过参数域网格相应的三角网格在b样条参数场中的位置得到。最后，考虑到不同网格参数化方法得到的曲面展开结果不同，而网格参数化结果的好坏又会影响到b样条参数化方法的优化结果，本方法制定了一套基于网格扭曲度的使用准则，通过网格参数化后扭曲度的大小判断该展平结果是否满足b样条参数化方法的要求。如果扭曲度超过某一标准，则需要用割缝法对三维曲面进行切割，将切割后的网格参数化结果的扭曲度再一次与标准值进行比较，直至网格参数化结果的扭曲度满足b样条参数化要求为止。通过以上的网格参数化与b样条参数化方法的结合，可以实现复杂曲面的b样条参数化加筋结构优化设计，不需要额外过滤及光顺便可以得到光滑且高阶连续的、高度渐变的三维曲面上的加筋布局。
[0055]
为了使本领域技术人员更好地理解本发明，下面结合具体实施例对本发明进行详细说明。
[0056]
示例一：本发明用于受扭载荷下的圆柱筒结构加筋优化设计。参考图2，圆柱筒尺寸为半径12mm长度60mm，许可加筋高度h
max
为0.5mm，圆柱筒底端固定，上端施加两对对称点载荷模拟扭转载荷，点载荷大小为f＝1n。材料的杨氏模量和泊松比分别为e＝2.1
×
105pa,ν＝0.3。
[0057]
步骤1、对圆柱筒曲面用三角网格进行划分，对所得的三角曲面网格进行网格参数化。首先采用arap网格参数化算法，得到圆柱筒展开结果如图3。
[0058]
步骤2、计算展平后网格的扭曲度，得d
t
的最大值最小值分别为7.34、0.0013，明显超出给定的许可范围。因而三维网格需要进一步的切割展平，采用割缝法后，曲面的割缝位置及的网格参数化结果如图4。
[0059]
步骤3、对圆柱筒切割后进行网格参数化，切割位置及网格参数化结果如图3，计算网格扭曲度得到，d
t
最大值最小值均趋近于1，即网格几乎没有被拉伸压缩。因此该网格参数化后的结果可以应用于b样条参数化方法。
[0060]
步骤4、选取b样条阶数为p＝q＝3，在平面参数域上选取100
×
100个b样条控制点建立b样条参数场，并保证b样条参数场能完全包含网格参数化结果，设计变量取值范围为
[0061]
步骤5、在三维曲面上，三角形网格进行法向拉伸，得到总高度为0.5mm的5层三棱柱实体单元。每一个三角网格对应的加筋高度通过该三角网格三个顶点在b样条参数场的值取平均值得到h
tc
＝(h
t1
h
c2
h
c3
)/3。每一列三棱柱各层的实体单元通过加筋高度与实体单元位置关系近似，根据公式6，实体单元高度0.1，β＝59；可以得到进而可以计算每个三棱柱单元的伪密度。对有限元模型施加固定约束及载荷，即可得到分析模型。
[0062]
步骤6、以结构柔顺度为目标函数，优化柔顺度最小，以整体结构体积分数比0.3为体积约束条件，建立优化模型：
[0063]
find:
[0064]
min:
[0065]
s.t.
[0066]
步骤7、计算目标函数及体积约束对设计变量的灵敏度，将优化模型带入gcmma算法进行优化设计。优化迭代曲线如图5所示。表1展示优化前后的柔顺度与体分比对比。优化后参数域控制点分布及加筋布局如图6。
[0067]
表1
[0068] 结构柔顺度(kj)结构总体分比优化前59.680.387优化后23.950.3
[0069]
示例二：本发明用飞机前段模型的加筋优化。参照图7，模型为包含窗口的机头，工况条件为模型后端固定、整体受气动载荷。根据工程应用中的载荷施加方法、气动载荷可以等效为线框载荷。加筋许可高度为80mm；材料的杨氏模量和泊松比分别为e＝9
×
10
10
pa,ν＝
0.3。
[0070]
步骤1、对机头曲面用三角网格进行划分，对所得的三角曲面网格进行网格参数化。采用arap网格参数化算法，得到圆柱筒展开结果如图8。
[0071]
步骤2、计算展平后网格的扭曲度，得d
t
的最大值最小值分别为1.4、0.615，满足扭曲度的要求。则说明该网格参数化结果可以用于b样条参数化方法。
[0072]
步骤3、取b样条阶数为p＝q＝5，选取200
×
200个b样条控制点建立b样条参数场，并包含网格参数化结果，设计变量取值范围为
[0073]
步骤4、在机头网格上进行法向拉伸、得到总高度为80mm的8层三棱柱实体单元。加筋高度仍通过相应的参数域三角单元坐标的代入公式5得到。单元的伪密度可以表示为根据图7施加边界条件及截面载荷，得到有限元分析模型。
[0074]
步骤5、以结构柔顺度为目标函数，优化柔顺度最小，以整体结构体积分数比0.1为体积约束条件，建立优化模型：
[0075]
find:
[0076]
min:
[0077]
s.t.
[0078]
步骤6、利用gcmma算法，代入优化模型，计算目标函数及约束灵敏度，进行优化设计，优化迭代曲线图9所示。优化后参数域控制点分布及加筋布局如图10所示，这里加筋应位移机头的内表面，为了便于观察将其画于外表面。表2展示优化前后的柔顺度与体分比对比。
[0079]
表2
[0080] 结构柔顺度(kj)结构总体分比优化前1187.560.11优化后98.120.1
[0081]
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明公开的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。