CN113268910A - 一种重力驱动的自然对流异型热沉结构拓扑优化方法 - Google Patents

Abstract

一种重力驱动的自然对流异型热沉结构拓扑优化方法，先建立重力驱动自然对流传热问题物理模型并简化，然后采用有限元方法对自然对流传热简化模型进行分析，以强化散热弱度、降低平均温度等提升结构散热性能的手段为优化目标，建立变密度法的优化模型，通过移动渐近线方法(MMA)对设计变量进行更新，将得到的设计变量分配到网格点上，得到网格点密度，通过构造水平集函数，实现散热翅片结构的显式边界表达；最后对优化后的散热翅片结构进行光滑圆整处理；本发明既能保证优化结构的形状和拓扑，又能实现优化结果的显式边界表达，有利于生产制造，得到能强化自然对流换热效果的优化结构，达到最佳冷却效果。

Description

一种重力驱动的自然对流异型热沉结构拓扑优化方法

技术领域

本发明属于散热翅片结构设计技术领域，特别是一种重力驱动的自然对流异型热沉结构拓扑优化方法。

背景技术

随着现代设备技术的发展，设备正向高集成度、高密度组装、高性能、高可靠性方向发展；设备的功耗越来越大，相应的热流密度也越来越大，高温使得设备的性能降低、可靠性下降；自然对流散热因其可靠性高、成本低、无噪声等优点被视为设备系统理想的散热方式；基于自然对流传热效应的翅片散热结构在工程中应用广泛，例如电子设备的冷却、大型变压器的冷却、反应堆冷却、航空航天设计等等；特别是对于需要室外长期稳定工作的特殊环境下这种散热方式无可替代；因此，研究重力驱动自然对流散热系统的翅片散热成为了设备进一步发展的关键。

重力驱动自然对流散热系统的翅片散热常规尺寸及形状优化方法对翅片结构设计通常受到既定拓扑形态的限制，设计自由度受限，拓扑优化能突破自由度的局限性，找到非直觉和非预期的设计；目前发展比较成熟的拓扑优化方法是变密度方法；虽然变密度法理论模型相对简单，易于操作，但是优化结果边界仍存在中间密度单元，不能生成光滑边界，不利于生产制造。

发明内容

为了克服上述现有技术的缺点，本发明的目的在于提供了一种重力驱动的自然对流异型热沉结构拓扑优化方法，在变密度法的基础上，对边界实现光滑化的显式边界优化结构。

为了达到上述目的，本发明采取的技术方案为：

一种重力驱动的自然对流异型热沉结构拓扑优化方法，包括以下步骤：

1)建立重力驱动自然对流传热问题物理模型以及模型简化；

2)建立有限元分析模型；

3)建立自然对流传热问题优化模型；

4)迭代更新驱使设计结构收敛。

所述的步骤1)具体为：

1.1)建立重力驱动自然对流传热问题物理模型：

根据实际自然对流传热物理问题，提取问题的几何模型以及边界条件；整个计算区域被分为设计域和流体域，其中设计域包含流体和固体，整个设计域被流体域包围，流体流动会影响设计域中流体和固体的分布，设计域中固体对应的区域即为散热翅片结构；在整个计算区域中建立自然对流传热问题控制方程组，通过系数来控制区域的性态，其动量方程、质量方程和能量方程分别为：

其中u是流体的速度场，p表示压力场，T表示温度场；ρ₀是参考密度；T₀是对应的参考温度；β为体积膨胀系数；μ是流体的动力学粘度；g为重力加速度向量；c_p为比热容；Q(x)为空间中的体积热源项；k(x)为空间变化的有效导热率，定义为：

其中k_f是流体的有效导热率；k_s是固体的有效导热率；

整个计算区域受重力作用的影响；边界上有指定的温度边界、热绝缘边界以及热流密度边界；所有的边界均为无滑移边界条件；考虑的均是定常的、不可压层流；

1.2)重力驱动自然对流传热问题物理模型简化：

假设对流项

被忽略，并且将μΔu表示为速度的线性表达式，重力驱动自然对流传热问题物理模型简化为：

其中

是新引入的材料参数，定义为：

其中

分别是流体和固体的材料参数；

在固体中设置新的材料参数

取∞，即

通过与商业软件的自然对流传热全阶模型控制方程组进行对比，将温度的最小二乘误差作为因变量，通过调整流体中材料参数

的值，观察简化模型的逼近程度；最小的温度最小二乘误差所对应的参数即为流体的材料参数

所述的步骤2)具体为：

对计算区域用四边形网格进行离散，采用有限元方法进行分析；对压力和温度采用一阶多项式插值，对控制方程的两边乘以测试函数w，在计算区域进行积分得到自然对流传热问题物理模型的变分形式为：

对测试函数w进行扰动，即

其中稳定化参数τ为：

其中h表示四边形单元沿x、y方向尺寸的最大值；

采用牛顿迭代的方法进行求解，得到自然对流传热问题的温度场、压力场以及速度场。

所述的步骤3)具体为：

3.1)优化模型建立：

针对步骤2)所得到的自然对流传热问题的温度场，考虑的优化目标是使施加热流密度处的平均温度最小，满足的约束条件为自然对流传热方程组、设计域内固体材料所占的比例、设计变量的取值范围；建立自然对流传热问题优化模型：

其中N_d是指计算区域中施加了热流密度处所对应的所有节点总数；N_e是设计区域内单元总数；

为设计变量的集合；x_i是设计区域内第i个单元上的伪密度；v_i是设计区域内第i个单元的体积；f是指设计区域内所能容许的最大固体体积分数；

是指设计区域Ω_d的体积；等式约束为自然对流传热控制方程组的矩阵形式，其中k₁₁,k₁₂,k₂₂分别为系数矩阵，f₁,f₂分别为右端项；

3.2)材料参数连续化处理：

当设计变量x_i＝0时表示流体；当设计变量x_i＝1时表示固体，即散热翅片结构；x_i∈[0，1]表示流体和固体的过渡阶段；设计变量的变化，反映了散热翅片结构的变化；将新的材料参数和导热率表示为插值形式：

其中

p_k分别表示材料参数和导热率的惩罚因子，随着迭代步数的增加而增加；

3.3)敏度分析：

根据稳态耦合非线性系统的伴随方法，自然对流传热问题的敏度分析被实施如下；增广泛函

由拉格朗日乘子和余量的点积定义：

增广泛函关于设计变量微分的表达式为：

通过定义满足如下关系式的拉格朗日因子λ_t,λ_p，使得在敏度表达式中

的系数为0：

因此敏度被简化为：

3.4)迭代优化以及结果显式表达：

将设计变量、目标函数、约束条件以及敏度作为输入参数，通过移动渐近线方法(MMA)对设计变量进行更新，在每一步优化中，会对设计变量及敏度进行过滤处理；

在单元上均匀分布N个网格点，通过类似于过滤的方式将设计变量分配到节点上，得到节点密度，其过滤方式为：

其中ρ_n表示节点密度；

表示过滤后的设计变量；权重因子w_en为：

w_en＝max(0,r_nmin-Δ(e,n)) (20)

其中r_nmin表示单元和节点的过滤半径；Δ(e,n)表示单元e和节点n之间的距离；

将节点密度进行插值，得到网格点密度；通过实施Heaviside光滑函数得到0-1二元设计；

其中β表示陡度参数，每次迭代的增量为Λ；迭代过程中

由二分法来确定；

定义网格点上的水平集函数：

当相邻迭代设计变量总的改变量小于指定精度ε，式(23)或者循环次数k大于循环最大次数k_max时，迭代终止；

最后，通过每个单元上的节点密度来更新设计变量：

所述的步骤4)具体为：迭代更新驱使设计结构收敛，优化完成后得到设计变量的0-1二元分布，设计变量为1所对应的区域集合即为散热翅片结构；依据制造性工艺要求，对散热翅片结构进行人工圆整，形成重力驱动下自然对流传热问题最终散热翅片优化设计结构。

本发明与现有技术相比，其显著优点在于：

(1)本发明将自然对流传热问题的模型进行简化，减少了求解的未知量的个数，减少了系统资源消耗，提高计算效率；

(2)本发明收敛准则采用的是设计变量总的改变量来度量，在收敛性能方面有更好的性能；

(3)本发明新的材料参数和导热率采用变密度方法，能够使得中间密度快速惩罚到0-1二元设计，能够保证优化结构的形状和拓扑；

(4)本发明在实现变密度方法的前提下，将设计变量分配到网格点密度上，通过构造水平集函数，实现显式边界的优化结果。

附图说明

图1为本发明的流程图。

图2为本发明实施例中优化问题几何模型示意图。

图3为本发明实施例中单元上网格点密度分布示意图。

图4为本发明实施例中自然对流传热模型拓扑优化后散热翅片结构示意图。

图5为本发明实施例中自然对流传热模型拓扑优化后的温度分布示意图。

具体实施方式

下面结合附图和实施例对本发明作进一步详细描述。

参照图1，一种重力驱动的自然对流异型热沉结构拓扑优化方法，包括以下步骤：

1)建立重力驱动自然对流传热问题物理模型以及模型简化：

1.1)建立重力驱动自然对流传热问题物理模型：

根据实际自然对流传热物理问题，提取问题的几何模型以及边界条件；整个计算区域被分为设计域和流体域，其中设计域包含流体和固体，整个设计域被流体域包围，流体流动会影响设计域中流体和固体的分布，设计域中固体对应的区域即为散热翅片结构；在整个计算区域中建立自然对流传热问题控制方程组，通过系数来控制区域的性态，其动量方程、质量方程和能量方程分别为：

其中u是流体的速度场，p表示压力场，T表示温度场；ρ₀是参考密度，取值为1；T₀是对应的参考温度，取值为0；β为体积膨胀系数，取值为100；μ是流体的动力学粘度，取值为1；g为重力加速度向量，g＝(0,-1)^T；c_p为比热容，取值为1；Q(x)为空间中的体积热源项，取值为0；k(x)为空间变化的有效导热率，被定义为：

其中k_f是流体的有效导热率，取值为1；k_s是固体的有效导热率，取值为100；

如图2所示，整个计算区域是一个封闭区域，包括灰色翅片设计区域和白色流体流动区域；区域尺寸分别为L_d＝7，L_b＝4，H_d＝4,H_b＝2.5，L_q＝0.2；在左、右、上边界上温度为0，下边界中央边界L_q处有向内的热流密度，大小为q_h＝110，下边界其他部分为热绝缘；所有的边界均为无滑移边界条件，考虑的均是定常的、不可压层流；给定右上角的参考压力为p＝0；

1.2)重力驱动自然对流传热问题物理模型简化：

假设对流项

被忽略，并且将μΔu表示为速度的线性表达式，重力驱动自然对流传热问题物理模型简化为：

其中

是新引入的材料参数，定义为：

其中

分别是流体和固体的材料参数；

在固体中设置新的材料参数

取∞，即

通过与商业软件COMSOLMultiphysics 5.5的自然对流传热全阶模型控制方程组进行对比，将温度的最小二乘误差作为因变量，通过调整流体中材料参数

的值，观察简化模型的逼近程度；当

时取得最小的温度的最小二乘误差，自然对流传热简化模型达到很好的拟合效果；

2)建立有限元分析模型：

对计算区域用四边形网格进行离散，采用有限元方法进行分析；整个计算区域被离散为140*80的四边形单元，总体的计算节点数为141*81＝11421；对压力和温度采用一阶多项式插值，对控制方程的两边乘以测试函数w，在计算区域进行积分得到自然对流传热问题物理模型的变分形式为：

对测试函数w进行扰动，即

其中稳定化参数τ为：

其中h表示四边形单元沿x、y方向尺寸的最大值，取值为0.05；

采用牛顿迭代的方法进行求解，得到自然对流传热问题的温度场、压力场以及速度场；

3)建立自然对流传热问题优化模型：

3.1)优化模型建立：

针对步骤2)所得到的自然对流传热问题的温度场，考虑的优化目标是使施加热流密度处的平均温度最小，满足的约束条件为自然对流传热方程组、设计域内固体材料所占的比例、设计变量的取值范围；建立自然对流传热问题优化模型：

其中N_d是指计算区域中施加了热流密度处所对应的所有节点总数N_d＝5；N_e是设计区域内单元总数，N_e＝4000；

为设计变量的集合；x_i是设计区域内第i个单元上的伪密度；v_i是设计区域内第i个单元的体积，v_i＝0.0025；f是指设计区域内所能容许的最大固体体积分数，f＝0.5；

是指设计区域Ω_d的体积，

等式约束为自然对流传热控制方程组的矩阵形式，其中k₁₁,k₁₂,k₂₂分别为系数矩阵，f₁,f₂分别为右端项；

3.2)材料参数连续化处理：

当设计变量x_i＝0时表示流体；当设计变量x_i＝1时表示固体，即散热翅片结构；x_i∈[0，1]表示流体和固体的过渡阶段；设计变量的变化，反映了散热翅片结构的变化；将新的材料参数和导热率表示为插值形式：

其中

p_k分别表示材料参数和导热率的惩罚因子，取值分别为p_μ＝[8,8,8,20]、p_k＝[2,8，16,16]，当相邻迭代设计变量的改变量|x_k+1-x_k|＜0.01或者循环次数是50的整数倍时，惩罚因子取下一个值；

3.3)敏度分析：

根据稳态耦合非线性系统的伴随方法，自然对流传热问题的敏度分析被实施如下；增广泛函

由拉格朗日乘子和余量的点积定义：

增广泛函关于设计变量微分的表达式为：

通过定义满足如下关系式的拉格朗日因子λ_t,λ_p，使得在敏度表达式中

的系数为0：

因此敏度被简化为：

3.4)迭代优化以及结果显式表达：

将设计变量、目标函数、约束条件以及敏度作为输入参数，通过移动渐近线方法(MMA)对设计变量进行更新，在每一步优化中，会对设计变量及敏度进行过滤处理，过滤半径被设置为1.5个单元尺寸；

实际工程中的设计结构需要得到非0即1的二元设计，密度惩罚法的拓扑优化问题会在边界上产生中间密度，不利于实际制造，在单元上均匀分布9个网格点，单元上网格点密度分布示意图如图3所示，网格处黑色圆点表示网格点密度为1，白色小圆点表示网格点密度为0；通过类似于过滤的方式将设计变量分配到节点上，得到节点密度，其过滤方式为：

其中ρ_n表示节点密度；

表示过滤后的设计变量；权重因子w_en为：

w_en＝max(0,r_nmin-Δ(e,n)) (20)

其中r_nmin表示单元和节点的过滤半径，取值为1；Δ(e,n)表示单元e和节点n之间的距离；

将节点密度进行插值，得到网格点密度；通过实施Heaviside光滑函数得到0-1二元设计；

其中β表示陡度参数，初始取值为0.5；每次迭代的增量为Λ，取值为0.5；迭代过程中

由二分法来确定；

定义网格点上的水平集函数：

当相邻迭代设计变量总的改变量小于指定精度ε，式(23)或者循环次数k大于循环最大次数k_max时，迭代终止；本实施例中相邻迭代设计变量总的改变量指定精度设置为ε＝10^(-4)，循环的最大次数设置为k_max＝1000次；本实施例在迭代第241步收敛，优化后的自然对流散热翅片的结构如图4所示；图5给出了自然对流传热模型拓扑优化后温度分布示意图；

最后，通过每个单元上的节点密度来更新设计变量：

其中N表示每个单元上网格点的数量，取值为9；

4)迭代更新驱使设计结构收敛，优化完成后得到设计变量的0-1二元分布，设计变量为1所对应的区域集合即为散热翅片结构；依据制造性工艺要求，对散热翅片结构进行人工圆整，形成重力驱动下自然对流传热问题最终散热翅片优化设计结构。

Claims (5)

1.一种重力驱动的自然对流异型热沉结构拓扑优化方法，其特征在于，包括以下步骤：

1)建立重力驱动自然对流传热问题物理模型以及模型简化；

2)建立有限元分析模型；

3)建立自然对流传热问题优化模型；

4)迭代更新驱使设计结构收敛。

2.根据权利要求1所述的一种重力驱动的自然对流异型热沉结构拓扑优化方法，其特征在于，所述的步骤1)具体为：

1.1)建立重力驱动自然对流传热问题物理模型：

根据实际自然对流传热物理问题，提取问题的几何模型以及边界条件；整个计算区域被分为设计域和流体域，其中设计域包含流体和固体，整个设计域被流体域包围，流体流动会影响设计域中流体和固体的分布，设计域中固体对应的区域即为散热翅片结构；在整个计算区域中建立自然对流传热问题控制方程组，通过系数来控制区域的性态，其动量方程、质量方程和能量方程分别为：

ρ₀(u·▽)u-μΔu+▽p＝ρ₀(-β(T-T₀))g (1)

▽·u＝0 (2)

ρ₀c_p▽·(uT)-▽·(k(x)▽T)＝Q(x) (3)

其中u是流体的速度场，p表示压力场，T表示温度场；ρ₀是参考密度；T₀是对应的参考温度；β为体积膨胀系数；μ是流体的动力学粘度；g为重力加速度向量；c_p为比热容；Q(x)为空间中的体积热源项；k(x)为空间变化的有效导热率，定义为：

其中k_f是流体的有效导热率；k_s是固体的有效导热率；

整个计算区域受重力作用的影响；边界上有指定的温度边界、热绝缘边界以及热流密度边界；所有的边界均为无滑移边界条件；考虑的均是定常的、不可压层流；

1.2)重力驱动自然对流传热问题物理模型简化：

假设对流项ρ(u·▽u)被忽略，并且将μΔu表示为速度的线性表达式，重力驱动自然对流传热问题物理模型简化为：

其中

是新引入的材料参数，定义为：

其中

分别是流体和固体的材料参数；

在固体中设置新的材料参数

取∞，即

通过与商业软件的自然对流传热全阶模型控制方程组进行对比，将温度的最小二乘误差作为因变量，通过调整流体中材料参数

的值，观察简化模型的逼近程度；最小的温度最小二乘误差所对应的参数即为流体的材料参数

3.根据权利要求2所述的一种重力驱动的自然对流异型热沉结构拓扑优化方法，其特征在于，所述的步骤2)具体为：

对计算区域用四边形网格进行离散，采用有限元方法进行分析；对压力和温度采用一阶多项式插值，对控制方程的两边乘以测试函数w，在计算区域进行积分得到自然对流传热问题物理模型的变分形式为：

对测试函数w进行扰动，即w^*＝w+τ(u·▽)w，其中稳定化参数τ为：

其中h表示四边形单元沿x、y方向尺寸的最大值；

采用牛顿迭代的方法进行求解，得到自然对流传热问题的温度场、压力场以及速度场。

4.根据权利要求3所述的一种重力驱动的自然对流异型热沉结构拓扑优化方法，其特征在于，所述的步骤3)具体为：

3.1)优化模型建立：

针对步骤2)所得到的自然对流传热问题的温度场，考虑的优化目标是使施加热流密度处的平均温度最小，满足的约束条件为自然对流传热方程组、设计域内固体材料所占的比例、设计变量的取值范围；建立自然对流传热问题优化模型：

其中N_d是指计算区域中施加了热流密度处所对应的所有节点总数；N_e是设计区域内单元总数；

为设计变量的集合；x_i是设计区域内第i个单元上的伪密度；v_i是设计区域内第i个单元的体积；f是指设计区域内所能容许的最大固体体积分数；

是指设计区域Ω_d的体积；等式约束为自然对流传热控制方程组的矩阵形式，其中k₁₁,k₁₂,k₂₂分别为系数矩阵，f₁,f₂分别为右端项；

3.2)材料参数连续化处理：

当设计变量x_i＝0时表示流体；当设计变量x_i＝1时表示固体，即散热翅片结构；x_i∈[0，1]表示流体和固体的过渡阶段；设计变量的变化，反映了散热翅片结构的变化；将新的材料参数和导热率表示为插值形式：

其中

p_k分别表示材料参数和导热率的惩罚因子，随着迭代步数的增加而增加；

3.3)敏度分析：

根据稳态耦合非线性系统的伴随方法，自然对流传热问题的敏度分析被实施如下；增广泛函

由拉格朗日乘子和余量的点积定义：

增广泛函关于设计变量微分的表达式为：

通过定义满足如下关系式的拉格朗日因子λ_t,λ_p，使得在敏度表达式中

的系数为0：

因此敏度被简化为：

3.4)迭代优化以及结果显式表达：

将设计变量、目标函数、约束条件以及敏度作为输入参数，通过移动渐近线方法(MMA)对设计变量进行更新，在每一步优化中，会对设计变量及敏度进行过滤处理；

在单元上均匀分布N个网格点，通过类似于过滤的方式将设计变量分配到节点上，得到节点密度，其过滤方式为：

其中ρ_n表示节点密度；

表示过滤后的设计变量；权重因子w_en为：

w_en＝max(0,r_{n min}-Δ(e,n)) (20)

其中r_{n min}表示单元和节点的过滤半径；Δ(e,n)表示单元e和节点n之间的距离；

将节点密度进行插值，得到网格点密度；通过实施Heaviside光滑函数得到0-1二元设计；

其中β表示陡度参数，每次迭代的增量为Λ；迭代过程中

由二分法来确定；

定义网格点上的水平集函数：

当相邻迭代设计变量总的改变量小于指定精度ε，式(23)或者循环次数k大于循环最大次数k_max时，迭代终止；

最后，通过每个单元上的节点密度来更新设计变量：

5.根据权利要求4所述的一种重力驱动的自然对流异型热沉结构拓扑优化方法，其特征在于，所述的步骤4)具体为：迭代更新驱使设计结构收敛，优化完成后得到设计变量的0-1二元分布，设计变量为1所对应的区域集合即为散热翅片结构；依据制造性工艺要求，对散热翅片结构进行人工圆整，形成重力驱动下自然对流传热问题最终散热翅片优化设计结构。

Images