经典局部理论模型在分析低温，热冲击以及材料非均匀性较大的传热问题时，热响应结果与实际结果相比偏差较大，采用非局部模型可有效提高精度。J. F. Luciani[1]、G. D. Mahan[2]、S. L. Sobolev[3]及M. Grmela[4]使用非局部热传导模型，得到更精准的低温以及大温度梯度问题下的传热结果。P. Furmanski[5-6]采用总体平均方法给出非局部热传导模型并研究了材料微观结构对热响应结果的影响。近场动力学是一种新型非局部力学理论，其基本方程是一个积分-微分方程，与I. A. Kunin[7]提出的非局部理论相似，基本方程中不存在位移变量对空间坐标的导数，避免了经典局部理论的奇异性。该理论具有天然的无网格数值特征，特别方便进行数值求解[8]。
虽然PD理论已经能够应用到传热学，但还有许多问题需要探索。F. Bobaru和M. Duangpanya[9-10]依据能量平衡法则最先构建出包含裂纹的基于“联结”(bond-based) 的PD热传导方程，并得到PD热流密度函数，通过与解析解对比证明了PD解的正确性，分析了近场范围 (peridynamic horizon) 半径δ的取值对PD解收敛性的影响，计算中考虑了材料点自传导 (self-transfer) 和体积修正 (volume correction)，温度边界代入采用直接对最靠近边界的材料点进行温度赋值。这种边界温度代入方法在离散密度很高会逼近真实情况，网格密度较低时误差较大。A. Agwai[11]给出3种不同的PD热传导内核函数以及材料对应的微热导率计算公式，并发现在PD数值计算时需要对内核函数中的微材料参数进行修正，并给出表面修正因子完善了PD数值计算中的误差。随后，S. Oterkus和E. Madenci等[12]依据Euler-Lagrange方程导出基于状态 (state-based) 的PD热传导方程，并提出虚拟层温度边界条件的代入方法，但文中算例仍是基于“联结”的均匀材料热传导计算。Ziguang Chen和F. Bobaru[13]采用“形状因子”将3种热传导内核函数进行整合，研究“形状因子”对PD解的δ收敛性影响，分析温度边界代入方法、自传导项、表面修正因子等因素对均匀材料热传导PD解的影响。
当前PD方法大多是针对均匀材料进行的计算和分析，深入研究PD非均匀材料的传热已是亟待解决的问题。等效热导率比例系数、表面修正因子以及形状因子是PD法的内核参数中重要参数。本文将集中研究这几个内核参数对PD数值计算过程和非均匀材料瞬态温度响应的影响，可为采用PD理论分析非均匀材料热传导问题时提供参数选择参考和正确性验证依据。
1 非均匀材料的PD热传导理论近场动力学与经典局部理论的不同点在于PD考虑的是材料点近场范围内所有其他材料点与中心材料点间的相互作用，而经典局部理论只考虑与材料点临近接触点间的相互作用。近场范围指以材料点x为中心，有限大小半径δ的邻域，如图 1所示。平均划分区域网格，任一离散域的面积为Δ2，且其中只包含一个材料点。
|
| 图 1 PD空间离散和近场范围 |
在热传导中，材料点的相互作用就是以热流密度表达的热能交换。PD理论的非均匀材料热传导方程与均匀材料热传导方程形式一致[5]
|
(1) |
式中, ρ(x) 和c (x) 分别表示材料的密度和比热, hs(x, t) 表示单位体积上的热生成量。fq(x′, x, t) 表示x点近场范围内所有材料点x′与点x间的相互作用, 它本质上是一个双热流密度函数, 也被称作热响应函数或者“内核”函数。基于联结下的内核函数表达如下[6]
|
(2) |
式中, κ(x′, x) 表示材料的微热导率 (micro-conductivity), 该材料系数为PD理论中的热导率。整数n被称为形状因子, 该参数通常取为0、1和2。ξ是点x′和x间的向量差。τ(x′, x, t) 是t时刻点x′和x温差, 可表示为
|
(3) |
二维热传导常数形式的κ(x′, x) 表达式如下
|
(4) |
式中, h是材料点x在坐标z方法向上的厚度, 一般取为Δ。k为经典理论热导率, 它是近场范围内2个材料点热导率的函数, 可以写为
|
(5) |
式中, φ为等效热导率比例系数, 用来描述材料点物性值所占比重。ρ(x)、c (x) 以及k (x′) 和k (x) 由非均匀材料物性值函数决定。
将时间微分项替换为向前差分格式进行, 对空间积分项进行离散并采用Gauss中点数值积分, PD热传导方程变为
|
(6) |
式中, 下标i表示第i个材料点, j表示材料点i近场范围内的第j个材料点。N为材料点i近场范围内的材料点总数。上标m表示第m个离散时刻。Δt为时间步长, V(j)表示材料点j的体积。g(i)(j)称为表面修正因子, v(j)为体积修正因子[5]。
以图 2所示的非均匀材料矩形板为算例, 应用分离变量法与PD数值方法求解二维热传导问题。长度和宽度为l=b=20 mm, 四周作用边界温度T1=100 K, 初始温度为T0=0 K, 假定材料热物性值为指数函数:k=Decx+dy, ρc=Eecx+dy, 其中D、E以及c和d为材料非均匀参数, 取为D=E=1.0, c=d=0.1。其分离变量解为
|
(7) |
|
| 图 2 算例 |
式中
|
(8) |
PD算法中矩形板x方向和y方向上离散材料点总数均为80, Δ=0.25 mm, 近场范围δ=3.015 Δ, 时间步长Δt=10-4 s, 取形状参数n=1。对比t=1 s、5 s和9 s时x=10.125 mm板宽度方向上虚线的温度分布, 如图 3所示。对比2种方法的计算结果可知, 本文PD数值求解程序正确可靠。
|
| 图 3 矩形板x=10.125 mm宽度方向上温度验证 |
为了说明PD理论相比传统局部理论在传热学中的优势和应用价值, 设置图 2算例边界为低温1 K, 材料函数取法不变。给出PD方法, 传统FEM以及分离变量法在1.0 s时的结果对比, 表 1为图 2中从边界到中心实线上离散点的温度结果, 靠近矩形板中心位置的离散点温度接近为零。分析表 1结果可知, 在非均匀材料低温导热时, 采用PD方法的计算结果整体相比传统FEM结果误差要小。因此, PD非局部热传导模型更适用于低温以及材料非均匀性较大时的传热问题。
| PD/K | FEM/K | SVM/K | Error: PD | Error: FEM |
| 0.930 | 0.921 | 0.924 | 0.7% | 0.3% |
| 0.777 | 0.768 | 0.776 | 0.2% | 1.0% |
| 0.638 | 0.627 | 0.638 | 0.0% | 1.7% |
| 0.511 | 0.500 | 0.513 | 0.3% | 2.5% |
| 0.400 | 0.390 | 0.403 | 0.6% | 3.2% |
| 0.306 | 0.297 | 0.309 | 0.7% | 3.6% |
| 0.229 | 0.222 | 0.231 | 0.8% | 3.7% |
| 0.167 | 0.162 | 0.168 | 0.7% | 3.4% |
| 0.119 | 0.117 | 0.119 | 0.2% | 2.3% |
| 0.083 | 0.082 | 0.083 | 0.6% | 0.6% |
| 0.057 | 0.057 | 0.056 | 1.8% | 2.2% |
| 0.038 | 0.039 | 0.036 | 3.6% | 6.3% |
| 0.025 | 0.026 | 0.023 | 6.1% | 11.7% |
| 0.016 | 0.017 | 0.014 | 9.1% | 18.9% |
| 注：SVM表示由分离变量法得到的解析解。 | ||||
均匀材料的微热导率κ可由PD热势能与经典力学热势能等价条件得到。在计算κ时 (4) 式中近场范围内任意材料点的热导率是相等的, 而非均匀材料的热导率则需要考虑近场中心材料点与其它材料点的比重。为了研究比例系数φ对解的影响, 取φ为0.0、0.3、0.5、0.7以及1.0, 对比在不同形状因子n下的计算结果, 计算参数不变。如图 4所示为t=5.0 s时, 材料点 (10.125, 4.125) 的温度值与解析解的对比。
|
| 图 4 微热导率模型比例系数对温度的影响 (时间t=5.0 s) |
由图 4知, 不同比例系数φ下的温度响应基本成线性分布, 当系数φ取0.5时, PD解与解析解符合最好, φ离0.5越远PD解误差越大, 微热导率中的k取两点平均值是可行的。形状因子n对这一规律没有影响。
2.2 表面修正因子对非均匀材料PD解的影响为分析表面修正因子对PD解的影响, 计算不同形状因子n时考虑与不考虑表面效应2种情况的瞬态PD温度响应, 表 2给出材料点 (10.125, 4.125) 在t=5 s时的温度。
| g(i)(j) | T(n=0)/K | T(n=1)/K | T(n=2)/K |
| True | 15.898 | 15.865 | 15.834 |
| False | 12.487 | 13.030 | 13.165 |
| 注：True表示考虑表面修正因子g(i)(j)的情况, False为不考虑表面修正因子g(i)(j)的情况。 | |||
由此可知, 非均匀材料在温度边界下考虑表面修正因子的PD解误差更小 (分离变量解为15.822 K)。另外结合图 5, 还可以看出当微热导率取两点平均值时, 形状因子n取2时PD解与解析解最接近, 但影响并不明显。
|
| 图 5 计算近场范围对热响应的影响 |
材料点离散密度不变, 近场范围可由δ=mΔ控制。取m分别等于9.015、7.015、5.015、3.015以及1.015, PD解与分离变量解的对比如图 5所示。由图 5可知, 当m→0时 (即δ→0), 非均匀材料瞬态热传导的PD解与其他m值相比更趋于经典解, 但并不存在m值越大误差越大的单调现象。
3 结论采用近场动力学理论计算非均匀材料瞬态热传导响应时, 内核参数不同会影响温度响应。经研究, 各内核参数影响如下:
1) 不同比例系数φ下的温度响应可认为成线性分布, 当系数φ取0.5时, PD解最接近解析解, φ离0.5越远PD解误差越大, 微热导率取近场范围两点平均值最为可靠, “内核”形状系数n对这一规律没有影响。
2) 非均匀材料在温度边界下考虑表面修正因子的PD解更接近解析解。另外, 当微热导率取两点平均值形状系数n取2时, PD解收敛性更好。
3) 当m→0时, 非均匀材料瞬态热传导的PD解误差最小, 更趋于经典解, 采用PD理论求解非均匀材料热响应时, 仍需将δ取为材料点间距Δ的3倍。
| [1] | Luciani J F, Mora P, Virmont J. Nonlocal Heat Transport Due to Steep Temperature Gradients[J]. Physics Review Letters, 1983, 51(18): 1664–1667. DOI:10.1103/PhysRevLett.51.1664 |
| [2] | Mahan G D, Claro F. Nonlocal Theory of Thermal Conductivity[J]. Physics Review B, 1988, 38(3): 1963–1969. DOI:10.1103/PhysRevB.38.1963 |
| [3] | Sobolev S L. Equations of Transfer in Non-Local Media[J]. International Journal of Heat and Mass Transfer, 1994, 37(14): 2175–2182. DOI:10.1016/0017-9310(94)90319-0 |
| [4] | Grmela M, Lebon G. Finite-Speed Propagation of Heat: a Nonlocal and Nonlinear Approach[J]. Physica A, 1998, 248: 428–441. DOI:10.1016/S0378-4371(97)00552-9 |
| [5] | Furmanski P. Effective Macroscopic Description for Heat Conduction in Heterogeneous Materials[J]. International Journal of Heat and Mass Transfer, 1992, 35(11): 3047–3058. DOI:10.1016/0017-9310(92)90324-L |
| [6] | Furmanski P. A Mixture Theory for Heat Conduction in Heterogeneous Materials[J]. International Journal of Heat and Mass Transfer, 1994, 37(18): 2993–3002. DOI:10.1016/0017-9310(94)90353-0 |
| [7] | Kunin I A. Elastic Media with Microstructure I: One-Dimensional Model[M]. Springer Press, 1982. |
| [8] | Silling S A. Reformulation of Elasticity Theory for Discontinuities and Long-Range Forces[J]. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209. DOI:10.1016/S0022-5096(99)00029-0 |
| [9] | Bobaru F, Duangpanya M. The Peridynamic Formulation for Transient Heat Conduction[J]. International Journal of Heat and Mass Transfer, 2010, 53(19): 4047–4059. |
| [10] | Bobaru F, Duangpanya M. A Peridynamic Formulation for Transient Heat Conduction in Bodies with Evolving Discontinuities[J]. Journal of Computational Physics, 2012, 231(7): 2764–2785. DOI:10.1016/j.jcp.2011.12.017 |
| [11] | Abigail Agwai. A Peridynamic Approach for Coupled Fields[D]. Tucson: University of Arizona, 2011 |
| [12] | Oterkus S, Madenci E, Agwai A. Peridynamic Thermal Diffusion[J]. Journal of Computational Physics, 2014, 265(10): 71–96. |
| [13] | Chen Ziguang, Bobaru Florin. Selecting the Kernel in a Peridynamic Formulation: A Study for Transient Heat Diffusion[J]. Computer Physics Communications, 2015, 197: 51–60. DOI:10.1016/j.cpc.2015.08.006 |
| [14] | Cheng Zhanqi, Zhang Guanfeng, Wang Yenan, Bobaru Florin. A Peridynamic Model for Dynamic Fracture in Functionally Graded Materials[J]. Composite Structures, 2015, 133: 529–546. DOI:10.1016/j.compstruct.2015.07.047 |
| [15] |
黄丹, 章青, 乔丕忠, 等.
近场动力学方法及其应用[J]. 力学进展, 2010, 40 (4): 448–459.
Huang Dan, Zhang Qing, Qiao Pizhong, et al. A Review on Peridynamics (PD) Method and Its Applications[J]. Advances in Mechanics, 2010, 40(4): 448–459. DOI:10.6052/1000-0992-2010-4-J2010-002 (in Chinese) |
| [16] |
章青, 顾鑫, 郁杨天.
冲击载荷作用下颗粒材料动态力学响应的近场动力学模拟[J]. 力学学报, 2016, 48 (1): 56–63.
Zhang Qing, Gu Xin, Yu Yangtian. Peridynamics Simulation for Dynamic Response of Granular Materials Under Impact Loading[J]. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(1): 56–63. DOI:10.6052/0459-1879-15-291 (in Chinese) |
