渗流问题的拓扑优化

周向阳，范良志

（武汉纺织大学 机械工程与自动化学院，湖北 武汉 430073）

渗流问题的拓扑优化

周向阳，范良志

（武汉纺织大学 机械工程与自动化学院，湖北 武汉 430073）

分析拓扑优化中的密度惩罚函数插值法SIMP（Solid isotropic material with penalization model），将结构力学中的拓扑优化方法应用到渗流问题的拓扑优化设计中，用有限元的方法建立了简单的理想状态下的渗流问题的拓扑优化数学模型，并采用基于剃度法的数值解法－优化准则法（OC系列算法），以设计具有最小能量损耗情况下的流体最佳流动路径的拓扑分布。以一组二维渗流问题为例，说明了该模型的有效性，为渗流问题的优化设计提供了一种有效的新思路和方法。

拓扑优化；渗流；最小能量损耗；最佳流动路径

结构的拓扑优化设计是结构的尺寸优化设计和形状优化设计以后，在结构优化领域出现的一种新型的富有挑战性的研究方向[1]。从有关文献看，连续体结构的拓扑优化设计技术，特别是SIMP方法[1,2]，已经成功的用于机械设计[3,4]、MEMS系统[5,6]、材料设计[7]等方面。在航空航天、汽车制造等固体结构领域，拓扑优化技术不仅仅是学术上的应用，而且成为一种实用的设计工具被广泛采用。把拓扑优化理论用于流体场流体流动设计是拓扑优化技术应用的新的热点研究方向之一。其方法是：在给定边界条件的设计域Ω中，确定哪部分是流体哪部分是非流体，使得满足规定流体部分体积比的某个目标最小化。

由于渗流是流体中比较简单的情况，在实际工程中也比较常见，比如闸坝的渗流，农田地下排水系统，油气地下渗流，裂隙岩体渗流以及裂隙排水等，作者首先尝试将用于固体结构的刚度拓扑优化技术用于渗流问题的优化设计。本文主要考虑流体在多孔介质中的流动，例如水在地坝中的流动，在管道中或者围绕固体的流动，并且只考虑理想流体在稳太、无旋（即流体质点只是平动）、不可压缩（等质量密度）和无粘（没有粘性）状态下的情况，流体与表面之间的摩擦也被忽略，流体也并不渗透到周围物体或并不与物体的表面分开。取介质相对密度为设计变量，流体在整体设计域的体积比为约束，以系统的最小能量损耗为目标函数。要说明的是：目的在于探讨拓扑优化技术在流体方面的应用的可行性和方法，由于本文选取理想渗流，所以要使模型符合实际的应用还有待进一步的研究。

1 SIMP方法回顾

SIMP模型主要通过引入惩罚因子，在材料的弹性模量和单元相对密度之间建立起一种显示的非线性对应关系。它的作用是当设计变量的值在(0,1)之间时，对中间密度值进行惩罚，使中间密度值逐渐向0-1两端聚集，这样可以使连续变量的拓扑优化模型能较好地逼近原来0-1离散变量的优化模型。

SIMP材料模型的数学表达形式：

E表示插值以后的弹性模量，0E为初始弹性模量表示单元j的设计变量即j单元的相对密度，为了避免刚度矩阵奇异，给xj一个大于0的下限值表示j单元初始刚度矩阵，表示第j单元优化后的刚度矩阵。 为两数学模型中对中间密度材料的惩罚因子。为有效压缩中间密度材料，

以结构的最小柔度设计问题为例，其拓扑优化模型可表示为：

X表示设计变量， C表示结构的柔度，V表示优化后的有限元单元体积列向量， V*表示优化的目标体积即体积约束，F表示力矢量，K表示结构的刚度矩阵，U表示位移矢量。

这里只有一个约束条件，一般情况下采用优化准则法。它是由目标函数和约束条件构成的拉格朗日函数，在满足Kuhn-Tucker条件下推导出相应的迭代求解公式。在综合考虑设计变量上下限的情况下，可得问题（3）式的优化准则法求解公式如下：

式中 为阻尼系数，引入 的目的是为了确保数值计算的稳定性和收敛性，有关 取值范围的讨论详见文献[8]。

2 渗流问题拓扑优化模型

类似地，可以建立渗流问题的拓扑优化模型。考虑在多孔介质中的二维流体流动，根据质量守恒和达西定律，对于常数渗透系数流体的流动微分方程为：

式中，Kxx和Kyy分别为多孔介质在x和y方向的渗透系数，Φ是流体的水头或者速度势函数，是单位体积的体积流动率。

式中Cx和Cy为表面S2的单位法线矢量的方向余弦，与图1中所示的相同。式（8）表明，在边界表面S1上已知边界流体水头或者速度势ΦB，而式（9）则表明，在垂直于表面S2上的势的剃度或者速度为已知。在不渗透的边界上速度或剃度等于0。

一般来说，基本单元可以有内部源或汇，例如来自泵的源或汇，或来自河或江的表面边界流动率。为了包括这些影响，单元外载荷包括作用在整个单元上的均匀内源Q和作用在表面上的均匀表面流动率源*q，则力矩阵的项分别是：

由控制方程、边界条件和初始条件，构造渗流场的范函表达式，然后变分求解。在求解区域上离散单元，选取插值函数，经过单元分析和组装，即可得到渗流结构的有限元表达式

类似的可以建立渗流问题的拓扑优化数学模型为：

图1 流体流动的边界条件

式中[B]为联系水力剃度与节点势的矩阵，与固体结构有限元中的单元应变矩阵相似，其求法也与单元应变矩阵一样。[D]为材料性质矩阵，对于各向同性渗透流体中有则

3 算法和二维算例

目前适用于拓扑优化的优化算法主要包括两类：优化准则法（又称OC系列算法）[9]、序列线性规划或序列二次规划算法，后者较为典型的是移动渐进算法（又称MMA算法，the method of moving asymptotes）[10]。

优化准则法是一种间接的优化方法，它不直接优化目标函数，而是基于Kuhn—Tucker条件，通过构造Lagrange函数来形成设计变量的更新方案，一般适用于大量设计变量、单目标、单约束条件下问题的优化。对于OC系列方法的详细介绍见文献[9]。

MMA方法是一种更高级的数学规划算法，能够适用于单约束情况和多目标、多约束情况下问题的求解，但其计算收敛性不够理想。在实际工程应用中，应根据不同情况选用不同的优化求解算法。关于MMA系列方法的详细介绍参考Bruyneel等人的文献[10]。

这里的问题属于单目标单约束的优化问题，所以选择优化准则法（OC算法），并在matlab上实现。

以二维平面问题为例，设计出在不同的初始边界条件时，使系统具有最小能量损耗的最佳流体流动路径。

为了消除拓扑优化数值计算中的棋盘格和网格依赖现象，在算法中采用了惩罚项，取惩罚因子p=3.0，加入了阻力项，取阻力系数η= 0.5，采用了局部敏度过滤算法[11]。模型离散为四节点四边形单元。渗透系数均为1。拓扑分布图中，不可渗透材料部分用黑色表示，流体流动部分用白色表示。

算例一：区域为30×30单元，左右边界分别有两点水力，速度为1，水流进为正，流出为负。结构尺寸如图2（a），体积约束比为0.3，优化结果如图2（b），迭代过程如图2（c）。

算例二：区域为30×30单元，左上角下右角两点水力，速度为1。结构尺寸如图3（a），体积约束比为0.25，优化结果如图3（b），迭代过程如图3（c）。

算例三：区域为30×30单元，左边界均匀分布水力，速度为1，右边界中间节点势为0。结构尺寸如图4（a）。体积约束比 f为0.5，优化结果如图4（b），迭代过程如图4（c）。体积约束比 f为0.4，优化结果如图4（d），迭代过程如图4（e）。

算例四：区域为50×50单元，左右边界分别有一点水力，速度为1。结构尺寸如图5（a），在中间偏下位置有一个直径为15的圆形空洞。体积约束比 f为0.55，优化结果如图5（b），迭代过程如图5（c）。体积约束比 f为0.3，优化结果如图5（d），迭代过程如图5（e）。

图2（a）算例一设计域

图2（b）算例一优化结果

图2（c）算例一迭代过程

图3（a）算例二的设计域

图3（b）算例二优化结果

图3（c）算例二迭代过程

图4（a）算例三的设计域

图4（b）体积比0.5的优化结果

图4（c）体积比0.5的迭代过程

图4（d）体积比0.4的优化结果

图4（e）体积比0.4的迭代过程

图5（a）算例四的设计域

图5（b）体积比0.55的优化结果

图5（c）体积比0.55的迭代过程

图5（d）体积比0.3的优化结果

图5（e）体积比0.3的迭代过程

4 结论

从以上算例我们看出，拓扑优化方法在渗流问题的拓扑优化设计中的应用，没有出现象在固体结构设计中出现的数值不稳定现象，比如多孔材料、棋盘格、网格依赖性、局部极值等。在算例三和四中出现了一些小岛问题只与约束体积比有关，当系统采用的约束体积比太小时才会出现。算例表明了拓扑优化方法在解渗流问题时的可行性和有效性。

将拓扑优化方法应用于渗流问题的拓扑优化设计是一种新的思想和方法。该方法能够设计出具有最小总体能量损耗情况下的流体最佳流动路径，为渗流结构设计提供了一种新的方法和思路。

[1] Bendsøe M, Sigmund O. Topology optimization: Theory, Methods, and Applications [M]. New York: Springer, 2003.

[2] Bendsoe M, Sigmund O. Material interpolations in topology optimization [J]. Archive of Applied Mechanics, 1999,69: 635～654.

[3] 石琴，姚成，马恒永.集装箱半挂车车架结构拓扑优化设计[J].农业机械学报,2005,36(1):10～12，17.

[4] 范文杰，徐进永，张子达.基于双向渐进结构优化法的装载机动臂拓扑优化[J].农业机械学报,2006,37(11):24～27.

[5] 左孔天, 陈立平, 王书亭, 等. 用拓扑优化方法进行微型柔性机构的设计研究[J]. 中国机械工程, 2004, 15(21): 1886-1890.

[6] 罗震, 蒙永立, 郭文德, 等．分布式柔性机构拓扑优化设计的理论和算法[J]. 机械工程学报，2006，42(10): 27～36

[7] Sigmund O. A new class of extremal composites [J]. J Mech Phys Solid, 2000, 48(2)：397-428.

[8] 左孔天, 陈立平, 王书亭，等. 基于人工材料密度的新型拓扑优化理论和算法研究[J]. 机械工程学报, 2004, 40(12):23～29.

[9] Hassani B, Hinton E. Homogenization and Structural Topology Optimization Theory, Practice and Software [M]. London: Springer,1999.

[10] Bruyneel M, Duysinx P, Fleury C. A family of MMA approximations for structural optimization [J]. Structural and Multidiscipline Optimization, 2002, 24: 263～276.

[11] Sigmund O, Petersson J. Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards,mesh-dependencies and local minima [J]. Structural Optimization, 1998,16: 68～75.

Topology Optimization of Seepage Flow Problems

ZHOU Xiang-yang, FAN Liang-zhi
(College of Mechanical Engineering and Automation, Wuhan Textile University, Wuhan Hubei 430073, China)

The SIMP method used in topology optimization is analyzed. The topology optimization method is applied into the topology optimization design of seepage flow. The finite element method is used to model the simple flow problems under the ideal conditions. And the optimization problem is solved with a gradient-based math-programming algorithm (OC) that is driven by analytical sensitivities. The topology optimization model for seepage flow is established, which is used to minimize the energy dissipation in the system and find the optimal flow route. Several two-dimensional examples were tested. The results show the model is effective. A new approach for the design of seepage structure is presented.

Topology Optimization; Seepage Flow; Energy Dissipation Minimization; Optimal Flow Route

TH11

A

1009－5160(2011)06－0079－06

周向阳（1977-），女，讲师，博士，研究方向：数字化设计与制造.

国家自然科学基金（50805109）.