張晨群，鮑益東，安魯陵，宦 蕾

（南京航空航天大學(xué)機(jī)電學(xué)院，南京 210016）

熱壓罐工藝是當(dāng)前制造大型航空用復(fù)合材料的主要工藝[1]。在熱壓罐成型工藝中，模具工裝型面的溫度場分布是影響制件質(zhì)量的關(guān)鍵因素之一。成型過程中模具表面溫度分布不均會導(dǎo)致制件內(nèi)部存在溫度梯度，以至于制件固化不同步、在結(jié)構(gòu)內(nèi)部產(chǎn)生殘余應(yīng)力和殘余應(yīng)變，最終會造成制件出現(xiàn)內(nèi)部缺陷、發(fā)生初始破壞，嚴(yán)重時會影響制件的質(zhì)量和使用壽命[2–4]。因此，分析和研究熱壓罐成型過程時模具與復(fù)合材料構(gòu)件接觸的表面溫度分布特點對改善復(fù)合材料構(gòu)件最終的成型質(zhì)量具有重要意義。

當(dāng)前大多數(shù)學(xué)者在模擬熱壓罐成型過程使用傳統(tǒng)的以網(wǎng)格為基礎(chǔ)的計算機(jī)流體力學(xué)（Computational Fluid Dynamics，CFD）軟件，如CFX、Fluent 等[3,5–7]，在前處理時都需花費大量時間進(jìn)行網(wǎng)格劃分等工作。而且對于熱壓罐成型這類流固耦合問題，傳統(tǒng)的網(wǎng)格劃分技術(shù)難以滿足實際工程上的精度要求。因此需要尋找一種能夠減少前處理時間、提升計算精度的方法來解決目前的問題。

本文將使用新一代流體力學(xué)軟件XFlow，建立框架式模具在熱壓罐中強(qiáng)迫對流換熱的有限元模型，并分析模具表面溫度場溫度分布規(guī)律及其影響因素。

XFlow 軟件及其原理

1 XFlow概述

傳統(tǒng)的CFD 軟件大都采用有限體積法[8]，如STAR–CCM、Phoenics、ANSYS Fluent，此外，CFX 混合使用有限體積法和有限元法。研究人員使用時需消耗大量的時間劃分網(wǎng)格，且計算結(jié)果的可靠性高度依賴網(wǎng)格質(zhì)量。但計算精度不是隨著網(wǎng)格尺寸的減小而不斷增加，且網(wǎng)格尺寸減小會增加計算資源的消耗和計算時間。XFlow 是采用格子玻爾茲曼技術(shù)的新一代CFD 軟件。它的自動點陣生成和自適應(yīng)優(yōu)化功能避免了人工劃分網(wǎng)格，有效縮短了前處理的時間。

2 格子玻爾茲曼方法

XFlow 所使用的LBM 方法是一種介于宏觀流體連續(xù)性假設(shè)與微觀分子動力學(xué)之間的介觀模擬方法，同時具有微觀方法的適應(yīng)性廣和宏觀方法的不關(guān)注分子運動細(xì)節(jié)的特點，精度和計算量上均有較大優(yōu)勢[9]。在LBM 模型中，假設(shè)宏觀流體由大量虛擬流體粒子構(gòu)成，粒子僅可做兩種基本運動，如圖1 所示的對流和碰撞。微觀顆粒不斷進(jìn)行“對流–碰撞”，構(gòu)成了流體的宏觀運動。LBM 方法的基本思想是不研究每個微觀粒子的具體運動，而是利用概率密度函數(shù)研究大量粒子在格子方向上的運動概率[10–11]。

3 工作流程

XFlow 的前處理、求解器和后處理完全集成在同一用戶界面環(huán)境中，其建立數(shù)值模型與分析流程如圖2所示。

基于XFlow的熱壓罐成型過程的模擬方法

1 控制方程

由于熱壓罐內(nèi)的熱對流和熱傳導(dǎo)的換熱方式，所以需要使用流體流動與熱交換中的質(zhì)量、動能及能量守恒方程作為模擬的基本控制方程，進(jìn)行溫度場的模擬。所需要求解的方程在直角坐標(biāo)系下的控制方程如下所示：

質(zhì)量守恒方程：

式中，ρ為流體密度；t為時間；U為流體速度。

動量守恒（牛頓運動定律）方程：

能量守恒（熱力學(xué)第一定律）方程：

式中，h=h(p，T)，與流體壓強(qiáng)p和流體溫度T有關(guān)；λ為流體的導(dǎo)熱系數(shù)；Sh為流體的內(nèi)熱源；φ為耗散函數(shù)。

狀態(tài)方程，使得方程組封閉：

對于固體區(qū)域，其能量方程為：

式中，ρs為固體密度；Cs為固體比熱；Ts為固體溫度；QT為內(nèi)部熱源項。

2 模型的建立與導(dǎo)入

熱壓罐中為了使溫度場均勻利用風(fēng)扇加快壓罐內(nèi)空氣的流動，同時通過控制系統(tǒng)來控制熱壓罐內(nèi)的溫度與壓強(qiáng)。按實際情況模擬難以實現(xiàn)，需對模型進(jìn)行簡化。對熱壓罐裝置的簡化過程如下：僅模擬工作內(nèi)腔，將其簡化成圓柱體模型，一端為進(jìn)口，另一端為出口。模具上的一些輔助裝置對溫度場影響很小，對其進(jìn)行簡化：僅保留型板和支撐結(jié)構(gòu)。利用三維建模軟件CATIA建立尺寸為φ2500mm×7000mm的簡化熱壓罐模型及尺寸為1500mm×1500mm×400mm 的簡化框架式模具模型。圖3 為復(fù)合材料成型模具。

圖1 粒子在空間格子上的兩種基本動態(tài)運動Fig.1 Two basic dynamic motions of particles on a lattice

圖2 XFlow求解流程Fig.2 XFlow solution flow

把建好的模具模型與熱壓罐模型分別導(dǎo)入到XFlow 并調(diào)整相對位置，如圖4 所示。

3 定義求解類型

流體模型選擇Single phase，因其只涉及整個流體域中的一個連續(xù)流體相，流體領(lǐng)域的任何一點都是相同的流體物質(zhì)，符合熱壓罐成型過程中罐內(nèi)流體的情況。

熱模型選擇Couple energy，該解算器可用于解釋當(dāng)氣體高度壓縮/膨脹時的壓力/溫度變化，適用于熱壓罐成型中氣體溫度變化且伴隨加壓的情況。由于熱壓罐內(nèi)溫度在200℃以下，可忽略熱輻射，因此把熱輻射方式設(shè)置為禁用。

接下來就是定義湍流模型。熱壓罐內(nèi)復(fù)合材料模具的升溫主要依靠熱對流傳熱方式實現(xiàn)，因此還需研究分析模擬過程中與模具進(jìn)行熱交換的流體區(qū)域的流動狀態(tài)。流體的流動狀態(tài)一般分為層流與湍流，使用雷諾系數(shù)Re 的大小來判斷。

當(dāng)空氣在熱壓罐內(nèi)流動時，相當(dāng)于氣體在管內(nèi)流動，則其雷諾數(shù)為

式中，ρ為氣體密度，u為流體流動速率，μ為氣體的動力黏度，D為熱壓罐直徑。

熱壓罐內(nèi)流體為空氣，密度為1.225kg·m–3，動力黏性系數(shù)取值1.82×10–5Pa?s。風(fēng)速1.5m/s，熱壓罐直徑為2.5m。

通常雷諾數(shù)與流動狀態(tài)的對應(yīng)關(guān)系如下：

（1）Re ≤2300，為層流；

（2）2300 ≤Re ≤8000，為層流與湍流的過渡區(qū)；

（3）Re ≥8000，為湍流。

根據(jù)計算出的雷諾數(shù)可知罐內(nèi)氣體流動類型為湍流，選擇自動適應(yīng)壁面的局部渦流模型。

4 定義初始條件與材料屬性

重力對溫度場影響極小，為減少計算量把重力加速度設(shè)為0。初始溫度設(shè)為288K（室溫），初始壓強(qiáng)設(shè)為標(biāo)準(zhǔn)大氣壓。接下來定義材料屬性：罐內(nèi)流體為空氣，按空氣的屬性定義流體參數(shù)。熱壓罐成型過程包含復(fù)雜的流固耦合問題，定義模具為共軛傳熱屬性，讓模具熱傳導(dǎo)與流體對流同時進(jìn)行；模具的型面和支撐結(jié)構(gòu)均采用Q235?？諝馀c模具的熱屬性如表1 所示。

5 定義邊界條件與仿真參數(shù)

（1）入口設(shè)置。

圓柱模型一個端面作為進(jìn)口，邊界條件設(shè)置為Inlet–Velocity，設(shè)定入口速度為風(fēng)速，熱邊界條件為氣體溫度–時間變化函數(shù)。

（2）出口及壁面設(shè)置。

設(shè)置出口類型為Convective outlet，即出口處的靜壓和速度都是從內(nèi)部域推斷出來的。同時設(shè)置熱壓罐圓筒壁面為Wall，熱邊界條件類型為絕熱類型。

（3）仿真參數(shù)設(shè)置。

設(shè)置仿真時間，XFlow 無需劃分網(wǎng)格，只需定義求解尺寸即粒子大小。

模型精度驗證及模具型面溫度場分析

1 熱壓罐成型模具傳熱分析模型精度驗證

為了驗證XFlow 應(yīng)用在熱壓罐溫度場模擬上的可行性，需要分析仿真結(jié)果的精度。將仿真得到的模具表面的溫度數(shù)據(jù)和實驗室中利用熱電偶檢測得到的數(shù)據(jù)[5]進(jìn)行對比，溫度監(jiān)測點如圖5 所示。以0.025K/s的升溫速率升溫至第4800s，此時模具表面溫度分布如圖6 所示。

模具是對稱結(jié)構(gòu)，選擇沿對角線上的數(shù)據(jù)。選取X01、X09、X13、X17、X25 這5 個坐標(biāo)點的數(shù)據(jù)進(jìn)行對比，見表2[5]和表3。

計算得到的平均相對誤差為1.83%，最大相對誤差為4.44%。誤差來源主要有：簡化了熱壓罐模型和其他裝置；文獻(xiàn)[5]中的模具為較為平緩的曲面，但本文所用模具模型將其近化為平面；XFlow 設(shè)置的傳感器是試驗監(jiān)測模擬點的近似取值?？傊抡娴玫降臄?shù)據(jù)與試驗數(shù)據(jù)很接近，平均誤差較小，可以使用XFlow 軟件分析熱壓罐成型的模具溫度場。

表1 材料的熱屬性Table 1 Thermal properties of materials

圖3 框架式模具簡化結(jié)構(gòu)Fig.3 Simplified structure of frame die

圖4 框架式模具和熱壓罐模型Fig.4 Frame type mold and simplified autoclave model

圖5 試驗監(jiān)測點分布圖Fig5 Distribution map of experimental monitoring points

圖6 模具表面的溫度分布Fig.6 Temperature distribution of type plane mold

2 模具型面溫度場分析

在如圖7 所示的簡化工藝曲線下，設(shè)定風(fēng)速為1.5m/s，初始加壓加6個大氣壓，對熱壓罐成型過程模擬。

研究模具溫度場的最終目的是改善復(fù)合材料構(gòu)件在固化過程中的外部溫度場分布，模具上表面與復(fù)合材料構(gòu)件接觸，因此本節(jié)主要討論模具上表面的溫度場分布情況。升溫、保溫及降溫各階段結(jié)束時，模具上表面的溫度分布情況如圖8 所示。

在進(jìn)行實際試驗時，不方便測量模具上大量點的溫度值，一般通過放置熱電偶監(jiān)測個別點在試驗過程中的溫度。因此大多數(shù)學(xué)者使用最大溫差值表征模具的溫均性。但最大溫差值僅利用了溫度最大值和最小值兩個數(shù)據(jù)，不能反映模具整個上表面的溫均性。在有限元模擬時可以采集大量點在成型過程中的溫度變化，因此本文引入了標(biāo)準(zhǔn)差這一統(tǒng)計量，標(biāo)準(zhǔn)差將更多的數(shù)據(jù)引入計算，更加能反映一組數(shù)據(jù)的均勻性，溫度方程能更好地體現(xiàn)出模具上表面的溫均性。標(biāo)準(zhǔn)差是在幾何表面A上積分標(biāo)量與平均值的偏差，計算公式為：

表2 監(jiān)測點熱電偶試驗數(shù)據(jù)Table 2 Monitoring thermocouple test data

表3 監(jiān)測點XFlow模擬數(shù)據(jù)Table 3 Monitoring XFlow simulation data

通過XFlow 后處理中的表面積分功能可以計算模具上表面溫度標(biāo)準(zhǔn)差分布，將結(jié)果整理得到圖9。

結(jié)果表明：升溫階段，模具迎風(fēng)端的溫度高，最高點在迎風(fēng)端的邊角處，溫度最低處集中在背風(fēng)端，溫度呈階梯狀分布，溫差逐漸增大；保溫階段，溫度分布比較均勻，溫差隨時間持續(xù)減小；降溫階段，趨勢與升溫過程相反，背風(fēng)端溫度最高，迎風(fēng)端溫度最低，溫度呈階梯狀分布，溫差不斷增大。升溫時溫度分布呈階梯狀的原因主要是：模具型面的熱量主要來自于流經(jīng)上表面流體的對流換熱和底部支撐框格的熱傳導(dǎo)[10]。氣體流過模具表面可看作是氣體外掠平板，緊貼模具的氣體會形成熱邊界層，流體的黏性損耗了動能，因此邊界層中速度沿風(fēng)向降低，且邊界層厚度變厚，如圖10 所示邊界層越厚，熱阻越大，傳熱效率越低。底部支撐框格的熱量由流體的射流沖擊換熱產(chǎn)生，迎風(fēng)端框格阻礙氣體流向背風(fēng)端，迎風(fēng)端框格溫度高于背風(fēng)端與模具的傳熱更多。此外可以看出，模具上表面的前端和側(cè)面的溫度均高于內(nèi)側(cè)的溫度，原因是模具周邊和高溫氣流之間的對流傳熱更強(qiáng)。

工藝因素對模具型面溫度場的影響規(guī)律

熱壓罐工藝中可能對模具溫度場產(chǎn)生影響的工藝因素有：罐內(nèi)風(fēng)速、升溫速率及降溫速率等[12]。由圖9 可知，保溫階段罐內(nèi)溫差逐漸減小，溫差出現(xiàn)在升、降溫階段，因此在研究升降溫速率對模具表面的溫度影響規(guī)律時，本文僅對升溫階段和降溫階段進(jìn)行模擬。

1 風(fēng)速

在XFlow 中將風(fēng)速分別設(shè)置為：0.5m/s、1.5m/s、3m/s 和5m/s，入口溫度按圖7 設(shè)置。整理不同風(fēng)速下模具型面溫度標(biāo)準(zhǔn)差，如圖11、圖12所示。

結(jié)果表明：在升溫階段，模具型面的最大溫度標(biāo)準(zhǔn)差出現(xiàn)在升溫結(jié)束時；在升溫后期，風(fēng)速越大，模具型面的溫差最大值越小，溫差的增加越緩慢。在保溫階段，風(fēng)速越大，模具型面溫度趨于均勻所需的時間越短。在降溫階段，模具型面的最大溫度標(biāo)準(zhǔn)差出現(xiàn)降溫過程中，溫差在降溫后期停止增加，風(fēng)速越大，溫差最大值越小，出現(xiàn)時間越早。綜上，提高風(fēng)速可以有效減小模具型面溫差。

圖7 簡化工藝曲線Fig.7 Simplified cured temperature curve

圖8 模具上表面的溫度分布Fig.8 Temperature distribution of type plane mold

圖9 模具上表面溫度標(biāo)準(zhǔn)差分布Fig.9 Temperature standard deviation distribution of type plane mold

圖10 邊界層示意圖Fig.10 Boundary layer of plane

圖11 不同風(fēng)速下的溫度標(biāo)準(zhǔn)差對比Fig.11 Comparison of temperature standard deviation under different wind speeds

圖12 不同風(fēng)速下的最大溫度標(biāo)準(zhǔn)差Fig.12 Maximum standard deviation of temperature at different wind speeds

2 升溫速率

將升溫速率分別設(shè)為：1K/min、1.5K/min、2K/min 和3K/min，模具的初始溫度設(shè)為288K，經(jīng)不同升溫速率的升溫過程，溫度升高至468K 結(jié)束，風(fēng)速設(shè)為1.5m/s。整理不同升溫速率下模具型面溫度標(biāo)準(zhǔn)差，如圖13、圖14 所示。

結(jié)果表明：模具型面最大溫度標(biāo)準(zhǔn)差隨著升溫速率的增大而明顯增大，最大溫度標(biāo)準(zhǔn)差的值隨升溫速率的增長基本呈線性增長。

3 降溫速率

將降溫速率分別設(shè)為：0.5K/min、1K/min、1.5K/min 和2K/min，模具的初始溫度設(shè)為468K，經(jīng)不同降溫速率的降溫過程，溫度降低至288K 結(jié)束。整理不同降溫速率下模具型面溫度標(biāo)準(zhǔn)差，如圖15 和16所示。

結(jié)果表明：模具型面最大溫度標(biāo)準(zhǔn)差隨著降溫速率的增大而明顯增大，最大溫度標(biāo)準(zhǔn)差的值隨升溫速率的增長基本呈線性增長。

結(jié)論

（1）基于XFlow 軟件建立了熱壓罐成型模具溫度場仿真的數(shù)值模型，將模擬數(shù)據(jù)與試驗結(jié)果對比，平均誤差為1.83%，最大相對誤差為4.44%。

（2）分析了熱壓罐成型過程中模具表面溫度分布情況，引入了溫度標(biāo)準(zhǔn)差作為模具表面溫度均勻性的表征值，結(jié)果表明：溫度標(biāo)準(zhǔn)差在升、降溫階段不斷增大，保溫階段不斷減小。

圖13 不同升溫速率下的溫度標(biāo)準(zhǔn)差對比Fig.13 Comparison of temperature standard deviation at different heating rates

圖14 不同升溫速率下的最大溫度標(biāo)準(zhǔn)差Fig.14 Maximum standard deviation of temperature at different heating rates

圖15 不同降溫速率下的溫度標(biāo)準(zhǔn)差對比Fig.15 Comparison of temperature standard deviation at different cooling rates

圖16 不同降溫速率下的最大溫度標(biāo)準(zhǔn)差Fig.16 Maximum standard deviation of temperature at different cooling rates

（3）研究了風(fēng)速、升溫速率和降溫速率對模具型面溫度場均勻性的影響。結(jié)果表明：增大風(fēng)速、減小升降溫速率均可以有效降低模具型面溫度標(biāo)準(zhǔn)差。

（4）采用基于LBM方法的XFlow 流體力學(xué)軟件分析框架式模具在熱壓罐中模具表面溫度場分布，可以大幅減少前處理時間，并能夠提高計算精度。