段献葆,李飞飞,秦新强

(西安理工大学理学院，西安 710048)

1 引言

大量的工程实际问题表明，初始设计阶段选择一个好的拓扑或形状可以使得材料等有很大的节省.从Pironneau学者[1,2]开始，在过去的几十年中，许多数学家和工程师致力于这一领域的研究，并且取得了非常好的成效[3-6].由于流体流动及其数学描述的复杂性，直到最近，特别是随着计算技术的发展，流体中的形状最优控制问题才得到很好的发展.但是，对流体力学中的拓扑优化问题进行系统的研究则只有十几年的时间[7].

迄今为止，已经有很多用来求解形状最优控制问题的算法.其中最出名的，也用得非常成功的算法是由Bendsoe和Kikuchi[8]提出的均质化方法(homogenization approach)，该算法把整个计算区域离散成足够多的微小区域，然后根据所求解的物理问题对这些微小区域进行分析，并按照一定的优化算法删除或增加某些微小区域，用最终保留下来的微小区域来描述最优拓扑或形状.随后，Bendsoe[9]对均质化方法加以改进，提出了SIMP(solid isotropic microstructure with penalization)方法.SIMP方法属于变密度法的一种.均质化方法是把计算区域中的离散点假设为相对密度在0到1之间的可变材料，而SIMP方法则是用连续函数来表达区域离散点处的相对密度与材料物理属性之间的对应关系，通过引入惩罚因子，使中间的密度函数值趋于0或1，从而取得与均质化方法一样的效果.该方法基于假设的各向同性材料，以每个单元的相对密度作为设计变量，不需要引入微结构和均质化过程.由于其程序实现简单，计算高效稳定，并且较好的消除了均质化方法的棋盘格现象，SIMP方法在形状最优控制领域取得了非常广泛的应用.

不管是均质化方法还是SIMP方法，一般都先要与其他优化算法相结合，如最优化准则方法(OC)[3]，然后进行迭代求解.当所求解的问题比较简单，计算区域不太复杂，方程组规模较小的情况下，计算时间会比较短，是一个极为有效的方法.由于流体流动的复杂性，采用这两种方法时总计算量都非常大.其中很大一方面的原因是在实际优化的过程中，一般都是对整个计算区域进行均匀剖分.而在形状最优控制问题中，我们更关心的是计算区域的内部边界处形状的改变.为此，我们提出了一种基于材料分布的自适应网格方法，并将该方法用于求解流体力学形状最优控制问题.材料分布的信息在使用SIMP方法时已经得到，不需要额外的计算.所提算法的计算主要集中在流体的内边界区域，从而在达到同样分辨率的情况下可以有效的减少计算量.

2 最优形状设计问题

在流体力学最优形状设计问题中，一般考虑在控制方程约束下使得流体流动的能量耗散达到最小.设D为包含所有可能形状Ω的工作区域，是R2中具有Lipschitz连续边界的开区域和p分别表示流体的速度和压力，f表示体积力.不失一般性，我们假设流体是在理想介质中流动，并且符合Darcy定律[7]，即是在点x的局部渗透性常数的倒数，可以表示为如下凸函数

其中是需要优化的设计变量，也是用来控制介质在每个点处的渗透性.ϱ=0表示固体材料，流体无法通过，而ϱ=1表示流体流动的区域.θ是一个可以用来控制α(ϱ)的正实数.在现实问题中，对于固体材料比如说墙壁，需要但从数值计算的角度只能选择一个尽可能大的数，在本文中最小值为状态方程为Stokes问题的最优控制问题的一般形式为：求u和ϱ，使得目标泛函

在满足

为求解最优控制问题(2)–(8)，需要目标泛函(2)关于区域的灵敏度分析结果.这可以用共轭方法得到，即，如果状态方程是Stokes问题(3)–(6)，则目标泛函(2)关于Ω的形状导数为[6]

其中ε＞0是一个小的正实数，(P,q)是下面问题的解

在本文中，使用

作为最优化准则方法中的下降方向.

3 最优化准则及自适应网格方法

最优化准则方法在上世纪50年代就已经被用于工程结构设计，后来数学家和工程师把Kuhn-Tucker条件作为最优结构满足的准则，从而使该方法有了数学理论的保障，通用性也得到提高.该方法的最大优点是对设计变量修改较大，所以收敛速度快，迭代次数与结构大小及复杂程度无关，并且原理简单、直观、易为工程设计人员接受与掌握.最优化准则方法已经被用于求解很大一类最优形状设计问题.本文所提的算法同样可以与移动渐进方法(MMA)[10]或Level set方法[11]等优化算法相耦合，用以解决更为复杂的问题，但为了简便且不失一般性，这里采用标准的最优化准则方法[3].优化过程中的设计变量采用下面的更新准则

其中为最小相对密度向量(为了避免奇性，没有取为0)，ρ＞0是移动的范围，n=1/2是控制灵敏度影响的系数.

本文所提算法的目标是在降低计算量的前提使得最优控制问题所得到的形状有更精确的边界，所以在网格自适应的过程中需要一个指示函数.这可以从计算区域的材料分布信息得到，而这些信息在计算流体问题的时候就已经求出.我们将用如下函数作为自适应网格的指示函数

其中在每个单元上的平均值，h为所有单元的最大边长.由(15)可以看出，粗网格单元被标记为1或0与该单元与边界的距离有关，确切的说，只有在边界附近的单元才被标记为1，其他单元被标记为0.

4 数值算例

4.1 算法

在上面理论分析的基础上，我们可以得到如下算法：

步骤1给定一个初始的设计变量并在整个计算区域上划分粗网格.令迭代次数k=1；

步骤2进行迭代，直到满足体积约束条件.第k次迭代由以下几步构成：

步骤2.1所有的网格根据指示函数(15)标记为1或者0，对其中标记为1的网格进行加密；

步骤2.2根据上一步得到的设计变量首先求解Stokes问题(3)–(6)得到速度场；

步骤2.3由(9)求出目标泛函的形状导数；

步骤2.4由(14)更新设计变量得到新的设计变量同时得到流体流动的新区域Ωk+1.相应的得到计算区域任一点的局部渗透性常数α(ϱ(x))；

步骤2.5如果需要的话，输出计算结果.

为了节省计算量，如果给出的初始设计变量过于简单或与可能的最优拓扑或形状差别较大，建议先对优化问题求解几次，然后进行网格自适应的求解.

4.2 算例

本小节给出了两个经典算例用以验证所提算法的可靠性和有效性，这两个算例也被很多研究者作为验证算例[6,7,12-16].

在两个算例中，固体区域都是假设不可渗透的，没有考虑液体本身的重力，没有液体流过的边界满足无滑移条件.在得到的最优化结果中，黑色区域表示ϱ=0(固体，液体不可渗透)，白色区域表示ϱ=1(液体).

算例1在第一个算例中，我们考虑改变沉浸在液体中物体的形状，使得流体流动的阻力(能量耗散)达到最小.据我们所知，这是最经典，也是最早被用于流体力学形状优化的算例[1].同样的问题曾经出现在很多流体力学形状或拓扑最优化的文献中[13-16].

问题的设计区域如图1所示.考虑Dirchlet问题，区域外边界上速度在x方向的分量为1，在y方向的分量为0，沉浸物体的边界是无滑移条件.经过我们的测试，最终结果对设计变量的初始值不敏感，所以我们假设在整个计算区域上的初始值都为0.4.得到最优形状时流体占整个计算区域的94%，雷诺数Re=1.这些假设与其他人工作中的设定是一致的.

图1:最小阻力问题的设计区域及边界条件

最优化结果如图2所示.两个图形分别为得到最优形状时对应的自适应网格(左)和最优形状(右).从结果可以看出，最优形状是一个橄榄球形，与文献[1]中的分析结果完全一致.

图2:算例1的计算结果

算例2在第二个算例中，我们考虑有两个入口和两个出口的情形，类似的问题同样被很多学者研究过[7,12,13,15,16].计算区域和边界条件如图3所示，有对称的两个入口和出口.在入口处速度在x方向的分量的大小为抛物状，出口处的压力为0，其它边界上面满足无滑移条件.区域的纵横比为1:1.优化的目标是使得流体在整个计算区域上面的能量耗散最小，同时满足最终流体所占区域为整个计算区域的40%.本算例的最终结果同样对设计变量的初始值不敏感，所以我们和上一算例一样假设在整个计算区域上的初始值都为0.4.

最终的结果如图4所示.左边为最优化所得拓扑(形状)相对应的自适应网格，右边最优形状.可以看出本算法所得结果要优于以前文献中所得到的结果.同样，本算例也说明了流体力学形状最优控制的重要性，因为我们在设计初期是无法确定最终优化结果的拓扑或形状的.

图3:算例2的设计区域及边界条件

图4:算例2的计算结果

5 结论

本文给出了一种求解流体力学形状最优控制问题的新算法，该算法非常容易实施.并且由于运算主要集中在流体流动的内部边界附近，从而在达到同样精度的情况下可以很大的节省计算量.尽管最优化算法采用的是最优性准则方法，并且只是以Stokes问题为例进行讨论，本算法完全可以采用其他更好的优化算法用以求解更为复杂的流体力学最优形状控制问题.提供的数值算例说明了算法的高效性和稳定性.

参考文献：

[1]Pironneau O.On optimal profiles in Stokes flow[J].Journal of Fluid Mechanics,1973,59(1):117-128

[2]Pironneau O.On optimum design in fluid mechanics[J].Journal of Fluid Mechanics,1974,64(1):97-110

[3]Bendsøe M P,Sigmund O.Topology Optimization,Theory,Methods,and Applications[M].New York:Springer-Verlag,2003

[4]Thévenin D,Janiga G.Optimization and Computational Fluid Dynamics[M].Berlin:Springer-Verlag,2008

[5]Plotnikov P,Sokolowski J.Compressible Navier-Stokes Equations:Theory and Shape Optimization[M].Birkhuser:Springer Basel,2012

[6]Mohammadi B,Pironneau O.Applied Shape Optimization for Fluids(2nd Edition)[M].Oxford:Oxford University Press,2010

[7]Borrvall T,Petersson J.Topology optimization of fluid in Stokes flow[J].International Journal for Numerical Methods in Fluid,2003,41(1):77-107

[8]Bendsøe M P,Kikuchi N.Generating optimal topologies in structural design using a homogenization method[J].Computer Methods in Applied Mechanics and Engineering,1988,71(2):197-224

[9]Bendsøe M P.Optimal shape design as a material distribution problem[J].Structural and Multidisciplinary Optimization,1989,1(4):193-202

[10]Svanberg K.The method of moving asymptotes—a new method for structural optimization[J].International Journal of Numerical Methods in Engineering,1987,24(2):359-373

[11]Osher S,Sethian J A.Level Set Methods and Dynamic Implicit Surface[M].New York:Springer-Verlag,2002

[12]Olesen L H,Okkels F,Bruus H.A high-level programming-language implementation of optimization applied to steady-state Navier-Stokes flow[J].International Journal of Numerical Methods in Engineering,2006,65(7):975-1001

[13]Challis V J,Guest J K.Level set topology optimization of fluids in Stokes flow[J].International Journal of Numerical Methods in Engineering,2009,79(10):1284-1308

[14]Duan X B,Ma Y C,Zhang R.Optimal shape control of fluid flow using variational level set method[J].Physics Letters A,2008,372(9):1374-1379

[15]Guest J K,Prevost J H.Optimization of creeping fluid flows using a Darcy-Stokes finite element[J].International Journal of Numerical Methods in Engineering,2006,66(3):461-484

[16]Pingen G,Evgrafov A,Maute K.Topology optimization of flow domains using the lattice Boltzmann method[J].Structural and Multidisciplinary Optimization,2007,34(6):507-524