基于参数化水平集方法的流体管道拓扑优化设计方法

摘要

本发明涉及一种基于参数化水平集方法的流体管道拓扑优化设计方法，属于流体管道拓扑优化设计领域。解决传统基于材料分布模型的流体管道拓扑优化方法无法直接准确获取管道边界信息，且设计变量较多的技术问题。所述方法是使用径向基函数对初始水平集函数进行插值，得到的插值系数则为优化问题的设计变量。再利用Heaviside函数将经插值得到的水平集函数映射成设计区域各个有限单元上的间接设计变量，从而将水平集模型转化为了材料分布模型，然后以流量驱动流通过设计区域的耗散能最小作为目标，以流体管道体积作为约束进行设计。该方法能够有效地减少设计变量的个数，且同时提供更加精确的管道边界描述。

说明书

技术领域

本发明涉及流体管道拓扑优化设计领域，具体涉及一种基于参数化水平集方法的流体管道拓扑优化设计方法。

背景技术

现有流体管道拓扑优化设计方法主要有两种，一种方法是文献1提出的基于伪密度的材料分布模型拓扑优化方法，该方法在满足特定体积约束的情况下，以流体耗散能极小为目标，使用有限单元方法，通过伴随敏度分析并应用梯度优化算法更新迭代得到流体管道的最优拓扑(T.Borrvall,J.Petersson,Topology optimization of fluid in Stokes flow,Int.J.Numer.Methods Fluids 41(2003)77-107.)。上述方法设计变量为每个有限单元对应的伪密度，该方法不依赖初始给定结构，能够有效地寻找出最优拓扑，且经过十几年的发展，已经形成一整套成熟算法。不过该方法不能直接获取管道边界信息，设计者只能通过其他方法间接得到，具有很大的人为因素。且网格划分较密集时，设计变量数量较大。

另一种方法是文献2提出的将管道边界用高一维度的水平集函数的零水平集隐式表达，并通过求解汉密尔顿-雅克比方程来演化管道边界最终找到较优拓扑结构的方法，该方法也是将流体体积作为约束条件，以高一维度的水平集函数为设计变量，通过演化水平集函数找到最优拓扑，在优化迭代过程中具有较为清晰准确的边界信息(S.Zhou,Q.Li,Avariational level set method for the topology optimization of steady-state Navier-Stokes flow,J.Comput.Phys.227(2008)10178-10195.)。但是该方法由于有限差分方法的局限性，收敛较慢，且最终结果较依赖初始给定结构，边界信息也不能由解析的数学式表达。

发明内容

本发明要解决传统基于材料分布模型的流体管道拓扑优化方法无法直接准确获取管道边界信息，且设计变量较多的技术问题，提供一种基于参数化水平集方法的流体管道拓扑优化设计方法。

为了解决上述技术问题，本发明的技术方案具体如下：

一种基于参数化水平集方法的流体管道拓扑优化设计方法，包括以下步骤：

步骤一、确定管道设计区域Ω和管道的出入口位置，以及紧支撑径向基函数的数量N及在设计区域中的分布利用紧支撑径向基函数对初始水平集函数进行拟合近似，得到初始设计变量α＝(α₁,α₂,...,α_N)；

步骤二、定义流体拓扑优化模型，将设计区域Ω离散，定义Φ(α)为优化目标函数，在本问题里为流体在管道中的耗散能极小，约束条件为管道体积小于或等于

findα＝(α₁,α₂,...,α_N)

minΦ(α)

u＝u_D,onΓ_D

α_i,min≤α_i≤α_i,max,i＝1,...,N

其中u为流体速度，p为流体压强，η为流体粘滞系数，ρ为流体密度，Γ_D为Dirichlet边界，Γ_N为Neumann边界；u_D为定义在Γ_D上的速度分布，g为定义在Γ_N上的应力分布；n为边界上的单位外法向量，I为单位对角张量；α_i,min与α_i,max分别为第i个设计变量α_i的下边界及上边界；β代表材料的不渗透性，它的值可以由下式给出：

其中q为调节β凸性的参数，β_max通常为一个较大的值，β_min通常取为0；H是经光滑化处理的Heaviside函数；

步骤三、应用有限单元法分析计算流场状态变量u和p，根据得到的状态变量计算流体耗散能；

步骤四、利用伴随变量法及有限单元法计算目标函数对于设计变量的灵敏度；

步骤五、根据上述求得的灵敏度进行优化，选取梯度优化算法进行优化迭代，最终得到优化结果。

所述基于参数化水平集方法的流体管道拓扑优化设计方法的一种具体方案为：

步骤一、管道拓扑优化设计区域Ω长度与宽度均为1，管道入口与出口呈90度，紧支撑径向基函数在水平方向与竖直方向等间距0.1分布在设计区域内，得到紧支撑径向基函数在设计区域的分布初始水平集函数的值统一设为0.2，经插值得到初始设计变量α＝(α₁,α₂,…,α₁₂₁)；在本例中所使用的紧支撑径向基函数为C4连续的Wendland径向基函数，任何其他连续性不小于C2连续的径向基函数都可以使用：

其中(1-r)₊＝max{0,1-r}，r是支撑半径，在二维欧几里得空间中的定义如下：

其中x_i,y_i定义第i个径向基函数的中心坐标，R_i表示该径向基函数的影响范围，在本例中取为2.5倍的相邻径向基函数之间的距离，即0.25；

步骤二、定义流体拓扑优化模型，将设计区域Ω离散为100×100的正方形网格，Φ(α)为优化目标函数，在本问题里为流体在管道中的耗散能极小，约束条件为管道体积小于或等于0.08π：

findα＝(α₁,α₂,...,α₁₂₁)

minΦ(α)

u＝u_D,onΓ_D

V(α)＝∫_Ω(1-H(Φ))dΩ≤0.08π

α_i,min≤α_i≤α_i,max,i＝1,...,N

其中u为流体速度，p为流体压强，η为流体粘滞系数，均取为1；ρ为流体密度，取为1；Γ_D为Dirichlet边界，Γ_N为Neumann边界；u_D为定义在Γ_D上的速度分布：

u_D＝-20(y-0.7)(0.9-y)n

g为定于在Γ_N上的应力分布，取为0；n为边界上的单位外法向量，I为单位对角张量；α_i,min与α_i,max分别为第i个设计变量α_i的下边界及上边界，对于第k个迭代步，设计变量的下边界及上边界均由下式给出：

β代表材料的不渗透性，它的值可以由下式给出：

其中q为调节β凸性的参数，在优化的初始阶段取为0.1，迭代的后半段取为1；β_max取为10000，β_min取为0；H是经光滑化处理的Heaviside函数，在本例中取为：

其中h为一个正的常数，表示光滑区域的宽度，在本例中取为0.5。

步骤三、应用有限单元法分析计算流场状态变量u和p，根据得到的状态变量计算流体耗散能；

步骤四、利用伴随变量法及有限单元法计算目标函数对于设计变量的灵敏度；

步骤五、根据上述求得的灵敏度进行优化，选取移动渐近线方法(The Method of Moving Asymptotes)进行优化迭代，最终得到优化结果。

本发明的有益效果是：

本发明提供的基于参数化水平集方法的流体管道拓扑优化设计方法采用紧支撑径向基函数对隐式表达的水平集函数插值拟合，使其可以解析表达，并将插值所得到的水平集函数通过Heaviside函数映射成各个单元上的间接设计变量，从而转化为传统的材料分布模型。在优化过程中以径向基函数的插值系数作为设计变量，以流体耗散能极小为目标函数，应用伴随方法求解目标函数对于设计变量的灵敏度，使用有限单元方法求解，并使用成熟的梯度优化算法更新迭代，最终得到设计结果。该方法能够有效地利用材料分布方法较易产生新拓扑的优点，并且通过结合参数化的水平集函数，得到更加精确的管道边界信息。且设计变量的个数不依赖于有限单元的划分，可以很大程度上减少设计变量的个数。本发明经过实施例120步迭代后得到设计结果。比较上述文献1的结果，基本相同，设计变量的个数减少了98％以上，且边界具有水平集方法的光滑性，并可以解析表达。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细说明。

图1为本发明方法实施例中设计区域几何模型图。

图2为本发明方法实施例的设计初始水平集函数分布图。

图3为本发明方法实施例的设计初始水平集函数对应的材料分布图。

图4为本发明方法实施例的设计结果水平集函数分布图。

图5为本发明方法实施例的设计结果水平集函数对应的材料分布图。

具体实施方式

下面结合附图对本发明做以详细说明。

参照附图1-5：本发明提供的基于参数化水平集方法的流体管道拓扑优化设计方法具体步骤如下：

以二维的流量驱动流通过一段出口与入口呈90度的弯管为实例说明本发明。

步骤一、管道拓扑优化设计区域Ω长度与宽度均为1，管道入口与出口呈90度，紧支撑径向基函数在水平方向与竖直方向等间距0.1分布在设计区域内，得到紧支撑径向基函数在设计区域的分布初始水平集函数的值统一设为0.2，经插值得到初始设计变量α＝(α₁,α₂,…,α₁₂₁)；在本例中所使用的紧支撑径向基函数为C4连续的Wendland径向基函数，任何其他连续性不小于C2连续的径向基函数都可以使用：

其中(1-r)₊＝max{0,1-r}，r是支撑半径，在二维欧几里得空间中的定义如下：

其中x_i,y_i定义第i个径向基函数的中心坐标，R_i表示该径向基函数的影响范围，在本例中取为2.5倍的相邻径向基函数之间的距离，即0.25；

步骤二、定义流体拓扑优化模型，将设计区域Ω离散为100×100的正方形网格，Φ(α)为优化目标函数，在本问题里为流体在管道中的耗散能极小，约束条件为管道体积小于或等于0.08π：

findα＝(α₁,α₂,...,α₁₂₁)

minΦ(α)

u＝u_D,onΓ_D

V(α)＝∫_Ω(1-H(Φ))dΩ≤0.08π

α_i,min≤α_i≤α_i,max,i＝1,...,N

其中u为流体速度，p为流体压强，η为流体粘滞系数，均取为1；ρ为流体密度，取为1；Γ_D为Dirichlet边界，Γ_N为Neumann边界；u_D为定义在Γ_D上的速度分布：

u_D＝-20(y-0.7)(0.9-y)n

g为定于在Γ_N上的应力分布，取为0；n为边界上的单位外法向量，I为单位对角张量；α_i,min与α_i,max分别为第i个设计变量α_i的下边界及上边界，对于第k个迭代步，设计变量的下边界及上边界均由下式给出：

β代表材料的不渗透性，它的值可以由下式给出：

其中q为调节β凸性的参数，在优化的初始阶段取为0.1，迭代的后半段取为1；β_max取为10000，β_min取为0；H是经光滑化处理的Heaviside函数，在本例中取为：

其中h为一个正的常数，表示光滑区域的宽度，在本例中取为0.5。

步骤三、应用有限单元法分析计算流场状态变量u和p，根据得到的状态变量计算流体耗散能；

步骤四、利用伴随变量法及有限单元法计算目标函数对于设计变量的灵敏度；

步骤五、根据上述求得的灵敏度进行优化，选取移动渐近线方法(The Method of Moving Asymptotes)进行优化迭代，最终得到优化结果。

本发明方法经过实施例120步迭代后得到设计结果，经与文献1中实例对比，管道边界基本相同。设计变量的个数从文献1的10000个下降到121个，下降了98.79％。且边界更加光滑，并可由插值函数解析表达出边界信息。

显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本发明创造的保护范围之中。