流体拓扑优化的参数选择与分析

2018-08-14董馨刘小民

西安交通大学学报 2018年8期

董馨, 刘小民

(西安交通大学能源与动力工程学院, 710049, 西安)

水平集方法具有优秀的运动界面捕捉能力[1]。任晓辉和封建湖提出了一种不用求解复杂Hamilton-Jacobi方程的水平集方法,在对Heaviside函数正则化处理中,基于获得的形状导数和拓扑导数的相关信息加快收敛速度[2];段献葆和秦新强结合了传统的形状灵敏度分析与改进的水平集方法,发展了一种无需网格重新剖分的水平集方法,减少了运算工作量[3];张彬和刘小民提出了一种可以保持界面位置不动的新型隐式重新初始化的水平集方法,这种新型的水平集方法只需使时间步长满足原始的CFL条件,保证了界面附近水平集函数节点值符号不变,从而实现高效的优化[4];Dunning等基于水平集拓扑优化方法发展了适用于非结构化三维网格的快速行进法和迎风格式,增加水平集方法的鲁棒性和高效性[5];Duan等提出了一种自适应网格方法来求解不可压缩Navier-Stokes的拓扑优化问题,将隐式表示的材料分布信息作为设计变量,并将计算量集中在界面附近,从而显著降低了计算代价[6]。Shintaro等提出了一种基于水平集方法的自由灰度拓扑优化方法,在结构边界处运用可适应网格,它通过移动欧拉网格的节点来保持水平集函数,从而改善网格质量[7]。但是他们没有对水平集函数演化中的参数和时间步长参数进行选择分析和研究,在保证优化精度的前提下,合理的参数可使改进的水平集方法实现高效计算,充分发挥该方法的优点。

通过对改进的水平集方法[8]中的参数进行选择和分析,研究优化过程中迭代和收敛参数选择对优化进程的影响,使计算能够快速、稳定收敛,达到最优,同时保证优化结果的精度。改进的水平集方法采用了保持界面位置不动的新型隐式重新初始化水平集方法,结合无需样条参数化网格重新划分的水平集优化方法。

1 数学模型

Ω代表流体域,它是一个具有连续边界Γ=∂Ω的开区域,建立二维不可压缩Stokes方程如下

式中:u为流体的流动速度;p为流体压强;f为体积力;ν为动力黏性系数。ΓS是流体与固体间的交界面;ΓD是除ΓS外的Dirichlet边界,有Γ=ΓS∪ΓD,u0是边界ΓD上已知的速度分布。

本文研究的优化问题可描述为

(2)

式中:u满足状态式(1),D是优化问题的工作区域。

水平集函数φ(x,t)满足

水平集函数对运动界面进行追踪的基本思想如下:用高一维空间中的水平集函数去描述低一维空间中的运动边界,那么水平集函数的零等值线即为所追踪的运动边界。

水平集函数演化过程中区域界面始终表达为

φ(x,t)=0

(4)

对上式两端同时对t求物质导数,可得

(5)

取Vn=x′(t)为界面处的水平集法向速度,水平集演化方程为

(6)

2 数值方法

2.1 优化步骤

基于改进水平集方法的流体拓扑优化执行过程如下,其中n代表第n次优化迭代。

步骤1初始化原始流体区域,在设计区域内将水平集函数初始化为符号距离函数;

步骤2对流体区域进行网格重新划分;

步骤3求解流动控制方程和优化共轭方程;

步骤4求解水平集法向速度Vn;

步骤5将Vn由界面扩展到整个设计区域;

步骤6执行面积约束;

步骤7求解水平集方程,实现水平集函数演化;

步骤8将水平集函数重新初始化为符号距离函数。

2.2 时间步长

在水平集函数求解过程中,为了防止数值振荡,本文方程的离散采用一阶Essentially Non-Oscillatory(ENO)格式[9-10]。下面给出上述二维问题的离散形式

式中:x和y分别代表网格的横向和纵向坐标。下标i代表x方向第i个网格节点,下标j代表y方向第j个网格节点。

(8)

其中x+=max(x,0)

式中:S(φ)是符号函数;ε是光滑参数,(εi,j)2采用中心差分进行数值离散,离散形式为

其中时间步长满足下式

2.3 水平集函数重新初始化

隐式重新初始化方法是通过迭代求解下式的偏微分方程来实现的[11-12]。

φ(x,0)=φ0(x)

(15)

为了数值求解的需要,对S(φ0)光滑化

ε的取值为

(17)

3 数值算例

图1 拓扑优化的工作区域与边界条件

图2给出了翼型优化达到最优时的流线分布情况,从图2中可以看出本文结果与文献[15]结果吻合良好。

(a)γ=0.8 (b)γ=0.85

(c)γ=0.9 (d)γ=0.95图2 最优形状及其流线分布

3.1 水平集函数进化次数的分析

在水平集函数演化过程中,每次优化迭代时进行多次水平集函数进化会带来一定的偏差,因此进行水平集函数重新初始化,以消除水平集函数在演化过程中偏离符号距离函数的现象。

引入参数m定义为每次迭代时水平集函数的进化次数。分别选取m值为10、8、6、4、2和1这6组数据进行研究,如图3、图4所示,分别为流体区面积比和其波动值随着迭代步数n的变化。为了实现在优化前期加速计算,在优化需要精确结果的地方提高准确度,使计算程序能够达到自动选择合理参数的目的,将水平集函数进化次数m与迭代步数n的关系分别采用分段线性函数(式(18))、三次函数(式(19))和幂函数(式(20))描述如下

(a)m取单一值

(b)m取函数图3 不同m取值流体区面积比随迭代步数的变化

图4 不同m取值时流体区面积比的波动值随迭代步数的变化

由图3可看出,流体区面积由初始面积逐渐上升至约束面积,之后基本保持不变直到迭代收敛。当m取值分别为10和8时,迭代收敛后在面积比为0.95附近呈现锯齿状的波动,原因是较大的m值导致结构参数(面积比和目标函数)变化较大,面积增加速度太快而导致流场面积大于目标面积,较大面积对应着较小目标函数,约束的施加使得面积反向变化,面积减小,目标函数增大。继续演化,面积又增大到大于目标面积,如此反复,形成了锯齿状演化曲线。参数m分别取值为6、4、2和1时都可以达到稳定收敛,收敛速度依次下降,收敛时的迭代步数与计算时间见表1。

表1 m取单一值时的收敛指标

图4中面积比的波动值,即迭代第n步与n-1步的面积比差值的绝对值。m取10和8时的波动值较大,且收敛之后的波动也较大。而其他参数的波动值都随着迭代步数大体呈下降趋势,尤其在收敛时有近乎竖直的下降梯度,收敛之后波动值接近0并呈现稳定趋势。

目标函数随迭代步数的收敛如图5所示,目标函数值随着迭代次数的增加而下降,达到0.85时几乎不再随着迭代次数发生变化,此时翼型绕流能量耗散值最小,翼型结构为最优拓扑结构。当m取值为10和8时,收敛到最优结果时产生了波动,收敛情况不理想。线性函数、三次函数与幂函数均在第13步时达到收敛,其中三次函数在优化过程中收敛速度最快,这是由于它在收敛前期有较大的m,收敛后m逐渐减小为1,使得收敛趋于稳定。

(a)m取单一值

(b)m取函数图5 不同m取值下目标函数随迭代步数的变化

3.2 时间步长的分析

时间步长的选择对于程序的收敛速度和计算速度起着至关重要的作用。时间步长过大会导致不满足courant-Friedrichs-Lewy(CFL)条件,程序计算不稳定,最后达不到收敛指标;而时间步长过小又会导致计算耗时、收敛速度太慢。因此,选取合适的时间步长尤为重要。

本文计算中时间步长设置为:Δt=kmin(Δx,Δy),其中k为常数,为满足CFL条件,k≤0.5。k分别选取0.5、0.2和0.1时进行计算研究,流体区面积比和波动值随迭代步数的关系分别如图6和图7所示,迭代收敛步数和计算时间见表2。为了实现优化前期快速达到最优值,优化后期稳定收敛,使计算程序能够自动选择合理参数,将描述时间步长的参数k与迭代步数n的关系分别表达为线性函数和指数函数,公式如下

k=-0.006n+0.55,n≤99

(21)

k=0.833e-0.03n

(22)

由图6可看出,流体区面积由初始面积逐渐上升,之后由于加入了面积约束,保持面积比基本不变直到迭代收敛。图7中的波动值大体呈现下降趋势,最后接近0稳定收敛,随着k取值的减小,波动更平缓但收敛更慢。

(a)k取单一值

(b)k取函数图6 不同k取值流体区面积比随迭代步数的变化

图7 不同k取值时流体区面积比的波动值随迭代步数的变化

图8为目标函数随迭代步数的变化,目标函数随着迭代次数的增加而下降,最后逐渐稳定保持目标函数不变,能量耗散数为0.85,此时可得到最优拓扑结构;线性函数与指数函数分别在第18步与第14步达到收敛,由于它们随着迭代步数变化,因此在计算前期可以达到快速收敛的目的,在后期防止因为时间步长过大引起的数值不稳定现象。

(a)k取单一值

(b)k取函数图8 不同k取值下目标函数值随迭代步数的变化

k迭代步数迭代时间/s0.5185 4150.24314 7100.18227 860

(a)参数改进后的流体区面积比

(b)参数改进后的目标函数图9 参数改进后的面积比及目标函数随迭代步数的变化

综合选取两个参数对优化过程的影响规律,流体区面积比及目标函数随迭代步数的变化情况如图9所示。从图9a中可以看出,流体区面积随迭代步数的增加逐渐上升至面积约束值,达到稳定收敛,收敛速度最慢的为只改进k的方案,其次为只改进m的方案,最快的为同时改进m和k的选取方式。由图9b可得,目标函数在优化前期下降速度最快,接近收敛时目标函数变化逐渐减小,直至稳定收敛。同时改进m和k的选取方式是在优化之初不需要精确的收敛数值的情况下,增加每一次迭代时水平集函数的进化次数从而提高计算速度,同时时间步长增加使迭代速度加快,随着计算进行,由于需要获取精确的收敛数值,迭代速度减慢。改进前与改进后的收敛指标见表3。参数选取方式如下:改进前m=1,k=0.1;只改进m的选取方式为三次函数,k=0.1不变;只改进k取指数函数,m=1;改进m选取为三次函数,k选取指数函数。经过合理选择m和k,收敛速度提高96%。