一种基于自动微分的可压缩湍流流体拓扑优化方法

1.本发明属于计算流体力学中流体拓扑优化技术领域，具体是涉及一种基于自动微分的可压缩湍流流体拓扑优化方法。


背景技术：

2.可压缩性是所有流体的固有特性之一，通常在水和其他液体中，会忽略这一特性，而类似空气这类的气体，在高速高压情况下，可压缩性是至关重要的考虑因素之一。高速可压缩流体流动广泛存在于航空航天领域，在诸如发动机流道部件、飞行器外形、机翼型面、涡轮机械等气动设计中，都必须考虑流体的可压缩性。同时，随着计算机技术的发展，针对可压缩流的优化方法接连出现，cummings et al.([1]cummings r m,yang h t,oh y h.supersonic,turbulent flow computation and drag optimization for axisymmetric afterbodies[j].computers&fluids,1995,24(4):487-507)提出了针对可压缩流体的尺寸优化，burgreen and baysal([2]burgreen g w.three-dimensional aerodynamic shape optimization using discrete sensitivity analysis.[d].old dominion university.1994)则提出了形状优化。
[0003]
流体拓扑优化方法有别于尺寸优化和形状优化，其可实现优化对象的拓扑构型改变，在初步概念设计阶段具有很大的优势。borrval和petersson([3]borrvall t,petersson j.topology optimization of fluids in stokes flow[j].international journal for numerical methods in fluids,2003,41(1):77-107)首次提出了流体拓扑优化方法，并将其应用于流体通道的拓扑优化问题中。自此，流体拓扑优化方法被广泛应用于各种流动。而目前流体拓扑优化主要集中在不可压缩流动领域。最近，lfn s
á
et al.([4]lfn s
á
,okubo c m,silva e.topology optimization of subsonic compressible flows[j].structural and multidisciplinary optimization,2021(6))实现了亚声速可压缩流体拓扑优化，但未对自动微分实现离散伴随具体过程进行说明，且未考虑湍流的影响。而众多科研及工程领域中，机翼翼型、涡轮叶珊、发动机内流通道等工作环境都是高速可压缩湍流流动。因此，亟需发展考虑湍流模型的可压缩流体拓扑优化方法。
[0004]
流体拓扑优化中灵敏度求解通常采用连续或离散伴随方法，a.carnariuset al.([5]carnarius a,thiele f,oezkaya e,et al.adjoint approaches for optimal flow control[c]//5th flow control conference.2010)介绍了这两种方法的优缺点。对于复杂的、采用分离和迭代求解方法的cfd求解器，手动推导离散伴随方程是一项重要且容易出错的任务。自动微分方法是一种利用计算机程序对函数梯度进行精确求解的技术，具有灵活性强，计算精确等特点。r.gieringet al.([6]r.giering,t.kaminski,recipes for adjoint code construction,acm trans.math.softw.24(4)(1998)437
–
474)详细介绍自动微分方法。cetin b.dilgen et al.([7]dilgen cb,dilgen s b,fuhrman d r,et al.topology optimization of turbulent flows[j].computer methods in applied mechanics and engineering,2018,331(apr.1):363-393)提出了基于自动微分法的不可
压湍流拓扑优化。利用自动微分方法来求解可压缩流体伴随方程是本发明要讨论的，基于此，本发明提出了一种基于自动微分且考虑湍流模型的可压缩流体拓扑优化方法。


技术实现要素：

[0005]
本发明的目的在于提供可实现针对可压缩流体的流体拓扑优化，并考虑湍流模型对优化的影响，通过自动微分方法构建伴随方程，利用自动微分法的特性将雅可比矩阵高效且精准显示生成，避免手动推导可压缩流体伴随方程灵活度不足、过程繁琐易出错等问题的一种基于自动微分的可压缩湍流流体拓扑优化方法。
[0006]
本发明包括以下步骤：
[0007]
步骤1：根据优化对象的实际需求，建立对应流道拓扑优化的几何模型，确定优化设计区域ω和流道的出入口位置；
[0008]
步骤2：获得拓扑优化对象的基本参数，基于可压缩流体控制方程构建cfd拓扑优化模型；
[0009]
步骤3：根据基本参数将优化目标与流体控制方程转换成弱形式并构建目标函数，构成非线性规划问题；
[0010]
步骤4：基于伴随方法推导灵敏度方程；
[0011]
步骤5：设定合适的设计变量，中间变量初始值，求解流体控制方程得到流场结果和目标函数的值，输出新流场结果；
[0012]
步骤6：为求解伴随方程，利用自动微分技术组建雅可比矩阵和梯度向量；
[0013]
步骤7：将步骤6得到的雅可比矩阵和梯度向量构建矩阵形式伴随方程，得到伴随乘子；
[0014]
步骤8：将步骤5得到的新流场结果与步骤7得到的伴随乘子带入步骤4的灵敏度方程，求解灵敏度；
[0015]
步骤9：使用mma数值优化方法，结合已知拓扑优化的非线性规划问题搭建数学模型；
[0016]
步骤10：使用mma算法对流道进行优化求解，更新设计变量，若迭代已收敛，即得到最优解，输出最优二维拓扑构型；否则跳转至步骤5，并将流场结果设置为新初始值。
[0017]
在步骤2中，所述获得拓扑优化对象的基本参数，主要包括：优化目标、流体属性、来流条件、边界条件、热传导参数；所述基于可压缩流体控制方程构建cfd拓扑优化模型的具体步骤可为：
[0018]
(1)基于变密度拓扑优化方法指定材料密度为设计变量γ，在流体控制方程添加源项f，表示多孔介质所带来的阻力，其与设计变量的函数如下：
[0019]
f＝-au
[0020]
α(γ)＝α
min
+(α
max-α
min
)q(1-γ)/(q+γ)
[0021]
其中，α为渗透系数，且α值越大，流体所受多孔介质造成的阻力越大；u为流体速度场，q为惩罚参数，γ为材料密度。
[0022]
(2)基于可压缩流体控制方程构建cfd拓扑优化模型：
[0023]
min:j(u,p,t,γ)；
[0024][0025][0026][0027]
其中：
[0028][0029]
0≤γ≤1
[0030]
其中，e为能量密度，ρ为质量密度，u为流体速度场，p为流体压力场，t为流体温度场，k为导热系数，为粘性张量，μ为动力粘度，μ
t
为湍流模型求解得到的湍流粘度，v
θ
为体积分数上限。
[0031]
在步骤3中，所述构建目标函数，构成非线性规划问题的具体步骤可为：根据基本参数将优化目标与流体控制方程转换成弱形式并构建目标函数将体积分数限制转换成约束条件其中，设计变量γ即函数自变量x，然后构成下述非线性规划问题：
[0032]
minimize f0(x)
[0033]
subject to fi(x)≤0，i＝1,
…
,m
[0034]
在步骤4中，所述基于伴随方法推导灵敏度方程如下：
[0035][0036][0037][0038]
其中，u为中间变量，包括u,p,t；
[0039]
令伴随方程得到灵敏度方程：
[0040]
在步骤6中，所述得雅可比矩阵和梯度向量的具体方法可为：
[0041]
(1)通过自动微分法组装具有高度稀疏结构的雅可比矩阵：其展开形式为：
[0042][0043]
其中：m为中间变量的个数，n为设计变量维度。
[0044]
(2)通过自动微分法组装梯度向量l，其值为目标函数对中间变量的偏导：
[0045][0046]
在步骤7中，所述矩阵形式伴随方程如下：
[0047][0048]
求解方程后得到伴随乘子的值，输出伴随乘子λ。
[0049]
在步骤8中，所述求解灵敏度如下式；
[0050][0051]
在步骤9中，所述数学模型的关键参数包括流场结果u、p、t，伴随乘子λ，灵敏度以及目标函数与约束条件的一阶、二阶微分表达式等。
[0052]
在步骤10中，判断迭代是否已收敛，以更新前后每个网格单元最大的相对设计变量改变值与设定迭代终止条件对比作为收敛依据。
[0053]
本发明针对可压缩流体进行流体拓扑优化，在流体控制方程中考虑了湍流模型，并采用自动微分方法求解灵敏度分析中的伴随方程。本发明可压缩navier-stokes方程与基于变密度材料的惩罚方程耦合；并根据设定目标函数及体积限制等约束条件进行流体拓扑优化，其中优化算法中的梯度通过自动微分实现离散伴随，进而求解。本发明的优势在于利用了自动微分法实现并简化考虑湍流模型下可压缩流体拓扑优化的灵敏度求解，降低软件开发和维护者的负担，提高流体拓扑优化应用程序的可维护性和灵活性。这种基于自动微分的可压缩湍流流体拓扑优化方法为航空航天等工程领域的流体拓扑优化应用提供求解手段。本发明对可压缩层流、湍流的流体拓扑优化均适用。
附图说明
[0054]
图1是一种基于自动微分求解伴随灵敏度的可压缩流体拓扑优化方法具体实施流程图。
[0055]
图2是本发明实施例提供的可压缩湍流二维直通流道的拓扑优化设计示意图。
[0056]
图3是本发明实施例优化过程中目标函数和体积约束收敛图。其中，(a)为目标函
数收敛图；(b)为体积分数收敛图。
[0057]
图4是本发明实施例求解后的最优流道拓扑构型图。
具体实施方式
[0058]
下面将结合附图对本发明实施例的具体实施方式展开详细描述，出于解释而非限制性的目的，下面描述中，阐述了具体细节，以帮助全面地理解本发明。
[0059]
以可压缩湍流在简易二维直通道中优化为实施例，基于自动微分求解伴随灵敏度的可压缩流体拓扑优化方法主要构建流程详见图1。本实施例主要方法如下：
[0060]
1、根据直通道实施例要求，构建一个方型设计域ω，如图2。左侧为入口，右侧为出口，令入口宽度为出口的两倍；
[0061]
2、设定拓扑优化对象的基本参数，主要包括：基于变密度拓扑优化方法指定材料密度为设计变量的γ、优化目标、流体属性、来流条件、边界条件等。其中优化目标为流动损耗，流体属性为空气，边界条件为无滑移和热绝缘边界：
[0062]
设计变量：材料密度γ
[0063]
中间变量：u，p，t
[0064]
渗透系数：
[0065]
雷诺数：
[0066]
达西数：
[0067]
目标函数：其中，μ为动力粘度，μ
t
为湍流模型求解得到的湍流粘度；
[0068]
约束条件：其中，v
θ
是体积分数的人为限制值。
[0069]
在本实施例中雷诺数re＝2.6
×
106，入口速度为抛物线速度，且u
max
＝167m/s，马赫数ma＝0.5，体积约束上限v
θ
＝0.8。
[0070]
3、流体拓扑优化的约束条件分为等式约束条件和不等式约束条件，其中，等式约束条件包括添加源项f＝-au后的流体控制方程：
[0071][0072][0073]
[0074]
其中：
[0075]
不等式约束条件包括体积分数上限：综上所述，将优化目标与流体控制方程转换成弱形式并构建目标函数将体积分数限制转换成约束条件其中，设计变量γ即函数自变量x，然后构成下述非线性规划问题：
[0076]
minimize f0(x)
[0077]
subject to fi(x)≤0，i＝1,
…
,m
[0078]
4、基于伴随方法推导灵敏度方程：
[0079][0080][0081][0082]
其中，u为中间变量，包括u,p,t；
[0083]
令伴随方程得到灵敏度方程：
[0084]
5、设定合适的设计变量，中间变量初始值，求解流体控制方程得到流场结果和目标函数的值，输出流体速度场u，p，j(u,p,t,γ)；
[0085]
6、为求解伴随方程，使用自动微分技术简化的求解：
[0086]
a)通过自动微分法组装具有高度稀疏结构的雅可比矩阵：其展开形
[0087]
式为：
[0088][0089]
其中：m为中间变量的个数，n为设计变量维度。
[0090]
b)通过自动微分法组装梯度向量l，其值为目标函数对中间变量的偏导：
[0091][0092]
7、将第6步得到的雅可比矩阵k和梯度向量l组建矩阵形式伴随方程：
[0093][0094]
求解方程得到伴随乘子的值，输出λ；
[0095]
8、将第4步得到的流场结果u与第7步得到的伴随乘子λ带入第5步的灵敏度方程，求解灵敏度；
[0096][0097]
9、使用mma数值优化方法，结合已知拓扑优化的非线性规划问题搭建数学模型。数学模型关键参数包括上述步骤求解得到的流场结果：u、p、t，伴随乘子：λ，灵敏度：以及目标函数与约束条件的一阶、二阶微分表达式等。
[0098]
10、通过mma算法对直通流道进行优化求解，不断更新设计变量，以更新前后每个网格单元最大的相对设计变量改变值与设定迭代终止条件对比作为收敛依据，若迭代未收敛，则跳转至第5步，并将流场结果设置为新初始值。直至迭代收敛，得到最优解，停止计算，输出最低能量损耗的二维直通流道拓扑构型。优化过程中，惩罚参数q与α
max
决定了优化结果的好坏，在q与α
max
选取不当时，会出现流固界限不清淅的问题。图3是本发明实施例q＝0.1，α
max
＝10000参数设置下的目标函数和体积约束收敛图。从图中可以看出，优化过程中前10步的变化幅度较大，从第10步到第60步逐渐收敛。其中，目标函数收敛图(a)展示了不同迭代步数下的拓扑构型。图4是最终优化收敛后的拓扑构型结果。
[0099]
本发明提出一种基于自动微分的可压缩流体拓扑优化方法，其关键在于利用自动微分技术求解灵敏度分析中的伴随方程。本发明不仅实现了考虑湍流影响的可压缩流体拓扑优化，并且其自动微分实现离散伴随的特点使得无需根据不同研究对象的目标函数和边界条件重复推导，相比繁琐的手动推导伴随方程方法具有显著的优越性。
[0100]
显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定。这里无需也无法对所有的实例予以穷举。任何依据本发明思路下的修改和改进，均视为本发明保护范围。