摘要:利用有限变形理论, 以无限长压电圆杆为研究对象, 考虑了在横向惯性、等效泊松比效应以及在热电弹耦合共同作用下, 基于Hamilton原理, 并引入Euler方程推导出压电圆杆的纵向波动方程. 采用Jacobi椭圆函数展开法, 求解压电圆杆的波动方程和对应的解. 最后, 通过Matlab软件得到不同波速比下的色散曲线, 以及温度场对压电圆杆的波形、波幅和波数的影响曲线. 数值分析结果表明: 随着温度的升高, 波速逐渐降低, 温度场的改变可影响和控制孤立波的传播特性.
关键词: 压电圆杆/
Hamilton原理/
温度效应/
孤波

English Abstract

全文HTML

--> --> -->

随着现代科学技术的飞速发展, 压电材料因其独特的性能而被广泛应用^[1]. 在航空航天、智能结构等领域, 装备工作环境恶劣, 差异大, 特别是温度变化比较大, 对实施控制有着非常大的影响. 因此在实际工况条件的精确建模过程中, 需要考虑电场和温度场的耦合作用. 而杆作为常用的构件, 也吸引了不少****的关注. 刘延柱等^[2]和He等^[3]利用广义热-弹耦合理论求解并研究了半无限压电杆的边值问题. 目前对于波动的问题, 主要采用多尺度法、齐次平衡法等展开了对压电圆杆波动的研究^[4,5]. 冯依虎^[6]利用泛函分析变分迭代的方法求出各次孤子波近似解, 进而研究了强非线性波动方程的行波解. Guo等^[7]利用Hamilton变分原理, 根据有限变形理论的拉格朗日描述, 推导出弹性细杆的非线性波动方程, 利用多尺度法得到了稳定的行波解. 李敏等^[8]通过对薛定谔方程的相平面分析, 约化得到其同异宿轨道, 并在相应条件下得到方程的明、暗孤立波解.
然而上述求解方法具有一定的局限性, 只能求出波动方程的冲击波解、孤波解和初等函数的周期解^[9-11]. 但采用Jacobi椭圆函数法便可求出波动方程的广义周期解和对应的孤立波解^[12,13], 刘志芳和张善元^[14-16]利用Jacobi椭圆函数展开法求得了无限长圆杆的非线性扭转波解、孤波解以及非圆截面杆的行波解、周期解.
由于压电结构在工程领域具有广泛应用, 压电介质中波的传播吸引了很多****的关注. 邓庆田等^[17]用逐步近似法对位移函数进行了假设并通过变动参数法求解, 对压电层和圆杆中的几何非线性波进行了研究. Seadawy和Manafian^[18]通过扩展尝试方程法和积分方法, 推导出了磁电弹圆杆纵波方程的暗孤子、亮孤子、孤波、周期孤波、有理函数解和椭圆函数解等不同形式的新的显式精确解. Baskonus等^[19,20]利用Sin-gordon展开法对磁电弹性圆杆的纵波方程的解析解进行了研究, 得到了更多新的解析解, 给出了所有解的数值模拟, 很好地解释了一些实际物理问题. Wang^[21]研究了压电耦合圆柱壳结构中波的传播, 从理论上得到了双模壳模型的频散曲线, 推导出波数极限情况下的截止频率和相速度. Xue和Pan^[22]考虑几何非线性以及横向泊松比引起的弥散效应对无限长磁电弹圆杆进行了研究, 建立了纵波方程, 并通过Jacobi椭圆函数法对其进行了求解. Samsonov^[23]首先报道了杆中存在孤波的实验研究, 利用聚苯乙烯的弹性介质, 设计了一套通过光学原理构造和记录孤波的实验方法, 用全息照相法记录下了孤波轨迹, 从事实上证实了弹性固体中孤波的存在. 2013年, Toffoli等^[24]在一个大的定向波池实验中探究了平面波对斜摄动的调制和有限水深下异常波动的产生, 并采集了实验数据, 对流体中的波动进行了实验探究.
综上所述, 由于非线性波的激发和观测是非常困难的, 导致实验上工作比较少. 所以本文采用建立模型和数值分析的方法, 研究不同温度场下的非线性波动问题. 通过Hamilton变分原理, 引入Euler方程, 采用Jacobi椭圆函数展开法, 推导出压电圆杆的波动方程和对应的解, 并讨论温度变化对压电介质的波形、波幅以及波速等的影响.

图1为无限长压电圆杆示意图, 建立圆柱坐标系$\left( {r, \theta, z} \right)$, z是沿着杆的轴向, 也是波传播的方向. θ = [0, 2π], 0 ≤ r ≤ R, 其中R为压电圆杆半径. 为了研究方便, 假设: 1)变形后, 板中初始状态垂直于中心平面的点仍然垂直于中心平面; 2)杆的截面是轴对称的, 即${U_\theta } = 0$和${\partial / {\partial \theta = {\rm{0}}}}$, 其中${U_\theta }$为θ方向位移; 3)考虑泊松比效应, 纵向位移U和径向位移${U_r}$之间满足${U_r} = {v_{{\rm {eff}}}}r{{\partial U} / {\partial z}}$, 其中${v_{{\rm {eff}}}}$是有效泊松比.
图 1 压电圆杆示意图
Figure1. Schematic diagram of piezoelectric rod.

对于横观各向同性的压电材料圆杆, 在考虑温度效应时的本构方程如下^[25]:

$\left\{ \begin{aligned}& {\sigma _r} = {c_{11}}{\varepsilon _r} + {c_{12}}{\varepsilon _\theta } + {c_{13}}{\varepsilon _z} - {e_{31}}{E_z} - {\lambda _{11}}\varTheta, \\& {\sigma _\theta } = {c_{12}}{\varepsilon _r} + {c_{11}}{\varepsilon _\theta } + {c_{13}}{\varepsilon _z} - {e_{31}}{E_z} - {\lambda _{11}}\varTheta, \\ &{\sigma _z} = {c_{13}}{\varepsilon _r} + {c_{13}}{\varepsilon _\theta } + {c_{33}}{\varepsilon _z} - {e_{33}}{E_z} - {\lambda _{33}}\varTheta, \\ &{\tau _{rz}} = {c_{44}}{\gamma _{rz}} - {e_{15}}{E_r}, \\ &{\tau _{\theta z}} = {c_{44}}{\gamma_{\theta z}} - {e_{15}}{E_\theta }, \\ &{\tau _{r\theta }} = {c_{66}}{\gamma _{r\theta }}, \\ &{D_r} = {e_{15}}{\gamma _{rz}} + {\varepsilon _{11}}{E_r} + {d_1}\varTheta, \\& {D_\theta } = {e_{15}}{\gamma _{\theta z}} + {\varepsilon _{11}}{E_\theta } + {d_1}\varTheta, \\ &{D_z} = {e_{31}}{\varepsilon _r} + {e_{31}}{\varepsilon _\theta } + {e_{33}}{\varepsilon _z} + {\varepsilon _{33}}{E_z} + {d_3}\varTheta,\end{aligned} \right.$

其中σ_i是法向应力; τ_ij是切向应力; ε_i是法向应变; γ_ij是切向应变; E_i是电场; D_i是电位移; c_ij是弹性常数; ε_ij是介电常数; e_ij是压电耦合系数; d_i是热电耦合系数; λ_ii是热机耦合系数^[26,27], ${\lambda _{{\rm{11}}}} = ( {c_{11}} + {c_{12}} + {c_{13}} ){\alpha _{\rm{1}}}$, ${\lambda _{33}} = \left( {2{c_{13}} + {c_{33}}} \right){\alpha _{3}}$; α_i为热膨胀系数; Θ为相对于初始温度T₀的温度增量.
有限(非线性)弹性应变位移关系为

$\begin{split}&{\varepsilon _r} = \!\frac{{\partial {U_r}}}{{\partial r}},\; \;{\varepsilon _\theta } = \!\frac{{\partial {U_\theta }}}{{r\partial \theta }} + \frac{{{U_r}}}{r},\; \;{\varepsilon _z} =\! \frac{{\partial U}}{{\partial z}}\! +\! \frac{1}{2}{\left( {\frac{{\partial U}}{{\partial z}}} \right)^2},\\&{\gamma _{r\theta }} = \frac{{\partial {U_r}}}{{r\partial \theta }} + \frac{{\partial {U_\theta }}}{{\partial r}} - \frac{{{U_\theta }}}{r},\; \;{\gamma _{\theta z}} = \frac{{\partial U}}{{r\partial \theta }} + \frac{{\partial {U_\theta }}}{{\partial z}},\; \;\;\\&{\gamma _{rz}} = \frac{{\partial U}}{{\partial r}} + \frac{{\partial {U_r}}}{{\partial z}}.\\[-15pt]\end{split} $

由于这是一维问题, 杆的横向边界的牵引力应该为零. 因此可以得出${\sigma _r} = 0$, ${\tau _{rz}} = 0$, ${\tau _{r\theta }} = 0$, ${D_r} = 0$, 从中可以得到如下关系:

${D_\theta } = {\rm{0}},\;\;\;{{\gamma _{rz}}} = {\gamma _{\theta z}} = 0,\;\;\;{{E_r} = {E_\theta } = 0} ,$

${\varepsilon _r} = \frac{{{e_{31}}{E_z} + {\lambda _{11}}\varTheta - {c_{12}}{\varepsilon _\theta } - {c_{13}}{\varepsilon _z}}}{{{c_{11}}}}.$

根据广义Hamilton变分原理可得

$\delta \int_{{t_0}}^t {\int_{{z_0}}^z {L{\rm{d}}z{\rm{d}}t = \delta } } \int_{{t_0}}^t {\int_{{z_0}}^z {(T - {E_{\rm{P}}}} } + {W_{\rm{e}}}){\rm{d}}z{\rm{d}}t = 0,$

式中L为Lagrange密度函数, T为系统的动能, E_P为系统的势能, W_e为系统的电能, 具体表达式分别为

$T = \frac{{{1}}}{{{2}}}\int\nolimits_V \rho {\rm{d}}V{\left(\frac{{\partial U}}{{\partial t}}\right)^2} + \frac{1}{2}\int\nolimits_V \rho {\rm{d}}V{\left(\frac{{\partial {U_r}}}{{\partial t}}\right)^2},$

${E_{\rm{P}}} = \frac{1}{2}\int\nolimits_V {{{{S}}^{\rm{T}}}{{Q}}{\rm{d}}V}, $

${W_{\rm{e}}} = \frac{1}{2}\int\nolimits_V {{{{E}}^{\rm{T}}}{{D}}{\rm{d}}} V,$

这里$\rho $为压电材料的密度, V为压电材料的体积, ${{S}} = {\left\{ {{\varepsilon _r}, {\varepsilon _\theta }, {\gamma _{r\theta }}} \right\}^{\rm{T}}}$和${{Q}} = {\left\{ {{\sigma _\theta }, {\sigma _z}, {\tau _{\theta z}}} \right\}^{\rm{T}}}$分别表示材料的应变向量和应力向量, ${{E}} = {\left\{ {0, 0, {E_z}} \right\}^{\rm{T}}}$表示电场向量, ${{D}} = {\left\{ {0, 0, {D_z}} \right\}^{\rm{T}}}$表示电位移向量.
根据Euler方程, 若$L = L( U, {U_z}, {U_t}, {U_{zz}}, {U_{tt}}, {U_{zt}}, \cdot \cdot \cdot )$, 则

$\begin{split}&\frac{{\partial L}}{{\partial U}} - \frac{\partial }{{\partial z}}\frac{{\partial L}}{{\partial {U_z}}} - \frac{\partial }{{\partial t}}\frac{{\partial L}}{{\partial {U_t}}} + \frac{{{\partial ^2}}}{{\partial {z^2}}}\frac{{\partial L}}{{\partial {U_{zz}}}} \\&+ \frac{{{\partial ^2}}}{{\partial {t^2}}}\frac{{\partial L}}{{\partial {U_{tt}}}} + \frac{{{\partial ^2}}}{{\partial z\partial t}}\frac{{\partial L}}{{\partial {U_{zt}}}} - \cdot \cdot \cdot = 0,\end{split}$

此处, ${U_z} \!=\! \dfrac{{\partial U}}{{\partial z}}, {U_t} \!=\! \dfrac{{\partial U}}{{\partial t}}, {U_{zt}} \!= \! \dfrac{{{\partial ^2}U}}{{\partial z\partial t}} \cdots $, 其余类同.
根据以上关系可以得到如下表达式:

$\left\{ \begin{aligned} &\frac{{\partial L}}{{\partial U}} = 0, \\ &\frac{\partial }{{\partial z}}\frac{{\partial L}}{{\partial {U_z}}} \\ =& - \frac{1}{2}V\bigg[ \left( {2{k_1}v_{\rm{eff}}^2 + 2{k_2}v_{\rm{eff}}^2 + 2{k_3} + {k_4} + {k_5}} \right) \Big({\frac{{\partial U}}{{\partial z}}} \Big) \\ & + 3\Big( {\frac{1}{2}{k_2}{v_{\rm{eff}}} + {k_5}} \Big){{\Big( {\frac{{\partial U}}{{\partial z}}} \Big)}^2} + {k_5}{{\Big( {\frac{{\partial U}}{{\partial z}}} \Big)}^3} \bigg], \\& \frac{\partial }{{\partial t}}\frac{{\partial L}}{{\partial {U_t}}} = \rho V\frac{{{\partial ^2}U}}{{\partial {t^2}}},~~\frac{{{\partial ^2}}}{{\partial {z^2}}}\frac{{\partial L}}{{\partial {U_{zz}}}} = 0,\\ &\frac{{{\partial ^2}}}{{\partial {t^2}}}\frac{{\partial L}}{{\partial {U_{tt}}}} = \rho V\frac{{{\partial ^2}}}{{\partial {t^2}}}\left( {\frac{{{\partial ^2}U}}{{\partial {t^2}}}} \right), \\ &\frac{{{\partial ^2}}}{{\partial z\partial t}}\frac{{\partial L}}{{\partial {U_{zt}}}} = \rho Vv_{\rm{eff}}^2\frac{{{\partial ^2}}}{{\partial z\partial t}}\left( {\frac{{{\partial ^2}U}}{{\partial z\partial t}}} \right),\end{aligned} \right.$

其中

$\begin{split}& {k_1} = \left( {{c_{11}} - \frac{{c_{12}^2}}{{{c_{11}}}}} \right),\;\;\;{{k_2} = 2\left( {{c_{13}} - \frac{{{c_{12}}{c_{13}}}}{{{c_{11}}}}} \right)} ,\\&{{k_3} = \left( {\frac{{{c_{13}}{e_{31}}}}{{{c_{11}}}} - {e_{33}}} \right){E_z}} ,\\& {k_4} = \left( {\frac{{{c_{13}}}}{{{c_{11}}}}{\lambda _{11}} - {\lambda _{33}}} \right)\varTheta ,\;\;\;{{k_5} = \left( {{c_{33}} - \frac{{c_{13}^2}}{{{c_{11}}}}} \right)} .\end{split} $

令位移梯度$u = \dfrac{{\partial U}}{{\partial z}}$, 由方程(9)可得出压电圆杆波导的纵波运动方程

$\frac{{{\partial ^{\rm{2}}}u}}{{\partial {t^2}}} - {c_0}^2\frac{{{\partial ^2}u}}{{\partial {z^2}}} = \frac{{{\partial ^2}}}{{\partial {z^2}}}\left( {\alpha {u^2} + \beta \frac{{{\partial ^2}u}}{{\partial {t^2}}}} \right),$

$\begin{split}&c_0^2 = \frac{{{k_1}v_{\rm{eff}}^2 + {k_3} + \dfrac{1}{2}{k_4} + {k_5}}}{\rho },\;\;\;\\&{\alpha = \frac{{3{k_2}{v_{\rm{eff}}} + 6{k_5}}}{{4\rho }}},\;\;\;{\beta = v_{\rm{eff}}^2{r^2}}, \end{split}$

其中, c₀是压电圆杆的线性纵向波速; α为耗散系数, β为弥散系数, 两者都由材料的性质和几何参数决定. 值得说明的是, 如果忽略电场、热-电耦合, 考虑纯弹性杆, 结果与文献[16]一致.
在推导方程(12)时, 使用了有效泊松比v_eff^[28]. 可以得到

$\begin{split}\;& {v_{\rm{eff}}} =- {{{\varepsilon _r}}}/{{{\varepsilon _z}}} \\=\;& \frac{{{\varepsilon _{33}}\left( {{c_{12}}{c_{13}} + {c_{11}}{c_{13}}} \right) - {e_{31}}{e_{33}}({c_{11}} + {c_{12}})}}{{2c_{12}^2{\varepsilon _{33}} - 2{c_{11}}e_{31}^2 + {c_{12}}e_{31}^2 - c_{11}^2{\varepsilon _{33}}}}.\end{split}$

假设方程(12)的行波解为

$u = u\left( \xi \right),\;\;\;{\xi = k\left( {z - ct} \right)},$

将方程(15)代入方程(12)可以转化为一个常微分方程:

$\beta {c^2}{k^2}\frac{{{{\rm{d}}^4}u}}{{{\rm{d}}{\xi ^4}}} - \left( {{c^2} - c_0^2} \right)\frac{{{{\rm{d}}^2}u}}{{{\rm{d}}{\xi ^2}}} + \alpha \frac{{{{\rm{d}}^2}{u^2}}}{{{\rm{d}}{\xi ^2}}} = 0,$

其中k和c分别是波数和波速. 用下列Jacobi椭圆函数表示方程(16)的解:

$u\left( \xi \right) = \sum\limits_{j = 0}^n {{a_j}} {{\rm{cn}}^j}\xi .$

可以看出$u\left( \xi \right)$的最高阶数是n, 即

$O\left( {u\left( \xi \right)} \right) = n.$

为了讨论方便, 列出三种Jacobi椭圆函数之间的关系以及渐近值

${{\rm{sn}}^2}\xi + {{\rm{cn}}^2}\xi = 1,\;\;\;{{{\rm{dn}}^2}} \xi + {m^2}{{\rm{sn}}^2}\xi = 1,$

其中${\rm{cn}}\xi $为Jacobi椭圆余弦函数, ${\rm{dn}}\xi $为第三类Jacobi椭圆函数, m为模数$\left( {0 \leqslant m \leqslant 1} \right)$.

$\begin{split}&\frac{\rm{d}}{{{\rm{d}}\xi }}{\rm{sn}}\xi = {\rm{cn}}\xi {\rm{dn}}\xi,\;\;\;{\frac{\rm{d}}{{{\rm{d}}\xi }}{\rm{cn}}\xi = - {\rm{sn}}\xi {\rm{dn}}\xi },\;\;\;\\&{\frac{\rm{d}}{{{\rm{d}}\xi }}{\rm{dn}}\xi = - {m^2}{\rm{sn}}\xi {\rm{dn}}\xi }. \end{split}$

根据上述微分关系, 很容易有

$O\left( {\frac{{{\rm{d}}u}}{{{\rm{d}}\xi }}} \right) = n + 1.$

类似地, 可以推出:

$O\left( {{u^2}} \right) = 2n,\;\;\;\;{O\left( {\frac{{{{\rm{d}}^2}u}}{{{\rm{d}}{\xi ^2}}}} \right)} = n + 2.$

因此, 将方程(16)进一步化为

${k^2}\frac{{{{\rm{d}}^4}u}}{{{\rm{d}}{\xi ^4}}} - {N_1}\frac{{{{\rm{d}}^2}u}}{{{\rm{d}}{\xi ^2}}} + {N_2}\frac{{{{\rm{d}}^2}{u^2}}}{{{\rm{d}}{\xi ^2}}} = 0,$

其中

${N_1} = \frac{{{c^2} - c_0^2}}{{\beta {c^2}}},\;\;\;{{N_2} = \frac{\alpha }{{\beta {c^2}}}} .$

对ξ积分两次, 为计算方便, 令积分常数为零. 可以得到

${k^2}{u_{\xi \xi }} - {N_1}u + {N_2}\left( {{u^2}} \right) = 0.$

通过谐波平衡法, 使方程(25)中的非线性项次数和微分项最高阶数相等, 结合方程(22)可以确定方程(18)中的最高阶数n = 2.
根据Jacobi椭圆余弦函数展开法, 方程(17)的解有如下表达形式:

$u\left( \xi \right) = {a_0} + {a_1}{\rm{cn}}\xi + {a_2}{{\rm{cn}}^2}\xi .$

方程(26)对ξ微分两次可求得

$\begin{split}\frac{{{{\rm{d}}^2}u}}{{{\rm{d}}{\xi ^2}}} =\;& 2{a_2}(1 - {m^2}) - {a_1}(1 - 2{m^2}){\rm{cn}}\xi \\&+ 4{a_2}(2{m^2} - 1){{\rm{cn}}^2}\xi - 2{a_1}{m^2}{{\rm{cn}}^3}\xi \\&- 6{a_2}{m^2}{{\rm{cn}}^4}\xi .\end{split}$

将方程(26)和方程(27)代入方程(25), 并比较${\rm{cn}}\xi $相同次幂的系数, 可得

${a_0} = \frac{{{N_1} - 4{k^2}\left( {2{m^2} - 1} \right)}}{{2{N_2}}},\;\;\;{{a_1}}=0,\;\;\;{{a_2}} = \frac{{6{k^2}{m^2}}}{{{N_2}}}.$

因此方程(23)的精确周期解为

$u\left( \xi \right) = \frac{{{N_1} - 4{k^2}(2{m^2} - 1)}}{{2{N_2}}} + \frac{{6{k^2}{m^2}}}{{{N_2}}}{{\rm{cn}}^2}(\xi,m),$

其中m为模数(0 < m < 1). 应该指出, 当$m \to 1$时, ${\rm{cn}}\xi \to{\rm{sech}}\xi $, 于是非线性波动方程(29)的孤波解可以写成:

$u\left( \xi \right) = \frac{{{N_1} - 4{k^2}}}{{2{N_2}}} + \frac{{6{k^2}}}{{{N_2}}}{{\rm{sech}}^2}\xi .$

由于此处只讨论波的特性, 令常数项为零. 即

${k^2} = \frac{{{c^2} - c_0^2}}{{4\beta {c^2}}}.$

将方程(31)代入方程(30), 得到方程(25)的一个标准孤波解:

$u\left( \xi \right) = A{{\rm{sech}}^2}\frac{{z - ct}}{\varLambda },$

式中A为波幅, Λ是波长,

$A = \frac{{3\left( {{c^2} - c_0^2} \right)}}{{2{\alpha }}},\;\;{\varLambda = \frac{{2{\rm{\pi}} }}{k}} = 4{\rm{\pi}} \sqrt {\frac{{{c^2}\beta }}{{{c^2} - c_0^2}}}, $

其中c > c₀是孤立波存在的条件.
利用Jacobi椭圆正弦函数展开法, 方程(17)的解有如下表达形式:

$u\left( \xi \right) = {a_0} + {a_1}{\rm{sn}}\xi + {a_2}{{\rm{sn}}^2}\xi .$

将方程(34)对ξ微分两次求得

$\begin{split}\frac{{{{\rm{d}}^2}u}}{{{\rm{d}}{\xi ^2}}} =\;& 2{a_2} - {a_1}(1 + 2{m^2}){\rm{sn}}\xi - 4{a_2}(1 + {m^2}){{\rm{sn}}^2}\xi \\&+ 2{a_2}{m^2}{{\rm{sn}}^3}\xi + 6{a_2}{m^2}{{\rm{sn}}^4}\xi.\\[-10pt] \end{split}$

将方程(34)和方程(35)代入方程(25), 并比较${\rm{sn}}\xi $相同次幂的系数, 可得

${a_0} = \frac{{4{k^2}\left( {1 + {m^2}} \right) + {N_1}}}{{2{N_2}}},\;\;\;{{a_1} = 0},\;\;\;{{a_2} = - \frac{{6{k^2}{m^2}}}{{{N_2}}}} .$

因此方程(23)的精确周期解为

$u\left( \xi \right) = \frac{{4{k^2}(1 + {m^2}) + {N_1}}}{{2{N_2}}} + \frac{{ - 6{k^2}{m^2}}}{{{N_2}}}{{\rm{sn}}^2}\left( {\xi,m} \right).$

将方程(19)代入到方程(37)求得

$u\left( \xi \right) = \frac{{4{k^2}\left( {1 - 2{m^2}} \right) + {N_1}}}{{2{N_2}}} + \frac{{6{k^2}{m^2}}}{{{N_2}}}{{\rm{cn}}^2}(\xi,m),$

这就是方程(23)的另一个精确周期解, 可以观察到取m = 1时, 方程(38)退化为与方程(29)同样的形式.

本文通过Matlab软件, 对如图1所示的模型进行数值模拟, 取无限长压电圆杆为BaTiO₃材料, 初始温度为T₀ = 20 ℃, 杆的半径分别为R₁ = 0.025 m, R₂ = 0.05 m和R₃ = 0.075 m. 将这些参数代入相关方程进行计算, 通过改变温度及波速比值, 得到相应的模拟结果.
根据表1所列参数, 通过改变温差Θ的大小, 可以计算出不同温度下波速c₀的大小. 如表2所列, 可以看出, 当温差Θ = 10 ℃时, 波速c₀最大, 随着温度的升高, 波速c₀逐渐减小, 但也可以看出波速c₀并没有存在大幅度的衰减.

参数	参数值		参数	参数值
c11/(N·m–2)	166 × 109		e15/(C·m–2)	11.6
c12/(N·m–2)	77 × 109		ε11/(C2·N–2·m–1)	11.2 × 10–9
c13/(N·m–2)	78 × 109		ε33/(C2·N-2·m–1)	12.9 × 10–9
c33/(N·m–2)	162 × 109		d1/(C·km–2)	–5.4831 × 10–6
c44/(N·m–2)	43 × 109		d3/(C·km–2)	–5.4831 × 10–6
e31/(C·m–2)	–4.4		α1 = 0.5α3	10 × 10–6
e33/(C·m–2)	18.6		β1/m2	0.5278 × 10–4
β2/m2	2.112 × 10–4		β3/m2	4.7502 × 10–4
α/(N·m·kg–1)	3.5561 × 107		veff	0.2906
ρ/(kg·m-3)	5.8 × 103

表1钛酸钡材料参数^[29]
Table1.Material parameters of BaTiO₃^[29].

Θ/℃	c0/(103 m·s–1)
10	5.4482
50	5.4367
90	5.4251

表2不同温度下波速比较
Table2.Comparison of the wave velocities at different temperature.

通过改变波速比c/c₀的大小, 可求得在不同比值下的波幅A、波长Λ和波数k. 如表3所列, 可以看出, 随着波速比的增大, 波幅A和波数k增大, 而波长Λ随之减小. 同样地, 如表4所列, 固定波速比c/c₀ = 1.1, 改变压电圆杆的半径, 分别取R₁ = 0.025 m, R₂ = 0.05 m和R₃ = 0.075 m. 能够发现随着压电圆杆半径R的增大, 波长Λ随之增大, 而波数k随之减小.

c/c0	A/m	Λ/m	k
1.1	0.2567	0.5429	11.5675
1.2	0.5378	0.4092	15.3470
1.3	0.8434	0.3540	17.740

表3R = 0.05 m时不同波速比下参数比较
Table3.Comparison of parameters under different wave velocity ratios when R = 0.05 m.

R/m	Λ/m	k
0.025	0.2705	23.2162
0.050	0.5429	11.5675
0.075	0.8116	7.7375

表4波速比c/c₀ = 1.1时不同半径下参数比较
Table4.Comparison of parameters under different radii when c/c₀ = 1.1.

图2给出了在不考虑温度影响时不同波速比c/c₀下, 孤波波速u和变量ξ的关系. 可以看出, 当ξ = 0时, u达到最大值, 并且幅值关于ξ = 0对称, 同时可以观察到, 随着c/c₀的增大, 孤波幅值增大, 波长减小. 简而言之, 孤波振幅越大, 波长越小, 这体现了非线性孤波的弥散特性.
图 2 不同波速比c/c₀下孤波波速u与变量ξ的关系
Figure2. Relationship between solitary wave u and variable ξ under different wave velocity ratio c/c₀ values.

图3给出了在波速比c/c₀ = 1.3, 时间固定在t = 0.01 s时, 当Θ分别取10, 50, 90 ℃三种不同温度下的波形^[30], 观察到随着温度Θ的改变, 孤立波在传播过程中波形并没有发生改变, 体现了其稳定性. 同时可以看出, 当波速比c/c₀一定时它们的波幅非常接近, 温度的改变对波幅的影响并不是很明显, 但是随着温度的逐渐升高, 波速却逐渐降低, 这一点也与表2的数据相符合.
图 3 波速比c/c₀ = 1.3时三种不同温度下的波形
Figure3. Three waveforms at different temperatures when the velocity ratio of c/c₀ = 1.3.

图4给出了当温差Θ = 50 ℃时, 波速比c/c₀分别取为1.1, 1.2和1.3时的波形, 可知, 孤立波的能量主要集中在中间有效区域, 并且没有因为波速比的改变而扩散, 体现了孤立波的稳定性. 而当温度一定时, 随着波速比c/c₀的升高, 波幅增大, 波长减小. 这一点与图2表示的孤波特性相一致.
图 4 当Θ = 50 ℃时不同波速比c/c₀下的波形
Figure4. Waveforms under different wave velocity ratio c/c₀ values when Θ = 50 ℃.

图5给出了孤立波的三维曲面, 图6给出了压电圆杆的孤波特性. 可以看出, 当给定某一时间t和z之后, 孤立子就会出现, 这说明了孤立波不是单独关于时间或空间的单一变量, 而是以时间和空间为组合的变量, 随着时间的变化, 波在传播. 并且相同的波形在给定t和z的组合后会重复出现, 且能量比较集中, 这与孤波的稳定性相符合. 因此它在当代通信技术、缺陷检测等许多方面有很大的研究意义和应用潜力.
图 5 当波速比c/c₀ = 1.3, Θ = 50 ℃时三维曲面图
Figure5. Three-dimensional surface figure when the wave ratio c/c₀ = 1.3, Θ = 50 ℃.

图 6 压电圆杆的孤波特性
Figure6. Solitariness of piezoelectric rod.

图7给出了在三种不同温度下, 压电圆杆的波速c和波数k的关系图. 可以看出, 当温度一定时, 随着波数k的增加, 波速c也呈增加趋势, 同样地, 当波数k一定时, 随着温度的升高波速c反而呈减小趋势, 这与表2以及图3和图4的模拟结果相符合.
图 7 三种不同温度下波速c和波数k的关系图
Figure7. Graph of wave velocity c and wave number k at three different temperatures.

孤立波在数学上是一类非线性偏微分方程的局部行波解, 这类方程可表示为运动项+色散项+非线性项+(耗散项) = 0. 从物理本质上讲, 非线性效应使得波形在传播过程中出现陡突(能量聚集). 孤立子就是由非线性场激发的、能量不弥散的、形态上稳定的粒子. 本文研究了温度效应下无限长压电圆杆的孤波问题, 由于结构的有限变形(如轴向压缩等)可引起非线性效应, 而二次运动和变形(如横向泊松效应等)可分散这些效应. 相互作用在非线性效应和色散效应之间, 在一定条件下就产生了孤波这一稳定传播的行波. 利用Hamilton变分及Euler方程推导出圆杆的波动方程, 并采用Jacobi椭圆余弦函数展开法和椭圆正弦函数展开法对推导出的波动方程进行求解. 通过Matlab软件进行数值模拟, 发现压电圆杆中不仅有孤波存在, 而且在不同温度下具有不同的性质. 可以得到如下结论:
1)由计算结果可分析得到, 当固定波速比时, 温度的改变对波速的影响比较明显, 随着温度的升高, 波速逐渐降低.
2)当固定温度时, 波速比的改变对孤波的幅值影响比较明显, 随着波速比的增大, 波幅逐渐升高, 这也是孤波的特点之一.
3)温度的改变虽然对孤立波有一定的影响, 但在传播的过程中, 孤立波仍然是关于ξ对称的钟型波, 这也体现了非线性和色散效应共同作用下孤立波的稳定特性. 另外, 利用Jacobi椭圆函数求解得到压电圆杆的精确周期解, 周期解可退化为孤波解, 从理论上也证明了压电圆杆中可能有稳定传播的孤波. 目前, 波动理论在结构的无损检测和提高信息传输质量等方面得到了较为广泛的应用. 因此, 将波动理论运用到压电材料中具有现实工程意义和理论研究价值.