一种适配风扇特性曲线的电子设备异型风道设计方法与流程

1.本发明属于电子设备散热设计技术领域，具体涉及一种适配风扇特性曲线的电子设备异型风道设计方法。


背景技术：

2.随着半导体技术的发展，各类芯片的晶体管密度日益增加，性能愈发强大，但芯片的散热设计功耗也越来越高。为了维持芯片的正常工作温度，散热器的体积、重量和成本也越来越高。当前如手机、平板、笔记本电脑等电子设备的发展趋势是体积越来越小，重量越来越轻，功能越来越强，其主要矛盾在于有限的空间内实现机体内高效的散热设计。
3.常见的电子设备风冷应用场景，如笔记本电脑，其散热模组以风扇为动力源，以空气为主要冷却工质，使用热管或均热板将主要的热源芯片的热量转移到鳍片，再通过风扇实现与鳍片的对流热交换，将热量排出机体外，以实现降温目的。影响风冷效果的不仅有散热模组的性能，还与机体腔内风道对冷空气的组织相关。合理的风道设计能提高机体腔内冷空气利用效率，避免热空气的回流，降低主板和键盘温度，提升用户体验。此外，由于风扇本身的工作曲线决定其所能提供的风量与其风压负相关，因此风道设计还能降低沿程阻力，提升风扇的工作效率。
4.目前对于风道的设计，包括对进出风口的位置、形状和尺寸设计、内部风冷风道结构的尺寸和形状设计，较依赖于工程师的主观经验；此外，在风道结构设计过程中，由于风扇工作特性曲线的存在，其工作点随内部风阻的变化而实时改变，导致较强的非线性，增加了散热问题的设计难度。
5.风道结构的另一种可行的设计方法是拓扑优化方法。由于空气工质的动力粘度较小，且风冷应用场景中往往风速较高，这导致空气的流态以湍流流动为主；此外，对于大部分应用场景，风速一般小于0.3马赫，可以按照不可压缩流体处理，因此常见的风冷导流问题本质是不可压缩湍流流固共轭传热问题。对于湍流流动的模拟，目前工业上最为流行的方法是基于求解雷诺平均的纳维-斯托克斯方程(rans)，但由于rans方程本身求解变量相比层流流动更多，且偏微分方程存在很强的非线性，因此在数值求解过程需要占用大量内存以及对方程迭代求解，计算成本较高；且拓扑优化过程通常需要成百上千次的迭代寻优，因此直接以rans方程进行拓扑优化计算成本高昂，时间成本难以接受；此外，为适配风扇特性曲线，相比常用的速度边界条件和压力边界条件，进一步增加了散热问题的复杂性。
6.由于风扇特性曲线带来的设计问题的复杂性，电子设备散热设往往先完成散热模组的设计，再根据尺寸和性能的约束，进行风扇的定制和选型。然而厂家提供的风扇型号往往不是连续的，为保证散热性能，往往要求风扇的风量和风压尽量大，这造成了风扇性能不能得到充分的发挥，产生了额外的成本。


技术实现要素：

7.为了克服现有上述现有技术的缺点，本发明的目的在于提供一种适配风扇特性曲
线的电子设备异型风道设计方法，在系统压降损失允许的范围内，最小化热源均温，同时可实现风扇工作点的精确控制，适用于风扇作为动力源的电子设备散热场景。
8.为了达到上述目的，本发明采取的技术方案为：
9.一种适配风扇特性曲线的电子设备异型风道设计方法，使用达西简化模型近似湍流流场，使用变密度法以渗透率描述流体和固体；建立风扇特性曲线同系统流场及温度场耦合的简化模型，实现风扇工作点与当前工况实时匹配；使用移动渐近线优化算法，在系统压降约束和体积约束条件下，最小化热源均温，实现电子设备内风道的设计。
10.一种适配风扇特性曲线的电子设备异型风道设计方法，包括以下步骤：
11.1)建立达西流模型描述下的系统控制方程，包括质量守恒方程、动量守恒方程和能量守恒方程，并建立对应的有限元刚度方程；
12.2)配置变密度法中流体材料和固体材料对应的伪密度值x，将其映射到0-1的区间，其中x＝0表示固体，x＝1表示流体；确定流体材料体积的上限确定初始的流体材料和固体材料的分布；
13.3)确定系统所使用风扇特性曲线，用风压表征强制对流过程中系统阻力；
14.4)配置具有简化风扇模型的边界条件，并在对应的有限元单元和节点施加载荷和约束；
15.5)迭代准备，初始化移动渐近线mma优化算法的参数配置，启动优化迭代过程；
16.6)建立网格过滤器，有限元网格单元的伪密度为其自身伪密度与其邻近单元的伪密度的加权平均，更新单元伪密度；根据当前伪密度计算单元的材料属性；
17.7)计算各单元的压力单元刚度矩阵，并组装得到压力场求解的整体刚度矩阵；两步迭代计算当前系统阻抗匹配的风扇工作点；
18.8)根据当前系统各单元的热导率以及迭代后的压力场，计算各有限元单元的单元刚度方程，并组装得到温度场计算的整体刚度矩阵；根据体热源及热流密度边界条件，计算温度场计算的载荷向量；风扇进出风口边界温度耦合，修正温度场计算的整体刚度矩阵和载荷向量，并计算得各有限元节点温度；
19.9)计算目标函数、压降约束函数及体积约束函数，并进行灵敏度分析；将上述灵敏度输入到移动渐近线优化算法中，优化获得新的设计变量，即各单元的伪密度值；
20.10)收敛条件判断，满足以下任一条件，停止迭代；否则返回步骤6)；其一，当前迭代步数大于最大迭代步数，强制停止优化迭代；其二，当相邻两迭代步设计变量的最大变化值小于收敛标准。
21.所述的步骤1)中达西流模型描述下的系统控制方程为：
[0022][0023]
式中，u为速度场，p为压力场，t为温度场，q为体热源，α为渗透率，μ为动力粘度，ρ为密度，c
p
为比热容，k为热导率；
[0024]
达西流模型描述下的有限元刚度方程为：
[0025][0026]
式中，k
p
、f
p
、p分别为压力场求解的总体刚度矩阵、载荷向量、解向量，k
t
、c(p)、f
t
、t分别为温度场求解的热传导总体刚度矩阵、热对流总体刚度矩阵、载荷向量、解向量，n为有限元形函数。
[0027]
所述的步骤2)中的变密度法描述拓扑时，拓扑优化问题公式表达为：
[0028][0029]
式中，xi为各单元的伪密度，x为所有设计变量构成的向量，c为目标函数，nh为目标函数所关注的节点数量，ti为对应节点温度，vc和为当前流体体积分数及其上限，p
fin
和p
fout
分别为风扇进风口的压力和出风口的压力。
[0030]
所述的步骤3)中风扇特性曲线为风量-风压特性曲线，描述的风量与风压应为负相关。
[0031]
所述的步骤4)中简化风扇模型是指用不考虑风扇内部因扇叶转动产生的流场，而是将风扇进风口和出风口作为边界条件；边界条件包括压力边界条件、速度边界条件、温度边界条件、体热源、边界热流密度；
[0032]
确定流体从外界流入系统的进口和出口在达西流模型下均设置为狄利克雷边界条件，即设置系统进口节点处的压强和出口节点处的压强；设置系统进口在温度控制方程中为狄利克雷边界条件，即设置系统进口节点处的温度；
[0033]
风扇模型区域内部不分配网格单元和节点，风扇模型内部的流动和传热现象作为黑箱处理；风扇的进风口和出风口在达西流模型下均设置为诺伊曼边界条件，即设置进口单元面的法向速度和出口单元面的法向速度，满足“进风口面积
×
进风速度＝出风口面积
×
出风口速度”；
[0034]
风扇出风口温度和进风口温度耦合，设置风扇出风口节点温度为进风口节点温度均值与风扇自身发热功率相关的温度增量的加和；其表达式为：
[0035][0036]
式中，t
fout
和t
fin
分别为风扇出风口节点温度和进风口节点温度，n
fin
为风扇进风口节点数量，δt为风扇自身热损耗产生的温差。
[0037]
所述的步骤6)中网格过滤器的配置包括网格的过滤半径和单元的权重，其形式为：
[0038][0039]
式中，r
min
为网格过滤半径，和x分别表示过滤后和过滤前的单元伪密度，e单元应位于过滤半径内的中心位置，dist为两单元之间的距离，h
e,i
为i单元对e单元的过滤权重；
[0040]
网格过滤后单元的材料属性表达式为：
[0041][0042]
式中，α、k、ρ、c
p
分别表示单元的渗透率、热导率、密度、比热，f和s分别用于表示流体和固体，p
α
、pk、p
ρ
、分别表示单元的渗透率、热导率、密度、比热的插值惩罚因子。
[0043]
所述的步骤7)中两步迭代计算当前系统阻抗匹配的风扇工作点，第一步需配置初始假设风扇出风口风速，并根据进出口风量相等，计算对应风扇进风口风速，并计算当前压力场，获取当前风扇进出风口的压差；第二步需根据达西流流场特性，系统阻抗与风扇风速存在线性关系，即为过原点斜率为正的直线，计算得到该直线与风扇的风量-风压曲线的交点，即为当前风扇工作点；取该工作点的风量，对应计算风扇进出风口的速度，作为边界条件，再次计算压力场，并计算风扇进出风口的压强差，即为该风扇匹配当前系统阻抗的风压；
[0044]
两步迭代计算当前系统阻抗匹配的风扇工作点，具体表达式如下：
[0045][0046]
式中，u1、f
p1
、p1分别表示初始假设的风扇出风口速度及对应计算所得载荷向量和压力场解向量，a表示风扇出风口的面积，p1和q1表示初始假设计算所得风压和风量，f(p,q)表示风扇风量风压特性曲线，p2和q2表示匹配当前系统阻抗的风扇工作点对应的风压和风量；u2、f
p2
、p2分别表示工作点对应的风扇出风口速度及载荷向量和压力场解向量。
[0047]
所述的步骤9)中对目标函数、压降约束函数及体积约束函数进行灵敏度计算，使用伴随法分析；某函数对设计变量的伴随敏度表达形式为：
[0048][0049][0050][0051][0052][0053][0054][0055][0056]
式中，φ和分别为设计变量相关的某函数(目标函数或约束函数)及其对应的伴
随形式，λ
tt
和为伴随向量，为迎风修正后的形函数，b为形函数n对空间的偏导；
[0057]
目标函数即热源均温对各设计变量的偏导，结合式(7)表达为：
[0058][0059]
式中，ωh表示目标函数所关注的节点的集合；
[0060]
压降约束函数即风扇风压对各设计变量的偏导，结合式(7)表达为：
[0061][0062]
式中，ω
fout
和ω
fin
分别表示风扇出风口节点集合和进风口节点集合，n
fout
和n
fin
分别表示风扇出风口节点数量和进风口节点数量；在优化过程，通过初值和压降约束函数的设置，获得在某预设风扇工作点下的最优风道结构；如若在优化过程中，随着目标函数的下降，压降约束函数呈下降趋势，则通过约束风扇进出风口压降大于某设定值，控制风道在该预设风扇风压值下进行优化；若在优化过程中，压降约束函数随目标函数的下降呈上升趋势，则约束风扇进出风口压降小于某设定值，实现风道结构优化过程中风扇工作点的精确控制；
[0063]
体积约束函数即流体体积占比对各设计变量的偏导，结合式(7)表达为：
[0064][0065]
本发明的有益效果为：
[0066]
1)本发明使用简化湍流计算模型和简化风扇模型，实现湍流状态下流固共轭传热问题与风扇工作特性曲线的高效耦合，仅需两步迭代即可获取适配当前风道结构的风扇工作点，具有更高的计算效率和更低的计算成本；
[0067]
2)本发明采用拓扑优化的设计方法，将风扇模型作为边界条件耦合其中，风扇进出风口的风量和风压随拓扑结构的迭代优化实时改变，可实现风扇先选型，后优化风道，并获取适配当前风扇特性曲线的最优风道以及最优风扇工作点；相比恒流量作为边界条件的
优化，更符合工程实际，可获得满足设计需求的更优的风道构型；
[0068]
3)本发明根据实际工况约束，可优化风道到任意指定风扇工作点，如风扇能量转换效率最高工作点、风扇寿命最长工作点、指定系统风压工作点，满足定制化设计需求。
附图说明
[0069]
图1是本发明方法的流程图。
[0070]
图2是本发明两步迭代计算当前系统阻抗匹配的风扇工作点的原理图。
[0071]
图3是本发明适配风扇特性曲线的风道设计域(图(a))和设计结果(图(b))。
[0072]
图4是本发明优化前与优化后对应的风扇工作点。
具体实施方式
[0073]
下面结合实施例和附图对本发明做进一步详细描述。
[0074]
如图1所示，一种适配风扇特性曲线的电子设备异型风道设计方法，包括以下步骤：
[0075]
1)建立达西流模型描述下的系统控制方程，包括质量守恒方程，动量守恒方程和能量守恒方程，并建立对应的有限元刚度方程：
[0076]
达西流模型描述下的系统控制方程为：
[0077][0078]
式中，u为速度场，p为压力场，t为温度场，q为体热源，α为渗透率，μ为动力粘度，ρ为密度，c
p
为比热容，k为热导率；达西流模型可通过渗透率α，对材料物性以渗透率进行插值，进而实现对固体单元、流体单元及变密度法中的灰度单元的统一描述。
[0079]
达西流模型描述下的有限元刚度方程为：
[0080][0081]
式中，k
p
、f
p
、p分别为压力场求解的总体刚度矩阵、载荷向量和解向量，k
t
、c(p)、f
t
、t分别为温度场求解的热传导总体刚度矩阵、热对流总体刚度矩阵、载荷向量和解向量，n为有限元形函数；对于图3(a)所述的二维问题，整体尺寸360
×
260mm，模拟系统内单风扇用于五个热源的降温，如何设计风道可最大化地发挥风扇性能，使热源温度最低；本实施例采用四节点四边形单元划分结构化网格，单元尺寸2mm，单元内形函数采用四高斯积分点形式；
[0082]
2)配置变密度法中流体材料和固体材料对应的伪密度值x，将其映射到0-1的区间，其中x＝0表示固体，x＝1表示流体；确定初始的流体材料和固体材料的分布，初始设置所有单元的伪密度均为0.8；利用变密度法描述拓扑时，拓扑优化问题公式表达为：
[0083][0084]
式中，xi为各单元的伪密度，x为所有设计变量构成的向量，c为目标函数，nh为目标函数所关注的节点数量，ti为对应节点温度，vc和为当前流体体积分数及其上限，p
fin
和p
fout
分别为风扇进风口的压力和出风口的压力；这里以图3(a)所示的设计域内的网格单元的伪密度值为设计变量x；以最小化图3(a)所示的五个体热源的平均温度为目标函数，确定流体材料体积的上限可全部为流体；g(p
fout-p
fin
)＝p
fout-p
fin-15≤0，即约束风扇风压小于15pa；
[0085]
3)确定系统所使用风扇的风量-风压特性曲线，风量-风压特性曲线为第一象限斜率为负的直线，如图2所示，最大风压为100pa，最大风速(由风量转换而来)为8m/s；用风压表征强制对流过程中系统阻抗；该风扇参数是在笔记本电脑、台式电脑等电子设备的风冷散热场景中可接受的风扇参数；
[0086]
4)配置具有简化风扇模型的边界条件，并在对应的有限元单元和节点施加载荷和约束；如图3(a)所示，本实施例在360
×
260mm的二维区域，有1个风扇、5个体热源、系统的进出口、风扇的进出风口，其位置和尺寸如附图所示；
[0087]
简化风扇模型是指用不考虑风扇内部因扇叶转动产生的流场，而是将风扇进风口和出风口作为边界条件；边界条件包括压力边界条件、速度边界条件、温度边界条件、体热源；
[0088]
确定流体从外界流入系统的进口和出口在达西流模型下均设置为狄利克雷边界条件，即设置系统进口节点处的压强和出口节点处的压强，这里设置为相对大气压为0pa；设置系统进口在温度控制方程中为狄利克雷边界条件，即设置系统进口节点处的温度，这里可设置温度为环境温度300k；体热源设置为105w/m3；
[0089]
风扇模型区域内部不分配网格单元和节点，风扇模型内部的流动和传热现象作为黑箱处理；风扇的进风口和出风口在达西流模型下均设置为诺伊曼边界条件，即设置进口单元面的法向速度和出口单元面的法向速度；为保证风扇进风口风量和出风口风量，由于本实施例风扇进出风口面积相同，因此风扇进出风口风速相同；
[0090]
风扇出风口温度和进风口温度耦合，设置风扇出风口节点温度为进风口节点温度均值与风扇自身发热功率相关的温度增量的加和；其表达式为：
[0091]
[0092]
式中，t
fout
和t
fin
分别为风扇出风口节点温度和进风口节点温度，n
fin
为风扇进风口节点数量，δt为风扇自身热损耗产生的温差；由于风扇通常自身功率相比电子设备整体功率为小值，通常可忽略不计，因此本实施例不考虑风扇自身热损耗，取δt＝0；
[0093]
5)迭代准备，初始化移动渐近线(mma)优化算法的参数配置，启动优化迭代过程；
[0094]
6)建立网格过滤器，以避免变密度法(simp)中的“棋盘格”现象；有限元网格单元的伪密度为其自身伪密度与其邻近单元的伪密度的加权平均，更新单元伪密度；根据当前伪密度计算单元的材料属性；
[0095]
网格过滤器的配置包括网格的过滤半径和单元的权重，其形式为：
[0096][0097]
式中，r
min
为网格过滤半径；和x分别表示过滤后和过滤前的单元伪密度，e单元应位于过滤半径内的中心位置，dist为两单元之间的距离，h
e,i
为i单元对e单元的过滤权重；本实施例根据实际网格数量，调整参数后选择取r
min
＝1.5
[0098]
网格过滤后单元的材料属性表达式为：
[0099][0100]
式中，α、k、ρ、c
p
分别表示单元的渗透率、热导率、密度、比热，f和s分别用于表示流体和固体，p
α
、pk、p
ρ
、分别表示单元的渗透率、热导率、密度、比热的插值惩罚因子；按照国际标准单位，这里流体的材料参数按照空气取值αf＝1.05e-6，kf＝0.024，ρf＝1.225，导流板固体材料仅考虑其导流作用，因此热导率尽量小，取值αs＝1.05e-16，ks＝0.00242，ρs＝2719，c
ps
＝871；惩罚因子取值p
α
＝3，pk＝3，p
ρ
＝1，此外，体热源区域的网格不参与设计，材料属性参考硅，热导率为148；
[0101]
7)计算各单元的压力单元刚度矩阵，并组装得到压力场求解的整体刚度矩阵；两步迭代计算当前系统阻抗匹配的风扇工作点，过程如图2所示；
[0102]
第一步需配置初始假设风扇出风口风速，该假定风速可任意取值，这里取5m/s；根据进出口风量相等，计算对应风扇进风口风速；计算当前压力场，获取当前风扇进出风口的压差；
[0103]
第二步需根据达西流流场特性，系统阻抗与风扇风速存在线性关系，即为过原点斜率为正的直线，可计算得到该直线与风扇的风量-风压曲线的交点，即为当前风扇工作点；取该工作点的风量，对应计算风扇进出风口的速度，作为边界条件，再次计算压力场，并计算风扇进出风口的压强差，即为该风扇匹配当前系统阻抗的风压；
[0104]
两步迭代计算当前系统阻抗匹配的风扇工作点，具体表达式如下：
[0105][0106]
式中，u1、f
p1
、p1分别表示初始假设的风扇出风口速度及对应计算所得载荷向量和压力场解向量，a表示风扇出风口的面积，p1和q1表示初始假设计算所得风压和风量，f(p,q)表示风扇风量风压特性曲线，p2和q2表示匹配当前系统阻抗的风扇工作点对应的风压和风量；u2、f
p2
、p2分别表示工作点对应的风扇出风口速度及载荷向量和压力场解向量；
[0107]
8)根据当前系统各单元的热导率以及迭代后的压力场，计算各有限元单元的单元刚度方程，并组装得到温度场计算的整体刚度矩阵；根据体热源及热流密度边界条件，计算温度场计算的载荷向量；风扇进出风口边界温度耦合，修正温度场计算的整体刚度矩阵和载荷向量，并计算得各有限元节点温度；
[0108]
9)计算目标函数、压降约束函数及体积约束函数，并进行灵敏度分析；将上述灵敏度输入到移动渐近线优化算法中，优化获得新的设计变量，即各单元的伪密度值；
[0109]
对目标函数、压降约束函数及体积约束函数进行灵敏度计算，使用伴随法分析；某函数对设计变量的伴随敏度表达形式为：
[0110][0111][0112][0113][0114][0115]
[0116][0117][0118]
式中，φ和分别为设计变量相关的某函数(目标函数或约束函数)及其对应的伴随形式，λ
tt
和为伴随向量，为迎风修正后的形函数，b为形函数n对空间的偏导；
[0119]
目标函数即热源均温对各设计变量的偏导，可表达为：
[0120][0121]
式中，ωh表示目标函数所关注的节点的集合；本实施例中，ωh即为图3(a)所示的5个体热源包含的节点，以体热源节点的平均温度作为目标函数；
[0122]
压降约束函数即风扇风压对各设计变量的偏导，可表达为：
[0123][0124]
式中，ω
fout
和ω
fin
分别表示风扇出风口节点集合和进风口节点集合，n
fout
和n
fin
分别表示风扇出风口节点数量和进风口节点数量；
[0125]
体积约束函数即流体体积占比对各设计变量的偏导，可表达为：
[0126][0127]
10)收敛条件判断，满足以下任一条件，停止迭代；否则返回步骤6)；其一，当前迭代步数大于最大迭代步数，强制停止优化迭代，这里取最大迭代步数1000；其二，当相邻两迭代步设计变量的最大变化值小于收敛标准，这里取max(|x
k 1-xk|)＜0.0001。
[0128]
本实施例优化得到的异形风道结构如图3(b)所示；图4为优化前(设计域全为流体且未设计风道结构)和优化后风扇工作点的变化。在k-ω湍流模型验证下，对比优化前后温
度场，结果如表一所示，结温从417k下降到330k，平均温度(包括所有固体和流体的体积平均)从316k下降到304k；无风道时风扇的风量并不能充分流经热源区域，而异型风道可以使风扇性能得到充分发挥，合理分配风量到各个热源，有效降低热源的结温和整体的均温。
[0129]
表一
[0130]