基于蒙特卡罗法的C/SiC波纹点阵结构传热性能分析

Analysis of heat transfer performance of C/SiC corrugated lattice structure based on Monte Carlo method

摘要

为系统了解C/SiC材料的传热机理，以 C/SiC波纹点阵结构复合材料为研究对象，基于蒙特卡罗法对波纹点阵结构进行传热性能的理论分析。考虑能束发射位置、发射方向、遮挡判定这3个因素建立了波纹点阵结构的传热理论模型，分析了结构几何参数、固体发射率及外载施加温度对波纹点阵结构传热性能的影响。结合蒙特卡罗法预测了C/SiC波纹点阵结构的热传递系数，样本数量精确到4.0×10⁶，热传递系数收敛趋于0.554 9 \mathrm{W}/(\mathrm{m}\cdot \mathrm{K})。结果表明，芯子高度、宽度和倾斜角度对热传递系数的影响都呈正相关，而芯子厚度对热传递系数的影响呈负相关；热传递系数随着固体发射率和上面板所受温度载荷的增加而增加，且增加的趋势越来越大；对航空航天复合材料隔热性能研究具有一定的工程应用价值。

Abstract

In order to gain a theoretical understanding of the heat transfer mechanism of C/SiC materials，Monte Carlo method was applied to analyze heat transfer of C/SiC corrugated lattice structure composite material.The heat transfer mechanism of the corrugated lattice structure was studied from three aspects： the beam emission position，the beam emission direction and the judgment of the beam occlusion.At the same time，the influencing factors of the heat transfer coefficient of the lattice structure were studied by changing the structural geometric parameters，the emissivity of the solid surface and the temperature applied by the upper panel.The number of samples can be accurate to 4.0×10⁶ and the heat transfer coefficient converges to 0.554 9 \mathrm{W}/(\mathrm{m}\cdot \mathrm{K}) by using the heat transfer coefficient derived from the formula of Monte Carlo method.The results show that，the linear relationship between the heat transfer coefficient，the geometric parameters of the corrugated lattice and the reasons for its change are given，the optimal thermal insulation model parameters of the corrugated lattice structure are obtained.The theoretical research provides certain engineering application value for thermal insulation of aerospace material and composite material research.

关键词

C/SiC / 蒙特卡罗法 / 波纹点阵结构 / 传热性能 / 复合材料

Key words

C/SiC / Monte Carlo method / corrugated lattice structure / heat transfer performace / composite material

隔热结构是保护高速飞行器安全工作和完成各项任务的关键因素，有效的隔热结构可以满足高速飞行器隔热的基本需求，使其免受极端气动热环境的伤害^［1］。C/SiC复合材料展现出优异的机械性能，如轻量化、高比强度、高比刚度及卓越的耐高温和耐腐蚀能力^［2］。其优异的力学性能和抗氧化性能可作为高温复合材料部件的首选^［3］，在高效轻质防热一体化、航空航天及高端汽车领域被广泛应用^［4］。研究热端部件传热性能的理论方法有很多，作为随机抽样的数值模拟计算方法，蒙特卡罗法被广泛应用于求解复杂问题的实际计算中，该方法与理论公式相结合可准确预测实验结果^［5］。为了获取最佳隔热结构设计，国内外学者通过实验、理论分析和数值仿真等方法，作出了大量研究。

Dai等 ^［6］通过解析模型的方法计算C/SiC金字塔结构的夹芯板等效导热系数，并针对极限温度因素对金字塔夹层结构热传递影响的研究进行了深入分析和讨论。Yang等^［7］采用蒙特卡罗法建立了一种简化的二阶模型，将复杂的网格结构简化为三维杆系结构，并基于有限元原理对热传导、对流和辐射进行了简化的二阶模拟。随着杆长比的增大，模型的计算误差减小，计算效率提高1 000倍。Yao等^［8］通过实验方法测定了不同温度下C/SiC复合材料波纹夹芯板的等效导热系数，并对其传热特性进行了深入分析和阐释。Chan等^［9］开发了一个理论模型，用于确定充满填料的石墨烯涂层金属泡沫的有效导热率。Yang等^［10］采用分析建模和数值模拟相结合的方法，对高孔隙率、流体饱和金属丝桁架点阵结构的有效导热系数进行了估算。Kovtanyuk等^［11］提出了一种改进的蒙特卡罗法，用于由2个反射和辐射平面限定的介质中的非线性耦合辐射-传导传热问题的数值处理。Lin等^［12］将凹面三棱柱泡沫韧带引入到四面体泡孔结构中，并考虑了韧带取向和填充介质的影响，提出一种新的有效导热系数预测模型。吴一昊^［13］研究了金属蜂窝结构的热传方式和理论，建立了有限元传热模型，将理论和数值模拟相结合研究传热原理。Wei等^［14］针对点阵结构在工程应用中热膨胀的影响，提出了一种零热膨胀系数（coefficient of thermal expansion，CTE）的三维桁架结构设计方法。

这些研究工作在点阵结构隔热方面作出了大量研究，但采用蒙特卡罗法对点阵结构传热性能理论的研究较少。因此，本文基于蒙特卡罗法对波纹点阵结构传热性能进行理论分析。

1 理论模型

1.1 热传输机理

波纹点阵结构由中心的波纹芯子结构与上、下复合材料层合板共同构成，这种结构设计能够在保证可靠性的同时降低部件质量，以此来减少结构件质量。

图1为波纹点阵结构传热物理模型图，上面板的温度为T _u，下面板的温度为T _d，T _u大于T _d。由于上、下面板厚度与中间波纹点阵芯子高度相比非常小，因此在考虑传热问题时，忽略上、下面板厚度对热传递系数的影响，同时认为在上面板施加的温度载荷是均匀分布的。将波纹点阵结构的传热方式分为以下4部分：

1）C/SiC复合材料固体之间的热传导。

2）波纹点阵结构芯子中气体的热传导。

3）传至下面板后气体之间的热对流。

4）芯子外表面与上下面板内表面间的热辐射。

1.2 热传导模型

C/SiC波纹点阵结构复合材料的主要热传递机制包括固体间的热传导和波纹芯子表面与上、下复合材料层合板之间的辐射换热。通过经典传热公式推导得出一个简洁而有效的式（1），用以表达C/SiC波纹点阵结构复合材料的等效导热系数，该式能更准确地理解和预测其热传导性能。

k_{\mathrm{e}\mathrm{q}\mathrm{u}} = k_{\mathrm{c}} + k_{\mathrm{r}} （1）

式中： k_{\mathrm{e}\mathrm{q}\mathrm{u}}为C/SiC波纹点阵结构复合材料的等效导热系数； k_{\mathrm{c}}为波纹点阵结构导热系数； k_{\mathrm{r}}为C/SiC复合材料热辐射系数。

根据傅立叶热传导定律，可以推导出波纹点阵结构导热系数如式（2）所示。

k_{\mathrm{c}} = k_{\mathrm{s}} \overline{\rho } \mathrm{s} \mathrm{i} {\mathrm{n}}^{\mathrm{2}} \theta （2）

式中： k_{\mathrm{s}}为C/SiC复合材料的母体热导率； \overline{\rho }为波纹点阵结构复合材料的无量纲化相对密度，由式（3）表示。式（3）中，单胞长度l和宽度w分别由式（4）、（5）表示。

\overline{\rho } = \frac{\mathrm{4}ad}{lh} + \frac{\mathrm{4}bt\left(h-\mathrm{2}d\right)}{wlh\mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{c}\mathrm{o}\mathrm{s}\omega } （3）

l = \mathrm{4} a + \frac{\mathrm{2}\left(h-\mathrm{2}d\right)}{\mathrm{t}\mathrm{a}\mathrm{n}\theta } （4）

w = \mathrm{2} \left(\frac{\left(h-\mathrm{2}d\right)\mathrm{t}\mathrm{a}\mathrm{n}\omega }{\mathrm{s}\mathrm{i}\mathrm{n}\theta }+\frac{b}{\mathrm{s}\mathrm{i}\mathrm{n}\omega }\right) （5）

C/SiC复合材料波纹点阵结构的几何参数如表1所示，波纹点阵结构整体示意图如图2所示，芯子单胞示意图如图3所示。

1.3 等效导热系数

为了探究C/SiC波纹点阵结构复合材料的传热机理，本文对传热问题进行深度讨论后作出以下假设：1）忽略光的衍射效果对波纹点阵结构传热所造成的影响；2）由于C/SiC复合材料上、下表面极为粗糙，因此镜面反射对传热的影响可以忽略不计；3）C/SiC复合材料为灰体，因此光谱无法对波纹点阵结构产生透射。

在上述假设成立的条件下，可将模型的结构热传递系数表示为式（6）。

k_{\mathrm{r}} = \frac{q_{\mathrm{r},\mathrm{n}\mathrm{e}\mathrm{t}}}{\left(T_{\mathrm{u}}-T_{\mathrm{d}}\right)h} （6）

综上，结合式（2）、（3）和（6）整理化简得到波纹点阵结构模型的等效热传递系数表达式，如式（7）所示。

k_{\mathrm{e}\mathrm{q}\mathrm{u}} = \frac{k_{\mathrm{s}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta }{\mathrm{2}h\left(h-\mathrm{2}d\right)} + \frac{\mathrm{4}bt{k}_{\mathrm{s}}\mathrm{s}\mathrm{i}\mathrm{n}\theta \left(h-\mathrm{2}d\right)}{h\mathrm{c}\mathrm{o}\mathrm{s}\omega } + \epsilon F_{\mathrm{1}-\mathrm{2}} \sigma h \left(T_{\mathrm{u}}^{\mathrm{2}}+T_{\mathrm{d}}^{\mathrm{2}}\right) \left(T_{\mathrm{u}}+T_{\mathrm{d}}\right) （7）

式中： \epsilon为C/SiC复合材料表面的固体发射率； F_{\mathrm{1}-\mathrm{2}}为C/SiC复合材料层合板上、下面板的配置因子； \sigma为常数，其值为 \mathrm{5.669} \times {\mathrm{10}}^{-\mathrm{8}}； h为芯子高度。

2 蒙特卡罗法

作为统计模拟法，蒙特卡罗法是将概率现象作为研究对象的模拟方法，本文通过将蒙特卡罗法的特性与微观传热相结合的等效方式来研究C/SiC波纹点阵结构模型的热传递系数。 F_{i-j}表示上下面板配置因子，其表达式如式（8）所示。

F_{i-j} = \frac{m}{N} （8）

式中：N为从波纹点阵结构上表面i所发射的能束数量；m为被波纹点阵结构下表面j所吸收到的能束数量。

2.1 能束发射位置的概率分布计算

定义上面板表面i的所有能束均携带相等的能量 \frac{E_j}{N}，若用 P \left(x\right)表示在上面板x和y之间区域的能束发射位置分布概率方程，则 R_x表示在坐标a和x之间区域能束发射出辐射能与上面板总辐射能的比值，能束在上面板表面发射位置平面概率分布如图4所示， R_x表达式如式（9）所示。

R_x = {\int }_a^x P\left(x\right) \mathrm{d} x = \frac{\mathrm{1}}{E_j} {\int }_{x=a}^{x=b} E_j \mathrm{d} x （9）

考虑到C/SiC波纹点阵结构模型中，上、下面板均呈现7×7的二维矩形平面，且这些平面的材料属性和所施加的温度载荷均匀分布，可以认为能束的发射位置是在上面板均匀分布的。因此，发射点的概率分布公式如式（10）、（11）所示。

R_x = \left(X_{\mathrm{m}\mathrm{a}\mathrm{x}}-X_{\mathrm{m}\mathrm{i}\mathrm{n}}\right) + X_{\mathrm{m}\mathrm{i}\mathrm{n}} （10）

R_y = \left(Y_{\mathrm{m}\mathrm{a}\mathrm{x}}-Y_{\mathrm{m}\mathrm{i}\mathrm{n}}\right) + Y_{\mathrm{m}\mathrm{i}\mathrm{n}} （11）

式中：自变量x的取值范围为 \left[X_{\mathrm{m}\mathrm{i}\mathrm{n}},X_{\mathrm{m}\mathrm{a}\mathrm{x}}\right]；y的取值范围为 \left[Y_{\mathrm{m}\mathrm{i}\mathrm{n}},Y_{\mathrm{m}\mathrm{a}\mathrm{x}}\right]，且R_x 和R_y 均遵循［0，1］区间内的均匀分布。

2.2 能束发射方向的概率分布计算

该模型辐射发射方向的概率分布用球坐标方程可表示为

E_{\lambda } = \frac{\mathrm{1}}{\mathrm{\pi }} {\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }} {\int }_{\mathrm{0}}^{\frac{\mathrm{\pi }}{\mathrm{2}}}\epsilon \lambda \mathrm{c}\mathrm{o}\mathrm{s}\phi \mathrm{s}\mathrm{i}\mathrm{n}\phi \mathrm{d} \phi \mathrm{d} \theta （12）

能束发射方向的概率分布仅与上面板表面的发射率有关，与能束的具体发射位置无关。将能束点理解为半球结构，则其上面板能束发射方向示意图如图5所示。

2.3 能束遮挡的判定

定义 l_{ij}为长方体i的质心最靠近长方体j之间的距离；同理 l_{ji}为长方体j的质心最靠近长方体i之间的距离，其表达式如式（13）、（14）所示。

l_{ij} = \frac{p_i\left[\left(R_j-R_i\right)\left(\mathrm{1}-p_i^{\mathrm{2}}\right)\right]}{\mathrm{1}-{\left(p_i{p}_j\right)}^{\mathrm{2}}} （13）

l_{ji} = \frac{p_j\left[\left(R_j-R_i\right)\left(p_i^{\mathrm{2}}-\mathrm{1}\right)\right]}{\mathrm{1}-{\left(p_i{p}_j\right)}^{\mathrm{2}}} （14）

当且仅当 \left|l_{ij}\right| \le \frac{\mathrm{1}}{\mathrm{2}} ‖P_i‖和 \left|l_{ji}\right| \le \frac{\mathrm{1}}{\mathrm{2}} ‖P_j‖两个条件同时满足时，可以得出长方体i与长方体j之间存在三维相交，从而推断能束被芯子杆遮挡。如若上述条件有一条不满足，则芯子杆并未与能束相交。

蒙特卡罗法选取样本的数量会直接影响结果的可靠性。图6为不同能束基数对波纹点阵结构等效热传递系数的影响，选取的单胞高度为12 mm，加强台厚度为1.5 mm。当能束的数量较少时，热传递系数并未收敛。随着能束数量的增加，热传递系数呈现递减趋势，当能束数量超过4.0×10⁶时，热传递系数开始趋于稳定，最终稳定于0.554 9 \mathrm{W}/(\mathrm{m}\cdot \mathrm{K})。

3 结果与讨论

3.1 上面板温度对C/SiC材料的影响

为验证本文基于蒙特卡罗法建立的波纹点阵结构传热模型的预测有效性，对比了理论推导的结构热传递系数、基于蒙特卡罗法预测的结构等效导热系数及Wei等^［15］实验测得的结构导热系数。

图7为C/SiC波纹点阵结构等效导热系数值对比图。从图7可以看出，基于蒙特卡罗法的热传递系数预测值位于理论值和实验值之间，其大小更贴合于热传递系数的理论计算值。在样本数量过少时，热传递系数过于发散，无法预测出精确的数据。由于样本基数大于4.0×10⁶时模型处于收敛状态，热传递系数趋于0.554 9 \mathrm{W}/(\mathrm{m}\cdot \mathrm{K})。计算出来的结果更趋近于热传递系数的真实值。

3.2 结构几何参数对热传递系数的影响

为研究结构几何参数对波纹点阵结构的等效导热系数的影响，在保证波纹点阵结构夹芯板体积分数不变的条件下，通过变化加强台面长度参数以维持所有构型的芯子相对密度基本保持恒定。图8为不同几何参数对热传递系数的影响。由图8可知，芯子高度、宽度和倾斜角。对热传递系数的影响都呈正相关，而芯子厚度对热传递系数的影响呈负相关。随着芯子厚度的增加，能束在上、下面板穿透过程中损失更多的能量，从而使能量束在下面板接收位置较芯子厚度小的夹芯板会减少很多。为了提高不同几何参数对热传递系数的影响因子可靠性，考虑了相关分散性对热传递系数的影响，如图8所示。在热传递系数从4.48 \mathrm{W}/(\mathrm{m}\cdot \mathrm{K})上升至4.60 \mathrm{W}/(\mathrm{m}\cdot \mathrm{K})中，芯子厚度的变化趋势极小可忽略，芯子倾斜角的上升趋势非常明显。

3.3 固体发射率对热传递系数的影响

在相同粗糙度条件下，改变固体发射率即改变上、下面板及材料内部的固体发射率属性来研究固体发射率对热传递系数的影响，如图9所示。由图9可知，C/SiC复合材料波纹点阵结构的热传递系数随着固体发射率的增加而增加，且增加的趋势越来越大。

研究结果表明，固体发射率增大，导致其反射率减小，由上面板发射出的能束更多抵达下面板，少部分能束发生反射，因此发射在总辐射中所占的比例较大。

3.4 温度对热传递系数的影响

图10为温度对热传递系数的影响，根据图10中上面板所受不同温度条件下对热传递系数的影响，得出波纹点阵结构的热传递系数随着上面板温度的增加而增大的结论。这种增大呈现非线性增长的趋势，且随着温度升高，非线性趋势更加明显。导致这种现象的原因是SiC基体材料的导热系数随温度升高而下降，由于下降趋势不明显，因此对波纹点阵结构的等效导热系数的影响较小。

4 结论

本文构建了综合考虑热传导和热辐射效应的C/SiC复合材料波纹点阵结构的热传递模型。基于蒙特卡罗法对能束发射的位置、方向及遮挡判定进行准确描述，可将样本数量精确到4.0×10⁶，热传递系数收敛趋于0.554 9 \mathrm{W}/(\mathrm{m}\cdot \mathrm{K})。给出了细观结构尺寸、温度、材料发射率等因素对波纹点阵夹层结构等效导热系数的影响规律，对比分析其隔热性能。得出以下结论：

1）研究C/SiC复合材料波纹点阵结构的几何参数（芯子高度、芯子宽度和芯子倾斜角）对热传递系数的影响。造成影响的原因是在中间芯子板穿透过程中损失部分能量束，从而使得能量束在下面板接收位置较上面板发射位置会减少很多。

2）研究了固体发射率对热传递系数的影响。研究结果表明，随着固体发射率的增加，C/SiC复合材料波纹点阵结构的热传递系数也随之增加，且增加的趋势越来越大。这是由于固体发射率越大，其反射率越弱。

3）研究了温度梯度对波纹点阵结构热传递系数的影响，其热传递系数呈现非线性增长的趋势，且随着温度的升高，非线性趋势更加明显。

参考文献

[1] Daryabeigi K.Heat transfer in adhesively bonded honeycomb core panels［J］.Journal of Thermophysics and Heat Transfer，2002，16（2）：217-221.

[2] Zhao C Y， Tassou S A， Lu T J.Analytical considerations of thermal radiation in cellular metal foams with open cells［J］.International Journal of Heat and Mass Transfer，2008，51（3/4）：929-940.

[3] Fu J W， Akbarzadeh A H， Chen Z T，et al.Non-fourier heat conduction in a sandwich panel with a cracked foam core［J］.International Journal of Thermal Sciences，2016，102（11）：263-273.

[4] Wei K， Peng Y， Qu Z L，et al.High temperature mechanical behaviors of lightweight ceramic corrugated core sandwich panel［J］.Composite Structures，2017，176（5）：379-387.

[5] Boomsma K， Poulikakos D.On the effective thermal conductivity of a three-dimensionally structured fluid-saturated metal foam［J］.International Journal of Heat and Mass Transfer，2001，44（4）：827-836.

[6] Dai Z， Nawaz K， Park Y G，et al.Correcting and extending the boomsma-poulikakos effective thermal conductivity model for three-dimensional，fluid-saturated metal foams［J］.International Communications in Heat and Mass Transfer，2010，37（6）：575-580.

[7] Yang X H， Kuang J J， Lu T J，et al.A simplistic analytical unit cell based model for the effective thermal conductivity of high porosity open-cell metal foams［J］.Journal of Physics D：Applied Physics，2013，46（25）：255302.

[8] Yao Y P， Wu H Y， Liu Z Y.A new prediction model for the effective thermal conductivity of high porosity open-cell metal foams［J］.International Journal of Thermal Sciences，2015，97（6）：56-67.

[9] Chan K C， Tso C Y， Hussain A，et al.A theoretical model for the effective thermal conductivity of graphene coated metal foams［J］.Applied Thermal Engineering，2019，161：114112.

[10] Yang X H， Bai J X， Kang K J，et al.Effective thermal conductivity of wire-woven bulk Kagome sandwich panels［J］.Theoretical and Applied Mechanics Letters，2014，4（5）：051010.

[11] Kovtanyuk A E， Botkin N D， Hoffmann K H.Numerical simulations of a coupled radiative-conductive heat transfer model using a modified Monte Carlo method［J］.International Journal of Heat and Mass Transfer，2012，55（4）：649-654.

[12] Lin X， Zheng C， Huang W，et al.A novel heliostat aiming optimization framework via differentiable Monte Carlo ray tracing for solar power tower systems［J］.Applied Energy，2025，388：125640.

[13] 吴一昊.陶瓷梯度点阵结构热机械损伤行为研究［D］.哈尔滨：哈尔滨工业大学，2016.

[14] Wei K， Peng Y， Wang K Y，et al.Three dimensional lightweight lattice structures with large positive，zero and negative thermal expansion［J］.Composite Structures，2018，188（1）：287-296.

[15] Wei K， He R J， Cheng X M，et al.Fabrication and heat transfer characteristics of C/SiC pyramidal core lattice sandwich panel［J］.Applied Thermal Engineering，2015，81（2）：10-17.

基金资助

国家自然科学基金(11972141)

辽宁省教育厅系列项目(JYT2020121)

航空科学基金(2020Z006054002)

国防重点学科实验室项目(SHSYS202208)