段献葆， 党 妍， 秦 玲

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

由于工业领域的迫切需要, 自Pironneau[1]开创性的工作以来，流体中的形状优化问题就吸引了很多工程师和数学家的关注。大量的实例表明，在工业设计的早期阶段，选择合适的形状或拓扑结构可以大大地节省成本。在过去的几十年里，流体动力学中的形状优化方法已取得了很大进展，许多方法已经被引入并应用到工业设计问题中[2-7]。流体形状优化的细节和更多内容参见文献[8]。

Osher 等[9]提出的水平集方法(level set method)最初用来数值跟踪或捕捉流体动力学中的动态界面和自由边界。水平集方法不直接跟踪界面，可以很容易地处理演化过程中的拓扑变化。水平集方法目前已经被广泛地应用到与形状优化相关的领域，如材料科学、数字图像处理、几何光学和生物学等[10-15]。在经典的水平集方法中，水平集函数在演化的过程中需要保持为距离函数，但这一性质是很难保持的。通常的处理方法是在演化的过程中对水平集函数进行重新初始化，而这一过程往往是非常复杂费时的，且对边界有很强的依赖性。为了解决这个问题，Chan 等[16]提出了变分水平集方法，该方法在水平集函数演化的过程中不再需要重新初始化，从而可以进一步提高计算效率。

在经典的水平集方法中，水平集函数是定义在整个求解区域、采用均匀网格的有限差分格式实现，这使得求得较高精度数值解的时候会增加计算量。而很多时候人们关心的只是水平集函数的零水平集，为了克服这一局限性，提出了窄带法[17]和局部水平集等[18]方法。此外，人们也提出了自适应方法来提高求解效率和精度[19]。

本工作的主要目的是研究流体流动引起的形状优化问题。为了提高优化过程的计算效率，提出了一种几何自适应网格方法。该方法具有灵活性、通用性等优点，易于与现有的有限元方法软件集成，非常容易实施。

1 控制方程和形状优化问题

设计算区域Ω是R2中一个具有Lipschitz 连续边界Γ:=∂Ω的有界连通区域。以偏微分方程为状态方程的形状优化问题的一般形式如下：

使得

式中：γ(x)(0 ≤γ(x)≤1)是优化变量，用来控制当地介质的渗透性；VD是设计领域的体积；β是流体可能占据的最大区域。

式中：φ(x)是水平集函数；γ=0表示固体材料，γ=1表示液体材料。方程(1)中目标泛函J是关于状态变量u和设计变量γ的函数，对于阻力最小化问题，具有如下形式。

式中：ν=1/Re(Re是雷诺数)是运动黏性系数。

式(2)～(3)构成了形状优化问题中的Dirichlet 型状态约束。若状态方程为Stokes 方程，则对应的Dirich-let 边值问题为

式中：函数u=(u1,u2)；g、f和标量函数p分别表示流体的速度、边界上规定的流速、外力和压力。假设在理想多孔介质中流动的流体满足Darcy 定律，即外力f和流体速度成正比。因此，f=P(x)u，P(x)是渗透介质在位置x的渗透率。基于不可压缩性条件，利用散度定理容易看出函数g必须满足相容性条件，

式中：n为单位外法向量。

与Borrvall 等[20]的研究相同，本工作选择一个凸插值函数表示流体的渗透性参数P，

式中：0 ≤γ≤1,q是一个实的正参数，用于调整P(γ)的形状。对于实际问题来说，不能渗透的固壁需要Pmax=∞, 但在实际计算中只能选择一个有限值。本工作取Pmax=104，Pmin=0。

2 灵敏度分析和变分水平集方法

2.1 灵敏度分析

为了对水平集函数进行演化，需要得到优化问题中目标函数(式(1))的灵敏度分析信息。设Ωα ⊂R2是区域Ω的一个扰动，Ωα=(Id+α)(Ω):={x+α(x),其中α是正的小参数，使得x ∈Ω}，δΩ=ΩαΩ，边界为Γα={x+α(x)n(x),∀x ∈Γ:=∂Ω}。

设w=u(Ωα)-u(Ω)≡uα-u,是u在区域Ωα的零延拓。对于向量α, 开集Ω(α)是区域Ω小的扰动，(Id+α)是R2中一个微分同胚。

目标泛函J(Ω)关于区域Ω的形状导数为W1,∞(R2,R2)中的Fr´echet 导数，

式中：DJ(Ω)是一个连续线性函数。

引理1[8]考虑Sα关于边界S的一个小扰动

式中：α是关于x(s)∈S的函数(s曲线横坐标)；λ是一个趋向于0 的正数。设Ωα=Ω(Sα)，则对于S附近的任意连续函数f有

定理1如果约束是Stokes 问题(式(6)～(8))，则目标函数J(u(γ),γ)关于区域Ω的形状导数为

式中：n为局部单位外法向量。

证明 对于目标泛函

有

由引理1，可得

因此，

若∇u是光滑的，则有

将式(6)乘以试验函数w，并利用Green 公式，可得

因此，

由泰勒展开，可得

而u｜Γ=0，所以有

设s表示切线分量，则有

把式(11)和(12)代入式(10)，可得

因此，本工作使用

作为水平集函数的演化速度。

2.2 变分水平集方法

水平集方法的基本思想是将低维界面表示为一个函数在高维空间中的零水平集。对于开集Ω(x,t)，x ∈R2，t ∈R+，可以通过

定义R2×R+上的函数φ(x,t)来确定区域Ω(x,t)。边界Γ(x,t) :=∂Ω(x,t)为零水平集，

在区域Ω(x,t)外部有φ(x,t)＞0，即水平集函数φ(x,t)是满足

的函数D为字体区域。

对方程(15)两边关于时间进行微分，并应用链式法则，可得到Hamilton-Jacobi 型方程，

式中：V= dφ(x,t)/dt是移动边界的法向速度。式(16)即为水平集函数的运动方程。在水平集方法中，单位外法向量n和表面曲率κ分别为

在经典的水平集方法中，水平集函数在演化过程中需要保持为符号距离函数的形式，而重新初始化的过程相当复杂。水平集函数φ(x,t)是否是一个距离函数当且仅当它满足｜∇φ(x,t)｜= 1。所以在文献[16]中积分

被用来描述水平集函数φ(x,t)在Ω ∈R2中与距离函数的密切程度。假设水平集函数φ(x,t)的法向导数在边界上为0，

则φ(x,t)在ψ方向的Gateaux 导数为

在本算法中，水平集函数将由基于形状灵敏度的Hamilton-Jacobi 方程求解，

式中：λ是平衡(Δφ(x,t)-κ)影响的一个正参数。方程(18)的第二项能使水平集函数的零水平集趋向区域的边界。与经典的水平集方法相比，使用方程(18)求解水平集函数时，重新初始化的过程在优化过程中可以忽略。

3 自适应网格方法

本工作通过一个数值阻尼项将参考域中的无滑移条件正则化，该数值阻尼项对于固体或流体速度较慢的材料来说很大的，而对于流体材料是接近于0。目标函数被重新表述为流体通过多孔介质时的功率耗散。解用0-1 拓扑表示，不包含中间情形，即为0 的区域是流体，为1的区域是固体，通过固体材料的流动是由Darcy 定律控制的。

基于上述分析，下面给出自适应网格法的主要步骤。

（1）对整个计算区域进行均匀剖分，在该网格上求解水平集函数。粗网格水平集函数的取值是对粗网格进行细化的关键。使用以下水平集函数作为细化指示子，

式中：φ(x)是水平集函数；h是到边界距离的支撑（水平集函数的零水平集）。可以看出，粗网格单元是否标记为1 还是0 取决于粗网格单元到界面的距离。

（2）对细化指标标记为1 的粗网格进一步划分为均匀的细网格。优化过程中，粗网格内均匀细网格的生成是自动完成的。

（3）对优化问题进行求解。

自适应水平集方法的具体步骤如下。

（1）给出初始形状Ω0和初始化水平集函数φ0(x)构建统一的粗网格的基础上最初的形状。设迭代次数k=1。

（2）进行迭代，直到达到体积约束。第k次迭代包括以下步骤：①所有粗网格按细化准则(式(19))进行分类，标记为1 的网格进一步划分为均匀细网格，细化后重新排列单元和节点编号；②由所得数据可以得到γk，也即关于材料的分布，从而可以得到计算区域，利用细化网格求解Stokes 问题(式(6)～(8))，得到速度场；③由式(9)求出灵敏度分析结果；④通过求解粗网格上的水平集方程(18)更新水平集函数，由式(5)可以得到新的区域Ωk+1；⑤若需要，则进行可视化处理。

4 数值算例

本工作采用一个经典的二维形状优化问题验证所提算法的可行性和有效性，该算例可以在很多流体形状优化的文献中找到[1-2,15,21-22]。假设固体区域是不能渗透的并且忽略外力项，固壁边界采用无滑移边界条件。在所有数据中，灰度图片显示最优结果，黑色对应的设计值γ= 0(固体)，白色对应γ= 1(流体)。考虑的是浸入外部Stokes 流的物体形状的优化问题，目标是使物体引起的阻力最小。最小阻力问题的设计域如图1 所示，其中状态约束是Stokes方程，上、下及入流边界上速度在x方向上是1，在y方向上是0。

图1 最小阻力问题的设计域[1]Fig.1 Design domain and boundary conditions for the minimum drag example[1]

在本算例中，Stokes 问题(式(6)～(8))的求解采用有限元方法，在得到灵敏度分析结果后，通过插值映射到求解Hamilton-Jacobi 方程(式(18))的有限差分网格上。而在求解方程(18)的过程中，需要满足CFL(Courant-Friedrichs-Lewy) 条件。本工作所采用的时间步长Δt需要满足Δt≤Δx/｜V ｜max，其中｜V ｜max是由式(13)所得到结果在整个求解区域内的最大值。在求解过程中，对于水平集函数的演化没有重新初始化的过程，这也相应地减少了很多求解过程中的迭代，进一步提升了计算效率。本工作使用多种不同的初始形状和体积约束进行优化，发现结果对初始形状和体积约束不敏感。本算例采用的初始形状是一个简单的圆形，流体体积为20%。最终的液体体积限制在80%和94%，结果在Re= 1 时得到。所有参数都与已有文献保持一致。优化结果如图2 所示，其中最优形状对应的自适应网格见图(a)和(b)，相应的水平集函数等值线图见(c)和(d)。水平集函数的零水平集对应内部实心物体的边界，最优形状如图2(e)和(f)所示。可以看出，该优化过程得到了类似橄榄球的形状，这符合文献[1]的理论分析得到的形状优化结果，与文献[15, 21-22]所得到的数值模拟结果也完全吻合。

图2 最小阻力问题的最优结果Fig.2 Optimal results of the minimum drag example

本算例的初始网格所用三角元1 332 个，节点数是714，求解水平集方程用的差分网格是20×20. 求解时间是516.24 s。为获得相同的分辨率，本工作也通过均匀的细化网格（所用三角元13 126 个，节点6 108 个，差分网格是75×75）求解了该问题，最终结果一致，但运行时间1 652.96 s，约为自适应网格方法的3.2 倍，表明本算法的计算效率得到非常大的提升。

5 结束语

本工作给出了一种求解流体力学中形状优化问题的自适应水平集方法。该方法计算量小，数值实现简单。使用这种自适应方法，计算量主要集中在界面(边界)附近，因此，与全区域均匀网格相比，达到了相同的分辨率时的计算成本大大降低。虽然本工作的讨论是以二维Stokes 问题为例，但是该方法可以推广到求解复杂流体流动的形状或拓扑优化问题。