CN108491574A - 基于光滑变形隐式曲线的结构形状设计方法 - Google Patents

Abstract

本发明提供一种固定网格下基于光滑变形隐式曲线的结构形状设计方法，具体按照下述步骤进行的：步骤一、采用光滑变形隐式曲线表达待优化结构的可设计边界；步骤二、基于R函数构建待优化结构的隐式模型Φ_Ω；步骤三、建立优化模型；步骤四、设置变量ε,P₂,P₃,...,P_n‑1的上下限，根据链式求导法则计算式(5)中目标函数和约束函数的灵敏度，进而采用基于梯度的优化算法进行优化设计求解。本发明在形状优化设计中采用变量少且变形能力强的光滑变形隐式曲线表达待优化结构的可设计边界，不仅有利于扩大优化设计空间以得到更优的解，还便于判断高斯积分点与结构区域的位置关系以提高优化设计效率。

Description

基于光滑变形隐式曲线的结构形状设计方法

技术领域

本发明属于结构优化设计领域，涉及一种结构形状优化设计方法，特别涉及一种固定网格下基于光滑变形隐式曲线的高精度形状优化设计方法。

背景技术

结构形状优化设计技术经过近半个世纪的发展，已广泛应用于航空航天、汽车、机械和土木工程等领域，解决的问题涉及到降低结构重量和应力水平、提高结构性能和安全寿命等多个方面。然而，传统形状优化设计方法在结构力学性能分析阶段，常采用基于贴体网格的有限元等数值分析方法，网格需随着结构边界的变化而不断地调整甚至重划分。这不仅使优化灵敏度无法求解或者在求解中涉及到繁琐的节点速度场计算环节，还容易造成网格畸变，从而迫使优化迭代提前终止。而近年来发展十分迅速的固定网格形状优化设计方法能够有效克服上述弊端，代表着形状优化设计技术未来的发展方向。

文献1“Cai SY,Zhang WH,Zhu JH,Gao T.Stress constrained shape andtopology optimization with fixed mesh:A B-spline finite cell method combinedwith level set function.Computer Methods in Applied Mechanics andEngineering,2014,278:361-387.”公开了一种固定网格下基于隐式三次样条曲线的高精度形状优化设计方法。该方法采用隐式三次样条曲线表达待优化结构的可设计边界，采用有限胞元方法(finite cell method)在固定网格下对结构进行高精度的应力分析，因此具有无需网格更新和灵敏度分析简便等优点。然而，该方法所采用的隐式三次样条曲线的变形能力有限且不便于表达闭合的结构边界。

文献2“ZhangWH,Huang QQ.Unification ofparametric and implicitmethodsfor shape sensitivity analysis and optimization with fixed mesh.InternationalJournal forNumerical Methods inEngineering,2016,109:326-344.”公开了一种固定网格下基于Bezier等自由曲线的高精度形状优化设计方法。该方法采用变形能力强的Bezier曲线表达待优化结构的可设计边界，采用有限胞元方法在固定网格下对结构进行高精度的应力分析，并基于Hamilton–Jacobi方程推导出解析灵敏度，使Bezier这类自由曲线能够和隐式曲线一样直接应用于固定网格形状优化设计中。然而，该方法在每步优化迭代中都要构建出逼近结构边界(由Bezier曲线描述)的多边形，以在固定网格分析时采用多边形内外点判别算法(The point-in-polygon algorithm)判断高斯积分点与结构区域的位置关系。

发明内容

为了克服现有固定网格形状优化设计方法在描述结构可设计边界方面的不足，本发明提供一种固定网格下基于光滑变形隐式曲线的结构形状设计方法。该方法将超二次曲线的半轴长扩展为极坐标中角度的函数，不仅具有较强的变形能力，还能在固定网格下简便快速地判断任意一点与结构区域的位置关系，有利于扩大优化设计空间和提高优化设计效率。

本发明的技术方案是这样实现的：一种基于光滑变形隐式曲线的结构形状设计方法，是按照下述步骤进行的：

步骤一、采用光滑变形隐式曲线表达待优化结构的可设计边界，光滑变形隐式曲线的表达式为：

式中，ε为形状指数，f(θ)为如下恒为正的连续函数：

式中，N_i,p(ξ)为B样条基函数，n是B样条基函数的个数，p是B样条基函数的次数，ξ为参数区域坐标且其取值范围为[0,1)，P_i为恒为正的控制参数且满足：P₁＝P_n＝(P₂+P_n-1)/2，θ为直角坐标(x,y)的函数：

式中，sgn(x)为x的符号函数；根据待优化结构的可设计边界即可设置光滑变形隐式曲线的n-1个参数(ε,P₂,P₃,...,P_n-1)的初始值；

步骤二、基于R函数构建待优化结构的隐式模型Φ_Ω，其中R函数系定义如下：

式中，Φ₁和Φ₂为两个几何体的隐函数，∧_α和∨_α分别为R-合取运算符和R-析取运算符，它们分别与布尔交操作∩和布尔并操作∪相对应，α为满足-1<α≤1的任意函数；

步骤三、建立优化模型，以给定体积上限为约束，以最小化结构最大vonMises应力为目标，可建立如下数学模型：

式中，V_lim为体积上限，H(·)为Heaviside函数，Ω为结构区域，D为计算区域，σ_von(x_Ω)为结构区域内任意一点x_Ω处的vonMises应力响应：

式中其中x_Ω处的应力分量计算如下：

(σ_x(x_Ω),σ_y(x_Ω),τ_xy(x_Ω))^T＝σ(x_Ω)＝DB(x_Ω)U (7)

式中，D为弹性矩阵，B为应变矩阵，U为由平衡方程组KU＝F求解得到的位移向量。根据有限胞元方法(finitecellmethod)，总体刚度矩阵K和总体载荷向量F的计算公式分别为：

K＝∫_DB^TDBH(Φ_Ω)dΩ (8)

式中，N为形函数矩阵，f为体积力向量，t为边界力向量，Γ_N为结构受力边界；

步骤四、设置变量ε,P₂,P₃,...,P_n-1的上下限，根据链式求导法则计算式(5)中目标函数和约束函数的灵敏度，进而采用基于梯度的优化算法进行优化设计求解。

本发明在形状优化设计中采用变量少且变形能力强的光滑变形隐式曲线表达待优化结构的可设计边界，不仅有利于扩大优化设计空间以得到更优的解，还便于判断高斯积分点与结构区域的位置关系以提高优化设计效率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是具体实施例的几何模型尺寸和载荷边界条件示意图。

图2是采用本发明方法得到的优化结果的几何模型示意图。

图3是具体实施例初始结构与优化结果的应力分布示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

一种基于光滑变形隐式曲线的结构形状优化设计方法，包括以下步骤：

步骤一、采用光滑变形隐式曲线表达待优化结构的可设计边界。光滑变形隐式曲线的表达式为：

式中，ε为形状指数，f(θ)为如下恒为正的连续函数：

式中，N_i,p(ξ)为B样条基函数，n是B样条基函数的个数，p是B样条基函数的次数，ξ为参数区域坐标且其取值范围为[0,1)，P_i为恒为正的控制参数且满足：P₁＝P_n＝(P₂+P_n-1)/2，θ为直角坐标(x,y)的函数：

式中，sgn(x)为x的符号函数。根据待优化结构的可设计边界即可设置光滑变形隐式曲线的n-1个参数(ε,P₂,P₃,...,P_n-1)的初始值。

步骤二、基于R函数构建待优化结构的隐式模型Φ_Ω。最常用的R函数系定义如下：

式中，Φ₁和Φ₂为两个几何体的隐函数，∧_α和∨_α分别为R-合取运算符和R-析取运算符，它们分别与布尔交操作∩和布尔并操作∪相对应，α为满足-1<α≤1的任意函数。将式(1)定义的光滑变形隐式曲线函数Φ_S和一些简单几何体(椭圆、矩形、三角形、梯形等)的隐函数通过上式进行一系列合取/析取运算，即能得到结构隐式模型。

步骤三、根据所考虑的优化问题建立起优化模型。以给定体积上限为约束，以最小化结构最大von Mises应力为目标，可建立如下数学模型：

式中，V_lim为体积上限，H(·)为Heaviside函数，Ω为结构区域，D为计算区域，σ_von(x_Ω)为结构区域内任意一点x_Ω处的vonMises应力响应：

式中其中x_Ω处的应力分量计算如下：

(σ_x(x_Ω),σ_y(x_Ω),τ_xy(x_Ω))^T＝σ(x_Ω)＝DB(x_Ω)U (7)

式中，D为弹性矩阵，B为应变矩阵，U为由平衡方程组KU＝F求解得到的位移向量。根据有限胞元方法(finite cell method)，总体刚度矩阵K和总体载荷向量F的计算公式分别为：

K＝∫_DB^TDBH(Φ_Ω)dΩ (8)

式中，N为形函数矩阵，f和t分别为体积力向量和边界力向量，Γ_N为结构受力边界。上面两式的积分计算方法与有限元方法基本相同：先将计算区域D离散成规则的胞元，再在胞元里面采用高斯积分法进行计算。而不同之处在于：在实施高斯积分时需根据结构隐式模型Φ_Ω判断高斯积分点是否属于结构区域Ω内部。

步骤四、设置变量ε,P₂,P₃,...,P_n-1的上下限，根据链式求导法则计算式(5)中目标函数和约束函数的灵敏度，进而采用基于梯度的优化算法进行优化设计求解。

本发明的方法在结构应力分析方面采用的是与背景技术文献1和2相同的有限胞元方法，但在结构可设计边界描述方面采用的是一种能够自由且光滑变形的隐式曲线。该曲线将超二次曲线的半轴长扩展为极坐标中角度的函数，不仅具有较强的变形能力，还能在固定网格下简便快速地判断任意一点与结构区域的位置关系，有利于扩大优化设计空间和提高优化设计效率。

下面结合具体实施方式说明本发明的应用：

参照图1～3，以开孔平板结构形状优化设计为例说明本发明。开孔平板结构的板厚为3mm，长宽均为40mm，中心圆孔半径为5mm，材料杨氏弹性模量和泊松比分别为207.4GPa和0.3，四条边受三角形分布拉应力且最大应力值均为50MPa。设计开孔平板结构的中心圆孔形状，在保证结构体积不增加的条件下使其最大von Mises应力最小，在建模时不失一般性地将圆孔圆心作为坐标原点。方法步骤如下：

步骤一、采用光滑变形隐式曲线表达开孔平板结构的孔洞边界。光滑变形隐式曲线的表达式为：

式中，ε为形状指数，f(θ)为如下恒为正的连续函数：

式中，N_i,2(ξ)为定义在均匀节点矢量{0,0,0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1,1,1}上的二次B样条基函数，ξ为参数区域坐标且其取值范围为[0,1)，P_i为恒为正的控制参数且满足：P₁＝P₂₂＝(P₂+P₂₁)/2，θ为直角坐标(x,y)的函数：

式中，sgn(x)为x的符号函数。

由于开孔平板结构的初始孔洞边界是半径为5mm的圆形边界，则将形状指数ε的初值设置为1，将控制参数P_i(i＝2,3,…,21)的初值均设置为5mm。此时，根据式(2)可知f(θ)恒为5mm，再根据式(1)可知Φ_S就是半径为5mm的圆形的隐函数。

步骤二、基于R函数构建开孔平板结构的隐式模型Φ_Ω。最常用的R函数系定义如下：

式中，Φ₁和Φ₂为两个几何体的隐函数，∧_α和∨_α分别为R-合取运算符和R-析取运算符，它们分别与布尔交操作∩和布尔并操作∪相对应，α为满足-1<α≤1的任意函数。

开孔平板结构的几何模型可通过对其外边界构成的矩形和中心圆孔模型进行布尔差操作而得到，其中矩形的隐函数可由四条边的隐函数构建为：Φ_R＝(20-x)∧₀(x+20)∧₀(20-y)∧₀(y+20),则开孔平板结构的隐式模型为：Φ_Ω＝Φ_R∧₀(-Φ_S)。

步骤三、根据所考虑的优化问题建立起优化模型。以给定体积上限为约束，以最小化结构最大von Mises应力为目标，可建立如下数学模型：

式中，V_lim选为初始结构的体积4564.3mm³，H(·)为Heaviside函数，Ω为结构区域，D为计算区域，本实施例将开孔平板结构外边界围成的矩形区域作为D，σ_von(x_Ω)为结构区域内任意一点x_Ω处的vonMises应力响应：

式中其中x_Ω处的应力分量计算如下：

(σ_x(x_Ω),σ_y(x_Ω),τ_xy(x_Ω))^T＝σ(x_Ω)＝DB(x_Ω)U (7)

式中，D为弹性矩阵，B为应变矩阵，U为由平衡方程组KU＝F求解得到的位移向量。根据有限胞元方法(finite cell method)，总体刚度矩阵K和总体载荷向量F的计算公式分别为：

K＝∫_DB^TDBH(Φ_Ω)dΩ (8)

式中，N为形函数矩阵，f为体积力向量(在本实施例中f为零向量)和t为边界力向量(边界力分布见图1)，Γ_N为结构外边界。上面两式的积分计算方法与有限元方法基本相同：先将计算区域D均匀离散成40×40个胞元，再在胞元里面采用高斯积分法进行计算。而不同之处在于：在实施高斯积分时需根据结构隐式模型Φ_Ω判断高斯积分点是否属于结构区域Ω内部。

步骤四、将变量ε的上下限设置为10和0.01，将变量P₂,P₃,...,P₂₁的上下限统一设置为20mm和0.01mm，根据链式求导法则计算式(5)中目标函数和约束函数的灵敏度，进而采用基于梯度的移动渐近线优化算法(The Method of Moving Asymptotes)进行优化设计求解。

由图2可见，优化结果的孔形虽然不规则但仍保持几何对称性，这是因为开孔平板结构所受的载荷是对称的。由图3可见，优化结果的孔边应力集中程度下降明显。表1列出了优化前后的结构最大von Mises应力和结构体积，表明优化结果在体积不增加的情况下，最大von Mises应力的降幅达45.7％。与背景技术中文献1的方法相比较，本发明方法由于采用了具有较强变形能力且便于表达闭合边界的光滑变形隐式曲线，因此具有较大的优化设计空间且便于对本实施例的孔洞边界进行形状优化设计；与背景技术中文献2的方法相比较，本发明方法由于能够通过隐式模型Φ_Ω的正负性快速判断高斯积分点与结构区域的位置关系，无需在每步优化迭代中都要构建出逼近结构边界的多边形，因此具有较高的优化设计效率，表2对比了这两种方法基于MATLAB平台对本实施例的初始结构进行一次应力分析所使用的时间。

表1

	最大vonMises应力	结构体积
			优化前	39.2MPa	4564.3mm3
优化后	21.3MPa	4563.8mm3

表2

	文献2的方法	本发明方法
			分析时间	52s	37s

对于开孔平板结构形状优化设计问题，本发明方法仅需一条光滑变形隐式曲线就可描述闭合的孔洞边界，无需网格更新操作，无需构建逼近结构边界的多边形，而且优化结果在体积不增加的情况下，其最大von Mises应力(21.3MPa)相对于初始结构最大vonMises应力(39.2MPa)的降幅达45.7％。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于光滑变形隐式曲线的结构形状设计方法，其特征在于是按照下述步骤进行的：

步骤一、采用光滑变形隐式曲线表达待优化结构的可设计边界，光滑变形隐式曲线的表达式为：

式中，ε为形状指数，f(θ)为如下恒为正的连续函数：

式中，N_i,p(ξ)为B样条基函数，n是B样条基函数的个数，p是B样条基函数的次数，ξ为参数区域坐标且其取值范围为[0,1)，P_i为恒为正的控制参数且满足：P₁＝P_n＝(P₂+P_n-1)/2，θ为直角坐标(x,y)的函数：

式中，sgn(x)为x的符号函数；根据待优化结构的可设计边界即可设置光滑变形隐式曲线的n-1个参数(ε,P₂,P₃,...,P_n-1)的初始值；

步骤二、基于R函数构建待优化结构的隐式模型Φ_Ω，其中R函数系定义如下：

式中，Φ₁和Φ₂为两个几何体的隐函数，∧_α和∨_α分别为R-合取运算符和R-析取运算符，它们分别与布尔交操作∩和布尔并操作∪相对应，α为满足-1<α≤1的任意函数；

步骤三、建立优化模型，以给定体积上限为约束，以最小化结构最大von Mises应力为目标，可建立如下数学模型：

Find ε,P₂,P₃,…,P_n-1

S.t.KU＝F

V＝∫_DH(Φ_Ω)dΩ≤V_lim

ε＞0,P_i＞0i＝2,3,…,n-1 (5)

式中，V_lim为体积上限，H(·)为Heaviside函数，Ω为结构区域，D为计算区域，σ_von(x_Ω)为结构区域内任意一点x_Ω处的vonMises应力响应：

式中其中x_Ω处的应力分量计算如下：

(σ_x(x_Ω),σ_y(x_Ω),τ_xy(x_Ω))^T＝σ(x_Ω)＝DB(x_Ω)U (7)

式中，D为弹性矩阵，B为应变矩阵，U为由平衡方程组KU＝F求解得到的位移向量。根据有限胞元方法(finite cell method)，总体刚度矩阵K和总体载荷向量F的计算公式分别为：

K＝∫_DB^TDBH(Φ_Ω)dΩ (8)

式中，N为形函数矩阵，f为体积力向量，t为边界力向量，Γ_N为结构受力边界。

步骤四、设置变量ε,P₂,P₃,...,P_n-1的上下限，根据链式求导法则计算式(5)中目标函数和约束函数的灵敏度，进而采用基于梯度的优化算法进行优化设计求解。