Thermo-Electric-Mechanical Coupling Bending Property and Strength Analyses of Thermoelectric Devices With the Negative Poisson's Ratio Architecture
-
摘要: 随着智能可穿戴设备的快速发展,对供电元件的续航时间、便捷性以及轻量化等提出了更高要求. 热电器件可以直接将人体新陈代谢释放的热能转换为电能,为可穿戴设备持续供电. 利用整体-局部、细观-宏观相结合的分析方法,该文研究了负Poisson比热电器件的热-电-力耦合弯曲行为及其强度失效问题. 首先,通过建立负Poisson比热电器件的均质化分析模型,获取了器件的宏观弯曲特征,并给出了应力最大的截面. 然后,建立热电蜂窝的受力分析模型,利用热力学强度理论导出了胞壁的细观强度失效临界荷载. 研究发现:热电蜂窝的应力水平随着内凹角增大呈现先减小后增加的趋势;对于负Poisson比热电蜂窝,强度失效首先发生在中间部位;对于传统的六边形热电蜂窝,端部比中间部位先发生强度破坏;热电器件发生断裂破坏时,中间和端部的临界裂纹长度近似相等,可以拟合为内凹角的指数函数.
-
关键词:
- 热电器件 /
- 负Poisson比结构 /
- 强度失效 /
- 多场耦合
Abstract: The rapid development of smart wearable devices makes a higher requirement for the power supply components, including endurance, convenience and lightweight and so on. The thermoelectric devices can directly convert the thermal energy released by human metabolism into electricity, which can be further used to continuously power the wearable devices. With the global-local and micro-macro combined analysis method, the thermo-electro-mechanical coupling bending behavior and strength failure of a negative Poisson's ratio thermoelectric device (NPR-TEG) were analyzed. Firstly, the macroscopic bending characteristics and the section with the largest stress were given through the establishment a homogeneous analysis model for the NPR-TEG. Then, the force analysis model for the thermoelectric honeycomb was built. The critical load for the strength failure of a mesoscopic cell wall was also derived with the thermodynamic strength theory. The results show that, the stress level of the thermoelectric honeycomb decreases first and then increases with the re-entrant angle. For the NPR-TEG, the strength failure occurred first in the middle part of the device. For the thermoelectric device with the traditional hexagonal honeycomb, the strength failure occurs at the end of the device rather than the middle part. With the fracture failure occurring in the thermoelectric device, the critical crack length of the middle fracture approximately equals that of the end fracture. The critical crack length could be fitted as an exponential function of the re-entrant angle. -
0. 引言
疫情防控需要和人口老龄化加剧使人们对主动健康监测的需求与日俱增,以智能可穿戴设备为代表的数字健康产业在医疗健康领域的应用场景也日趋丰富. 然而传统电池因体积庞大、使用寿命短、需要定期更换等缺陷,无法满足可穿戴设备轻型化、小型化和长期稳定的供电需求. 采用能量收集技术把人体释放的热能转换为电能,是实现智能可穿戴设备持续自供电的理想选择[1]. 因此,研究和开发自供电技术代替传统电池,对于人体生物医疗设备等不易更换电源或充电的装置具有重大意义.
热电器件可以直接将热能转化为电能,具有结构简单、无机械部件、可靠性高、使用寿命长、体积小等优点[2]. 高柔韧性和高延展性是智能可穿戴设备的基本要求,直接影响人体和设备之间的热交换效率以及穿戴舒适度. 与传统的无机热电材料相比,有机热电材料具有良好的变形能力[3]. 目前,最常用且制备工艺最成熟的柔性有机热电材料是PEDOT: PSS,已经有大量科研学者对其热、电性能进行了广泛研究[4]. 虽然柔性热电材料具有良好的变形能力,但是能量转换效率远低于纯无机热电材料,严重阻碍了有机柔性热电器件的商业化应用[5]. 因此,部分学者再次聚焦无机热电材料,期望通过结构设计提高变形能力,从而制备出高效的柔性可穿戴热电器件.
可穿戴热电器件主要承受弯曲荷载,因此必须具备优异的抗弯强度. 将无机热电材料制成薄膜结构,可以在一定程度上提高热电器件的机械性能. 例如,将蛇形热电薄膜固定在预拉伸的弹性基底上,释放拉应力后热电薄膜会自动弯曲成三维螺旋结构[6]. 该热电薄膜具有沿面内拉伸60%,沿垂直方向压缩30%的变形能力,但是经过200次单轴拉伸后电阻增大22%,并且输出功率仅有纳瓦级. 另一方面,Kong等在柔性基底上制备了平面型热电薄膜,将输出功率提高到数百纳瓦,且经过2 000次弯曲循环后电阻增大不足5%[7]. Zhao等利用碲化铋与纳米纤维复合物制备的热电薄膜,其电导率经过1 000次弯曲循环后仅下降8%,展现出稳定的力学和电学性能[8]. 在此基础上,Karthikeyan等研究了曲率半径和弯曲循环次数对热电薄膜性能的影响规律[9]. 虽然这些研究都证实了薄膜结构可以提高热电材料的变形能力,但是薄膜热电器件的缺陷也十分明显. 首先,薄膜热电器件在制备过程中,热电薄膜与基底之间微裂纹的萌生是不可避免的,随着裂纹的扩展会导致热电器件发生断裂、脱层和屈曲破坏[10]. 其次,热电薄膜的实际加热面积非常小,导致薄膜热电器件很难实现稳定的高功率输出[11]. 薄膜热电器件仅能满足微瓦级以下电子产品的供电需求,因此需要开发新型高效柔性热电器件.
传统热电器件由块体热电腿、电极、散热器和衬底组成. 热电腿是热电器件的核心元件,将其设计为螺旋结构可以提高柔韧性和延展性. 例如,将无机热电腿制备成双臂螺旋结构,可以实现100%以上的延展性[12]. Xu等以螺丝钉为模具成功制备出三维螺旋无机热电腿,可以承受100%的应变[13]. 另一方面,将电极和覆盖在电极表面的散热器设计成可拉伸结构,也可以提高无机热电器件的柔韧性. 例如,Feng等将铜电极设计成正弦蛇形结构,大幅提升了无机热电器件的变形能力,但经过200次弯曲变形后,热电器件的热电性能急剧下降[14]. Lee等利用相变材料和柔性聚合物制备散热器,使弯曲变形由柔性填充物承担,脆性的相变材料不发生弯曲形变,从而提高热电器件的变形能力[15]. 此外,科研学者们也尝试对衬底进行结构设计,从而改善热电器件的柔韧性. 例如,Fukuie等利用折纸技术把不可拉伸的衬底设计为折叠结构,从而使热电器件可以承受20%的拉伸变形[16]. Park等利用胶木制成了凹形支架,然后将无机热电块放入支架并通过柔性电极连接,从而制备出可全方位变形的柔性热电器件[17]. 但是,目前已有的无机块体热电器件仍然笨重、变形能力有限且抗弯强度较差,还不满足可穿戴设备轻质、高柔韧和高强度的使用需求. 为此,从结构上构筑一种轻质、高效、高强度的新型柔性热电器件,是智能可穿戴设备发展的迫切需求.
在单轴拉伸作用下,传统材料在垂直于荷载的方向上会发生收缩变形,而负Poisson比材料却沿垂直方向膨胀. 与传统结构相比,负Poisson比结构具有抗冲击/断裂性能好、吸能隔振以及可设计性强的优点,被广泛应用于航空航天、航海、机械自动化、生物医疗、国防军事、纺织工业等领域[18]. 因此,科研学者们开展了负Poisson比结构对热电器件性能影响的研究工作. 针对可穿戴热电器件使用过程中的弯曲断裂问题,Cui等建立了负Poisson比热电材料夹芯板的裂纹扩展分析模型,给出了热电夹芯板的疲劳寿命表达式,揭示了负Poisson比结构和人体穿戴部位对热电材料断裂性能的影响规律[19]. 考虑到多孔结构会阻碍热流的传导,Cui等建立了多孔热电材料与弹性基体的非Fourier热冲击断裂分析模型,发现与传统的正Poisson比蜂窝结构相比,采用拉胀蜂窝结构可大幅降低多孔热电材料的应力水平[20]. 此外,研究发现,当热电腿与电极之间夹角为60°,且受到20%的循环拉伸应变作用时,负Poisson比热电器件的疲劳寿命是传统块体热电器件的1 000倍[21]. 然而,负Poisson比热电器件在热-力耦合作用下的弯曲失效机理尚不明确. 为此,本文建立了负Poisson比热电器件的宏观弯曲变形分析模型,并基于此研究了热电蜂窝的细观强度失效机理. 在此基础上,给出了增强能量转换性能、柔韧性和抗弯强度的器件构筑方案. 本项目的相关研究成果可为新型负Poisson比热电器件的设计、制造及性能评价提供理论指导和技术支持.
1. 负Poisson比热电器件宏观弯曲特征
图 1(a)是一个典型的柔性可穿戴热电器件实物图,主要由有机热电薄膜(TE film)和柔性基底(PI substrate)组成[22]. 隔热衬底(fabric)放置于人体皮肤和柔性基底之间,使热电薄膜只有热端与人体皮肤接触,因此人体释放的热量无法直接沿厚度方向传入热电层,只能沿热电薄膜长度方向传导,如图 1(b)所示[22]. 本文研究的负Poisson比柔性热电器件由封装层、热电层和隔热衬底组成,热电层是多孔负Poisson比蜂窝结构,几何尺寸如图 1(c)、1(d)所示. 斜胞壁与垂直方向的夹角θ<0°表示内凹负Poisson比蜂窝,θ>0°代表传统正六边形蜂窝. 本文采用直角坐标系Oxyz,且xOy平面位于热电器件的中平面. 为了模拟任意边界位移条件,器件中平面的左右两端连接大量的线性弹簧以及扭转弹簧,其单位厚度的弹簧刚度系数分别表示为k1,k2,k3,k4和k5,k6. 与典型柔性可穿戴热电器件的工作机理相同,隔热衬底(insulating layer)使热电层只有热端与皮肤接触,阻止热电层除热端之外的其余部分吸收人体新陈代谢释放的热流量. 同时,封装层(packaging layer)可以隔绝热电层与外界环境发生热交换,使得热电层的内部热流只沿x方向传导. 热电器件热端(x =0)贴合人体皮肤并吸收人体释放的热量,其温度与皮肤温度相同,记作Th. 热电器件冷端(x=L)暴露在外界环境中,器件稳定工作时其冷端温度等于室温,记作Tc. 在温差荷载ΔT=Th-Tc作用下,热电层的冷热两端形成开路电压Δφ=λΔT并向外输出电流I=Δφ/(RI+RL),其中λ表示Seebeck系数,RI和RL分别是内电阻和外电阻.
在穿戴过程中,热电器件主要受弯曲荷载Mh和Mc. 在外力作用下,水平胞壁向外移动同时联动斜胞壁弯曲变形,从而导致器件膨胀产生负Poisson比效应. 接下来,我们将采用复合材料的均质化思想,研究负Poisson比多孔热电器件的宏观弯曲特征.
1.1 均质化模型
复合材料宏观力学表明,当蜂窝单胞的尺寸远小于整体结构时,可以采用均质化方法将多孔蜂窝结构等效为连续性结构. 针对图 1(b)所示周期性排列的负Poisson比蜂窝,可以利用Ohm定理给出热电层的等效电阻率. 假设电流由节点1输入从节点8输出,此时节点2→3→5→7组成的电阻R2-7与节点2→4→6→7组成的电阻R2-7*并联. 利用电阻的串并联关系,可以导出蜂窝单胞的总电阻R2=R1-2+R2-7R2-7*/(R2-7+R2-7*)+R7-8,其中,R1-2=R7-8=ρ20h/(2Bt),R2-7=R2-7*=ρ20(2l+h)/(Bt),ρ20表示致密热电材料的电阻率,B表示热电器件沿y方向的宽度. 由几何关系可知,蜂窝单胞的等效长度(x方向)和等效宽度(z方向)分别是h*=2h+2lsin θ和l*=2lcos θ,则等效蜂窝单胞的电阻R2*=ρ2h*/(l*B). 令R2=R2*, 给出负Poisson比热电层的等效电阻率:
$$\rho_2=\frac{2 l+3 h}{2 t} \frac{l^*}{h^*} \rho_{20}=\frac{2+3 h / l}{2 t / l} \frac{\cos \theta}{h / l+\sin \theta} \rho_{20} .$$ (1) 研究表明,多孔蜂窝结构的等效热导率可以写作[23]
$$K_2=\frac{t / l \cos \theta}{h / l+\sin \theta} K_{20}, $$ (2) 其中K20表示致密热电材料的热导率. 考虑负Poisson比蜂窝结构的剪切变形、弯曲变形以及轴向拉压变形,利用梁变形理论给出蜂窝结构的等效弹性模量E2和υ2,如下所示[23]:
$$E_2 =\frac{\cos \theta}{\frac{\cos }{h / l+\sin \theta}\left(\frac{\cos ^2 \theta}{K_{\mathrm{f}}}+\frac{\cos ^2 \theta}{K_{\mathrm{h}}}+\frac{2 h / l+\sin ^2 \theta}{K_{\mathrm{s}}}\right)},$$ (3) $$v_2 =-\frac{\sin \theta(h / l+\sin \theta)\left(\frac{1}{K_{\mathrm{f}}}+\frac{1}{K_{\mathrm{h}}}+\frac{1}{K_{\mathrm{s}}}\right)}{\frac{\cos ^2 \theta}{K_{\mathrm{f}}}+\frac{\cos ^2 \theta}{K_{\mathrm{h}}}+\frac{2 h / l+\sin ^2 \theta}{K_{\mathrm{s}}}},$$ (4) 其中,Kf=E20(t/l)3是弯曲刚度系数,Ks=E20(t/l)是拉伸刚度系数,Kh=0.5E20(t/l)/(1+υ20)是剪切刚度系数,E20和υ20分别表示致密热电材料的弹性模量和Poisson比. 至此,多孔负Poisson比热电层被转化为均质热电模型,随后可以继续使用经典连续介质力学分析热电器件的弯曲变形特征.
1.2 温度场
温度场是连接热电材料内部电流与应力的桥梁,为了分析热电器件的弯曲特征,需要先获取器件的温度场. 当无内热源的热电材料处于稳态时,其内部的电流密度和热流满足[24]J =-∇φ/ρ2-λ∇T/ρ2以及q=JλT-K2∇T,其中φ和λ分别表示电势和Seebeck系数. 同时,热电材料也满足能量守恒方程:∇J=0和∇q+ J∇φ=0. 为了提高热电器件的柔韧性,图 1中热电层的厚度H2远小于其长度L. 同时,考虑到热荷载沿热电层的轴向施加,可以将热流和电流在热电层内传导视为一维问题. 此时,由∇J=0可知等效热电层内部的电流密度是常数J=I/(l*B). 忽略材料性能的温度相关性,将电流密度和热流代入∇q+J∇φ=0,可以获得一维热电材料的稳态温度控制方程K2∇2T+J2ρ2=0. 结合热边界条件T(x =0)=Th和T(x=L)=Tc,可以直接给出图 1中热电层的温度场:
$$T=-\frac{J^{2} \rho_{2}}{2 K_{2}}\left(x^{2}-L x\right)+\frac{x}{L}\left(T_{\mathrm{c}}-T_{\mathrm{h}}\right)+T_{\mathrm{h}} . $$ (5) 1.3 热-电-力耦合应力场
在穿戴过程中,热电器件的弯曲变形由温度荷载和弯曲荷载共同产生. 由于热电器件沿y方向的挠度是均匀的,且器件的长度值大于厚度和宽度,因此可以将图 1中的负Poisson比热电器件视为夹芯复合梁. 在接下来的分析中,只需考虑x方向的轴向位移和z方向的挠度. 根据Timoshenko梁理论,热电器件的位移分量可表示为
$$u(x, z)=u_{0}+z \vartheta_{x}, w(x)=w_{0}, $$ (6) 式中,u(x, z)和w(x)分别表示热电器件的轴向位移和挠度,u0和w0是器件中平面的轴向位移和挠度,ϑx描述器件中平面绕y轴的旋转角度.
对于长厚比较大的夹芯梁,弹性力学表明垂直于中平面的应力分量远小于面内应力分量,因此分析热电器件的应力状态和变形时只需考虑其面内应力. 同时考虑温度荷载和弯曲荷载,满足平面应力的夹芯梁应力分量可以写作
$$\sigma_{x}=Q_{11}\left[\varepsilon_{x}-\alpha\left(T-T_{\mathrm{c}}\right)\right], \tau_{x z}=Q_{55} \gamma_{x z}, $$ (7) 式中,α是热膨胀系数,Q11=Ee/(1-υe2)和Q55=Ee/[2(1+υe2)]是刚度系数,σx和τxz分别表示轴向正应力和横向剪切应力. 对于封装层和隔热衬底,其温度变化可以忽略不计,即此时式(7)中α=0. 在日常生活中,人体与周围环境的温度差小于40 ℃,且热电材料的热膨胀系数约10-5. 由热弹性力学可知,可穿戴热电器件的热应变,其数量级仅有10-4. 另一方面,热电器件在佩戴过程中由机械荷载引起的应变可达0.3[21]. 由此可知,对于可穿戴热电器件而言,热应变相对于机械荷载引起的变形可以忽略. 换言之,虽然式(7)考虑了温度梯度对可穿戴热电器件的影响,但是在工程应用中可以令α=0从而简化分析过程. 在小变形假设下,热电器件的正应变和剪切应变写作
$$ \varepsilon_{x}=\frac{\partial u_{0}}{\partial x}+z \frac{\partial \vartheta_{x}}{\partial x}, \gamma_{x z}=\frac{\partial w_{0}}{\partial x}+\vartheta_{x} . $$ (8) 接下来,我们将通过能量法获取热电器件的弯曲变形特征. 根据弹性力学,热电器件中的弹性应变能为
$$U_{1}=\frac{B}{2} \int_{0}^{L} \int_{-H / 2}^{-H / 2}\left\{\sigma_{x}\left[\varepsilon_{x}-\alpha\left(T-T_{\mathrm{c}}\right)\right]+\tau_{x z} \gamma_{x z}\right\} \mathrm{d} z \mathrm{d} x, $$ (9) 式中,H=H1+H2+H3是热电器件的总厚度. 计算封装层和隔热衬底的弹性应变能时,可直接令式(9)中的热膨胀系数α=0. 当热电器件变形时,其左右两端的弹簧会同步发生变形,因此储存在弹簧中的弹性应变能可写作
$$U_{2}=\left.\frac{B}{2}\left(k_{1} u_{0}^{2}+k_{3} w_{0}^{2}+k_{5} \vartheta_{x}^{2}\right)\right|_{x=0}+\left.\frac{B}{2}\left(k_{2} u_{0}^{2}+k_{4} w_{0}^{2}+k_{6} \vartheta_{x}^{2}\right)\right|_{x=L} .$$ (10) 由两端的弯曲荷载Mh和Mc引起的外力功为
$$W=\left.M_{\mathrm{h}} \frac{\mathrm{d} w}{\mathrm{d} x}\right|_{x=0}-\left.M_{\mathrm{c}} \frac{\mathrm{d} w}{\mathrm{d} x}\right|_{x=L} .$$ (11) 根据定义,热电器件的Lagrange能量函数(Lagrangian energy function)为Π=U1+U2-W. 由式(7)—(11)可知,Π是关于位移u0, w0和ϑx的函数. 随后,我们将采用Ritz法(Ritz method)求热电器件的位移近似解. 基于Gram-Schmidt归一化的正交多项式,给出热电器件的如下位移容许函数[25]:
$$ \left\{\begin{array}{l} u_{0}(\zeta)=\sum\limits_{n=1}^{N} A_{n} \frac{\chi_{n}(\zeta)}{\sqrt{\int_{0}^{1}\left[\chi_{n}(\zeta)\right]^{2} \mathrm{d} \zeta}}, w_{0}(\zeta)=\sum\limits_{n=1}^{N} B_{n} \frac{\chi_{n}(\zeta)}{\sqrt{\int_{0}^{1}\left[X_{n}(\zeta)\right]^{2} \mathrm{d} \zeta}}, \\ \vartheta_{x}(\zeta)=\sum\limits_{n=1}^{N} C_{n} \frac{\chi_{n}(\zeta)}{\sqrt{\int_{0}^{1}\left[\chi_{n}(\zeta)\right]^{2} \mathrm{d} \zeta}}, \end{array}\right.$$ (12) 式中,ζ=x/L描述归一化的坐标,N表示正交多项式的截断项,An,Bn和Cn是待求参数,χn(ζ)是正交多项式,可表示为[25]
$$\chi_{1}(\zeta)=1, \chi_{2}(\zeta)=\left(\zeta-D_{1}\right) \chi_{1}(\zeta), \chi_{n+1}(\zeta)=\left(\zeta-D_{n}\right) \chi_{n}(\zeta)-I_{n} \chi_{n-1}(\zeta), \quad n \geqslant 2 , $$ (13) 其中,$ D_{n}=\int_{0}^{1} \zeta\left[\chi_{n}(\zeta)\right]^{2} \mathrm{d} \zeta / \int_{0}^{1}\left[\chi_{n}(\zeta)\right]^{2} \mathrm{d} \zeta, I_{n}=\int_{0}^{1} \zeta \chi_{n}(\zeta) \chi_{n-1}(\zeta) \mathrm{d} \zeta / \int_{0}^{1}\left[\chi_{n-1}(\zeta)\right]^{2} \mathrm{d} \zeta$.
将式(11)和(12)代入能量函数Π,得到一个关于待求参数An,Bn和Cn的函数表达式. 基于最小势能原理∂Π/∂Xn=0,形成大小为3N×3N的线性方程组M3N×3N{Xn}=0,其中Xn=A1,A2,…,AN,B1,B2,…,BN,C1,C2,…,CN. 结合具体的位移边界条件求解该线性方程组,可确定图 1中热电器件的位移场,简要计算流程如图 2所示.
2. 热电器件强度失效细观分析
热电材料是典型的脆性材料,当热电层的拉应力大于材料强度极限时,热电层发生强度失效. 因此,建立细观尺度下热电蜂窝的强度失效模型,并给出器件强度破坏的临界荷载,对于指导热电器件的制备和安全使用至关重要. 具体分析时,先基于第1节的宏观弯曲变形分析模型,得到热电器件的等效位移场和应力场;然后,将所得应力场施加于蜂窝单胞,建立热电蜂窝的细观强度失效分析模型.
2.1 热电蜂窝细观受力分析
建立图 3(a)所示的负Poisson比热电蜂窝受力分析示意图,其中σx和τxz分别表示均质热电模型中的正应力和剪切应力. 根据力平衡原理和叠加原理,可以获取热电蜂窝水平胞壁和斜胞壁的应力状态. 当图 3(b)所示的正应力σx单独作用于蜂窝时,任意截面上水平胞壁受到的荷载为Fσ=σxl*B. 对于单个斜胞壁,其总荷载0.5Fσ可以分解为垂直于截面的轴向力0.5Fσsin θ,以及平行于截面的剪切力0.5Fσcos θ. 当热电蜂窝仅受面内剪切荷载τxz时,建立图 3(c)所示的胞元受力分析模型. 对于水平胞壁,沿z方向由力平衡条件可得Fτ=τxzl*B. 由于水平胞壁与两个对称的斜胞壁组成受力平衡体,因此斜胞壁上的剪切力是水平胞壁的一半,即Fn=Fτ/2. 同时,斜胞壁为了满足x方向的力平衡条件,还受到荷载Fs的作用,且2Fs=τxzh*B.
当拉伸和剪切荷载同时作用于热电蜂窝时,根据叠加原理可得斜胞壁承受的轴向力和剪切力,分别为
$$\sigma_{\mathrm{in}} B t=-\left(0.5 F_{\sigma}-F_{\mathrm{s}}\right) \sin \theta+F_{\mathrm{n}} \cos \theta \Rightarrow \sigma_{\mathrm{in}}=-\frac{l^{*}}{2 t} \sigma_{x} \sin \theta+\left(\frac{h^{*}}{2 t} \sin \theta+\frac{l^{*}}{2 t} \cos \theta\right) \tau_{x z}, $$ (14) $$\tau_{\mathrm{in}} B t=-\left(0.5 F_{\sigma}-F_{\mathrm{s}}\right) \cos \theta-F_{\mathrm{n}} \sin \theta \Rightarrow \tau_{\mathrm{in}}=-\frac{l^{*}}{2 t} \sigma_{x} \cos \theta-\left(\frac{l^{*}}{2 t} \sin \theta-\frac{h^{*}}{2 t} \cos \theta\right) \tau_{x z} . $$ (15) 2.2 热电蜂窝细观断裂——热力学强度理论
在热-电-弯曲耦合荷载作用下,脆性无机热电层会断裂破坏,干扰热电器件内部热流和电流的传导,进而影响能量转换性能. 当热电蜂窝的裂纹较小时,裂纹对热传导和电传导的干扰较弱,热电器件的热流和电流仍然可以视为一维问题,即式(5)及由此导出的式(14)和(15)仍然适用. 在此基础上,建立热电蜂窝的细观断裂分析模型. 由于材料结构强度失效是非平衡且不可逆的,接下来我们将利用非平衡热力学理论研究断裂问题. 考虑熵、内能与外力功之间的关系,利用热力学定律可以给出Gibbs自由能及其微分形式[26]:
$$G=U-W-T S \Rightarrow \mathrm{d} G=\mathrm{d} U-{\rm{ \mathsf{ δ}}} W-S \mathrm{d} T-T \mathrm{d} S, $$ (16) 式中,U是系统内能,W表示外界对系统做功,S和T分别为熵和温度. 根据热力学定义,系统的熵变dS=dSe+dSi,其中dSe=δQ/T表示因系统与外界发生能量交换δQ而引发的熵变,dSi是由于系统内部产生不可逆过程(例如断裂破坏)而导致的熵变. 利用能量守恒定理(热力学第一定律):δQ=dU-δW,可以将式(16)改写为dG=-TdSi-SdT. 同时,式(16)表明,Gibbs自由能是状态变量Y={Y1, Y2, …, Yn}的函数. 因此,Gibbs自由能的增量表达式可以写作[27]
$${\rm{ \mathsf{ δ}}} G=G(Y+{\rm{ \mathsf{ δ}}} Y)-G(Y)=\frac{\partial G}{\partial Y_{i}} {\rm{ \mathsf{ δ}}} Y_{i}+\frac{1}{2} \frac{\partial^{2} G}{\partial Y_{i} \partial Y_{j}} {\rm{ \mathsf{ δ}}} Y_{i} {\rm{ \mathsf{ δ}}} Y_{j}+\cdots .$$ (17) 当系统处于平衡状态时,∂G/∂Yi=0,此时Gibbs自由能取极值. 当∂2G/(∂Yi∂Yj)>0时,Gibbs自由能的增量δG>0,虚功原理表明此时系统的平衡态是稳定的;反之,若∂2G/(∂Yi∂Yj)<0即δG<0,则系统的平衡态是不稳定的,此时系统处于强度失效的临界态.
统一考虑斜胞壁及作用于其上的荷载σin和τin,则细胞壁和荷载系组成一个封闭系统. 因此,式(16)对于图 1的热电蜂窝仍然适用. 式(16)中TS是系统的耗散能,针对本文研究的热电蜂窝断裂问题,-TS表示裂纹引起的表面能变化量. 假设斜胞壁沿壁厚方向含长度为2a的中心贯穿裂纹,如图 4所示,其中x1和y1是局部坐标.
根据定义,耗散项表示为-TS=4βat,其中β是表面能密度(单位面积对应的表面自由能). 虽然斜胞壁的长度和宽度有限,但远大于裂纹长度,因此图 4中的斜胞壁可视作无限远处受荷载的含裂纹薄板[28]. 图 1中热电蜂窝的内能等于弹性应变能,U=U0+Ua,其中U0表示无裂纹时的弹性应变能,Ua是裂纹生产导致的弹性应变能改变量. 此外,系统在恒定荷载作用下弹性应变能的改变量等于外力功的一半(即Ua=W/2),而剩余一半的能量供给裂纹扩展[28]. 基于此,给出Gibbs自由能关于裂纹长度的如下微分表达式:
$$\frac{\partial G}{\partial(2 a)}=\frac{\partial(U-W)}{\partial(2 a)}+\frac{\partial(-T S)}{\partial(2 a)}=-\frac{\partial U_{a}}{\partial(2 a)}+2 \beta t .$$ (18) 接下来,本文将以Westergard应力函数为基础,导出斜胞壁在轴向力σin和剪切力τin作用下,使长度为2a的中心裂纹闭合所需的外力功. 在局部坐标系下,取复变量z1=x1+iy1,其中i表示虚数单位. 研究指出,含裂纹板在拉伸和剪切荷载共同作用下的应力函数为[29]
$$f(Z)=\operatorname{Re} \bar{\bar{Z}}+y_{1} \operatorname{Im} \bar{Z}+a_{1} y_{1} \operatorname{Re} \bar{Z}+a_{2} y_{1}^{2}, $$ (19) 式中,$ Z=\sigma_{\mathrm{in}} z_{1} / \sqrt{z_{1}^{2}-a^{2}}$是复变解析函数,Z和$ \bar{\bar{Z}}$分别表示Z的一次积分和二次积分. 根据应力分量与应力函数之间的关系,可以给出[29]
$$\left\{\begin{array}{l} \sigma_{x_{1}}=\frac{\partial^{2} f(Z)}{\partial y_{1}^{2}}=\operatorname{Re} Z-y_{1} \operatorname{Im} Z^{\prime}-2 a_{1} \operatorname{Im} Z-a_{1} y_{1} \operatorname{Re} Z^{\prime}+2 a_{2}, \\ \sigma_{y_{1}}=\frac{\partial^{2} f(Z)}{\partial x_{1}^{2}}=\operatorname{Re} Z+y_{1} \operatorname{Im} Z^{\prime}-a_{1} y_{1} \operatorname{Re} Z^{\prime}, \\ \tau_{x_{y_{1}}}=-\frac{\partial^{2} f(Z)}{\partial x_{1} \partial y_{1}}=-y_{1} \operatorname{Re} Z^{\prime}-a_{1} \operatorname{Re} Z+a_{1} y \operatorname{Im} Z^{\prime}, \end{array}\right.$$ (20) 式中,上角标“′”表示导数. 结合力边界条件σx1(x1→∞,y1=0)=σin以及τx1y1(x1→∞,y1=0)=τin可以确定未知参数a1=-τin/σin,a2=-σin/2. 然后,利用弹性力学的几何方程和平面应力物理方程,可以确定裂纹面的水平位移和垂直位移,如下所示:
$$ \begin{align*} & u\left(z_{1}\right)=\int \varepsilon_{x_{1}} \mathrm{d} x_{1}=\frac{1}{E_{20}} \int\left(\sigma_{x_{1}}-v_{20} \sigma_{y_{1}}\right) \mathrm{d} x_{1}= \\ & \quad \quad\frac{1}{E_{20}} \int\left[\left(1-v_{20}\right) \operatorname{Re} Z-\left(1+v_{20}\right) y_{1} \operatorname{Im} Z^{\prime}+\left(1+v_{20}\right) \frac{\tau_{\text {in }}}{\sigma_{\text {in }}} y_{1} \operatorname{Re} Z^{\prime}+\frac{2 \tau_{\text {in }}}{\sigma_{\text {in }}} \operatorname{Im} Z-\sigma_{\text {in }}\right] \mathrm{d} x_{1}= \\ & \quad\quad \frac{1}{E_{20}}\left[\left(1-v_{20}\right) \operatorname{Re} \bar{Z}-\left(1+v_{20}\right) y_{1} \operatorname{Im} Z+\left(1+v_{20}\right) \frac{\tau_{\text {in }}}{\sigma_{\text {in }}} y_{1} \operatorname{Re} Z+\frac{2 \tau_{\text {in }}}{\sigma_{\text {in }}} \operatorname{Im} \bar{Z}-\sigma_{\text {in }} x_{1}\right], \end{align*}$$ (21) $$ \begin{align*} &v\left(z_{1}\right) = \int \varepsilon_{y_{1}} \mathrm{d} x_{1}=\frac{1}{E_{20}} \int\left(\sigma_{y_{1}}-v_{20} \sigma_{x_{1}}\right) \mathrm{d} y_{1}= \\ & \quad\quad\frac{1}{E_{20}} \int\left[\left(1-v_{20}\right) \operatorname{Re} Z+\left(1+v_{20}\right) y_{1} \operatorname{Im} Z^{\prime}-\left(1+v_{20}\right) \frac{\tau_{\text {in }}}{\sigma_{\text {in }}} y_{1} \operatorname{Re} Z^{\prime}-\right. \\ &\quad\quad \left. v_{20} \frac{2 \tau_{\text {in }}}{\sigma_{\text {in }}} \operatorname{Im} Z+v_{20} \sigma_{\text {in }}\right] \mathrm{d} y_{1}=\\ &\quad\quad \frac{1-v_{20}}{E_{20}} \int \operatorname{Re} Z \mathrm{d} y_{1}+\frac{1+v_{20}}{E_{20}} \int y_{1} \operatorname{Im} Z^{\prime} \mathrm{d} y_{1}- \\ & \quad\quad\frac{1+v_{20}}{E_{20}} \frac{\tau_{\text {in }}}{\sigma_{\text {in }}} \int y_{1} \operatorname{Re} Z^{\prime} \mathrm{d} y_{1}-v_{20} \frac{2 \tau_{\text {in }}}{\sigma_{\text {in }}} \int \operatorname{Im} Z \mathrm{d} y_{1}+\frac{v_{20} \sigma_{\text {in }}}{E_{20}} y_{1}= \\ &\quad\quad \frac{1}{E_{20}}\left[2 \operatorname{Im} \bar{Z}-\left(1-v_{20}\right) \frac{\tau_{\text {in }}}{\sigma_{\text {in }}} \operatorname{Re} \bar{Z}-\left(1+v_{20}\right) \frac{\tau_{\text {in }}}{\sigma_{\text {in }}} y_{1} \operatorname{Im} Z-\left(1+v_{20}\right) y_{1} \operatorname{Re} Z-v_{20} \sigma_{\text {in }} y_{1}\right] . \end{align*} $$ (22) 裂纹闭合时y1=0,即z1=x1且除裂纹尖端外|x1|<a,由此可得
$$ Z=\sigma_{\mathrm{in}} x_{1} / \sqrt{x_{1}^{2}-a^{2}}=-\mathrm{i} \sigma_{\mathrm{in}} x_{1} / \sqrt{a^{2}-x_{1}^{2}}, \bar{Z}=\int Z \mathrm{d} z_{1}=\sigma_{\mathrm{in}} \sqrt{z_{1}^{2}-a^{2}}=\mathrm{i} \sigma_{\mathrm{in}} \sqrt{a^{2}-x_{1}^{2}} . $$ 此时,式(21)和(22)简化为$ u\left(x_{1}\right)=(2 \tau_{\text {in }} \sqrt{a^{2}-x_{1}^{2}}-\sigma_{\mathrm{in}} x_{1}) / E_{20}$以及$ v\left(x_{1}\right)=2 \sigma_{\text {in }} \sqrt{a^{2}-x_{1}^{2}} / E_{20}$. 对于线弹性材料,假设应力σ和τ从0增大到σin和τin时,位移由0线性增加到u(x1)和v(x1)[28]. 在裂纹面的任意位置取长度为dx1的微裂纹,在其闭合过程中弹性应变能的改变量为
$$\mathrm{d} U_{a}=2 \mathrm{d} x_{1} t\left(\int_{0}^{v} \sigma \mathrm{d} \xi+\int_{0}^{u} \tau \mathrm{d} \xi\right)=2 \mathrm{d} x_{1} t\left[\frac{1}{2} \sigma_{\text {in }} v\left(x_{1}\right)+\frac{1}{2} \tau_{\mathrm{in}} u\left(x_{1}\right)\right].$$ (23) 当长度为2a的裂纹完全闭合时,斜胞壁弹性应变能的改变量为
$$ \begin{align*} U_{a}= & \int_{-a}^{a} \mathrm{d} U_{a}=t \int_{-a}^{a} \sigma_{\text {in }} v\left(x_{1}\right) \mathrm{d} x_{1}+t \int_{-a}^{a} \tau_{\text {in }} u\left(x_{1}\right) \mathrm{d} x_{1}= \\ & \sigma_{\text {in }}^{2} \frac{2 t}{E_{20}} \int_{-a}^{a} \sqrt{a^{2}-x_{1}^{2}} \mathrm{d} x_{1}+\frac{t}{E_{20}} \tau_{\text {in }} \int_{-a}^{a}\left(2 \tau_{\text {in }} \sqrt{a^{2}-x_{1}^{2}}-\sigma_{\text {in }} x_{1}\right) \mathrm{d} x_{1}= \\ & \left(\sigma_{\text {in }}^{2}+\tau_{\text {in }}^{2}\right) \frac{2 t}{E_{20}} \int_{-a}^{a} \sqrt{a^{2}-x_{1}^{2}} \mathrm{d} x_{1}= \\ & \left(\sigma_{\text {in }}^{2}+\tau_{\text {in }}^{2}\right) \frac{2 t}{E_{20}}\left[\frac{x_{1} \sqrt{a^{2}-x_{1}^{2}}}{2}+\frac{a^{2}}{2} \arcsin \left(\frac{x_{1}}{a}\right)\right]_{-a}^{a}=\frac{{\rm{ \mathsf{ π}}} a^{2} t}{E_{20}}\left(\sigma_{\text {in }}^{2}+\tau_{\text {in }}^{2}\right) . \end{align*}$$ (24) 将式(24)代入式(18),给出裂纹扩展的临界条件,如下所示:
$$\frac{\partial G}{\partial(2 a)}=-\frac{{\rm{ \mathsf{ π}}} a t}{E_{20}}\left(\sigma_{\mathrm{in}}^{2}+\tau_{\mathrm{in}}^{2}\right)+2 \beta t=0 \Rightarrow a_{\mathrm{c}}=\frac{2 \beta E_{20}}{{\rm{ \mathsf{ π}}}\left(\sigma_{\mathrm{in}}^{2}+\tau_{\mathrm{in}}^{2}\right)} . $$ (25) 式(7)、(14)和(15)表明,应力σin和τin是荷载(Mh,Mc,Th,Tc)以及几何结构尺寸(h,l,t,θ)的函数. 结合式(25),可以确定不同荷载及结构尺寸对应的临界裂纹长度,进而指导负Poisson比热电器件的结构优化和安全性设计.
3. 数值算例及讨论
碲化铋(Bi2Te3,bismuth telluride)是商业化热电材料,其在室温下具有高热电转换性能,因此本算例中采用碲化铋模拟热电层. 封装层和隔热衬底采用高分子聚合物PDMS(polydimethylsiloxane)薄膜,其材料属性和几何尺寸见表 1. 热电器件冷热两端的温度分别等于室温和人体温度,即Tc=298 K,Th=310 K. 在本算例中,我们假设沿z方向只有一个热电蜂窝,即l*+t=H2. 为了模拟穿戴过程,将热电器件视为不可移动简支梁(k1=k2=k3=k4=1×1012 N/m2,k5=k6=0 N/m2). 左右两端施加的参考弯矩,定义为Mh=Mc=0.1E2I/r=0.1E2B(H1+H2+H3)3/(12r),其中r表示热电器件的曲率半径. 为了获得最大输出功率,假设外电阻等于热电器件的内电阻,即RL=RI=ρ2L/(BH2).
表 1 可穿戴热电器件材料属性和几何尺寸Table 1. Material properties and geometric dimensions of the wearable thermoelectric deviceparameter value Bi2Te3 K20=1.1 W/(m·K)
E20=63.3 GPaρ20=9.4×10-6 Ω·m
υ20=0.23λ=194 μV/K
α=8.98×10-6 KPDMS E1=E3=10 MPa υ1=υ3=0.495 dimension L=43 mm B=39 mm H1=H3=0.5 mm H2=2 mm h/l=2 t=0.5 mm 3.1 有限元模拟
针对斜胞壁与z方向夹角θ=-π/12、曲率半径r=4.2 cm(成年人手腕半径的平均值)的情况,利用ABAQUS对热电器件进行有限元模拟,并与理论结果对比. 有限元模拟简要步骤如下:(a) 根据表 1中几何尺寸,建立热电器件的等效三维有限元模型. (b) 设置热电层的等效材料参数:K2=0.393 W/(m·K), r2=3.24×10-5 Ω·m,E2=8.53 GPa,υ2=0.307;封装层和隔热衬底的材料参数如表 1所示. (c) 根据不可移动简支条件,固定热电器件左右两端水平方向和垂直方向的位移,并在结构两端施加弯矩荷载Mh=Mc=1 782.54 N·mm. (d) 在预定义场中设置结构的初始温度为298 K,并沿长度方向施加温度荷载T=-5.717 2×10-5 x2-276.87×10-3x+310. (e) 利用线性六面体单元C3D8对几何模型进行网格划分,获得单元总数119 560、单元节点数127 224的有限元模型.
针对热电器件的上表面z=(H1+H2+H3)/2,图 5(a)对比了挠度的理论值与有限元模拟结果,其中插图是挠度的有限元云图. 在理论计算时,线性弹簧的刚度系数设置为1×1012 N/m2,虽然弹簧变形很小但仍然有部分外力功转化为弹簧势能. 然而建立有限元模型时,器件左右两端水平方向和垂直方向的位移被完全固定,即外力功全部被转化为器件的弹性势能. 因此,器件挠度的理论值必然小于有限元结果,且二者之间的误差无法消除. 对于最大挠度,理论值和有限元结果的相对误差为(1.81-1.68)/1.81=7%,这也说明了采用弹性边界条件模拟实际约束的重要性. 由图 5(a)也可以发现,虽然热-力荷载沿轴线方向是非对称的,但器件的挠度是对称变形的,这是因为热荷载远小于对称施加的弯矩荷载.
针对器件中间位置x=L/2,图 5(b)给出了热电器件的轴向力沿厚度方向的分布规律. 由于封装层、隔热衬底(PDMS)的刚度远小于热电层(Bi2Te3),因此器件发生弯曲变形时,封装层和隔热衬底的轴向应力相比于热电层近似为零. 根据弹性力学,热电层的上表面(z=H2/2)和下表面(z=-H2/2)有最大拉伸应力和压缩应力. 理论分析与有限元模拟的最大相对误差为(66.2-62.1)/66.2=6%,在工程应用的可接受范围内. 为了研究弹性边界对最大轴向应力的影响,表 2给出了轴向应力最大值随弹簧刚度系数的关系.
表 2 轴向应力最大值随弹簧刚度系数变化规律Table 2. Material properties and geometric dimensions of the wearable thermoelectric deviceparameter value k/(N/m2) 1×106 1×109 1×1012 1×1015 |σ|max/MPa 61.6 61.9 62.1 62.1 由表 2可知,随着弹簧刚度的系数增大,理论模型和有限元模拟结果和相对误差逐渐减小. 当刚度系数取1×1012 N/m2和1×1015 N/m2时,轴向应力最大值相同,说明刚度系数取1×1012 N/m2计算结果已经趋于稳定. 因此,本文数值部分刚度系数取1×1012 N/m2.
考虑到Bi2Te3是脆性材料,接下来将基于热电层上表面分析器件的强度失效问题,并给出热电器件强度失效的临界热-力耦合荷载.
3.2 内凹角对应力状态和强度失效的影响
针对不同大小的内凹角,图 6(a)给出了热电层上表面的轴向应力和剪切应力随内凹角的变化曲线. 当内凹角绝对值|θ|的增大相同时,负Poisson比热电蜂窝的应力水平低于传统六边形热电蜂窝. 从物理角度分析,负Poisson比结构的等效弹性模量低于传统蜂窝结构,如图 6(b)所示. 因此,使相同内凹角的负Poisson比热电蜂窝产生同等大小的曲率半径,所需要的外力矩必然小于传统蜂窝结构. 此外,研究发现负Poisson比热电蜂窝的应力水平随着内凹角增大呈现先减小后增加的趋势,如图 6(a)所示. 根据公式Mh=Mc=0.1E2I/r以及图 6(b)可知,热电蜂窝的等效弹性模量E2以及施加的外力矩M随内凹角绝对值的增大而降低,由此导致σx和τxz随|θ|的增加而减小. 当内凹角θ<-45°时,内凹角绝对值|θ|的增大会导致热电层端部的轴向应力σx增加. 这是因为随着|θ|的增加,热电层的相对密度χ=t/l*·(h*/l*+2)/[2cos θ·(h*/l*+sin θ)]逐渐提高. 由图 6(a)可以看出,虽然热电层中间部位的轴向力大于其端部,但其中间部位x=L/2的剪切应力为零. 弹性力学表明剪切应力关于对称中心呈现中心对称,因此剪切应力在热电层中间位置等于零. 此外,式(7)表明热电层上表面的轴向应力σx由总应变εx和热应变α(T-Tc)共同决定,其中εx是位置坐标x的函数,温度T为坐标x和内凹角θ的函数. 当θ=-52°或32°时,εx - α(T-Tc)在x=0和x=0.5L处恰好相等,因此图 6(a)中出现交点.
虽然热电层中间部位的剪切应力为零,但是该位置的轴向应力较大. 因此需要进一步分析斜胞壁的轴向应力和剪切应力,才能准确判断中间还是端部首先发生强度失效. 对于不同的内凹角,图 7(a)给出了中间(x=0.5L)以及端部(x=0)斜胞壁的轴向应力和剪切应力分布曲线. 从图中可以看出:对于负Poisson比热电蜂窝,中间部位的应力水平高于端部(σin2+τin2);对于传统六边形热电蜂窝,端部的应力水平反而高于中间部位. 因为轴向应力和剪切应力的数量级相同,由断裂力学可知斜胞壁发生Ⅰ型和Ⅱ型共存的混合断裂模式,相应的临界裂纹长度如图 7(b)所示. 从图中可以看出:负Poisson比热电蜂窝的临界裂纹长度大于六边形热电蜂窝,证实负Poisson比结构具有更好抗断裂性能;中间部位和端部发生断裂破坏时,二者的临界裂纹长度近似相等,可以拟合为内凹角的指数函数ac/β=2.265×10-10exp(-0.463 5θ)+4.198exp(-0.016 89θ), 单位μm/N. 该拟合公式适用于弯曲半径等于成人手臂、室温下的碲化铋热电器件,且蜂窝单胞需要满足h/l=2. 与图 6(a)相似,斜胞壁的轴向力σin由热电层等效应力σx和τxz以及内凹角θ决定,如式(14)所示. 当内凹及等于30°时,x=0和x=0.5L处斜胞壁的轴向力相同,因此图 7(a)中存在交点.
4. 结论
利用整体-局部、细观-宏观相结合的分析方法,本文研究了负Poisson比热电器件的热-电-力耦合弯曲行为及其强度失效问题. 研究发现:热电蜂窝的应力水平随着内凹角增大呈现先减小后增加的趋势;对于负Poisson比热电蜂窝,强度失效首先发生在中间部位;对于传统六边形热电蜂窝,端部比中间部位先发生强度破坏;热电器件发生断裂破坏时,中间和端部的临界裂纹长度近似相等,可以拟合为内凹角的指数函数.
-
表 1 可穿戴热电器件材料属性和几何尺寸
Table 1. Material properties and geometric dimensions of the wearable thermoelectric device
parameter value Bi2Te3 K20=1.1 W/(m·K)
E20=63.3 GPaρ20=9.4×10-6 Ω·m
υ20=0.23λ=194 μV/K
α=8.98×10-6 KPDMS E1=E3=10 MPa υ1=υ3=0.495 dimension L=43 mm B=39 mm H1=H3=0.5 mm H2=2 mm h/l=2 t=0.5 mm 表 2 轴向应力最大值随弹簧刚度系数变化规律
Table 2. Material properties and geometric dimensions of the wearable thermoelectric device
parameter value k/(N/m2) 1×106 1×109 1×1012 1×1015 |σ|max/MPa 61.6 61.9 62.1 62.1 -
[1] JIANG F, ZHOU X, LV J, et al. Stretchable, breathable, and stable lead-free perovskite/polymer nanofiber composite for hybrid triboelectric and piezoelectric energy harvesting[J]. Advanced Materials, 2022, 34(17): 2200042. doi: 10.1002/adma.202200042 [2] 崔有江, 王保林, 王开发. 多孔泡沫热电器件断裂及其对能量转化性能的影响规律研究[J]. 应用数学和力学, 2023, 44(11): 1291-1298. doi: 10.21656/1000-0887.440147 CUI Youjiang, WANG Baolin, WANG Kaifa. Evaluation of fracture and its effects on energy conversion performance of porous foam thermoelectric generators[J]. Applied Mathematics and Mechanics, 2023, 44(11): 1291-1298. (in Chinese) doi: 10.21656/1000-0887.440147 [3] SHI X L, SUN S, WU T, et al. Weavable thermoelectrics: advances, controversies, and future developments[J]. Materials Futures, 2024, 3(1): 012103. doi: 10.1088/2752-5724/ad0ca9 [4] YANG Y, DENG H, FU Q. Recent progress on PEDOT: PSS based polymer blends and composites for flexible electronics and thermoelectric devices[J]. Materials Chemistry Frontiers, 2020, 4(11): 3130-3152. doi: 10.1039/D0QM00308E [5] SUN T T, ZHOU B Y, ZHENG Q, et al. Stretchable fabric generates electric power from woven thermoelectric fibers[J]. Nature Communications, 2020, 11(1): 572. doi: 10.1038/s41467-020-14399-6 [6] NAN K W, KANG S D, LI K, et al. Compliant and stretchable thermoelectric coils for energy harvesting in miniature flexible devices[J]. Science Advances, 2018, 4(11): eaau5849. doi: 10.1126/sciadv.aau5849 [7] KONG D Y, ZHU W, GUO Z P, et al. High-performance flexible Bi2Te3 films based wearable thermoelectric generator for energy harvesting[J]. Energy, 2019, 175: 292-299. doi: 10.1016/j.energy.2019.03.060 [8] ZHAO X, ZHAO C S, JIANG Y F, et al. Flexible cellulose nanofiber/Bi2Te3 composite film for wearable thermoelectric devices[J]. Journal of Power Sources, 2020, 479: 229044. doi: 10.1016/j.jpowsour.2020.229044 [9] KARTHIKEYAN V, SURJADI J U, WONG J C K, et al. Wearable and flexible thin film thermoelectric module for multi-scale energy harvesting[J]. Journal of Power Sources, 2020, 455: 227983. doi: 10.1016/j.jpowsour.2020.227983 [10] CUI Y J, WANG B L, WANG P. Analysis of thermally induced delamination and buckling of thin-film thermoelectric generators made up of pn-junctions[J]. International Journal of Mechanical Sciences, 2018, 149: 393-401. doi: 10.1016/j.ijmecsci.2017.10.049 [11] KOGO G, XIAO B, DANQUAH S, et al. A thin film efficient pn-junction thermoelectric device fabricated by self-align shadow mask[J]. Scientific Reports, 2020, 10(1): 1067. doi: 10.1038/s41598-020-57991-y [12] ROJAS J P, SINGH D, CONCHOUSO D, et al. Stretchable helical architecture inorganic-organic hetero thermoelectric generator[J]. Nano Energy, 2016, 30: 691-699. doi: 10.1016/j.nanoen.2016.10.054 [13] XU X J, ZUO Y, CAI S, et al. Three-dimensional helical inorganic thermoelectric generators and photodetectors for stretchable and wearable electronic devices[J]. Journal of Materials Chemistry C, 2018, 6(18): 4866-4872. doi: 10.1039/C8TC01183D [14] FENG R, TANG F, ZHANG N, et al. Flexible, high-power density, wearable thermoelectric nanogenerator and self-powered temperature sensor[J]. ACS Applied Materials & Interfaces, 2019, 11(42): 38616-38624. [15] LEE G, KIM C S, KIM S, et al. Flexible heatsink based on a phase-change material for a wearable thermoelectric generator[J]. Energy, 2019, 179: 12-18. doi: 10.1016/j.energy.2019.05.018 [16] FUKUIE K, IWATA Y, IWASE E. Design of substrate stretchability using origami-like folding deformation for flexible thermoelectric generator[J]. Micromachines, 2018, 9(7): 315. doi: 10.3390/mi9070315 [17] PARK H, LEE D, KIM D, et al. High power output from body heat harvesting based on flexible thermoelectric system with low thermal contact resistance[J]. Journal of Physics D: Applied Physics, 2018, 51(36): 365501. doi: 10.1088/1361-6463/aad270 [18] 周世奇, 侯秀慧, 邓子辰. 一般宏观应力状态下凹角蜂窝结构的屈曲性能分析[J]. 应用数学和力学, 2023, 44(1): 12-24. doi: 10.21656/1000-0887.430202ZHOU Shiqi, HOU Xiuhui, DENG Zichen. Buckling analysis of re-entrant honeycomb structures under general macroscopic stress states[J]. Applied Mathematics and Mechanics, 2023, 44(1): 12-24. (in Chinese) doi: 10.21656/1000-0887.430202 [19] CUI Y J, LIU C, WANG K F, et al. Effect of negative Poisson's ratio architecture on fatigue life and output power of flexible wearable thermoelectric generators[J]. Engineering Fracture Mechanics, 2023, 281: 109142. doi: 10.1016/j.engfracmech.2023.109142 [20] CUI Y J, LI W J, WANG K F, et al. Thermal shock fracture of honeycomb-based porous thermoelectric materials with non-Fourier heat conduction[J]. Ceramics International, 2024, 50(1): 2151-2161. doi: 10.1016/j.ceramint.2023.10.328 [21] CUI Y J, WANG B L, WANG K F, et al. An analytical model to evaluate influence of negative Poisson's ratio architecture on fatigue life and energy conversion performance of wearable thermoelectric generator[J]. International Journal of Solids and Structures, 2022, 258: 112000. doi: 10.1016/j.ijsolstr.2022.112000 [22] WE J H, KIM S J, CHO B J. Hybrid composite of screen-printed inorganic thermoelectric film and organic conducting polymer for flexible thermoelectric power generator[J]. Energy, 2014, 73: 506-512. doi: 10.1016/j.energy.2014.06.047 [23] HU J S, WANG B L, HIRAKATA H, et al. Interfacial thermal damage and fatigue between auxetic honeycomb sandwich and underneath substrate[J]. International Journal of Solids and Structures, 2023, 279: 112364. doi: 10.1016/j.ijsolstr.2023.112364 [24] PENG J, LI D K, HUANG Z X, et al. Interfacial behavior of a thermoelectric film bonded to a graded substrate[J]. Applied Mathematics and Mechanics(English Edition), 2023, 44(11): 1853-1870. doi: 10.1007/s10483-023-3045-8 [25] MIAO X Y, LI C F, PAN Y C. Research on the dynamic characteristics of rotating metal-ceramic matrix DFG-CNTRC thin laminated shell with arbitrary boundary conditions[J]. Thin-Walled Structures, 2022, 179: 109475. doi: 10.1016/j.tws.2022.109475 [26] 王彪. 热力学强度理论[J]. 力学进展, 2023, 53(3): 693-712.WANG Biao. Thermodynamic strength theory (TST)[J]. Advances in Mechanics, 2023, 53(3): 693-712. (in Chinese) [27] WANG B. The principle of virtual energy for predicting the strength of material structures[J]. Engineering Fracture Mechanics, 2024, 300: 109997. doi: 10.1016/j.engfracmech.2024.109997 [28] JANSSEN M, ZUIDEMA J, WANHILL R J H. Fracture Mechanics[M]. 2nd ed. London: Spon Press, 2004: 83-106. [29] 蒋玉川, 蒲淳清. 用Westergaard应力函数求解Ⅰ-Ⅱ复合型平面裂纹问题的研讨[J]. 力学与实践, 2020, 42(4): 504-507.JIANG Yuchuan, PU Chunqing. The problem of Ⅰ-Ⅱ combined plane crack solved with Westergaard stress function[J]. Mechanics in Engineering, 2020, 42(4): 504-507. (in Chinese) -