Processing math: 8%

基于增量微分求积单元法的功能梯度材料夹层板非线性瞬态传热分析

张忠, 许家婧, 曹小建, 王艳超, 朱军, 姚潞

张忠, 许家婧, 曹小建, 等. 基于增量微分求积单元法的功能梯度材料夹层板非线性瞬态传热分析[J]. 复合材料学报, 2024, 41(11): 6284-6296. DOI: 10.13801/j.cnki.fhclxb.20240228.003
引用本文: 张忠, 许家婧, 曹小建, 等. 基于增量微分求积单元法的功能梯度材料夹层板非线性瞬态传热分析[J]. 复合材料学报, 2024, 41(11): 6284-6296. DOI: 10.13801/j.cnki.fhclxb.20240228.003
ZHANG Zhong, XU Jiajing, CAO Xiaojian, et al. Nonlinear transient heat transfer analysis of functionally graded material sandwich slabs by incremental differential quadrature element method[J]. Acta Materiae Compositae Sinica, 2024, 41(11): 6284-6296. DOI: 10.13801/j.cnki.fhclxb.20240228.003
Citation: ZHANG Zhong, XU Jiajing, CAO Xiaojian, et al. Nonlinear transient heat transfer analysis of functionally graded material sandwich slabs by incremental differential quadrature element method[J]. Acta Materiae Compositae Sinica, 2024, 41(11): 6284-6296. DOI: 10.13801/j.cnki.fhclxb.20240228.003

基于增量微分求积单元法的功能梯度材料夹层板非线性瞬态传热分析

基金项目: 国家自然科学基金(11802145;52308260);江苏省高等学校基础科学(自然科学)研究项目(23KJB560021);南通市科技项目(JC12022058);南通市社会民生科技计划面上研究项目(MS22022103)
详细信息
    通讯作者:

    姚潞,博士,讲师,研究方向为复合材料结构 E-mail: yaolu@ntu.edu.cn

  • 中图分类号: TK124;TB33

Nonlinear transient heat transfer analysis of functionally graded material sandwich slabs by incremental differential quadrature element method

Funds: National Natural Science Foundation of China (11802145; 52308260); Natural Science Foundation of the Jiangsu Higher Education Institutions of China (23KJB560021); Science and Technology Project of Nantong City (JC12022058); Nantong City Social Livelihood Science and Technology Project (MS22022103)
  • 摘要:

    作为首次尝试,采用增量微分求积单元法(IDQEM)开展了功能梯度材料(FGM)夹层板的一维非线性瞬态传热分析。夹层板组分材料的热工参数随空间位置变化,且具有温度依赖性。基于IDQEM,沿层界面将夹层板划分为3个空间子域,同时将整个受热过程划分为若干时间子域。采用微分求积技术对任一时间子域内的控制方程、初始条件、界面条件及边界条件进行离散处理。由于所获得的离散方程建立在不同区域的节点上,因此对方程进行修改并将其表示为矩阵形式,以便它们可以建立在同一区域中。采用Kronecker积将联立的矩阵方程转化为一系列代数方程组,并采用Newton-Raphson迭代法近似求解,即可获得单个时间子域内的温度解。由于每个时间子域的初始条件可由上一个时间子域最终时刻的温度分布决定,因此从第一个时间子域逐渐递推到最后一个子域,即可获得整个受热过程的温度分布。数值算例验证了本方法的快速收敛性,与已有文献的解析和数值结果的对比验证了本方法的正确性。最后,讨论了热工参数温度依赖性、体积分数指数及热边界条件对FGM夹层板温度分布的影响。

     

    Abstract:

    As a first attempt, the incremental differential quadrature element method (IDQEM) was adopted to perform the one-dimensional nonlinear transient heat transfer analysis of functionally graded material (FGM) sandwich slabs. The thermophysical properties of the slab were considered to be position- and temperature-dependent. To implement the IDQEM, the sandwich slab was divided into three spatial sub-domains along the layer interfaces, and the entire heating process was also divided into several temporal sub-domains. For each temporal sub-domain, the governing equations as well as the initial condition, interfacial condition, and boundary condition were discretized by the differential quadrature technique. Because the obtained discrete equations were built in different regions of grid points, a modification of the equations was proposed which were then expressed in the matrix forms so that they can be built in the same regions. Using the Kronecker product, the simultaneous matrix equations were transformed into a set of nonlinear algebraic equations, which were then solved by the Newton-Raphson iteration method to obtain the temperature profile for each temporal sub-domain. Because the initial condition of each temporal sub-domain was defined by the temperature results at the end of the previous sub-domain, the temperature profile of the slab during the entire heating process can be obtained by repeating the calculation procedure from the first temporal sub-domain to the last one. Numerical examples were carried out to verify the fast convergence of the present method. The correctness of the present method was verified through comparison with the analytical and numerical results reported in previous works. The effects of temperature-dependent thermophysical properties, volume fraction index, and thermal boundary on the temperature profile of the slab were discussed.

     

  • 夹层结构因其轻质、热力学性能优异等特点,近几十年来在飞行器、核反应堆、集成电路等领域得到广泛应用。与蜂窝[1]、泡沫[2-3]等夹层结构相比,功能梯度材料(FGM)夹层结构通过使组分材料在空间上连续分布,可消除相邻层界面的接触热阻和应力失配,因此FGM夹层结构常被用于高温工程领域。深入了解FGM夹层结构的传热过程,对其设计和制造具有重要意义。

    在FGM夹层结构的初步设计中,由于结构厚度较其他两个方向的尺寸小得多,因此常常将传热问题简化为沿厚度方向的一维问题[4]。Daikh等[5]给出了Dirichlet边界条件下FGM夹层板的一维稳态热传导问题的解析解。Dini等[6]考虑表面热对流和内部热源的作用,研究了旋转FGM夹层圆盘的径向稳态温度分布。由于高温会严重影响FGM的热工参数,因此传热分析需要考虑热工参数温度依赖性的影响,这将导致传热问题变成非线性问题。Shen和Li[7]考虑组分材料热工参数的温度依赖性,研究了FGM夹层板的非线性稳态热传导问题。

    以上工作仅研究了FGM夹层结构的稳态传热问题。然而,稳态分析得到的结构内温度分布处于平衡状态,这在整个加热过程中可能无法实现或者维持。因此,稳态分析可能会高估结构内温度分布,导致设计过于保守。为了能尽可能减少结构自重,有必要对非线性瞬态传热过程进行准确预测。Pandey等[8-9]采用有限差分法(FDM)研究了Dirichlet边界条件下FGM夹层板的一维非线性瞬态热传导问题。结果表明,数值结果与文献结果吻合良好。然而,Malekzadeh等[10]指出,传统数值方法如FDM[8-9]和有限元法[11]需要非常小的时间步长,以确保结果的收敛性和准确性,但这样会导致计算耗时。

    微分求积法(DQM)是一种求解微分方程的数值方法,该方法无需像有限元法那样推导弱形式的控制方程,因此可以大幅减少高阶近似中的公式化工作量[12]。相较于传统数值方法如FDM,DQM已被证明具有高精度和高计算效率的优势[10]。在已有工作中,学者主要采用DQM求解单域问题[13-15]。然而,DQM并不适用于求解由于荷载、材料和尺寸等不连续造成的多域问题。为了求解这类多域问题,Wang等[16]和Chen等[17]基于传统的DQM提出了微分求积单元法(DQEM)。基于DQEM,Malekzadeh等[18]分析了FGM夹层壳的稳态传热问题。Dai等[19]通过在空间域上应用DQEM并在时间域上应用Newmark法,给出了含FGM层的多层热防护系统的瞬态温度分布。然而Malekzadeh等[18]和Dai等[19]仅仅考虑了热传导系数的温度依赖性,却假定其他热工参数(如密度和比热容)为常量。在高温环境中,这种假设可能会造成一定的误差。

    本文首次尝试采用增量微分求积单元法(IDQEM)分析时变热边界条件下FGM夹层板的一维非线性瞬态传热问题。考虑了所有热工参数的温度依赖性及时变热边界条件,以提高板内温度场预测的准确性。为应用IDQEM,将时间域和空间域都划分为若干子域。采用微分求积技术对任一时间子域内的控制方程、初始条件、界面条件及边界条件进行离散处理,并基于矩阵原理对离散方程进行修改,最后采用Newton-Raphson迭代法近似求解,获得单个时间子域内的夹层板温度解。采用递推的方式,即可获得整个受热过程的温度分布。数值算例验证了本方法的快速收敛和高精度优势,并讨论了一些关键参数对FGM夹层板温度分布的影响。

    图1为一块FGM夹层板,其底层、芯层和顶层分别为金属单相材料、FGM材料和陶瓷单相材料。夹层板的总厚度为H,各单层厚度为Hs (s = 1、2、3)。考虑组分材料热工参数的温度依赖性,通常表述为[20]

    Pm(T)=Pm0(Pm1T1+1+Pm1T+Pm2T2+Pm3T3),Pc(T)=Pc0(Pc1T1+1+Pc1T+Pc2T2+Pc3T3) (1)

    其中:PmPc分别代表金属和陶瓷的有效热工参数;PmkPck(k=10123)分别为金属和陶瓷热工参数的温度系数;T为温度。令ˉPm0=Pm0ˉPc0=Pc0ˉPmj=PmjPm0ˉPcj=PcjPc0 (j=1123),则式(1)可以改写为

    Pm(T)=3k=1ˉPmkTk,Pc(T)=3k=1ˉPckTk (2)
    图  1  功能梯度材料(FGM)夹层板示意图
    Figure  1.  Schematic view of the functionally graded material (FGM) sandwich slab
    x—Global coordinate for the sandwich slab; x(s) (s = 1, 2, 3)—Local coordinate for the sth layer; H—Total thickness of the sandwich slab;HsThickness of the sth layer

    对于夹层板的任意第s (s = 1、2、3)层,可基于局部坐标x(s)(图1(b)),采用Voigt混合幂率模型来表征其有效热工参数P(s),即

    P(s)(x(s),T)=Pm(T)+[Pc(T)Pm(T)]V(s)(x(s)) (3)

    其中,V(s)为第s层陶瓷材料的体积分数。对于不同层,V(s)的取值分别为

    V(1)(x(1))=0, 0 (4)

    其中,η为FGM层的体积分数指数。

    本文FGM夹层板的传热分析基于以下两个假设:(1)上下表面的热边界条件不沿板面方向变化;(2)板厚相对于长度和宽度方向尺寸足够小。因此,板的传热问题可以简化为沿厚度方向的一维问题。基于局部坐标 {x^{(s)}} ,第s层的热传导方程为[8]

    \begin{split} & \left[\mathit{\lambda}^{(s)}(x^{(s)},T)T_{,x}^{(s)}(t,x^{(s)})\right]_{,x}= \\ & \quad\rho^{(s)}(x^{(s)},T)c^{(s)}(x^{(s)},T)\dot{T}^{(s)}(t,x^{(s)}) \end{split} (5)

    其中:变量顶部的点表示对时间的偏导;下标中的逗号表示对其后位置坐标的偏导;t为时间;T(s)λ(s)ρ(s)c(s)分别为第s层的温度、热传导系数、密度、比热容。

    假设初始时刻板内温度 T_0^{(s)} 任意分布,即

    {T^{(s)}}(0,{x^{(s)}}) = T_0^{(s)}({x^{(s)}}) (6)

    由相邻层的温度和热通量连续条件,得

    {T^{(s)}}(t,{H_s}) = {T^{(s + 1)}}(t,0) ,\quad s = 1, 2 (7)
    T_{,x}^{(s)}(t,{H_s}) = T_{,x}^{(s + 1)}(t,0) ,\quad s = 1, 2 (8)

    假设夹层板底面受自然对流作用,而上表面受气动加热和热辐射共同作用,则边界条件为

    - {\lambda ^{(1)}}(0,T)T_{,x}^{(1)}(t,0) = {h_{{\text{bot}}}}\left[ {{T^{(1)}}(t,0) - {T_{{\text{bot}}}}(t)} \right] (9)
    {\lambda ^{(3)}}({H_3},T)T_{,x}^{(3)}(t,{H_3}) =
    \qquad{q_{{\text{top}}}}(t) - \sigma \varepsilon \left\{ {{{\left[ {{T^{(3)}}(t,{H_3})} \right]}^4} - T_{{\text{top}}}^4(t)} \right\} (10)

    其中:Tbot为下部时变环境温度;hbot为下表面对流换热系数;Ttop为上部时变环境温度;qtop为上表面气动加热生成的时变热流;σ为Stephan-Boltzmann常数;ε为辐射率。需要说明的是,本方法还可将其他类型的热边界条件考虑在内,简单起见这里不再赘述。

    式(5)~(10)形成了FGM夹层板的瞬态传热定解问题。由于材料热工参数具有温度依赖性,同时辐射条件式(10)中含有温度的4次方项,因此该问题属于非线性问题,可采用IDQEM求解。为应用该方法,首先将整个加热过程划分为p (p ≥ 1)个时间子域,其中第r (r = 1、2、···、p)个子域的时长 \Delta {t_r} 由时变热边界条件确定;同时沿层界面将夹层板划分为3个空间子域,如图2所示。将第r个时间子域和第s个空间子域分别离散为 {M_r} {N_s} 个节点,则时间域和空间域上的总节点数分别为 M = \displaystyle\sum\nolimits_{r = 1}^p {{M_r}} N = \displaystyle\sum\nolimits_{s = 1}^3 {{N_s}} 个。节点分布规则依据Gauss-Lobatto-Chebyshev点[21]

    t_i^{(r)} = \frac{{\Delta {t_r}}}{2}\left[ {1 - \cos \frac{{\left( {i - 1} \right){\text{π}} }}{{{M_r} - 1}}} \right] , i = 1, 2, …, {M_r} ;
    x_j^{(s)} = \frac{{{H_s}}}{2}\left[ {1 - \cos \frac{{\left( {j - 1} \right){\text{π}} }}{{{N_s} - 1}}} \right] , j = 1, 2, …, {N_s} (11)
    图  2  r个时间子域下的节点分布
    Figure  2.  Grid point distribution in the rth temporal sub-domain
    Mr and NsGrid point numbers in the rth temporal sub-domain and the sth spatial sub-domain, respectively; t—Time variable

    则节点 (t_i^{(r)},x_j^{(s)}) 处的温度导数可以离散为[21]

    \frac{{{\partial ^{{k_t}}}}}{{\partial {{({t^{(r)}})}^{{k_t}}}}}\frac{{{\partial ^{{k_x}}}}}{{\partial {{({x^{(s)}})}^{{k_x}}}}}T_{i,j}^{(r,s)} = \sum\limits_{m = 1}^{{M_r}} {\sum\limits_{n = 1}^{{N_s}} {A_{i,m}^{({k_t},r)}B_{j,n}^{({k_x},s)}T_{m,n}^{(r,s)}} } (12)

    其中:T_{i,j}^{(r,s)} = T(t_i^{(r)},x_j^{(s)}) A_{i,m}^{({k_t},r)} 为待求温度在节点处对tkt阶偏导权系数;B_{j,n}^{({k_x},s)} 为待求温度在节点处对xkx阶偏导权系数。

    将式(12)代入式(5),则第r(r = 1、2、···、p)个时间子域的控制方程离散形式为

    \begin{split} \tilde \lambda _{i,j}^{(r,s)}&{\left( {\sum\limits_{n = 1}^{{N_s}} {B_{j,n}^{(1,s)}T_{i,n}^{(r,s)}} } \right)^2} + \hat \lambda _{i,j}^{(r,s)}\sum\limits_{n = 1}^{{N_s}} {B_{j,n}^{(1,s)}T_{i,n}^{(r,s)}}+\\ & \lambda _{i,j}^{(r,s)}\sum\limits_{n = 1}^{{N_s}} {B_{j,n}^{(2,s)}T_{i,n}^{(r,s)}} = \rho _{i,j}^{(r,s)}c_{i,j}^{(r,s)}\sum\limits_{m = 1}^{{M_r}} {A_{i,m}^{(1,r)}T_{m,j}^{(r,s)}} , \\ &i = 2, 3, …, {M_r} ; j = 2, 3, …, {N_s} - 1 ; s = 1, 2, 3 \end{split} (13)

    其中: \lambda _{i,j}^{(r,s)} \rho _{i,j}^{(r,s)} c_{i,j}^{(r,s)} 为相应热工参数在节点 (t_i^{(r)},x_j^{(s)}) 处的值; \tilde \lambda _{i,j}^{(r,s)} \hat \lambda _{i,j}^{(r,s)} 分别为 {\lambda ^{(r,s)}} 在节点处对Tx的偏导。由式(3), \lambda _{i,j}^{(r,s)} \rho _{i,j}^{(r,s)} c_{i,j}^{(r,s)} \tilde \lambda _{i,j}^{(r,s)} \hat \lambda _{i,j}^{(r,s)} 可以分别表示为

    \begin{gathered} P_{i,j}^{(r,s)} = \left[ {\sum\limits_{k = - 1}^3 {\left( {\bar P_k^{\text{c}} - \bar P_k^{\text{m}}} \right){{\left( {T_{i,j}^{(r,s)}} \right)}^k}} } \right]{V^{(s)}}(x_j^{(s)}) + \\ \qquad \sum\limits_{k = - 1}^3 {\bar P_k^{\text{m}}{{\left( {T_{i,j}^{(r,s)}} \right)}^k}} ,{\text{ }}P = \lambda ,{\text{ }}\rho ,{\text{ }}c, \\ \end{gathered}
    \begin{gathered} \tilde \lambda _{i,j}^{(r,s)} = \left[ {\sum\limits_{k = - 1}^3 {k\left( {\bar \lambda _k^{\text{c}} - \bar \lambda _k^{\text{m}}} \right){{\left( {T_{i,j}^{(s)}} \right)}^{k - 1}}} } \right]{V^{(s)}}(x_j^{(s)}) + \\ \qquad\sum\limits_{k = - 1}^3 {k\bar \lambda _k^{\text{m}}{{\left( {T_{i,j}^{(s)}} \right)}^{k - 1}}} , \\ \end{gathered}
    \hat \lambda _{i,j}^{(r,s)} = \left[ {\sum\limits_{k = - 1}^3 {k\left( {\bar \lambda _k^{\text{c}} - \bar \lambda _k^{\text{m}}} \right){{\left( {T_{i,j}^{(s)}} \right)}^{k - 1}}} } \right]{\hat V^{(s)}}(x_j^{(s)}) (14)

    其中, {\hat V^{(s)}}(x_j^{(s)}) {V^{(s)}} 在节点 x = x_j^{(s)} 处对x的导数。根据 {V^{(s)}} 的表达式(4),可得 {\hat V^{(s)}}(x_j^{(s)})

    {\hat V^{(1)}}(x_j^{(1)}) = 0 , {\hat V^{(2)}}(x_j^{(2)}) = \frac{\eta }{{{H_2}}}{\left( {\frac{{x_j^{(2)}}}{{{H_2}}}} \right)^{\eta - 1}} ,
    {\hat V^{(3)}}(x_j^{(3)}) = 0 (15)

    类似地,可将第r(r = 1、2、···、p)个时间子域的初始条件式(6)、界面条件式(7)和(8)及边界条件式(9)和(10)离散为

    式(6):T_{1,j}^{(r,s)} = T_0^{(r,s)}(x_j^{(s)}) , j = 1, 2, …, {N_s} ; s = 1, 2, 3 (16)
    式(7):T_{i,{N_s}}^{(r,s)} = T_{i,1}^{(r,s + 1)} , i = 2, 3, …, {M_r} ; s = 1, 2 (17)
    \begin{split} 式(8):& \sum\limits_{n = 1}^{{N_s}} {B_{{N_s},n}^{(1,s)}T_{i,n}^{(r,s)}} = \sum\limits_{n = 1}^{{N_{s + 1}}} {B_{1,n}^{(1,s + 1)}T_{i,n}^{(r,s + 1)}} ,\\ & i = 2, 3, …, {M_r} ; s = 1, 2 \end{split} (18)
    \begin{split} 式(9):&- \lambda _{i,1}^{(r,1)}\sum\limits_{n = 1}^{{N_1}} {B_{1,n}^{(1,1)}T_{i,n}^{(r,1)}} = {h_{{\text{bot}}}}\left[ {T_{i,1}^{(r,1)} - T_{{\text{bot}}}^{(r)}(t_i^{(r)})} \right] , \end{split}
    \qquad\quad i = 2, 3, …, {M_r} (19)
    \begin{split} 式(10):& \lambda _{i,{N_3}}^{(r,3)}\sum\limits_{n = 1}^{{N_3}} {B_{{N_3},n}^{(1,3)}T_{i,n}^{(r,3)}} = q_{{\text{top}}}^{(r)}(t_i^{(r)})-\\ & \sigma \varepsilon \left\{ {{{\left( {T_{i,{N_3}}^{(r,3)}} \right)}^4} - {{\left[ {T_{{\text{top}}}^{(r)}(t_i^{(r)})} \right]}^4}} \right\} , i = 2, 3, …, {M_r} \end{split} (20)

    图2可知,离散式(13)和式(16)~(20)建立在对应时间和空间子域的不同节点上(如式(13)中的待求节点温度位于圆点,而式(16)中的待求节点温度位于方点),因此直接求解这些离散式比较困难。基于矩阵理论,提出一种修正方法以使每个离散式都能建立在对应时间和空间子域的全部节点上。修正后的离散式(13)和式(16)~(20)可以写成如下矩阵方程的形式:

    式(13):{\boldsymbol{\tilde \lambda }}_{\text{d}}^{(r,s)} \circ {\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s)}} {{{{\bar {\boldsymbol{B}}}}}}_{:,2:{N_s} - 1}^{(1,s)}} \right)^{\circ 2}} +
    \qquad\quad\;\; {\boldsymbol{\hat \lambda }}_{\text{d}}^{(r,s)} \circ \left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s)}}{{\bar {\boldsymbol{B}}}}_{:,2:{N_s} - 1}^{(1,s)}} \right)+
    \qquad\quad\;\;{\boldsymbol{\lambda }}_{\text{d}}^{(r,s)} \circ \left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s)}}{{\bar {\boldsymbol{B}}}}_{:,2:{N_s} - 1}^{(2,s)}} \right) - {\boldsymbol{\rho }}_{\text{d}}^{(r,s)} \circ {\boldsymbol{c}}_{\text{d}}^{(r,s)} \circ
    \qquad\quad\;\; \begin{split} & \left( {{\boldsymbol{A}}_{2:{M_r},:}^{(1,r)}{{\boldsymbol{T}}^{(r,s)}}{\boldsymbol{\tilde I}}_{:,2:{N_s} - 1}^{(s)}} \right) = {{\boldsymbol{0}}}_{2:{M_r},2:{N_s} - 1}^{(r,s)} , \\ & s = 1, 2, 3 \end{split} (21)
    \begin{split} 式(16): {\boldsymbol{I}}_{1,:}^{(r)}{{\boldsymbol{T}}^{(r,s)}}{{\boldsymbol{\tilde I}}^{(s)}} = {{\boldsymbol{T}}^{({\text{ini}},r,s)}} ,s = 1, 2, 3 \end{split} (22)
    \begin{split} 式(17): {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s)}}{\boldsymbol{\tilde I}}_{:,{N_s}}^{(s)} = {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s + 1)}}{\boldsymbol{\tilde I}}_{:,1}^{(s + 1)} , s = 1, 2 \end{split} (23)
    \begin{split} 式(18):& {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s)}}{{\bar {\boldsymbol{B}}}}_{:,{N_s}}^{(1,s)} = {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s + 1)}}{{\bar {\boldsymbol{B}}}}_{:,1}^{(1,s + 1)} ,\\ & s = 1, 2 \end{split} (24)
    式(19): - {\boldsymbol{\lambda }}_{{\text{bot}}}^{(r,1)} \circ \left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,1)}}{{\bar {\boldsymbol{B}}}}_{:,1}^{(1,1)}} \right)=
    \qquad\qquad\quad {h_{{\text{bot}}}}\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,1)}}{\boldsymbol{\tilde I}}_{:,1}^{(1)} - {\boldsymbol{T}}_{2:{M_r}}^{({\text{bot}},r)}} \right) , (25)
    式(20): {\boldsymbol{\lambda }}_{{\text{top}}}^{(r,3)} \circ \left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,3)}}{{\bar {\boldsymbol{B}}}}_{:,{N_3}}^{(1,3)}} \right)=
    \qquad\quad\quad{\boldsymbol{q}}_{2:{M_r}}^{({\text{top}},r)} - \sigma \varepsilon \left[ {{{\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,3)}}{\boldsymbol{\tilde I}}_{:,{N_3}}^{(3)}} \right)}^{ \circ 4}} - {{\left( {{\boldsymbol{T}}_{2:{M_r}}^{({\text{top}},r)}} \right)}^{ \circ 4}}} \right] (26)

    其中: {{\boldsymbol{T}}^{(r,s)}} 为待求温度 T_{i,j}^{(r,s)} 形成的 {M_r} \times {N_s} 矩阵; {{\boldsymbol{T}}^{({\text{ini}},r,s)}} 为初始温度 T_0^{(r,s)}(x_j^{(s)}) 形成的 {N_s} \times 1 向量; {{\boldsymbol{T}}^{({\text{bot}},r)}} {{\boldsymbol{T}}^{({\text{top}},r)}} {{\boldsymbol{q}}^{({\text{top}},r)}} 分别为下部环境温度 T_{{\text{bot}}}^{(r)}(t_i^{(r)}) 、上部环境温度 T_{{\text{top}}}^{(r)}(t_i^{(r)}) 、上部热流 q_{{\text{top}}}^{(r)}(t_i^{(r)}) 形成的 {M_r} \times 1 向量; {{{\boldsymbol{0}}}^{(r,s)}} {M_r} \times {N_s} 零矩阵; {{\boldsymbol{I}}^{(r)}} {{\boldsymbol{\tilde I}}^{(s)}} 分别为 {M_r} 阶和 {N_s} 阶单位矩阵; {{\boldsymbol{A}}^{(1,r)}} {{\boldsymbol{B}}^{({k_x},s)}} ( {k_x} = 1,{\text{ 2}} )分别为权系数 A_{i,j}^{(1,r)} B_{i,j}^{({k_x},s)} 形成的 {M_r} 阶和 {N_s} 阶矩阵; {{{\bar {\boldsymbol{B}}}}^{({k_x},s)}} {{\boldsymbol{B}}^{({k_x},s)}} 的转置矩阵;其他表征热工参数的矩阵表达式见附录;矩阵或向量的下标表示索引范围,如 {\boldsymbol{I}}_{2:{M_r},:}^{(r)} 表示提取 {M_r} 阶矩阵 {{\boldsymbol{I}}^{(r)}} 的第二行到最后一行元素所形成的 ({M_r} - 1) \times {M_r} 矩阵;运算符\circ 表示Hadamard积,其定义为

    {\left( {{\boldsymbol{X}} \circ {\boldsymbol{Y}}} \right)_{i,j}} = {{\boldsymbol{X}}_{i,j}}{{\boldsymbol{Y}}_{i,j}} , {\left( {{{\boldsymbol{X}}^{ \circ k}}} \right)_{i,j}} = {\left( {{{\boldsymbol{X}}_{i,j}}} \right)^k} (27)

    其中,XY为同型矩阵。

    为求解联立矩阵方程式(21)~(26),可采用Kronecker积 \otimes [22]将其转化为如下代数方程:

    式(21): {\text{vec}}({\boldsymbol{\tilde \lambda }}_{\text{d}}^{(r,s)}) \circ {\left[ {{\boldsymbol{B}}_{2:{N_s} - 1,:}^{(1,s)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s)}})} \right]^{ \circ 2}}+
    \qquad\qquad\quad {\text{vec}}({\boldsymbol{\hat \lambda }}_{\text{d}}^{(r,s)}) \circ \left[ {{\boldsymbol{B}}_{2:{N_s} - 1,:}^{(1,s)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s)}})} \right]+
    \qquad\qquad\quad {\text{vec}}({\boldsymbol{\lambda }}_{\text{d}}^{(r,s)}) \circ \left[ {{\boldsymbol{B}}_{2:{N_s} - 1,:}^{(2,s)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s)}})} \right]-
    \qquad\qquad\quad \begin{split} & {\text{vec}}({\boldsymbol{\rho }}_{\text{d}}^{(r,s)}) \circ {\text{vec}}({\boldsymbol{c}}_{\text{d}}^{(r,s)}) \circ \left[ {{\boldsymbol{\tilde I}}_{2:{N_s} - 1,:}^{(s)} \otimes {\boldsymbol{A}}_{2:{M_r},:}^{(1,r)}} \right.\\ & \left.{{\text{vec}}({{\boldsymbol{T}}^{(r,s)}})} \right] = {\text{vec}}({{\boldsymbol{0}}}_{2:{M_r},2:{N_s} - 1}^{(r,s)}) ,\\ & s = 1, 2, 3 \end{split} (28)
    式(22): {{\boldsymbol{\tilde I}}^{(s)}} \otimes {\boldsymbol{I}}_{1,:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s)}}) = {\text{vec}}({{\boldsymbol{T}}^{({\text{ini}},r,s)}}) ,
    \qquad\qquad\quad s = 1, 2, 3 (29)
    式(23): {\boldsymbol{\tilde I}}_{{N_s},:}^{(s)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s)}})=
    \qquad\qquad\quad{\boldsymbol{\tilde I}}_{1,:}^{(s + 1)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s + 1)}}) , s = 1, 2 (30)
    式(24): {\boldsymbol{B}}_{{N_s},:}^{(1,s)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s)}})=
    \qquad\qquad\quad{\boldsymbol{B}}_{1,:}^{(1,s + 1)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s + 1)}}) , s = 1, 2 (31)
    式(25): - {\text{vec}}({\boldsymbol{\lambda }}_{{\text{bot}}}^{(r,1)}) \circ \left[ {{\boldsymbol{B}}_{1,:}^{(1,1)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,1)}})} \right]=
    \qquad\qquad\quad{h_{{\text{bot}}}}\left[ {{\boldsymbol{\tilde I}}_{1,:}^{(1)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,1)}}) - {\boldsymbol{T}}_{2:{M_r}}^{({\text{bot}},r)}} \right] (32)
    式(26): {\text{vec}}({\boldsymbol{\lambda }}_{{\text{top}}}^{(r,3)}) \circ \left[ {{\boldsymbol{B}}_{{N_3},:}^{(1,3)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,3)}})} \right] =
    \qquad\qquad\quad \begin{split} &{\boldsymbol{q}}_{2:{M_r}}^{({\text{top}},r)} - \sigma \varepsilon \left\{ {{\left[ {{\boldsymbol{\tilde I}}_{{N_3},:}^{(3)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,3)}})} \right]}^{ \circ 4}} -\right. \\ &\left.{{{\left( {{\boldsymbol{T}}_{2:{M_r}}^{({\text{top}},r)}} \right)}^{ \circ 4}}} \right\} \end{split} (33)

    其中,vec()为向量化算子;热工参数矩阵的向量化表达式见附录。

    由式(28)~(33)可得如下非线性代数方程组:

    \left[ {{\boldsymbol{E}}_{\text{L}}^{(r)} + {\boldsymbol{E}}_{{\text{NL}}}^{(r)}({\boldsymbol{T}}_{{\text{sum}}}^{(r)})} \right]{\boldsymbol{T}}_{{\text{sum}}}^{(r)} = {{\boldsymbol{f}}^{(r)}} (34)

    其中:向量 {{\boldsymbol{f}}^{(r)}} 包含了初始条件、界面条件及边界条件; {\boldsymbol{E}}_{\text{L}}^{(r)} 为线性刚度矩阵; {\boldsymbol{E}}_{{\text{NL}}}^{(r)} 为非线性刚度矩阵,由如下待求温度向量决定:

    {\boldsymbol{T}}_{{\text{sum}}}^{(r)} = {\left\{ {\begin{array}{*{20}{c}} {{\text{vec}}({{\boldsymbol{T}}^{(r,1)}})}&{{\text{vec}}({{\boldsymbol{T}}^{(r,2)}})}&{{\text{vec}}({{\boldsymbol{T}}^{(r,3)}})} \end{array}} \right\}^{\text{T}}} (35)

    采用Newton-Raphson迭代法求解式(34),即可得第r个时间子域内FGM夹层板的温度分布,该时间子域最终时刻的温度分布可作为下一个子域的初始温度条件。从第一个时间子域递推到最后一个子域,即可得到整个受热过程的温度分布。

    本文仅分析了FGM夹层板的一维非线性瞬态传热问题,但本方法仍可用于二、三维问题的计算。以二维问题为例,仅需在板长方向上布置节点,则二维瞬态传热控制方程仍然可以被离散处理,最后得到的非线性代数方程组也可以采用Newton-Raphson迭代法进行求解。这些内容将在后续工作中开展。

    考虑两种类型的金属-陶瓷复合FGM夹层板,第一类为不锈钢和氮化硅复合(记为SUS304/Si3N4),第二类为钛合金和氧化锆复合(记为Ti-6Al-4V/ZrO2)。表1给出了组分材料热工参数的温度系数[23-26]。如无特别说明,假设初始时刻夹层板处于均匀的温度场T0 = 300 K。

    表  1  组分材料热工参数的温度系数[23-26]
    Table  1.  Temperature coefficients of thermophysical properties for the component materials[23-26]
    Material Property P0 P−1 P1 P2 P3
    SUS304 λ/(W∙(m∙K)−1) 15.379 0 −1.264 × 10−3 2.092×10−6 −7.223×10−10
    c/(J∙(kg∙K)−1) 496.56 0 −1.151×10−3 1.636×10−6 −5.863×10−10
    ρ/(kg∙m−3) 8166.0 0 0 0 0
    Si3N4 λ/(W∙(m∙K)−1) 13.723 0 −1.032× 10−3 5.466 × 10−7 −7.876×10−11
    c/(J∙(kg∙K)−1) 555.11 0 1.016×10−3 2.920×10−7 −1.670×10−10
    ρ/(kg∙m−3) 2370.0 0 0 0 0
    Ti-6Al-4V λ/(W∙(m∙K)−1) 1.0000 0 1.704×10−2 0 0
    c/(J∙(kg∙K)−1) 625.30 0 −4.224×10−4 7.179×10−7 0
    ρ/(kg∙m−3) 4420.0 0 0 0 0
    ZrO2 λ/(W∙(m∙K)−1) 1.7000 0 1.276×10−4 6.648×10−8 0
    c/(J∙(kg∙K)−1) 487.34 0 3.049×10−4 −6.037×10−8 0
    ρ/(kg∙m−3) 5700.0 0 0 0 0
    Notes: λ, c, ρ—Thermal conductivity, specific heat, and density, respectively; Pk (k = −1, 0, 1, 2, 3)—Temperature coefficients; SUS304, Si3N4, Ti-6Al-4V, and ZrO2 denote stainless steel, silicon nitride, titanium alloy, and zirconia, respectively.
    下载: 导出CSV 
    | 显示表格

    由于本文基于IDQEM,因此需要研究数值结果的收敛性。考虑一块体积分数指数为η = 0.2、总厚度为H = 1 mm的1-2-1型SUS304/Si3N4夹层板,其中1-2-1表示底层、芯层、顶层的厚度之比,即芯层厚度为底层和顶层厚度的2倍。在后续算例中,仍然沿用这种格式表示各层厚度之比。夹层板上表面受热冲击作用,其上下表面温度分别为Tc = 400 K和Tm = 300 K。为应用IDQEM,假定底层、芯层、顶层对应的空间子域的节点数分别为N1 = 0.25NN2 = 0.5NN3 = 0.25N。由于热边界条件不随时间变化,因此无需进行时间域的划分(即p = 1)。表2给出了SUS304/Si3N4夹层板温度解的收敛结果及计算效率。可以看出,随着节点数MN的增加,数值结果快速收敛,但计算效率随之降低。当M × N = 20 × 28时,以5位有效数字表示的结果已不再变化,且此时计算效率也较高。此外,当仅采用12 × 20个节点时,精度已经令人满意。因此在后续算例中,如无特别说明,节点数取为20 × 28。

    表  2  体积分数指数η = 0.2时SUS304/Si3N4夹层板x = 0.75H处的温度结果(单位:K)
    Table  2.  Temperature results at x = 0.75H for the SUS304/Si3N4 sandwich slab with volume fraction index η = 0.2 (Unit: K)
    Time Temporal grid point number Spatial grid point number
    N = 12 N = 20 N = 28 N = 36
    t = 0.01 s M = 4 342.16 (0.0312 s) 342.96 (0.0625 s) 342.98 (0.0938 s) 342.99 (0.1154 s)
    M = 12 342.16 (0.1406 s) 343.01 (0.1719 s) 343.00 (0.2193 s) 343.01 (0.4688 s)
    M = 20 342.16 (0.1875 s) 343.01 (0.3281 s) 343.01 (0.5806 s) 343.01 (1.0312 s)
    M = 28 342.16 (0.2031 s) 343.01 (0.5156 s) 343.01 (1.1562 s) 343.01 (2.1094 s)
    t = 0.03 s M = 4 362.95 (0.0469 s) 362.80 (0.0712 s) 362.81 (0.1024 s) 362.80 (0.1193 s)
    M = 12 362.87 (0.1562 s) 362.69 (0.1736 s) 362.68 (0.2056 s) 362.68 (0.5021 s)
    M = 20 362.87 (0.1719 s) 362.70 (0.3598 s) 362.69 (0.6006 s) 362.69 (1.0156 s)
    M = 28 362.87 (0.2031 s) 362.70 (0.5469 s) 362.69 (1.0625 s) 362.69 (2.2188 s)
    Notes: M and N—Total numbers of grid points in the temporal and spatial domains, respectively; The content in parentheses represents the central processing unit (CPU) time.
    下载: 导出CSV 
    | 显示表格

    Pandey和Pradyumna[8-9]采用FDM开展了上述SUS304/Si3N4夹层板的非线性瞬态传热分析,其结果可用来验证本文IDQEM解的正确性。图3给出了η = 0.2和2两种情况下1-2-1型夹层板x = 0.5H处温度解的对比结果。图4给出考虑材料热工参数温度依赖性(TD)与不考虑温度依赖性(TI)两种状态下1-8-1型夹层板(η = 2)x = 0.5H处温度解的对比结果,其中TI状态下的热工参数取T = 0 K时的值。可以看出,IDQEM的预测结果与FDM结果具有很好的一致性。

    图  3  η = 0.2和2两种情况下SUS304/Si3N4夹层板x = 0.5H处温度变化
    Figure  3.  Temperature variations at x = 0.5H for the SUS304/Si3N4 sandwich slab with η = 0.2 and 2
    T—Temperature
    图  4  η = 2时SUS304/Si3N4夹层板x = 0.5H处温度变化
    Figure  4.  Temperature variations at x = 0.5H for the SUS304/Si3N4 sandwich slab with η = 2
    TI—Temperature independent; TD—Temperature dependent

    若将夹层板芯层材料由FGM换成均质材料,并忽略各组分材料热工参数的温度依赖性,则本文关于FGM夹层板的非线性问题将退化为普通夹层板的线性问题。对于此类问题,可采用一些解析方法求解,其结果可用于验证本文IDQEM解的正确性。Zhang等[27]采用Laplace变换法研究了对流条件下夹层板的瞬态温度解析解,其中材料热工参数、尺寸和热边界条件为λ(1) = λ(3) = 54.33 W∙(m∙K)−1λ(2) = 1.33 W∙(m∙K)−1c(1) = c(3) = 440 J∙(kg∙K)−1c(2) = 900 J∙(kg∙K)−1ρ(1) = ρ(3) =7850 kg∙m−3ρ(2) = 2300 kg∙m−3H1 = H3 = 20 mm,H2 = 160 mm,hbot = htop = 100 W∙m−2∙K−1T0 = Tbot = 293 K,Ttop = 573 K。图5给出温度分布对比结果。可以看出,本文数值解与解析解具有很好的一致性。

    图  5  不同时刻SUS304/Si3N4夹层板的温度分布
    Figure  5.  Temperature profiles at different time for the SUS304/Si3N4 sandwich slab

    Malekzadeh等[18]和Dai等[19]采用DQEM分析了含FGM层构件的传热问题,但都仅考虑了热传导系数的温度依赖性,而假定其他热工参数为常量。为分析这种假定对数值结果精度和计算效率的影响,采用本方法讨论两种情况下的结果:(1) 仅考虑热传导系数的温度依赖性,其他热工参数取T0 = 300 K时的值(Case 1);(2) 考虑全部热工参数的温度依赖性(Case 2)。仍以3.1节的夹层板为例,表3给出了两种情况下数值结果与CPU计算时间的对比。可以看出,Case 1的结果总是比Case 2的大,最大误差达到1.76%,但两种情况下计算效率并无显著差异。

    表  3  两种情况下数值结果和CPU计算时间的对比
    Table  3.  Comparison of numerical results and CPU time in the two cases
    Position Time/s Case 1/K Case 2/K Error CPU time of Case 1/s CPU time of Case 2/s
    x = 0.25H 0.01 302.43 301.40 0.34% 0.6094 0.6310
    0.03 314.28 312.23 0.66% 0.5938 0.5781
    x = 0.75H 0.01 349.05 343.01 1.76% 0.6094 0.6310
    0.03 365.93 362.69 0.89% 0.5938 0.5781
    Notes: In Case 1, only the thermal conductivity is considered to be TD; In Case 2, all the thermophysical properties are considered to be TD.
    下载: 导出CSV 
    | 显示表格

    考虑上述SUS304/Si3N4夹层板受3种不同上表面温度Tc = 400 K、500 K、600 K作用,而下表面温度保持为Tm = 300 K。为研究材料热工参数温度依赖性对瞬态传热的影响,图6给出了t = 0.01 s和0.02 s两种时刻下的温度分布。作为对比,同时绘制了TD和TI两种状态下结果,其中TI状态下的热工参数取T0 = 300 K时的值。可以看出,随着Tc的上升,TD和TI两种状态下的温度分布差异越来越大,这是由于两种状态下材料热工参数的差异随温度的上升而增大。同时可以发现,同一时刻同一表面温度作用下,TI状态下的内部温度总是比TD状态下的高。为解释这种现象,图7绘制了Tc = 600 K时夹层板的热扩散率α(α = λ/ρc)的时空分布。可以发现,对于陶瓷层和FGM层,TI状态下的α比TD下的大得多,导致TI状态具有更快的传热速度,进而具有更高的内部温度。因此在SUS304/Si3N4夹层板的设计中,忽略材料热工参数的温度依赖性将高估板内温度分布,从而导致设计过于保守。

    图  6  不同上表面温度Tc作用下SUS304/Si3N4夹层板的温度分布
    Figure  6.  Temperature profiles of the SUS304/Si3N4 sandwich slab under different top surface temperatures Tc
    图  7  Tc = 600 K时SUS304/Si3N4夹层板的热扩散率α时空分布
    Figure  7.  Thermal diffusivity α contour of the SUS304/Si3N4 sandwich slab under the surface temperature Tc = 600 K

    为进一步说明TD和TI两种状态下的温度差异,图8给出了t = 0.01 s、0.02 s、0.03 s这3种时刻下x = 0.5H处的温度T与上表面温度Tc上的关系。可以看出,TI状态下T随着Tc的上升线性增大,而TD状态下T随着Tc的上升非线性增大。当Tc达到600 K时,3种时刻下的温度差异分别达到了12.00%、12.21%、10.41%。因此在预测FGM夹层板的温度分布时,需要考虑材料热工参数的温度依赖性,特别是在高温工况下。在接下来的算例中,仅考虑TD状态下的温度结果。

    图  8  SUS304/Si3N4夹层板x = 0.5H处的温度与Tc的关系
    Figure  8.  Temperature at x = 0.5H versus Tc for the SUS304/Si3N4 sandwich slab

    考虑一块总厚度为H = 120 mm的1-2-1型Ti-6Al-4V/ZrO2夹层板。板的下表面受自然对流作用,上表面受气动加热和热辐射共同作用,其中hbot = 2 W∙m−2∙K−1Tbot = 300 K,qtop = 80 kW∙m−2Ttop = 4 K,σ = 5.67 × 10−8 W∙m−2∙K−4ε = 0.8。图9给出3种不同体积分数指数η = 0.2、1、5情况下夹层板下表面(x = 0)和上表面(x = H)的温度变化。可以看出,当t ≤ 250 s时,温度几乎不受η影响;随着加热时间的增长,η的影响越发显著。这是由于在早期阶段,传热主要发生在陶瓷层,因此传热不受η影响;随着时间的增长,整个夹层板都被加热,因此η对温度分布的影响越发显著。当t = 2000 s时,随着η从0.2增加到5,下表面温度从350 K增加到414 K,而上表面温度从1036 K减少到1021 K。为解释这种现象,图10绘制了不同体积分数指数η时夹层板的热扩散率α的时空分布。可以发现,随着η的增加,FGM层中陶瓷含量将会增加,导致FGM层整体具有更高的热扩散率,进而加快热传导。因此,体积分数指数对FGM夹层板的瞬态传热有着重要影响。

    图  9  η = 0.2、1、5这3种情况下Ti-6Al-4V/ZrO2夹层板的温度变化
    Figure  9.  Temperature variations of the Ti-6Al-4V/ZrO2 sandwich slab for η=0.2, 1, 5
    图  10  Ti-6Al-4V/ZrO2夹层板的热扩散率时空分布
    Figure  10.  Thermal diffusivity contour of the Ti-6Al-4V/ZrO2 sandwich slab

    考虑上述1-2-1型Ti-6Al-4V/ZrO2夹层板(η = 2)受时变热流作用[28],其余热边界条件与4.2节相同。图11给出了时变热流曲线,该曲线可用三次函数qtop = a3t3 + a2t2 + a1t + a0分段拟合,表4给出了分段区间及拟合系数。可以看出,整个加热过程可以划分为4个时间子域,取各子域节点数为M1 = 5、M2 = 4、M3 = 11、M4 = 8。

    图  11  Ti-6Al-4V/ZrO2夹层板上表面的热流
    Figure  11.  Heat flux on the top surface of the Ti-6Al-4V/ZrO2 sandwich slab
    qtop—Heat flux generated on the top surface
    表  4  时变热流的多项式拟合
    Table  4.  Polynomial fits for the time-dependent heat flux
    Time/s a0 a1 a2/10−3 a3/10−6
    0-431.1 0 0.003196 1.021 −1.383
    431.1-660.2 −228.7 1.595 −2.670 1.471
    660.2-1561 253.4 0.5961 0.6477 0.2043
    1561-2200 1257 2.307 −1.212 0.1929
    Note: ak (k = 0, 1, 2, 3)—Fitting coefficients.
    下载: 导出CSV 
    | 显示表格

    图12对比了时变热流(图11)和恒定热流qtop = 80 kW∙m−2作用下夹层板的温度时空分布。可以看出,靠近上部陶瓷层等温线更密集,表示温度梯度更大。图12(a)表明,当夹层板受恒定热流作用时,温度随着tx的增加单调上升。然而从图12(b)可以看出,由于外界热流先增后减,因此等温线在约1430 s后出现“U”型转折。此外,恒定热流情况下,整个加热过程中上下表面的最大温度分别约为1025 K和436 K;在时变热流情况下,整个加热过程中上下表面的最大温度分别约为1149 K和425 K。上下表面的巨大温差表明了FGM夹层板的优异隔热效果。

    图  12  Ti-6Al-4V/ZrO2夹层板的温度时空分布
    Figure  12.  Temperature contour of the Ti-6Al-4V/ZrO2 sandwich slab

    (1) 增量微分求积单元法具有快速收敛性,与已有文献的解析和数值结果的对比验证了本方法的正确性。

    (2) 对于SUS304/Si3N4夹层板,不考虑材料热工参数温度依赖性状态较考虑温度依赖性状态具有更快的热传导,且两种状态下的温度差异随着表面温度的上升而增加。当上表面温度为600 K时,时间t = 0.01 s、0.02 s、0.03 s的中部温度差异分别达到了12.00%、12.21%、10.41%。因此,材料热工参数的温度依赖性是准确预测板内温度分布的关键因素,尤其是在高温分析中。忽略材料热工参数的温度依赖性将高估板内温度分布,从而导致设计过于保守。

    (3) 在早期阶段,功能梯度材料(FGM)芯层的体积分数指数η对夹层板的温度分布几乎无影响,但随着加热时间的增长,η对温度分布的影响越发显著。当t =2000 s时,随着η从0.2增加到5,下表面温度从350 K增加到414 K,而上表面温度从1036 K减少到1021 K。

    (4) 当FGM夹层板受恒定热流作用时,其内部温度随着时间和位置的变化单调变化。然而当热流随时间先增后减时,等温线出现“U”型转折。

    附录:

    式(21)、(25)和(26)中,表征热工参数的矩阵的表达式为

    {\boldsymbol{P}}_{\text{d}}^{(r,s)} = \sum\limits_{k = - 1}^3 {\bar P_k^{\text{m}}{{\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s)}}{\boldsymbol{\tilde I}}_{:,2:{N_s} - 1}^{(s)}} \right)}^{ \circ k}}} + {\boldsymbol{V}}_{2:{M_r},2:{N_s} - 1}^{(s)} \circ \left[ {\sum\limits_{k = - 1}^3 {\left( {\bar P_k^{\text{c}} - \bar P_k^{\text{m}}} \right){{\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s)}}{\boldsymbol{\tilde I}}_{:,2:{N_s} - 1}^{(s)}} \right)}^{ \circ k}}} } \right],{\text{ }}P = \lambda ,{\text{ }}\rho {\text{, }}c,
    {\boldsymbol{\tilde \lambda }}_{\text{d}}^{(r,s)} = \sum\limits_{k = - 1}^3 {k\bar \lambda _k^{\text{m}}{{\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s)}}{\boldsymbol{\tilde I}}_{:,2:{N_s} - 1}^{(s)}} \right)}^{ \circ (k - 1)}}} + \left[ {\sum\limits_{k = - 1}^3 {k\left( {\bar \lambda _k^{\text{c}} - \bar \lambda _k^{\text{m}}} \right){{\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s)}}{\boldsymbol{\tilde I}}_{:,2:{N_s} - 1}^{(s)}} \right)}^{ \circ (k - 1)}}} } \right] \circ {\boldsymbol{V}}_{2:{M_r},2:{N_s} - 1}^{(s)} ,
    {\boldsymbol{\hat \lambda }}_{\text{d}}^{(r,s)} = \left[ {\sum\limits_{k = - 1}^3 {\left( {\bar \lambda _k^{\text{c}} - \bar \lambda _k^{\text{m}}} \right){{\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,s)}}{\boldsymbol{\tilde I}}_{:,2:{N_s} - 1}^{(s)}} \right)}^{ \circ k}}} } \right] \circ {\boldsymbol{\hat V}}_{2:{M_r},2:{N_s} - 1}^{(s)} ,
    {\boldsymbol{\lambda }}_{{\text{bot}}}^{(r,1)} = \sum\limits_{k = - 1}^3 {\bar \lambda _k^{\text{m}}{{\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,1)}}{\boldsymbol{\tilde I}}_{:,1}^{(1)}} \right)}^{ \circ k}}} + \left[ {\sum\limits_{k = - 1}^3 {\left( {\bar \lambda _k^{\text{c}} - \bar \lambda _k^{\text{m}}} \right){{\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,1)}}{\boldsymbol{\tilde I}}_{:,1}^{(1)}} \right)}^{ \circ k}}} } \right] \circ {\boldsymbol{V}}_{2:M,1}^{(1)} ,
    {\boldsymbol{\lambda }}_{{\text{top}}}^{(r,3)} = \sum\limits_{k = - 1}^3 {\bar \lambda _k^{\text{m}}{{\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,3)}}{\boldsymbol{\tilde I}}_{:,{N_3}}^{(3)}} \right)}^{ \circ k}}} + \left[ {\sum\limits_{k = - 1}^3 {\left( {\bar \lambda _k^{\text{c}} - \bar \lambda _k^{\text{m}}} \right){{\left( {{\boldsymbol{I}}_{2:{M_r},:}^{(r)}{{\boldsymbol{T}}^{(r,3)}}{\boldsymbol{\tilde I}}_{:,{N_3}}^{(3)}} \right)}^{ \circ k}}} } \right] \circ {\boldsymbol{V}}_{2:M,{N_3}}^{(3)}

    式(28)、(32)和(33)中,热工参数矩阵的向量化表达式为

    \begin{split} {\text{vec}}({\boldsymbol{P}}_{\text{d}}^{(r,s)}) =& \sum\limits_{k = - 1}^3 {\bar P_k^{\text{m}}{{\left[ {{\boldsymbol{\tilde I}}_{2:{N_s} - 1,:}^{(s)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s)}})} \right]}^{ \circ k}}} + \left\{ {\sum\limits_{k = - 1}^3 {\left( {\bar P_k^{\text{c}} - \bar P_k^{\text{m}}} \right){{\left[ {{\boldsymbol{\tilde I}}_{2:{N_s} - 1,:}^{(s)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s)}})} \right]}^{ \circ k}}} } \right\}\circ \\& {\text{vec}}({\boldsymbol{V}}_{2:{M_r},2:{N_s} - 1}^{(s)}) , P = \lambda ,{\text{ }}\rho {\text{, }}c , \end{split}
    \begin{split} {\text{vec}}({\boldsymbol{\tilde \lambda }}_{\text{d}}^{(r,s)}) =& \sum\limits_{k = - 1}^3 {k\bar \lambda _k^{\text{m}}{{\left[ {{\boldsymbol{\tilde I}}_{2:{N_s} - 1,:}^{(s)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s)}})} \right]}^{ \circ (k - 1)}}} + \left\{ {\sum\limits_{k = - 1}^3 {k\left( {\bar \lambda _k^{\text{c}} - \bar \lambda _k^{\text{m}}} \right){{\left[ {{\boldsymbol{\tilde I}}_{2:{N_s} - 1,:}^{(s)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s)}})} \right]}^{ \circ (k - 1)}}} } \right\} \circ\\ & {\text{vec}}({\boldsymbol{V}}_{2:{M_r},2:{N_s} - 1}^{(s)}) , \end{split}
    {\text{vec}}({\boldsymbol{\hat \lambda }}_{\text{d}}^{(r,s)}) = \left\{ {\sum\limits_{k = - 1}^3 {\left( {\bar \lambda _k^{\text{c}} - \bar \lambda _k^{\text{m}}} \right){{\left[ {{\boldsymbol{\tilde I}}_{2:{N_s} - 1,:}^{(s)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,s)}})} \right]}^{ \circ k}}} } \right\} \circ {\text{vec}}({\boldsymbol{\hat V}}_{2:{M_r},2:{N_s} - 1}^{(s)}) ,
    {\text{vec}}({\boldsymbol{\lambda }}_{{\text{bot}}}^{(r,1)}) = \sum\limits_{k = - 1}^3 {\bar \lambda _k^{\text{m}}{{\left[ {{\boldsymbol{\tilde I}}_{1,:}^{(1)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,1)}})} \right]}^{ \circ k}}} + \left\{ {\sum\limits_{k = - 1}^3 {\left( {\bar \lambda _k^{\text{c}} - \bar \lambda _k^{\text{m}}} \right){{\left[ {{\boldsymbol{\tilde I}}_{1,:}^{(1)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,1)}})} \right]}^{ \circ k}}} } \right\} \circ {\text{vec}}({\boldsymbol{V}}_{2:M,1}^{(1)}) ,
    {\text{vec}}({\boldsymbol{\lambda }}_{{\text{top}}}^{(r,3)}) = \sum\limits_{k = - 1}^3 {\bar \lambda _k^{\text{m}}{{\left[ {{\boldsymbol{\tilde I}}_{{N_3},:}^{(3)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,3)}})} \right]}^{ \circ k}}} + \left\{ {\sum\limits_{k = - 1}^3 {\left( {\bar \lambda _k^{\text{c}} - \bar \lambda _k^{\text{m}}} \right){{\left[ {{\boldsymbol{\tilde I}}_{{N_3},:}^{(3)} \otimes {\boldsymbol{I}}_{2:{M_r},:}^{(r)}{\text{vec}}({{\boldsymbol{T}}^{(r,3)}})} \right]}^{ \circ k}}} } \right\} \circ {\text{vec}}({\boldsymbol{V}}_{2:M,{N_3}}^{(3)})
  • 图  1   功能梯度材料(FGM)夹层板示意图

    Figure  1.   Schematic view of the functionally graded material (FGM) sandwich slab

    x—Global coordinate for the sandwich slab; x(s) (s = 1, 2, 3)—Local coordinate for the sth layer; H—Total thickness of the sandwich slab;HsThickness of the sth layer

    图  2   r个时间子域下的节点分布

    Figure  2.   Grid point distribution in the rth temporal sub-domain

    Mr and NsGrid point numbers in the rth temporal sub-domain and the sth spatial sub-domain, respectively; t—Time variable

    图  3   η = 0.2和2两种情况下SUS304/Si3N4夹层板x = 0.5H处温度变化

    Figure  3.   Temperature variations at x = 0.5H for the SUS304/Si3N4 sandwich slab with η = 0.2 and 2

    T—Temperature

    图  4   η = 2时SUS304/Si3N4夹层板x = 0.5H处温度变化

    Figure  4.   Temperature variations at x = 0.5H for the SUS304/Si3N4 sandwich slab with η = 2

    TI—Temperature independent; TD—Temperature dependent

    图  5   不同时刻SUS304/Si3N4夹层板的温度分布

    Figure  5.   Temperature profiles at different time for the SUS304/Si3N4 sandwich slab

    图  6   不同上表面温度Tc作用下SUS304/Si3N4夹层板的温度分布

    Figure  6.   Temperature profiles of the SUS304/Si3N4 sandwich slab under different top surface temperatures Tc

    图  7   Tc = 600 K时SUS304/Si3N4夹层板的热扩散率α时空分布

    Figure  7.   Thermal diffusivity α contour of the SUS304/Si3N4 sandwich slab under the surface temperature Tc = 600 K

    图  8   SUS304/Si3N4夹层板x = 0.5H处的温度与Tc的关系

    Figure  8.   Temperature at x = 0.5H versus Tc for the SUS304/Si3N4 sandwich slab

    图  9   η = 0.2、1、5这3种情况下Ti-6Al-4V/ZrO2夹层板的温度变化

    Figure  9.   Temperature variations of the Ti-6Al-4V/ZrO2 sandwich slab for η=0.2, 1, 5

    图  10   Ti-6Al-4V/ZrO2夹层板的热扩散率时空分布

    Figure  10.   Thermal diffusivity contour of the Ti-6Al-4V/ZrO2 sandwich slab

    图  11   Ti-6Al-4V/ZrO2夹层板上表面的热流

    Figure  11.   Heat flux on the top surface of the Ti-6Al-4V/ZrO2 sandwich slab

    qtop—Heat flux generated on the top surface

    图  12   Ti-6Al-4V/ZrO2夹层板的温度时空分布

    Figure  12.   Temperature contour of the Ti-6Al-4V/ZrO2 sandwich slab

    表  1   组分材料热工参数的温度系数[23-26]

    Table  1   Temperature coefficients of thermophysical properties for the component materials[23-26]

    Material Property P0 P−1 P1 P2 P3
    SUS304 λ/(W∙(m∙K)−1) 15.379 0 −1.264 × 10−3 2.092×10−6 −7.223×10−10
    c/(J∙(kg∙K)−1) 496.56 0 −1.151×10−3 1.636×10−6 −5.863×10−10
    ρ/(kg∙m−3) 8166.0 0 0 0 0
    Si3N4 λ/(W∙(m∙K)−1) 13.723 0 −1.032× 10−3 5.466 × 10−7 −7.876×10−11
    c/(J∙(kg∙K)−1) 555.11 0 1.016×10−3 2.920×10−7 −1.670×10−10
    ρ/(kg∙m−3) 2370.0 0 0 0 0
    Ti-6Al-4V λ/(W∙(m∙K)−1) 1.0000 0 1.704×10−2 0 0
    c/(J∙(kg∙K)−1) 625.30 0 −4.224×10−4 7.179×10−7 0
    ρ/(kg∙m−3) 4420.0 0 0 0 0
    ZrO2 λ/(W∙(m∙K)−1) 1.7000 0 1.276×10−4 6.648×10−8 0
    c/(J∙(kg∙K)−1) 487.34 0 3.049×10−4 −6.037×10−8 0
    ρ/(kg∙m−3) 5700.0 0 0 0 0
    Notes: λ, c, ρ—Thermal conductivity, specific heat, and density, respectively; Pk (k = −1, 0, 1, 2, 3)—Temperature coefficients; SUS304, Si3N4, Ti-6Al-4V, and ZrO2 denote stainless steel, silicon nitride, titanium alloy, and zirconia, respectively.
    下载: 导出CSV

    表  2   体积分数指数η = 0.2时SUS304/Si3N4夹层板x = 0.75H处的温度结果(单位:K)

    Table  2   Temperature results at x = 0.75H for the SUS304/Si3N4 sandwich slab with volume fraction index η = 0.2 (Unit: K)

    Time Temporal grid point number Spatial grid point number
    N = 12 N = 20 N = 28 N = 36
    t = 0.01 s M = 4 342.16 (0.0312 s) 342.96 (0.0625 s) 342.98 (0.0938 s) 342.99 (0.1154 s)
    M = 12 342.16 (0.1406 s) 343.01 (0.1719 s) 343.00 (0.2193 s) 343.01 (0.4688 s)
    M = 20 342.16 (0.1875 s) 343.01 (0.3281 s) 343.01 (0.5806 s) 343.01 (1.0312 s)
    M = 28 342.16 (0.2031 s) 343.01 (0.5156 s) 343.01 (1.1562 s) 343.01 (2.1094 s)
    t = 0.03 s M = 4 362.95 (0.0469 s) 362.80 (0.0712 s) 362.81 (0.1024 s) 362.80 (0.1193 s)
    M = 12 362.87 (0.1562 s) 362.69 (0.1736 s) 362.68 (0.2056 s) 362.68 (0.5021 s)
    M = 20 362.87 (0.1719 s) 362.70 (0.3598 s) 362.69 (0.6006 s) 362.69 (1.0156 s)
    M = 28 362.87 (0.2031 s) 362.70 (0.5469 s) 362.69 (1.0625 s) 362.69 (2.2188 s)
    Notes: M and N—Total numbers of grid points in the temporal and spatial domains, respectively; The content in parentheses represents the central processing unit (CPU) time.
    下载: 导出CSV

    表  3   两种情况下数值结果和CPU计算时间的对比

    Table  3   Comparison of numerical results and CPU time in the two cases

    Position Time/s Case 1/K Case 2/K Error CPU time of Case 1/s CPU time of Case 2/s
    x = 0.25H 0.01 302.43 301.40 0.34% 0.6094 0.6310
    0.03 314.28 312.23 0.66% 0.5938 0.5781
    x = 0.75H 0.01 349.05 343.01 1.76% 0.6094 0.6310
    0.03 365.93 362.69 0.89% 0.5938 0.5781
    Notes: In Case 1, only the thermal conductivity is considered to be TD; In Case 2, all the thermophysical properties are considered to be TD.
    下载: 导出CSV

    表  4   时变热流的多项式拟合

    Table  4   Polynomial fits for the time-dependent heat flux

    Time/s a0 a1 a2/10−3 a3/10−6
    0-431.1 0 0.003196 1.021 −1.383
    431.1-660.2 −228.7 1.595 −2.670 1.471
    660.2-1561 253.4 0.5961 0.6477 0.2043
    1561-2200 1257 2.307 −1.212 0.1929
    Note: ak (k = 0, 1, 2, 3)—Fitting coefficients.
    下载: 导出CSV
  • [1] 刘浩洋, 吕超雨, 石姗姗, 等. 芳纶纤维增韧碳纤维增强环氧树脂复合材料-铝蜂窝夹芯结构界面性能和增韧机制[J]. 复合材料学报, 2022, 39(2): 559-567.

    LIU Haoyang, LYU Chaoyu, SHI Shanshan, et al. Interfacial toughening and toughening mechanism of aramid staple fiber to carbon fiber reinforced epoxy resin composite-aluminum honeycomb sandwich structure[J]. Acta Materiae Compositae Sinica, 2022, 39(2): 559-567 (in Chinese).

    [2] 李华冠, 丁颖, 章月, 等. 玻璃纤维立体织物增强环氧树脂泡沫夹层复合材料的制备及力学性能[J]. 复合材料学报, 2023, 40(1): 601-612.

    LI Huaguan, DING Ying, ZHANG Yue, et al. Preparation and mechanical properties of glass fiber reinforced 3D fabric reinforced epoxy foam sandwich composites[J]. Acta Materiae Compositae Sinica, 2023, 40(1): 601-612 (in Chinese).

    [3]

    XIE H, LI W, FANG H, et al. Flexural behavior evaluation of a foam core curved sandwich beam[J]. Composite Structures, 2024, 328: 117729. DOI: 10.1016/j.compstruct.2023.117729

    [4]

    WANG X, WANG Z, ZENG T, et al. Exact analytical solution for steady-state heat transfer in functionally graded sandwich slabs with convective-radiative boundary conditions[J]. Composite Structures, 2018, 192: 379-386. DOI: 10.1016/j.compstruct.2018.03.006

    [5]

    DAIKH A A, MEGUENI A. Thermal buckling analysis of functionally graded sandwich plates[J]. Journal of Thermal Stresses, 2018, 41(2): 139-159. DOI: 10.1080/01495739.2017.1393644

    [6]

    DINI A, NEMATOLLAHI M A, HOSSEINI M. Analytical solution for magneto-thermo-elastic responses of an annular functionally graded sandwich disk by considering internal heat generation and convective boundary condition[J]. Journal of Sandwich Structures and Materials, 2021, 23(2): 542-567.

    [7]

    SHEN H S, LI S R. Postbuckling of sandwich plates with FGM face sheets and temperature-dependent properties[J]. Composites Part B: Engineering, 2008, 39(2): 332-344. DOI: 10.1016/j.compositesb.2007.01.004

    [8]

    PANDEY S, PRADYUMNA S. A finite element formulation for thermally induced vibrations of functionally graded material sandwich plates and shell panels[J]. Composite Structures, 2017, 160: 877-886. DOI: 10.1016/j.compstruct.2016.10.040

    [9]

    PANDEY S, PRADYUMNA S. Transient stress analysis of sandwich plate and shell panels with functionally graded material core under thermal shock[J]. Journal of Thermal Stresses, 2018, 41(5): 543-567. DOI: 10.1080/01495739.2017.1422999

    [10]

    MALEKZADEH P, RAHIDEH H. IDQ two-dimensional nonlinear transient heat transfer analysis of variable section annular fins[J]. Energy Conversion and Management, 2007, 48(1): 269-276. DOI: 10.1016/j.enconman.2006.04.005

    [11]

    CHEN Y, ZHANG L, HE C, et al. Thermal insulation performance and heat transfer mechanism of C/SiC corrugated lattice core sandwich panel[J]. Aerospace Science and Technology, 2021, 111: 106539. DOI: 10.1016/j.ast.2021.106539

    [12]

    STRIZ A, CHEN W, BERT C. Static analysis of structures by the quadrature element method[C]//35th Structures, Structural Dynamics, and Materials Conference. Reston, Virigina: American Institute of Aeronautics and Astronautics, 1994, 3(94): 1567–1574.

    [13] 沈景凤, 张翠, 刘方. 基于热弹性耦合理论的功能梯度材料薄壁旋转碟片力学性能[J]. 复合材料学报, 2019, 36(4): 1017-1028.

    SHEN Jingfeng, ZHANG Cui, LIU Fang. Thermally coupled analysis of thermoelastic field in a thin-walled rotating FGM circular disk[J]. Acta Materiae Compositae Sinica, 2019, 36(4): 1017-1028 (in Chinese).

    [14] 金晶, 邢誉峰, 廖选平, 等. 瞬态传热问题的微分求积和精细积分求解方法[J]. 北京航空航天大学学报, 2015, 41(8): 1526-1531.

    JIN Jing, XING Yufeng, LIAO Xuanping, et al. Application of differential quadrature and precise integration methods in analysis of transient heat transfer[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(8): 1526-1531 (in Chinese).

    [15] 林鹏程, 滕兆春. 热冲击下轴向运动FGM梁的自由振动分析[J]. 振动与冲击, 2020, 39(12): 249-256.

    LIN Pengcheng, TENG Zhaochun. Free vibration analysis of axially moving FGM beams under thermal shock[J]. Journal of Vibration and Shock, 2020, 39(12): 249-256 (in Chinese).

    [16]

    WANG X, GU H. Static analysis of frame structures by the differential quadrature element method[J]. International Journal for Numerical Methods in Engineering, 1997, 40(4): 759-772. DOI: 10.1002/(SICI)1097-0207(19970228)40:4<759::AID-NME87>3.0.CO;2-9

    [17]

    CHEN W, STRIZ A G, BERT C W. A new approach to the differential quadrature method for fourth-order equations[J]. International Journal for Numerical Methods in Engineering, 1997, 40(11): 1941-1956. DOI: 10.1002/(SICI)1097-0207(19970615)40:11<1941::AID-NME145>3.0.CO;2-V

    [18]

    MALEKZADEH P, HEYDARPOUR Y, HAGHIGHI M R G, et al. Transient response of rotating laminated functionally graded cylindrical shells in thermal environment[J]. International Journal of Pressure Vessels and Piping, 2012, 98: 43-56. DOI: 10.1016/j.ijpvp.2012.07.003

    [19]

    DAI T, LI B, TAO C, et al. Thermo-mechanical analysis of a multilayer hollow cylindrical thermal protection structure with functionally graded ultrahigh-temperature ceramic to be heat resistant layer[J]. Aerospace Science and Technology, 2022, 124: 107532. DOI: 10.1016/j.ast.2022.107532

    [20]

    SHEN H S, XIANG Y, REDDY J N. Assessment of the effect of negative Poisson's ratio on the thermal postbuckling of temperature dependent FG-GRMMC laminated cylindrical shells[J]. Computer Methods in Applied Mechanics and Engineering, 2021, 376: 113664. DOI: 10.1016/j.cma.2020.113664

    [21]

    WANG Y, FENG C, YANG J, et al. Nonlinear vibration of FG-GPLRC dielectric plate with active tuning using differential quadrature method[J]. Computer Methods in Applied Mechanics and Engineering, 2021, 379: 113761. DOI: 10.1016/j.cma.2021.113761

    [22]

    WANG X, YUAN Z. Accurate stress analysis of sandwich panels by the differential quadrature method[J]. Applied Mathematical Modelling, 2017, 43: 548-565. DOI: 10.1016/j.apm.2016.11.034

    [23]

    HAO Y X, CAO Z, ZHANG W, et al. Stability analysis for geometric nonlinear functionally graded sandwich shallow shell using a new developed displacement field[J]. Composite Structures, 2019, 210: 202-216. DOI: 10.1016/j.compstruct.2018.11.027

    [24]

    ZENKOUR A M, SOBHY M. Dynamic bending response of thermoelastic functionally graded plates resting on elastic foundations[J]. Aerospace Science and Technology, 2013, 29(1): 7-17. DOI: 10.1016/j.ast.2013.01.003

    [25]

    REDDY J N, CHIN C D. Thermomechanical analysis of functionally graded cylinders and plates[J]. Journal of Thermal Stresses, 1998, 21(6): 593-626. DOI: 10.1080/01495739808956165

    [26]

    WANG B L, TIAN Z H. Application of finite element-finite difference method to the determination of transient temperature field in functionally graded materials[J]. Finite Elements in Analysis and Design, 2005, 41(4): 335-349. DOI: 10.1016/j.finel.2004.07.001

    [27]

    ZHANG Z, ZHOU D, ZHANG J, et al. Transient analysis of layered beams subjected to steady heat supply and mechanical load[J]. Steel and Composite Structures, 2021, 40(1): 87-100.

    [28]

    UYANNA O, NAJAFI H, RAJENDRA B. An inverse method for real-time estimation of aerothermal heating for thermal protection systems of space vehicles[J]. International Journal of Heat and Mass Transfer, 2021, 177: 121482. DOI: 10.1016/j.ijheatmasstransfer.2021.121482

  • 目的 

    功能梯度材料(FGM)夹层结构是一种新型复合材料夹层结构。与蜂窝、泡沫等夹层结构相比,FGM夹层结构通过使组分材料在空间上连续分布,可消除相邻层界面的接触热阻和应力失配,因此FGM夹层结构被广泛运用于高温工程领域。为研究复杂高温工况下FGM夹层板沿厚度方向的传热过程,提出一种基于增量微分求积单元法(IDQEM)的数值分析模型,系统分析了一些关键参数对FGM夹层板传热的影响。

    方法 

    为应用IDQEM,沿层界面将FGM夹层板划分为三个空间子域,同时将整个受热过程划分为若干时间子域。采用微分求积技术对任一时间子域内的控制方程、初始条件、界面条件以及边界条件进行离散处理。由于所获得的离散方程建立在不同区域的节点上,因此对方程进行修改并将其表示为矩阵形式,以便它们可以建立在同一区域中。采用Kronecker积将联立的矩阵方程转化为一系列代数方程组,并采用Newton-Raphson迭代法近似求解,即可获得单个时间子域内的温度解。由于每个时间子域的初始条件可由上一个时间子域最终时刻的温度分布决定,因此从第一个时间子域逐渐递推到最后一个子域,即可获得整个受热过程的温度分布。通过数值算例讨论了本方法的收敛性与正确性。最后,讨论了热工参数温度依赖性、体积分数指数以及热边界条件对FGM夹层板温度分布的影响。

    结果 

    从数值结果可以看出IDQEM具有快速收敛性。与已有文献中基于有限差分法(FDM)的数值解和基于Laplace变换法的解析解的对比,可以验证本方法的正确性。参数化分析结果表明:①对于SUS304/Si3N4夹层板,不考虑材料热工参数温度依赖性状态较考虑温度依赖性状态具有更快的热传导,且两种状态下的温度差异随着表面温度的上升而增加。当上表面温度为600 K时, = 0.01 s、0.02 s、0.03 s三种时刻下的中部温度差异分别达到了12.00%、12.21%、10.41%。②在FGM夹层板受热初期,FGM芯层的体积分数指数对温度分布几乎无影响,但随着加热时间的增长,对温度分布的影响越发显著。当 = 2000 s时,随着从0.2增加到5,Ti-6Al-4V/ZrO夹层板下表面温度从350 K增加到414 K,而上表面温度从1036 K减少到1021 K。③当FGM夹层板受恒定热流作用时,其内部温度随着时间和位置的变化单调变化。然而当热流随时间先增后减时,等温线出现“U”型转折。

    结论 

    基于IDQEM的FGM夹层板非线性瞬态传热分析,考虑了组分材料所有热工参数的温度依赖性以及时变对流辐射条件,适用于复杂高温工况下FGM夹层板的温度场预测。相比传统的数值方法如FDM,本方法具有高精度和高计算效率的优势。参数化分析揭示了材料热工参数的温度依赖性是准确预测板内温度分布的关键因素,尤其是在高温工况中。

  • 功能梯度材料(FGM)夹层结构是一种新型复合材料夹层结构。通过使组分材料在某一方向连续分布,FGM夹层结构的热力学性质在空间内表现为连续变化。相较于传统的蜂窝、泡沫夹层结构,FGM夹层结构可以消除相邻层界面的接触热阻和应力失配,因此被广泛运用于高温工程领域,如飞行器和核反应堆。深入了解FGM夹层结构的传热过程,对其设计和制造具有重要意义。

    本文首次尝试采用增量微分求积单元法(IDQEM)开展了FGM夹层板的一维非线性瞬态传热分析。考虑了所有热工参数的温度依赖性以及时变热边界条件,以提高板内温度场预测的准确性。为应用IDQEM,沿层界面将夹层板划分为三个空间子域,同时将整个受热过程划分为若干时间子域。采用微分求积技术对任一时间子域内的控制方程、初始条件、界面条件以及边界条件进行离散处理,并基于矩阵原理对离散方程进行修改,接着采用Newton-Raphson迭代法近似求解,获得单个时间子域内的夹层板温度解。最后采用递推的方式,即可获得整个受热过程的温度分布。数值算例验证了本方法的快速收敛和高精度优势,并讨论了一些关键参数对FGM夹层板温度分布的影响。结果表明:材料热工参数的温度依赖性是准确预测FGM夹层板温度分布的关键因素,尤其是在高温分析中;在早期阶段,FGM芯层的体积分数指数η对夹层板的温度分布几乎无影响,但随着加热时间的增长,η对温度分布的影响越发显著;当FGM夹层板受恒定热流作用时,其内部温度随着时间和位置的变化单调变化,然而当热流随时间先增后减时,等温线出现“U”型转折。

    FGM夹层板的温度时空分布

图(12)  /  表(4)
计量
  • 文章访问数:  298
  • HTML全文浏览量:  142
  • PDF下载量:  22
  • 被引次数: 0
出版历程
  • 收稿日期:  2023-12-25
  • 修回日期:  2024-01-22
  • 录用日期:  2024-02-06
  • 网络出版日期:  2024-02-28
  • 发布日期:  2024-02-28
  • 刊出日期:  2024-11-14

目录

/

返回文章
返回