[發明專利]一種由零虧格網格曲面生成散熱拓撲的圖形學方法

有效

• 申請號: 202010632026.3
• 申請日: 2020-07-02
• 公開（公告）號: CN111832203B
• 公開（公告）日: 2022-12-09
• 發明（設計）人: 李寶童;張路寬;劉宏磊;唐文豪;洪軍
• 申請（專利權）人: 西安交通大學
• 主分類號: G06F30/23
• 分類號: G06F30/23;G06F111/04;G06F119/08
• 代理公司: 西安智大知識產權代理事務所 61215
• 代理人: 賀建斌

1.一種由零虧格網格曲面生成散熱拓撲的圖形學方法，其特征在于，包括以下步驟：

1)定義設計工況：

以二維方形平面結構為設計域，設計域尺寸參數為l×l，內部均布熱載荷為Q，設計域邊界均為絕熱邊界，流體入口、出口分別位于設計域的上下邊界中央，寬度均為l/10；

2)構建初始零虧格網格曲面：

取n個互不重合的點作為零虧格網格中幾何模型的頂點，構建方形網格，各個頂點坐標值分別為(x_i,y_i)，i＝1,2,…,n，參數化幾何將方形網格進一步劃分為三角形網格，形成網格曲面描述的初始組件布局，設置網格進行細化和分裂的閾值，分別為L_S和D_S，設置流道體積分數上限β₀，最大迭代步數Loop，變量容差t_v；

3)初始化等幾何分析網格：

3.1)利用非均勻有理B樣條(NURBS)曲線構建等幾何分析網格：

取m個控制點構建p次B樣條曲線，其控制向量為0～1之間的非遞減序列ξ＝{ξ₁，ξ₂，···，ξ_m+p+1}，B樣條基函數的遞推公式：

使用兩組B樣條基函數N_i,p(ξ)、N_j,p(ξ)和權重系數ω_i,j獲得二維NURBS曲面，作為底層等幾何分析網格，其基函數表達式

3.2)構建等幾何單元的單元剛度矩陣：

等幾何分析需要實現物理空間、參數空間到母空間之間的映射，求解熱流耦合問題的單元剛度為：

式中，表示第e號等幾何單元u方向上第i個和v方向上第j個的高斯積分點，J(u，v)為二維雅克比矩陣，R_(i,l)和R_(k,j)分別為u方向上第段和v方向上第段NURBS曲線方程的導數；

3.3)施加載荷和約束：

根據給定的邊界條件，在計算層上施加諾依曼邊界條件，均布的熱載荷Q，在流體入口、出口的等幾何網格控制點上施加狄利克雷邊界條件，分別為恒定溫度T_in，恒定壓力P_out；，給定流體入口速度U_in；

4)計算流道體積：

根據上層零虧格三角形網格的半邊數據結構，遍歷所有三角形單元獲得網格曲面的幾何邊界，將邊界上的網格頂點設為控制頂點，即優化變量，三角形網格體積V的計算方法如下：

式中ng是三角形單元的個數，l_ng是三角形單元的邊長，t為三角形單元厚度；

p_l的計算方法如下：

5)物理場求解：

5.1)計算控制點投影值：

以規定的投影規則將上層三角形網格描述的拓撲結構映射至下層等幾何分析層，計算等幾何單元控制點的投影值，計算公式為：

式中，N^Acg，N^Bcg以及N^Ccg的具體計算方法如下：

式中，其中ng為三角形網格曲面中的單元數量，e為歐拉數，b為一個極大的常數，此處取b＝100000；

5.2)材料屬性插值模型：

等幾何單元控制點導熱系數插值模型：

式中，k_e為單元e的導熱系數，k_s為固體材料的導熱系數，k_w為流體材料自身的導熱系數，是單元e第j個節點的坐標；

固體、流體滲透率插值模型：

式中，κ_e為單元e的滲透率，κ_s為固體材料的滲透率，κ_w為流體材料自身的滲透率；

固體、流體密度插值模型：

式中，ρ_e為單元e的密度，ρ_s為固體材料的密度，ρ_w為流體材料自身的密度；

固體、流體比熱容插值模型：

式中，c_pe為單元e的比熱容，c_ps為固體材料的比熱容，c_pw為流體材料自身的比熱容；

5.3)構建整體剛度矩陣并求解物理場：

不可壓縮流體在熱流耦合問題中的控制方程為：

K_pP＝f_p (12)

(K_t+C(P))T＝f_t (13)

式中，K_p為整體滲透矩陣，K_t為整體導熱矩陣，P為壓力場向量，T為溫度場向量，f_p為壓力載荷向量，f_t為熱流載荷向量，C(P)為整體對流矩陣，是壓力場P的函數；

整體滲透矩陣K_p、導熱矩陣K_t以及對流矩陣C(P)為：

式中，N_e為區域Ω離散后包含的有限元單元數量，為單元e的導熱矩陣，為單元e的滲透矩陣，c^e為單元e的對流矩陣；

6)確定優化模型：

以熱流結構最佳散熱性能為優化目標，將結構的平均溫度設為目標函數；設定散熱結構能量耗散不得超過許用能量耗散、材料用量不得超過許用材料用量，將流體入口與流體出口間的壓降、流道體積占比作為約束函數；

優化數學模型如下：

式中，d＝(d₁,d₂,d₃,…,d_ng)^T為設計變量的所有控制頂點的坐標向量，ng為優化問題中涉及的控制頂點的總數，I為優化問題的目標函數，ΔP^*為限定的壓降約束值，V₀為設計域總體積，β₀為設計域內流體冷卻通道的體積V的體積分數上限值，U_d是d中所有取值的集合；

7)敏度分析：

以任意設計變量a為代表，對散熱結構的目標函數、約束函數進行靈敏度分析；

7.1)目標函數敏度：

目標函數的敏度計算公式如下：

式中：

7.2)約束函數敏度：

壓降約束的敏度計算公式如下：

式中：

體積約束的敏度計算公式如下：

式中，ngp為第e個等幾何單元中的高斯點個數，ncp為第e個等幾何單元中的節點個數；

8)迭代優化：

8.1)移動漸近線優化算法更新設計變量：

將由等幾何分析網格計算得到的目標函數、約束函數值以及它們的靈敏度值代入移動漸近線優化算法(MMA)中，迭代更新控制頂點的坐標值；

8.2)更新上層網格曲面形狀：

8.2.1)網格操作判定：

對控制頂點更新后的坐標值進行網格操作判定；三角形網格邊長滿足L_maxL_s條件時進行網格細分；相鄰三角形網格未共享頂點之間距離滿足D_ijD_s條件時進行網格分裂；根據Greiner-Hormann算法對三角形網格進行自相交判定，采用Constrained DelaunayTriangulation算法對三角形網格進行重劃分；獲得上層網格邊界形狀；

8.2.2)更新非控制點坐標：

在獲得控制頂點更新坐標之后，利用保剛性算法計算出其他非控制頂點的更新坐標，非控制頂點的移動根據最小化能量函數獲得，計算公式如下：

8.3)循壞迭代獲得最優解：

循環迭代步驟4)至步驟8.2.2)過程，直至迭代步數k大于最大迭代步數Loop或者更新前后設計變量差值小于變量容差t_v為止，此時獲得滿足材料用量和能量耗散約束的最優流道散熱結構；

9)適應性處理：按照生產工藝要求圓整熱流結構布局，從而獲得最終流道分布。