A METHOD FOR COMPUTING THE EFFECTIVE IN-PLAN THERMOELASTIC PROPERTIES OF PERIODIC COMPOSITES BASED ON THE RADIAL INTEGRATION BOUNDARY ELEMENT METHOD
-
摘要: 具有人工设计结构的多材料夹杂的复合材料不仅具有可编辑的弹性模量和泊松比, 还可以实现热膨胀系数的调控. 这些复合材料的等效热弹性参数与其微观结构和材料分布密切相关. 文章提出了一种基于径向积分边界元法的周期性复合材料的均匀化热弹性参数的计算方法. 该方法基于代表性体积元均匀化方法建立材料的微观结构应力、应变与宏观等效热弹性参数之间的关联. 同时, 采用径向积分边界元法计算周期变形条件下代表性体积元的应力和应变场. 该方法无需内部网格和域积分, 仅依赖于边界的位移和面力等信息即可获得材料的等效热弹性参数, 具有容易实现参数化建模、容易考虑细小结构特征等优势. 另外, 在角点处采用非连续元处理, 大大简化了周期性边界条件的施加. 针对复合材料梁的弯曲问题, 通过与直接有限元模拟的结果进行对比, 验证了该方法的有效性. 然后, 还采用该方法对网状复合超材料的热膨胀系数进行了分析. 结果表明, 通过调整结构的尺寸, 可实现结构在热载荷下从正膨胀到负膨胀的转变.Abstract: One can program not only the elastic modulus and Poisson’s ratio but also the thermal expansion coefficients of the composites which consist of artificial microstructures and multiple materials. The macro effective thermoelastic properties of the composites are determined by their microstructures. In this paper, the authors proposed a novel method to evaluate the homogenized in-plan thermoelastic properties of the composites made of the periodic microstructures, which is based on the radial integration boundary element method. The proposed method employs the representative volume elements homogenized approach to calculate the effective properties based on the displacements and stresses of the micro unit cells, while utilizing the radial integral boundary element method for solving the displacements and stresses. The proposed method does not require internal grid and domain integration, and it only relies on the displacements and surface tractions on the boundary to obtain the effective thermoelastic parameters of materials, such as elastic modulus and thermal expansion coefficients. This method has advantages such as easy implementation of parameterized modeling, which leads to priorities in structure optimization. Besides, using the boundary element method makes it easy to analyze accurately the cells with tiny structural features. What’s more, the discontinuous elements are employed at the conners of structures in this paper to simplify the implementation of the periodic boundary conditions of the representative volume elements. The authors simulated the banding problem of a composite beam using both the homogenized thermoelastic properties and the direct finite element simulation. The results of the homogenized model coincide the direct simulation well, which indicates the effectiveness of the proposed method. Using the proposed method, the authors also conducted the parametric analysis of the coefficients of thermal expansion of a cellular microstructure. The results show that one can achieve the negative thermal expansion by changing the sizes of the cellular structure’s components.
-
引 言
太阳翼、天线及桅杆等航天器中的大型空间结构在剧烈的温度变化下, 结构稳定性受到了严重的影响[1-2]. 由于卫星的日出和日落过程时间极短, 在强烈的太阳辐射下, 失重、低温以及真空环境中的太阳翼等薄细结构会在短短数分钟内升高/降低200 K以上, 从而引发热致振动. 极端情况下, 快速的温度变化还会产生热颤振问题[3]. 随着增材制造[4-5]技术的发展, 可以通过改变材料的微观结构实现对材料热力性能的调控, 例如设计具有近零热膨胀系数[6-8]的超材料以削弱热致振动. 为了设计出符合特定需求的复合材料, 必须准确地评估设计参数对材料力学、热学行为的影响规律.
目前, 对材料的弹性模量、热膨胀系数及热导率等的试验测定方法已经相当成熟. 然而, 如果基于宏观层面的拉伸和扭转等试验来研究设计参数对结构性能的影响, 需要制备大量的样件. 而且, 随着设计参数的增加, 所需测试的样本量将会呈指数增长, 成本巨大. 因此, 理论手段或数值手段的材料等效物性分析也十分重要. 目前, 主流的超材料结构设计仍然以网状、点阵及蜂窝等形式为主. 这些结构由梁、杆和薄板等组成, 理论分析较为容易, 并且可以得到特定参数对材料物性影响规律的解析表达式[9-12]. 然而, 这些结构杆件等较细, 其承载力比常规材料小得多. 为了提高结构的承载力, 同时保证结构满足拉胀等特性, 其结构外形将会更加复杂[13], 很难简化为梁、杆的组合, 因而采用解析方法来分析其力学性能十分困难. 对于泡沫和蜂窝等具有周期性微结构的材料, 常采用代表性体积元(representative volume elements, RVE)均匀化方法评估其宏观力学性能[14-17].
相比于泊松比等力学性能, 目前针对材料的等效热膨胀系数的讨论较少. 但材料的热膨胀系数对承受高热载荷的航天器等的设计十分重要. Lim[18]基于解析方法, 对网状材料的热膨胀系数进行了分析, 并设计出了一种具有负热膨胀系数的材料. Hassanzadeh-Aghdam等[19]借助简化单晶胞模型对玻璃纤维增强的高分子复合材料的等效热膨胀系数进行了分析. Tian等[20-21]采用基于有限元法的RVE均匀化方法对复合材料的等效热力参数进行了分析, 并对边界条件的施加方式进行了探讨.
目前, RVE更多采用有限元(finite element method, FEM)求解. 但是, 有限元需要在求解域内部划分网格, 因而网格需求量一般较大, 而且对网格质量要求较高. 特别是对于焊缝和内凹角等细小特征, 难以建立高质量的有限元网格[22]. 边界元法(boundary element method, BEM)由于只需要在边界上划分网格, 在纤维材料多尺度分析[23]、弹性波、声学[24-25]和拓扑优化[26-27]等方面有着广泛的应用. 近些年, 基于边界元的RVE均匀化分析在计算材料的等效弹性常数中也得到了一定发展[13,17,28-30]. 然而, 当温度变化时, 热变形项导致的域积分大大增加了边界元法计算的复杂度. 截至目前, 尚未见到有关边界元法计算复合材料等效热弹性参数的讨论. 如今, 边界元法域积分的处理方法已经无需内部网格, 采用径向积分法[31]和双重互易法[32]等均可将其转化到边界上. 其中, 径向积分法适用于任意函数, 甚至是包含微分算子的项, 在结构热弹性问题分析中应用广泛[33].
本文拟基于多域边界元方法建立考虑热变形的复合材料RVE内部应力和应变场的计算方法, 以分析周期性复合材料的面内等效热力参数. 为了解决热应力导致的域积分问题, 将采用Gao等[31,34]提出的径向积分法将域积分转换为边界积分, 因而无需划分内部网格. 本文方法仅需边界网格和边界未知量, 具有容易处理复杂结构、容易实现参数化建模等优点. 因此, 本文方法可为后续的复合超材料的拓扑优化设计奠定基础.
1. 等效热力参数计算方法
1.1 热弹性问题控制方程
$$\qquad\qquad\quad \frac{{\partial {\sigma _{ij}}}}{{\partial {x_j}}} = {f_i} $$ (1a) $$\qquad\qquad\quad \frac{\partial }{{\partial {x_i}}}k\frac{{\partial T}}{{\partial {x_i}}} = Q $$ (1b) 其中, $ {f_i} $和$ Q $为体积力和热源. 材料参数并不受外载荷的影响, 因而在本文中无需特别讨论源项的处理. $k$为热导率, ${\sigma _{ij}}$为应力. 根据热弹性本构关系
$$ {\sigma _{ij}} = {D_{ijkl}}({\varepsilon _{kl}} - {\alpha _{kl}}\Delta T) $$ (2) 式中, $ {\alpha _{kl}} $为热膨胀系数, 对于各向同性材料${\alpha _{ij}} = {\delta _{ij}}{\alpha ^0}$.
可确定应力、温度和应变之间的关系. 本文考虑的组元材料均为各向同性材料, 弹性常数张量为
$$ {D_{ijkl}} = \mu \left( {\frac{{2\nu }}{{1 - 2\nu }}{\delta _{ij}}{\delta _{kl}} + {\delta _{ik}}{\delta _{jl}} + {\delta _{il}}{\delta _{jk}}} \right) $$ (3) 式中, $\mu $为剪切模量, $\nu $为泊松比. 应变由位移的梯度计算
$$ {\varepsilon _{ij}} = \frac{1}{2}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right) $$ (4) 1.2 径向积分边界元求解热弹性问题
引入位移基本解$ u_{ij}^* $和面力基本解$ t_{ij}^* $
$$ u_{ij}^* = \frac{1}{{8\text{π} (1 - \nu )\mu }}\left[ { - (3 - 4\nu ){\delta _{ij}}\ln r + {r_{,i}}{r_{,j}}} \right] $$ (5) $$\begin{split} &t_{ij}^* = \frac{{ - 1}}{{4\alpha \text{π} (1 - \nu ){r^\alpha }}}\frac{{\partial r}}{{\partial n}}\left\{ {\left[ {(1 - 2\nu ){\delta _{ij}} + (\alpha + 1){r_{,i}}{r_{,j}}} \right]}- \right.\\ &\qquad \left. { (1 - 2\nu )({r_{,i}}{n_j} - {r_{,j}}{n_i})} \right\} \end{split} $$ (6) 对于平面问题$\alpha {\text{ = }}1$, $r$为从源点$P$到场点$Q$的距离. 然后, 通过加权余量法可以得到热弹性问题的边界积分方程[31]
$$ {c_{ij}}{u_j} = \int_\varGamma {(u_{ij}^{\text{*}}{t_j} - t_{ij}^*{u_j})} {\text{d}}\varGamma + \int_\varOmega {u_{ij,j}^{\text{*}}} \tilde \alpha \Delta T{\text{d}}\varOmega $$ (7) 式中, 系数${c_{ij}}$由其所在位置决定, 在内部时为1. ${t_j} = {\sigma _{ij}}{n_j}$为边界面力, $ \tilde \alpha = \dfrac{{2\left( {1 + \nu } \right)}}{{1 - 2\nu }}\mu {\alpha ^0} $. 注意到, 式(7)右侧第2项为域积分, 因此, 引入变量
$$ {\varPsi _i} = u_{ij,j}^*\tilde \alpha = - \frac{{(1 + \nu ){\alpha ^0}}}{{2\alpha \text{π} (1 - \nu ){r^\alpha }}}{r_{,i}} $$ (8) 对于平面问题, 将域积分转化到以源点为中心的极坐标下可得
$$ {\text{d}}\varOmega = r{\text{d}}r{\text{d}}\varphi {\text{ = }}r{\text{d}}r\frac{1}{r}\frac{{\partial r}}{{\partial n}}{\text{d}}\varGamma $$ (9) 因此, 式(7)中的域积分转化为
$$\int_\varOmega {u_{ij,j}^{\text{*}}} \tilde \alpha \Delta T{\text{d}}\varOmega {\text{ = }}\int_\varOmega {{\varPsi _i}\Delta T{\text{d}}\varOmega } = \int_\varGamma \frac{{\partial r}}{{\partial n}}{\varPsi _i}\bar F(Q){\text{d}}\varGamma $$ (10) 式中, $ \bar F(Q){\text{ = }}\displaystyle\int_0^{r(Q)} \Delta T{\text{d}}r $. 经过式(10)的变换, 可将域积分项转化成边界积分项, 因而无需构建求解域内部的网格. 位移和面力在边界上的分布可由单元形函数插值获得
$$ \phi ({\xi ^e}){\text{ = }}{\phi ^\alpha }N_\alpha ^e $$ (11) 式中, $ {\xi ^e} $为单元局部坐标, $ N_\alpha ^e $为节点$\alpha $在单元$e$中的形函数. 将式(10)和式(11)代入式(7)可得离散形式下的弹性力学方程
$$ {\boldsymbol{Hu}} + {\boldsymbol{Gt}} = {\boldsymbol{f}} $$ (12) 其中, $ {\boldsymbol{u}} $和$ {\boldsymbol{t}} $为位移和面力构成的向量. 系数矩阵分量可计算为
$$\qquad H_{ij}^{\alpha \beta } = - c_{ij}^{\alpha \beta } - \sum\limits_e {\int_\varGamma {t_{ij}^{*\alpha }N_\beta ^e} {\text{d}}\varGamma } $$ (13a) $$\qquad G_{ij}^{\alpha \beta } = \sum\limits_e {\int_\varGamma {u_{ij}^{*\alpha }N_\beta ^e} {\text{d}}\varGamma }$$ (13b) $$\qquad f_i^\alpha = \int_\varGamma \frac{{\partial r}}{{\partial n}}\varPsi _i^\alpha \bar F{\text{d}}\varGamma $$ (13c) 式中, 上标$\alpha $表示由以点$\alpha $为源点的基本解所计算的值.
复合材料往往由多种具有不同物性的材料组成. 因此, 根据材料的不同, 可将整个求解域划分为多个子域. 在公共界面上, 两侧的子域满足位移和应力的连续性条件, 即
$$\qquad\qquad\qquad\qquad u_i^ + = u_i^ - $$ (14a) $$\qquad\qquad\qquad\qquad t_i^ + + t_i^ - {\text{ = }}0 $$ (14b) 类似地, 引入温度和热通量的基本解
$$\qquad\qquad\qquad {T^*} = \frac{1}{{2\text{π} }}\ln \left( {\frac{1}{r}} \right) $$ (15a) $$\qquad\qquad\qquad {q^{\text{*}}} = kT_{,i}^*{n_i} $$ (15b) 再采用加权余量法可得温度场的边界元离散方程
$$ {\boldsymbol{\tilde HT}} + {\boldsymbol{\tilde Gq}} = {\boldsymbol{0}} $$ (16) 1.3 等效热力参数的计算
如图1所示, 整个RVE所占据的区域为${V_0}$(包括轮廓线内的空白区域). 其平均应变为
![]()
图 1 某周期性复合材料的RVE, s和d代表一对周期边界的源面和目标面Figure 1. The RVE of a periodic composite, where s and d represent the source face and target face of a couple of periodic boundary$$ {\bar \varepsilon _{ij}} = \frac{1}{{{V_0}}}\int_{{V_0}} {{\varepsilon _{ij}}{\text{d}}\varOmega } $$ (17) 将式(4)代入式(17)并采用格林-高斯公式可得
$$ {\bar \varepsilon _{ij}} = \frac{1}{{2{V_0}}}\int_{\partial {V_0}} {{u_i}{n_j}{\text{ + }}{u_j}{n_i}{\text{d}}\varGamma } $$ (18) 同样地, 平均应力为
$$ {\bar \sigma _{ij}} = \frac{1}{{{V_0}}}\int_{{V_0}} {{\sigma _{ij}}{\text{d}}\varOmega } {\text{ = }}\frac{1}{{2{V_0}}}\int_{{V_0}} {{t_i}{x_j}{\text{ + }}{t_j}{x_i}{\text{d}}\varGamma } $$ (19) 根据本构关系, 等效的弹性常数张量应满足
$$ {\bar \sigma _{ij}} = {\bar D_{ijkl}}{\bar \varepsilon _{kl}} $$ (20) 根据应力和应变的对称性可知, ${\bar D_{ijkl}} = {\bar D_{jikl}} = {\bar D_{ijlk}}$. 由于是平面线性问题, 任意应变场都可以看作是如下应变的叠加
$$ {\bar {\boldsymbol{\varepsilon}} ^{11}} = \left[ {\begin{array}{*{20}{c}} 1&0 \\ 0&0 \end{array}} \right],\quad {\bar {\boldsymbol{\varepsilon}} ^{22}} = \left[ {\begin{array}{*{20}{c}} 0&0 \\ 0&1 \end{array}} \right],\quad {\bar {\boldsymbol{\varepsilon}} ^{12}} = \left[ {\begin{array}{*{20}{c}} 0&1 \\ 1&0 \end{array}} \right] $$ (21) 因此, 只需求得3种不同平均应变下对应的平均应力, 等效的弹性常数张量就可以唯一确定. 值得注意的是, 对于线性问题, 给定的应变大小并不会影响结论.
由于RVE的特征尺寸远小于宏观结构的特征尺寸, 单个RVE内部温度可以视作是均匀的. 此时, 平均意义下的本构关系可写为
$$ {\bar \sigma _{ij}} = {\bar D_{ijkl}}({\bar \varepsilon _{kl}} - {\bar \alpha _{kl}}\Delta \bar T ) $$ (22) 因此, 等效热膨胀系数为
$$ {\bar \alpha _{kl}}{\text{ = }}\frac{{{{\bar \varepsilon }_{kl}}{{ - }}{{\bar D}_{ijkl}}^{{{ - }}1}{{\bar \sigma }_{ij}}}}{{\Delta \bar T }} $$ (23) RVE的平均热流和平均温度梯度为
$$ {\bar q_i} = \frac{1}{{{V_0}}}\int_{{V_0}} {{q_n}{x_i}{\text{d}}\varGamma } $$ (24) $$ {\bar T_{,i}} = \frac{1}{{{V_0}}}\int_{{V_0}} {T{n_i}{\text{d}}\varGamma } $$ (25) 等效热导率应满足
$$ {\bar k_{ij}} = - {\bar q_i}/{\bar T_{,j}} $$ (26) 1.4 周期性边界条件
Hori等[35]研究表明, 采用周期边界条件可以获得比线性位移条件更可靠的结果. 在本文中, 位移和温度条件采用如下周期性边界条件(见Xia等[36] 的论文)
$$ {u_i}({{\boldsymbol{x}}^s}) - {u_i}({{\boldsymbol{x}}^d}) = {\bar \varepsilon _{ij}} \cdot (x_j^d - x_j^s) $$ (27a) $$ T({{\boldsymbol{x}}^s}) - T({{\boldsymbol{x}}^d}) = {\bar T_{,j}} \cdot (x_j^d - x_j^s) $$ (27b) 在边界元中, 面力和热流是基本变量, 还需要根据RVE的周期性和平衡条件补充如下关系
$$\qquad\qquad\quad {t_i}({{\boldsymbol{x}}^s}) + {t_i}({{\boldsymbol{x}}^d}) = 0 $$ (28a) $$\qquad\qquad\quad {q_n}({{\boldsymbol{x}}^s}) + {q_n}({{\boldsymbol{x}}^d}) = 0 $$ (28b) 在式(27)和式(28)中, $ {{\boldsymbol{x}}^s} $和$ {{\boldsymbol{x}}^d} $代表着一对周期性边界上的点(即图1所示的s边界和d边界). 对于周期边界上的角点, 情况略微复杂[21], 尤其是面力. 对于平面问题, RVE的周期边界上有4个角点. 如果是三维问题, 则有8个角点和12条边线, 情况较多. 边界元法对于位移和面力在边界上的连续性没有要求, 而基本解的奇异性使得即使两个点十分靠近, 也不会使离散系统的系数矩阵奇异, 因而在边界元法中可以使用非连续元. 为了规避角点识别和处理问题, 本文在尖锐的边界处采用非连续元. 非连续元上的节点不会被多个边界所共享, 极大简化了周期性边界条件的施加.
2. 计算结果及分析
2.1 变形场计算精度验证
考虑如图2(a)所示的两种场材料构成的复合结构, 选取如图2(b)所示的结构为RVE. 组元材料的物性如表1所示. 考虑式(21)所给的3种平均应变, 并按式(27)加载RVE的边界条件. 分别准备了两套疏密不同的边界元网格和有限元网格. 其中疏密边界元网格分别包含了80和160个线单元. 稀有限元网格为1790个线性三角形单元, 密有限元网格有
30276 个6节点三角形单元. 有限元和边界元的稀网格在图3中给出.表 1 图2复合材料组元的物性Table 1. The properties of the components of the composite shown in Fig. 2Parameter Mcterial A Mcterial B Young modulus/MPa $1 \times {10^4}$ $5 \times {10^4}$ Poisson's ratio 0.25 0.3 thermal expansion coefficient/K−1 $1 \times {10^{{{ - 5}}}}$ $1 \times {10^{{{ - }}6}}$ thermal conductivity/(mW·K−1·mm−1) 20 10 图4给出了3种不同平均应变对应的RVE内部von Mises应力. 由于材料B的刚度较大, 材料B的应力也更大. 从图5可以看出, 即使采用相当少的单元, 边界元法也能准确地计算出边界面力. 边界元法中位移与面力的精度相同, 因而即使只有80个线单元, 边界元法计算出的面力也与有限元法加密之后的网格几乎相同.
![]()
图 4 3种平均应变下径向积分边界元法计算出的变形(缩小10倍)和应力Figure 4. The deformations (scaled down to 10 times) and the von Mises stresses with the three averaged strains![]()
图 5 有限元和边界元法计算出的右侧面力大小Figure 5. The magnitudes of the traction on the right edge computed by FEM and BEM对于该案例中的结构, 由于其中的夹杂相区域为圆形, 可采用Zheng等[37]提出的高精度闭合单元计算, 该方法对类椭圆纤维夹杂的问题具有很好的效果. 本文考虑一般的构型, 在此不做过多讨论.
2.2 均匀化等效参数验证
根据式(20)、式(23)和式(26)可以分别获得材料的等效弹性常数张量、热膨胀系数和热导率如下
$$ {\bar D_{1111}} = {\bar D_{2222}} = 1.23 \times {10^{10}}{\text{ MPa}} $$ $$ {\bar D_{1122}} = {\bar D_{2211}} = 3.12 \times {10^9}{\text{ MPa}} $$ $$ {\bar D_{1212}} = {\bar D_{2121}} = 9.08 \times {10^9}{\text{ MPa}} $$ $$ {\bar D_{1112}} = {\bar D_{2212}} = {\bar D_{1222}} = {\bar D_{1211}} = 0 $$ $$ {\bar \alpha _{11}} = {\bar \alpha _{22}} = 8.52 \times {10^{ - 6}}{{\text{ K}}^{{{ - 1}}}} \text{, } \quad {\bar \alpha _{21}} = {\bar \alpha _{12}}{\text{ = }}0 $$ $$ {\bar k_{11}} = {\bar k_{22}} = 18.5{\text{ W/}}\left( {{\text{m}} \cdot {\text{K}}} \right) \text{, }\quad {\bar k_{12}} = {\bar k_{21}} = 0 $$ 得到热力参数之后, 分别采用均匀化的热力参数和直接有限元模拟计算图2(a)中的结构. 直接有限元模拟采用的网格如图6所示, 为确保结果可靠, 共采用了
713036 个6节点三角形单元. 考虑如下两种加载情况: (1)仅存在如图2(a)所示的压力载荷, 大小为1000 kPa; (2)压力载荷为0, 左侧和下侧固定300 K温度, 上侧施加100 kW/m2的热流.两种加载情况下横向位移如图7所示. 由于第2种加载情况上侧温度高, 而下侧温度低, 使得梁也产生了向下的弯曲. 如图8所示, 无论是热载荷还是压力载荷, 采用均匀化模型计算出的位移与直接采用有限元法对细观模型模拟的结果几乎一致, 可见本文所采用方法的有效性.
![]()
图 8 两种加载情况下中性面竖向位移分布Figure 8. The vertical displacements on the middle face under the two loads2.3 网状结构的热膨胀系数
继续讨论由两种材料所构成的网状结构的热膨胀系数(单个胞元如图9所示). 其中浅色和深色材料分别为表1中所列的材料A和B. 两种材料的热膨胀系数差别较大, 因此随着两种材料的组合方式不同会产生不同的热膨胀系数. 对于较薄的网格线, 可以基于梁模型采用解析方式求出[9]. 由于该网状结构网线厚度较大, 采用梁模型计算误差大.
为了验证径向积分边界元法的精度, 首先对该胞元(t = 0.05 mm)的4个连接面施加固定约束, 然后对整个结构施加100 K的温升. 分别采用有限元法和径向积分边界元法对材料热变形引起的应力场和变形场进行分析. 边界元法所采用网格(BEM mesh 1)如图10所示, 其包含了520个线单元. 另一个加密的网格(BEM mesh 2)则包含1040个线单元, 此处未画出. 而有限元网格(FEM mesh 1)如图11所示, 包含了4914个二次三角形单元, 并在角点处进行了加密, 共计
10581 个节点. 另外两个网格 (FEM mesh 2和 mesh 3)则分别包含了11580 个二次三角形单元和33686 个线性三角形单元, 其中FEM mesh 3并未对角点特别加密.![]()
图 10 边界元网格, 共计520个线单元及518个节点 (单位: mm)Figure 10. The mesh for BEM: There are 520 line elements and 520 nodes (unit: mm)![]()
图 11 有限元网格, 共计4914个6节点三角形单元及10581 个节点 (单位: mm)Figure 11. The mesh for BEM: There are 4914 six-node triangular elements and10581 nodes (unit: mm)如图12所示, 在结构的内凹角点处和材料A和B的连接处出现了明显的应力集中. 而由于此时4个连接面的位移全部被约束, 其附近也产生了很大的应力, 其角点也有应力集中现象. 图13则对比了不同网格和方法计算出的连接处法向面力. 可以发现, 虽然材料发生了热膨胀, 而该结构连接处的面力总和大于0, 也就是说该结构整体上倾向于收缩, 即负膨胀.
![]()
图 12 单胞受热载荷时的von Mises应力分布及变形 (变形放大100倍)Figure 12. The contours of von Mises stress and deformations of the unit cell under thermal loads (the deformation is magnified by a factor of 100)![]()
图 13 右侧连接面法向面力分布Figure 13. The normal tractions on the right surface, which are connected to other cells边界元法采用基本解作为权函数, 虽然带来了奇异性问题, 但其不仅降低了问题维度, 还使得其能更好地计算出应力集中现象. 从图13(a)中可以看出, 即使边界元采用图10中仅有518个节点的网格, 其角点上的面力值也与采用加密网格(mesh 2)的有限元结果仅相差约8%. 本文边界元法仅采用线性单元, 若采用高阶单元其结果可以进一步提高. 从图13(b)中可以看出, 采用线性单元的有限元法(mesh 3)即使采用了3万多单元, 其计算结果的精度也不如二次单元, 甚至不如仅有518个节点的边界元法的结果.
采用周期边界对该结构的热应力进行计算, 可得其等效热膨胀系数为−7.7 × 10−6. 此外, 作者还研究了外框的厚度$t$与材料的等效热膨胀系数之间的关系, 结果如图14所示. 由于中间部分材料热膨胀系数较小, 使得4个角的位移受到约束, 受到向内的约束力. 如图15(a)所示, 在$t$较小时, 由于外框的抗弯能力较弱, 结构整体向内回缩. 而随着厚度增加, 抗弯能力逐渐增强, 弯曲变形减小, 热变形引起的结构膨胀占优, 结构整体呈现膨胀趋势(如图15(b)所示).
![]()
图 14 热膨胀系数随着外框厚度的变化Figure 14. The coefficient of thermal expansion varies with the thickness of the outer frame![]()
图 15 100 K温升热载荷下单个胞元的变形 (在左侧和下侧施加对称约束其他边界自由, 变形放大100倍)Figure 15. The deformations of a unit cell under thermal loads of temperature rise of 100 K (the symmetric boundary conditions are applied on the left and bottom surface, and the others are free, and the deformation is magnified by a factor of 100)3. 结 论
边界元法可以将问题降低一个维度, 只需在边界上划分网格, 对于一些问题可有效降低分析的复杂度和计算量, 其在分析超材料等复合材料的等效热弹性性能方面具有很大潜力. 采用边界元法计算结构的等效热膨胀系数需要计算热应力, 从而产生了域积分项. 本文采用径向积分法将涉及到热膨胀的域积分转换到边界上, 避免了在内部划分网格, 使其能继续发挥边界元的优越性.
从本文方法与有限元法的对比研究可以发现, 边界元法即使采用很少的网格, 其计算出的边界面力依然可以获得较好的精度. 且对于有应力集中的问题, 边界元法仅需在边界上加密, 且面力作为未知量直接求出, 计算结果更好. 计算实例表明, 热力载荷作用下, 采用等效物性与直接有限元模拟的结果几乎一致, 可见本文方法的有效性. 采用本文方法对网状结构复合材料的热膨胀系数进行分析, 验证了本文分析复杂结构的能力. 由于本文方法只需要边界网格, 在超材料优化设计及结构参数辨识等问题中有着巨大的应用前景.
边界元法虽然可以降低问题维度, 但是其系数矩阵为满阵. 对于复杂的三维问题, 可通过进一步划分子域, 或与有限元法结合使用, 以避免求解一个过于稠密的系数矩阵构成的线性方程组.
-
![]()
图 1 某周期性复合材料的RVE, s和d代表一对周期边界的源面和目标面
Figure 1. The RVE of a periodic composite, where s and d represent the source face and target face of a couple of periodic boundary
![]()
图 4 3种平均应变下径向积分边界元法计算出的变形(缩小10倍)和应力
Figure 4. The deformations (scaled down to 10 times) and the von Mises stresses with the three averaged strains
![]()
图 5 有限元和边界元法计算出的右侧面力大小
Figure 5. The magnitudes of the traction on the right edge computed by FEM and BEM
![]()
图 8 两种加载情况下中性面竖向位移分布
Figure 8. The vertical displacements on the middle face under the two loads
![]()
图 10 边界元网格, 共计520个线单元及518个节点 (单位: mm)
Figure 10. The mesh for BEM: There are 520 line elements and 520 nodes (unit: mm)
![]()
图 11 有限元网格, 共计4914个6节点三角形单元及
10581 个节点 (单位: mm)Figure 11. The mesh for BEM: There are 4914 six-node triangular elements and
10581 nodes (unit: mm)![]()
图 12 单胞受热载荷时的von Mises应力分布及变形 (变形放大100倍)
Figure 12. The contours of von Mises stress and deformations of the unit cell under thermal loads (the deformation is magnified by a factor of 100)
![]()
图 13 右侧连接面法向面力分布
Figure 13. The normal tractions on the right surface, which are connected to other cells
![]()
图 14 热膨胀系数随着外框厚度的变化
Figure 14. The coefficient of thermal expansion varies with the thickness of the outer frame
![]()
图 15 100 K温升热载荷下单个胞元的变形 (在左侧和下侧施加对称约束其他边界自由, 变形放大100倍)
Figure 15. The deformations of a unit cell under thermal loads of temperature rise of 100 K (the symmetric boundary conditions are applied on the left and bottom surface, and the others are free, and the deformation is magnified by a factor of 100)
表 1 图2复合材料组元的物性
Table 1 The properties of the components of the composite shown in Fig. 2
Parameter Mcterial A Mcterial B Young modulus/MPa $1 \times {10^4}$ $5 \times {10^4}$ Poisson's ratio 0.25 0.3 thermal expansion coefficient/K−1 $1 \times {10^{{{ - 5}}}}$ $1 \times {10^{{{ - }}6}}$ thermal conductivity/(mW·K−1·mm−1) 20 10 
下载: 导出CSV
-
[1] 薛明德, 向志海. 大型空间结构热-动力学耦合分析方法综述. 力学学报, 2022, 54(9): 2361-2376 (Xue Mingde, Xiang Zhihai. Review of thermal-dynamical analysis methods for large space structures. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(9): 2361-2376 (in Chinese) doi: 10.6052/0459-1879-22-171 Xue Mingde, Xiang Zhihai. Review of thermal-dynamical analysis methods for large space structures. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(9): 2361-2376 (in Chinese) doi: 10.6052/0459-1879-22-171
[2] Fan LJ, Xiang ZH. Suppressing the thermally induced vibration of large-scale space structures via structural optimization. Journal of Thermal Stresses, 2015, 38(1): 1-21 doi: 10.1080/01495739.2014.950529
[3] Yuan X, Xiang Z. A thermal-flutter criterion for an open thin-walled circular cantilever beam subject to solar heating. Chinese Journal of Aeronautics, 2018, 31(9): 1902-1909 doi: 10.1016/j.cja.2018.07.002
[4] Askari M, Hutchins DA, Thomas PJ, et al. Additive manufacturing of metamaterials: A review. Additive Manufacturing, 2020, 36: 101562 doi: 10.1016/j.addma.2020.101562
[5] Abdulhameed O, Al-Ahmari A, Ameen W, et al. Additive manufacturing: Challenges, trends, and applications. Advances in Mechanical Engineering, 2019, 11(2): 1-27
[6] Lim TC. Composite metamaterial with sign-switchable coefficients of hygroscopic, thermal and pressure expansions. Advanced Composites and Hybrid Materials, 2019, 2(4): 657-669 doi: 10.1007/s42114-019-00118-3
[7] Yu C, Lin K, Jiang S, et al. Plastic and low-cost axial zero thermal expansion alloy by a natural dual-phase composite. Nature Communications, 2021, 12: 4701 doi: 10.1038/s41467-021-25036-1
[8] Chen J, Gao Q, Sanson A, et al. Tunable thermal expansion in framework materials through redox intercalation. Nature Communications, 2017, 8: 14441 doi: 10.1038/ncomms14441
[9] Prall D, Lakes RS. Properties of a chiral honeycomb with a Poisson’s ratio of — 1. International Journal of Mechanical Sciences, 1997, 39(3): 305-314
[10] Yang L, Harrysson O, West H, et al. Mechanical properties of 3D re-entrant honeycomb auxetic structures realized via additive manufacturing. International Journal of Solids and Structures, 2015, 69-70: 475-490 doi: 10.1016/j.ijsolstr.2015.05.005
[11] Wang Y, Wang L, Ma ZD, et al. Parametric analysis of a cylindrical negative Poisson’s ratio structure. Smart Materials and Structures, 2016, 25(3): 035038 doi: 10.1088/0964-1726/25/3/035038
[12] Yu B, Xu Z, Mu R, et al. Design of large-scale space lattice structure with near-zero thermal expansion metamaterials. Aerospace, 2023, 10(3): 294 doi: 10.3390/aerospace10030294
[13] Liu HY, Zheng YT, Gao XW, et al. The inverse design of auxetics using the boundary element method and the constrained conjugate gradient method. Engineering Analysis with Boundary Elements, 2024, 162: 17-27 doi: 10.1016/j.enganabound.2024.01.029
[14] Hollister SJ, Kikuchi N. A comparison of homogenization and standard mechanics analyses for periodic porous composites. Computational Mechanics, 1992, 10(2): 73-95 doi: 10.1007/BF00369853
[15] Gao YH, Xing YF. The multiscale asymptotic expansion method for three-dimensional static analyses of periodical composite structures. Composite Structures, 2017, 177: 187-195 doi: 10.1016/j.compstruct.2017.06.063
[16] Wang Y, Ye J, Liu L, et al. Evaluation the interface damage influences on mechanical properties of fiber-reinforced composites based on subdomain boundary element theory. Mechanics of Advanced Materials and Structures, 2023, 30(20): 4264-4279 doi: 10.1080/15376494.2022.2092793
[17] Yazdanparast R, Rafiee R. Determining in-plane material properties of square core cellular materials using computational homogenization technique. Engineering with Computers, 2023, 39(1): 373-386 doi: 10.1007/s00366-021-01562-w
[18] Lim TC. Anisotropic and negative thermal expansion behavior in a cellular microstructure. Journal of Materials Science, 2005, 40(12): 3275-3277 doi: 10.1007/s10853-005-2700-6
[19] Hassanzadeh-Aghdam MK, Ansari R, Darvizeh A. Micromechanical modeling of thermal expansion coefficients for unidirectional glass fiber-reinforced polyimide composites containing silica nanoparticles. Composites Part A: Applied Science and Manufacturing, 2017, 96: 110-121 doi: 10.1016/j.compositesa.2017.02.015
[20] Tian W, Qi L. Unified periodic boundary condition for homogenizing the thermo-mechanical properties of composites. Applied Mathematical Modelling, 2023, 121: 252-269 doi: 10.1016/j.apm.2023.04.024
[21] Tian W, Qi L, Chao X, et al. Periodic boundary condition and its numerical implementation algorithm for the evaluation of effective mechanical properties of the composites with complicated micro-structures. Composites Part B: Engineering, 2019, 162: 1-10 doi: 10.1016/j.compositesb.2018.10.053
[22] 谷金良. B样条边界面法及边界积分方程中的等几何方法研究. [博士论文]. 长沙: 湖南大学, 2013 (Gu Jinliang. Research of the B-spline boundary face method and the isogeometric analysis in the boundary integral equation. [PhD Thesis]. Changsha: Hunan University, 2012 (in Chinese) Gu Jinliang. Research of the B-spline boundary face method and the isogeometric analysis in the boundary integral equation. [PhD Thesis]. Changsha: Hunan University, 2012 (in Chinese)
[23] 高效伟, 刘健, 彭海峰. 集成单元边界元法及其在主动冷却热防护系统分析中的应用. 力学学报, 2016, 48(4): 994-1003 (Gao Xiaowei, Liu Jian, Peng Haifeng. Integrated unit BEM and its application in analysis of actively cooling TPS. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(4): 994-1003 (in Chinese) doi: 10.6052/0459-1879-15-437 Gao Xiaowei, Liu Jian, Peng Haifeng. Integrated unit BEM and its application in analysis of actively cooling TPS. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(4): 994-1003 (in Chinese) doi: 10.6052/0459-1879-15-437
[24] 郑昌军, 高海峰, 杜磊等. 边界元特征值分析及Burton-Miller法探究. 计算力学学报, 2016, 33(3): 335-342 (Zheng Changjun, Gao Haifeng, Du Lei, et al. BEM eigenvalue analysis and the investigation of the Burton-Miller method. Chinese Journal of Computational Mechanics, 2016, 33(3): 335-342 (in Chinese) Zheng Changjun, Gao Haifeng, Du Lei, et al. BEM eigenvalue analysis and the investigation of the Burton-Miller method. Chinese Journal of Computational Mechanics, 2016, 33(3): 335-342 (in Chinese)
[25] 余志强, 刘强, 郑昌军. 三维声学特征频率灵敏度分析的边界元法. 计算力学学报, 2023, 40(4): 634-640 (Yu Zhiqiang, Liu Qiang, Zheng Changjun. A boundary element method for the eigenfrequency sensitivity analysis of three-dimensional acoustic problems. Chinese Journal of Computational Mechanics, 2023, (4): 634-640 (in Chinese) Yu Zhiqiang, Liu Qiang, Zheng Changjun. A boundary element method for the eigenfrequency sensitivity analysis of three-dimensional acoustic problems. Chinese Journal of Computational Mechanics, 2023, (4): 634-640 (in Chinese)
[26] Zhang SQ, Xu BB, Gao ZH, et al. Topology optimization of heat transfer and elastic problems based on element differential method. Engineering Analysis with Boundary Elements, 2023, 149: 190-202 doi: 10.1016/j.enganabound.2023.01.026
[27] Zhang W, Jiang Q, Feng W, et al. Explicit structural topology optimization using boundary element method-based moving morphable void approach. International Journal for Numerical Methods in Engineering, 2021, 122(21): 6155-6179 doi: 10.1002/nme.6786
[28] Fernandes GR, Crozariol LHR, Furtado AS, et al. A 2D boundary element formulation to model the constitutive behavior of heterogeneous microstructures considering dissipative phenomena. Engineering Analysis with Boundary Elements, 2019, 99: 1-22 doi: 10.1016/j.enganabound.2018.10.018
[29] Lo Cascio M, Milazzo A, Benedetti I. A hybrid virtual–boundary element formulation for heterogeneous materials. International Journal of Mechanical Sciences, 2021, 199: 106404 doi: 10.1016/j.ijmecsci.2021.106404
[30] Yazdanparast R, Rafiee R. Developing a homogenization approach for estimation of in-plan effective elastic moduli of hexagonal honeycombs. Engineering Analysis with Boundary Elements, 2020, 117: 202-211 doi: 10.1016/j.enganabound.2020.04.012
[31] Gao XW. Boundary element analysis in thermoelasticity with and without internal cells. International Journal for Numerical Methods in Engineering, 2003, 57(7): 975-990 doi: 10.1002/nme.715
[32] 陈豪龙, 周焕林, 余波. 瞬态热传导问题的精细积分-双重互易边界元法. 应用力学学报, 2017, 34(5): 835-841 (Chen Haolong, Zhou Huanlin, Yu Bo. Precise integration DRBEM for solving transient heat conduction problems. Chinese Journal of Applied Mechanics, 2017, 34(5): 835-841 (in Chinese) Chen Haolong, Zhou Huanlin, Yu Bo. Precise integration DRBEM for solving transient heat conduction problems. Chinese Journal of Applied Mechanics, 2017, 34(5): 835-841 (in Chinese)
[33] Yang K, Wang J, Du JM, et al. Radial integration boundary element method for nonlinear heat conduction problems with temperature-dependent conductivity. International Journal of Heat and Mass Transfer, 2017, 104: 1145-1151 doi: 10.1016/j.ijheatmasstransfer.2016.09.015
[34] 高效伟, 杨恺. 功能梯度材料结构的热应力边界元分析. 力学学报, 2011, 43(1): 136-143 (Gao Xiaowei Yang Kai. Thermal stress analysis of functionally graded material structures using boundary element method. Chinese Journal of Theoretical and Applied Mechanics, 2011, 43(1): 136-143 (in Chinese) doi: 10.6052/0459-1879-2011-1-lxxb2009-585 Gao Xiaowei Yang Kai. Thermal stress analysis of functionally graded material structures using boundary element method. Chinese Journal of Theoretical and Applied Mechanics, 2011, 43(1): 136-143 (in Chinese) doi: 10.6052/0459-1879-2011-1-lxxb2009-585
[35] Hori M, Nemat-Nasser S. On two micromechanics theories for determining micro–macro relations in heterogeneous solids. Mechanics of Materials, 1999, 31(10): 667-682 doi: 10.1016/S0167-6636(99)00020-4
[36] Xia Z, Zhang Y, Ellyin F. A unified periodical boundary conditions for representative volume elements of composites and applications. International Journal of Solids and Structures, 2003, 40(8): 1907-1921 doi: 10.1016/S0020-7683(03)00024-6
[37] Zheng YT, Liu Y, Gao XW, et al. Improved hole and tube elements in BEM for elasticity problems. Engineering Analysis with Boundary Elements, 2024, 159: 17-35 doi: 10.1016/j.enganabound.2023.11.021
-
期刊类型引用(1)
1. 陈豪龙，唐欣越，王润华，周焕林，柳占立. 基于物理信息神经网络的多介质非线性瞬态热传导问题研究. 力学学报. 2025(01): 89-102 . 
本站查看
其他类型引用(0)
计量
- 文章访问数: 161
- HTML全文浏览量: 20
- PDF下载量: 35
- 被引次数: 1
