引言

充气薄膜结构是各类工程应用中的一种常见结构，如用于观测气象要素的气象气球、探空气球^[1-2]以及各类介电驱动器和传感装置^[3-5]等；类地行星和空间结构^[6-7]等也经常简化为充气薄膜结构进行分析；在生物力学领域，充气薄膜结构经常被用来模拟囊状动脉瘤^[8-10]和纺锤形动脉瘤^[11-13].

形状分岔是超弹性材料充气薄膜结构在内压载荷作用下的一个典型问题.以往的研究表明，受均匀内压作用的超弹性球形薄膜在膨胀过程中内压随体积的变化曲线为N形曲线，在曲线下降段的某个内压区间内，球膜不再保持均匀膨胀变形，而是一端局部变薄分岔为梨形形状^[14-16].薄壁圆管的典型特点为在临界内压作用下发生局部起鼓失稳现象^[17-21].也就是说，当内压达到某个临界值后，同一个内压下充气薄膜结构会存在多个平衡解.而实际上是否会出现，或者在实验中是否会观察到某个平衡解对应的形状，就需要对每个平衡解进行稳定性分析.

Fu等^[16]采用能量判据研究了轴对称扰动下，球膜在膨胀过程中出现的球形解和梨形解的稳定性.结果表明当球形解和梨形解共存时，球形解在质量控制和内压控制下均为不稳定解，而梨形解在质量控制下为稳定解，但在内压控制下为不稳定解，这与Alexander^[22]实验中观察到的现象相一致. Liang等^[23]分析了介电超弹性球膜在内压和电压共同作用下的形状分岔现象，发现在曲线下降段同样有梨形解的存在，并采用相同的能量判据分别研究了内压和电压控制、内压和电荷控制、质量和电压控制以及质量和电荷控制四种加载方式下曲线下降段分岔解的稳定性，指出在内压和电压控制时分岔解为不稳定解，但在其他三种加载方式下为稳定解.另外，Chen^[24]和Fu等^[25]同样利用能量判据分别分析了内压作用下薄壁圆管的均匀解和局部起鼓解的稳定性.

针对薄膜椭球结构在内压作用下的变形行为，部分学者也进行了相关研究. Feng^[26]根据对称性选取二分之一椭圆模型，模拟了黏弹性超弹薄膜椭球在内压作用下产生大变形时材料的蠕变和软化现象. Sagiv^[27]则研究了多种形状的薄膜半椭球膨胀过程中的应力分布规律，但并未涉及分岔现象.最近，Geng等^[28]将薄膜椭球分为橄榄形和南瓜形，分别研究了内压载荷作用下模型的形状分岔现象，指出橄榄形薄膜球的细长比存在一个临界值，当细长比大于该值时，橄榄形薄膜球在内压达到某载荷值时会发生梨形分岔；而当细长比小于该临界值时，将不会出现梨形分岔，但当细长比进一步减小时，会出现类似于圆管的局部分岔现象.另外，南瓜形薄膜球则无论圆扁，在内压达到某载荷值时都会出现梨形分岔.

本文将在Geng等^[28]研究薄膜椭球分岔的基础上，针对压强控制和质量控制两种加载方式，采用连续介质力学中广泛用于稳定性分析的方法^[29-30]--能量判据方法，分析不同几何形状的薄膜椭球在内压作用下的均匀解及分岔解的稳定性.

1 控制方程和分岔解

考察一个受均匀内压P作用的不可压超弹性薄膜椭球的膨胀问题.变形前的构型在柱坐标系下描述为

$ R(\theta ) = b\sin \theta, \;\;Z(\theta ) = a(1-\cos \theta ), \;\;0 \le \theta \le {\rm{\pi }} $

(1)

其中，R-Z平面为生成薄膜椭球的旋转平面，a和b分别为变形前椭圆的半轴长，椭圆绕Z轴旋转生成薄膜椭球.根据a和b的长短不同，旋转生成的薄膜椭球有两种不同形状：当a>b时，生成的薄膜椭球如图 1(b)所示，为橄榄状，即两个相等的半轴长度b比另外一个半轴长度a小；当a < b时，生成的薄膜椭球如图 1(c)所示，为南瓜状，即两个相等的半轴长度b比另外一个半轴长度a大.为方便起见，分别将橄榄状薄膜椭球中的a设为单位长度，即采用a=1，b < 1描述；将南瓜状薄膜椭球中的b设为单位长度，即采用b=1，a < 1描述.

变形后的构型描述为r=r(θ)，z=z(θ)，分别用dS和ds表示在参考构型和当前构型下角度增量dθ对应的弧长增量，可得

$ \left.\begin{array}{l} {\rm d}S = \sqrt {({\rm d}R)^2 + ({\rm d}Z)^2} = \chi (\theta ){\rm d}\theta \\ {\rm d}s = \sqrt {\left(\dfrac{{\rm d}r}{{\rm d}\theta }\right)^2 + \left(\dfrac{{\rm d}z}{{\rm d}\theta }\right)^2} {\rm d}\theta \\ \end{array}\right\} $

(2)

其中$ \chi (\theta ) = \sqrt {a^2\sin ^2\theta + b^2\cos^2\theta } $.主方向分别沿椭球的纬度、经度及与变形表面垂直方向，对应主方向上的主伸长分别为

$ \left.\begin{array}{lll}\lambda _1 = \dfrac{r}{R}, \quad \lambda _2 = \dfrac{{\rm d}s}{{\rm d}S} = \dfrac{\mu _2 }{\chi (\theta )} \equiv\\ \dfrac{1}{\chi (\theta )}\sqrt {\left(\dfrac{{\rm d}r}{{\rm d}\theta }\right)^2 + \left(\dfrac{{\rm d}z}{{\rm d}\theta }\right)^2}, \\ \lambda _3 = \dfrac{h}{H}\end{array}\right\} $

(3)

其中，H和h分别为变形前后薄膜椭球的厚度.

本文考虑不可压Ogden本构模型^[31]，其应变能函数为

$ \tilde {W} = \sum\limits_{r = 1}^3 {\tilde {\mu }_r } (\lambda _1^{\alpha _r } + \lambda _2^{\alpha _r } + \lambda _3^{\alpha _r } - 3) / \alpha _r $

(4)

该模型用于描述8%硫化天然橡胶时，上式中材料常数取为：α₁=1.3，α₂=5.0，α₃=-2.0，$ {{\tilde \mu }_1} = 1.491, {{\tilde \mu }_2} = 0.003, {{\tilde \mu }_3} =-0.023 $，并且这里的应变能函数$ \tilde{W} $为被剪切模量无量纲化之后的函数.

根据文献[28]，平衡方程可写为

$ \left. {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}{\lambda _1}}}{{{\rm{d}}\theta }} = \frac{{{\lambda _2}\chi (\theta )\sin \phi - b{\lambda _1}\cos \theta }}{{b\sin \theta }}}\\ {\frac{{{\rm{d}}{\lambda _2}}}{{{\rm{d}}\theta }} = \frac{{{W_1} - {\lambda _2}{W_{12}}}}{{b{W_{22}}}}\frac{{\chi (\theta )\sin \phi }}{{\sin \theta }} - \frac{{{W_2} - {\lambda _1}{W_{12}}}}{{{W_{22}}}}\cot \theta }\\ {\frac{{{\rm{d}}\phi }}{{{\rm{d}}\theta }} = \frac{{{W_1}}}{{{W_2}}}\frac{{\chi (\theta )\cos \phi }}{{b\sin \theta }} - \frac{{P{\lambda _1}{\lambda _2}\chi (\theta )}}{{{W_2}}}} \end{array}} \right\} $

(5)

其中，φ为变形后经度方向与z轴夹角，W(λ₁, λ₂)为初始构型下单位表面面积的应变能，与三维的应变能函数的关系为$ W(\lambda _1, \lambda _2 ) = \tilde {W}(\lambda _1, \lambda _2, \lambda _1^{ - 1} \lambda _2^{ - 1} ) $，另外，$ {W_1} = \partial W/\partial {\lambda _1}{\rm{, }}{W_{12}} = {\partial ^2}W/\partial {\lambda _1}\partial {\lambda _2}{\rm{, }}{W_{22}} = {\partial ^2}W/\partial \lambda _2^2 $，P为真实内压与薄膜厚度H的比值.

薄膜椭球的轴对称边界条件可写为

$ \left.\begin{array}{l} \lambda _1 (0) = \lambda _2 (0), ~\lambda _1 (\rm{\pi } ) = \lambda _2 (\rm{\pi })\\ \;\phi (0) = \rm{\pi } / 2, ~\phi (\rm{\pi } ) = - \rm{\pi } / 2 \\ \end{array}\right\} $

(6)

控制方程(5)在边界条件(6)下的解可用数值方法--打靶法求解^[28].结果表明，当薄膜椭球的初始几何形状为橄榄形时，如图 2所示，内压随体积变化曲线均为N形曲线，而在曲线下降段是否存在类似于球形薄膜膨胀过程中出现的梨形分岔区间(P_cr1, P_cr2)是由橄榄形薄膜的细长比决定的.取a=1时，短轴b存在临界值，当b大于该临界值时，存在梨形分岔区间(P_cr1, P_cr2)，即当载荷P取这个区间的值时，除了三个均匀解，还有两个关于θ=π/2平面对称的梨形解（这两个梨形解对应的变形后体积相同）存在，图 3为均匀解与梨形分岔解的变形模态示意图；随着b值减小，这个区间越来越小；当b小于临界值后，梨形分岔解消失.对于Ogden材料，b的临界值约取0.626.进一步减小b的取值，在下降段只有一个均匀解存在，直到约取b=0.25时，由于第一个上升段和下降段的两个均匀解对应的λ₁(0)非常接近，打靶法在此求解困难.进一步分析可得，当取a=1和b=0.25时下降段的均匀解类似于管状薄膜膨胀过程中的局部起鼓解.

与橄榄形薄膜不同的是，当初始几何形状为如图 1(c)所示的南瓜形时，无论a多小，在P-V曲线的下降段总会出现有分岔解存在的区间^[28].

2 压强控制下的稳定性分析

得到薄膜椭球在内压载荷作用下的均匀解和分岔解之后，现在考虑在压强控制下各个解对应的平衡状态的稳定性.薄膜椭球的总势能（除以2πH后）可以写为

$ \mathit{\Pi } = \int_0^{\rm{ \mathsf{ π} }} {} \left( {W({\lambda _1}, {\lambda _2})b\sin \theta \chi (\theta )-\frac{1}{2}P{r^2}z'} \right){\rm{d}}\theta $

(7)

为方便起见，上式可简写为

$ \mathit{\Pi } = \int_0^{\rm{ \mathsf{ π} }} {} L(\mathit{\boldsymbol{u}}, \mathit{\boldsymbol{u'}}){\rm{d}}\theta $

(8)

其中

$ L(\mathit{\boldsymbol{u}}, \mathit{\boldsymbol{u'}}) = W({\lambda _1}, {\lambda _2})b\sin \theta \chi (\theta )-P{r^2}z'/2 $

(9)

这里u=(r(θ), z(θ)).

假设u=(r, z)为考察的平衡状态，对应的内压为P，主伸长分别为λ₁和λ₂.若u=u对应平衡状态的能量为极小值，则该解为稳定解.这里考虑扰动解u=u+v.当扰动位移v足够小时，存在

$ \mathit{\Pi }(\mathit{\boldsymbol{\overline u}} + \mathit{\boldsymbol{v}})-\mathit{\Pi }(\mathit{\boldsymbol{\overline u}} ) = \frac{1}{2}{\delta ^2}\mathit{\Pi } + O(|\mathit{\boldsymbol{v}}{|^3}) $

(10)

此时，根据能量二阶变分δ²Π的正负号即可判断u=u的稳定性：若δ²Π > 0，说明扰动解的能量比平衡解的能量高，则平衡解为稳定解；相反，若δ²Π < 0，说明扰动解具有较平衡解更低的能量，则平衡解为不稳定解.

总能量表达式(8)的二阶变分可写为

$ {\delta ^2}\mathit{\Pi } = \int_0^{\rm{ \mathsf{ π} }} {} (\mathit{\boldsymbol{v}} \cdot \mathit{\boldsymbol{S}} \cdot \mathit{\boldsymbol{v}} + 2\mathit{\boldsymbol{v}} \cdot \mathit{\boldsymbol{R}} \cdot {\mathit{\boldsymbol{v}}^\prime } + {\mathit{\boldsymbol{v}}^\prime } \cdot \mathit{\boldsymbol{Q}} \cdot {\mathit{\boldsymbol{v}}^\prime }){\rm{ d}}\theta $

(11)

其中，S, R和Q均为矩阵函数，其分量可分别写为

$ \left. {\begin{array}{*{20}{c}} {{S_{ij}} = {{\left. {\frac{{{\partial ^2}L}}{{\partial {u_i}\partial {u_j}}}} \right|}_{\mathit{\boldsymbol{u}} = \mathit{\boldsymbol{\overline u}} }}, \quad {R_{ij}} = {{\left. {\frac{{{\partial ^2}L}}{{\partial {u_i}\partial {{u'}_j}}}} \right|}_{\mathit{\boldsymbol{u}} = \mathit{\boldsymbol{\overline u}} }}}\\ {{Q_{ij}} = {{\left. {\frac{{{\partial ^2}L}}{{\partial {{u'}_i}\partial {{u'}_j}}}} \right|}_{\mathit{\boldsymbol{u}} = \mathit{\boldsymbol{\overline u}} }}} \end{array}} \right\} $

(12)

假设当0 < θ < π，Q始终为正，二阶变分δ²Π在归一化条件^[24]

$ \int_0^{\rm{\pi }} {\left( {v_1^2 + v_2^{'2}} \right)} {\rm{d}}\theta = 1 $

(13)

限制下存在极小值，极小值的正负号可以用来判断u=u的稳定性.当极小值为正时，解u为稳定的；当该值为负时，解u则不稳定.

将L的表达式(9)代入式(12)，再代回式(11)中，经化简可得能量二阶变分形式

$ \begin{array}{*{20}{c}} {{\delta ^2}\mathit{\Pi } = \int_0^{\rm{ \mathsf{ π} }} {({a_1}v_1^2 + 2{a_2}{v_1}{{v'}_2} + } }\\ {{a_3}v_2^{'2} + {a_4}v_2^{'1} + 2{a_5}{{v'}_1}{{v'}_2}){\rm{ d}}\theta } \end{array} $

(14)

其中

$ \left.\begin{array}{l} a_1 = W_{11} \chi (\theta ) / R(\theta ) - (W_{12} {r}' / \mu _2 {)}' - \bar {P} \cdot {z}' \\ a_2 = W_{12} {z}' / \mu _2 - \bar {P}r \\ a_3 = R(\theta )(\lambda _2 W_{22} - W_2 )\mu _2^{ - 3} {z}'^2 + R(\theta )W_2 \mu _2 \\ a_4 = R(\theta )(\lambda _2 W_{22} - W_2 )\mu _2^{ - 3} {r}'^2 + R(\theta )W_2 \mu _2 \\ a_5 = R(\theta )(\lambda _2 W_{22} - W_2 )\mu _2^{ - 3} {r}'{z}' \\ \end{array}\right\} $

(15)

右面的表达式均为在平衡解u=u时的取值.此时，Q> 0的条件等价于a₃ > 0，a₄ > 0和a₃a₄-a₅²>0，经检验，所有要考察稳定性的平衡解均满足这些不等式.采用拉格朗日乘数法，能量二阶变分表达式(14)在归一化条件(13)下极小值问题的拉格朗日函数可写为

$ G = {\delta ^2}\mathit{\Pi }-\beta \left( {\int_0^{\rm{ \mathsf{ π} }} {} \left( {v_1^2 + v_2^{'2}} \right){\rm{ d}}\theta-1} \right) $

(16)

其中β为拉格朗日乘子，由δG=0可得

$ a_1 v_1 + a_2 {v}'_2 - (a_4 {v'}_1 + a_5 {v'}_2 )^\prime = \beta v_1 $

(17a)

$ a_2 v_1 + a_3 {v}'_2 + a_5 {v}'_1 = \beta {v}'_2 $

(17b)

式(17a)式两边同时乘以v₁，式(17b)式两边同时乘以v'₂，相加之后并积分，可得在归一化条件(13)下，能量的二阶变分δ²Π与β数值上相等.将上面两个表达式中的v'₂约去，可以得到一个简单的特征值问题

$ {[{b_2}{{v'}_1}]^\prime } + {b_0}{v_1} = 0, \;0 \le \theta \le {\rm{\pi }} $

(18)

$ {v_1}(0) = {v_1}({\rm{\pi }}) = 0 $

(19)

其中

$ \left. {\begin{array}{*{20}{l}} {{b_0} =-{a_1}-{{\left( {\frac{{{a_2}{a_5}}}{{{a_3}-\beta }}} \right)}^\prime } + \frac{{a_2^2}}{{{a_3} - \beta }} + \beta }\\ {{b_2} = {a_4} - \frac{{a_5^2}}{{{a_3} - \beta }}} \end{array}} \right\} $

(20)

由此，分析解u=u的稳定性可以转化为特征值问题(18)的求解，其最小特征值β的正负号可判断u=u是否为稳定解.当特征值问题(18)存在至少一个负的特征值时，可得δ²Π < 0，v取该特征值对应的特征函数，有Π(u+v)-Π(u) < 0，表明扰动解u=u+ v能量更小，因此u=u为不稳定解；反之，如果没有负特征值或零特征值存在，那么u=u为稳定解.所以在我们的数值计算中，只需寻找β的负数解是否存在即可判断u=u解的稳定性.另外注意到，式(20)有$ \frac{1}{{{a_3}-\beta }} $项，经验证a₃ > 0，所以至少在β < 0的情况下，式(20)中不会出现分母为零的情况.

注意到在θ=0, π处，R(0)=R(π)=0，使得系数a₁在这两点存在奇异性，那么特征值问题(18)的求解区间需要避开这两点，替换为θ ∈ [δ, π-δ]，其中δ为小量，在此取δ=0.01.

将v₁在θ=0附近的展开式

$ v_1 = {v}'_1 (0)\theta + \frac{1}{2}{v}''_1 (0)\theta ^2 + \frac{1}{6}v_1^{(3)} (0)\theta ^3 +\cdot \cdot \cdot $

(21)

代入式(18)左边，并令含θ各级次幂的系数为零，可得v"₁(0)及其他高阶导数关于v'₁(0)的表达式，比如

$ \frac{{v}''_1 (0)}{{v}'_1 (0)} = - \frac{2\beta }{3W_{11} } $

(22)

$ \begin{array}{*{20}{c}} {\frac{{v_1^{(3)}(0)}}{{{{v'}_1}(0)}} = \frac{1}{{8{b^2}\lambda _1^2{W_{11}}}}\left\{ {18{W_1}({a^2}{\lambda _1} + {b^2}{{\lambda ''}_2}- 3{b^2}{{\lambda ''}_1})} \right. + }\\ {\begin{array}{*{20}{c}} {2\lambda _1^2{b^2}{\beta ^2}W_{11}^{- 1} + \lambda _1^2{W_{11}}(30{a^2}- 8{b^2}) + }\\ {6{a^2}\lambda _1^2{W_{12}} - 6{b^2}{\lambda _1}(3{W_{11}} + {W_{12}})}\\ {(3{{\lambda ''}_1} - {{\lambda ''}_2}) + 3{b^2}\lambda _1^2}\\ {[{W_{111}}({{\lambda ''}_1}-3{{\lambda ''}_2})-{W_{112}}(5{{\lambda ''}_1} + {{\lambda ''}_2})-} \end{array}}\\ {\left. {2\bar P\sqrt {{a^2}\lambda _1^2 + {b^2}{\lambda _1}({{\lambda ''}_2} - 3{{\lambda ''}_1})}]} \right\}} \end{array} $

(23)

其中，式(22)和式(23)右边在平衡解u=u中θ=0处取值.类似的，可以得到在θ=π处的结果

$ \frac{{v}''_1 (\rm{\pi } )}{{v}'_1 (\rm{\pi } )} = \frac{2\beta }{3W_{11} } $

(24)

$v_1^{(3)} (\rm{\pi } ) / {v}'_1 (\rm{\pi }) $的表达式和式(23)形式相同，但在θ=π处取值.

另外，在θ=δ和θ=π-δ处分别有

$ \dfrac{v_1 (\delta )}{{v}'_1 (\delta )} = V_{\rm L} (\delta ) + O(\delta ^4) $

(25)

$ \frac{v_1 (\rm{\pi } - \delta )}{{v}'_1 (\rm{\pi } - \delta )} = V_{\rm R} (\delta ) + O(\delta ^4) $

(26)

其中

$ V_{\rm L} (\delta ) = \delta - \frac{{v}''_1 (0)}{2{v}'_1 (0)}\delta ^2 + \left[ {\frac{1}{2}\left(\frac{{v}''_1 (0)}{{v}'_1 (0)}\right)^2- \frac{1}{3}\frac{v_1^{(3)} (0)}{{v}'_1 (0)}} \right]\delta ^3 $

(27)

$ V_{\rm R} (\delta ) = - \delta - \frac{{v}''_1 (\rm{\pi } )}{2{v}'_1 (\rm{\pi } )}\delta ^2 + \left[{- \frac{1}{2}\left(\frac{{v}''_1 (\rm{\pi } )}{{v}'_1 (\rm{\pi } )}\right)^2 + \frac{1}{3}\frac{v_1^{(3)} (\rm{\pi } )}{{v}'_1 (\rm{\pi } )}} \right] $

(28)

求解特征值问题(18)时，边界条件(19)可替换为

$ v_1 (\delta ) = V_{\rm L} (\delta ), \quad \;{v}'_1 (\delta ) = 1 $

(29)

通过反复更新β的大小，使得靶心条件为

$ v_1 (\rm{\pi } - \delta ) - V_{\rm R} (\delta ){v}'_1 (\rm{\pi } - \delta ) = 0 $

(30)

在一定容差范围内满足，即可求解.计算中假设当$ \vert v_1 (\rm{\pi } - \delta ) - V_{\rm R} (\delta ){v}'_1 (\rm{\pi } - \delta )\vert < 10^{ - 7} $即满足要求.

2.1 橄榄形薄膜结果

针对薄膜椭球的两种轴对称形状--橄榄形和南瓜形，分别进行了计算.

首先分析橄榄状薄膜椭球中b> b_cr的情况，例如，当取a=1, b=0.8，P-V曲线的下降段除了存在均匀解外，在P∈ (0.8369, 1.1598)区间内还有梨形分岔解的存在^[28].

针对下降段的均匀解，特征值问题(18)所得特征值β随内压P的变化曲线如图 4所示，这里只画出了负的β存在的情况，其中，图中较小的特征值曲线是数值放大了10倍的结果，以便可以与较大的特征值曲线进行对照显示.

由图 4可见，对P∈ (P_min, P_max)区间的内压存在负的特征值β，说明在该区间内的均匀解是不稳定的.其中在P∈ (P_min, P_cr2)和P∈ (P_cr1, P_max)存在一个负特征值，而在P∈ (P_cr2, P_cr1)存在两个负特征值，说明在此区间均匀构型有向分岔构型转化的趋势，从而达到一个能量较低的状态. 图 5给出了特征值对应的特征函数曲线，当P=1.2时，负特征值只有一个，为-1.428，对应的特征函数曲线为关于θ=π/2的对称曲线；当P=1时，存在两个负特征值，分别为-1.444和-0.039，对应的特征曲线分别为关于θ=π/2的对称曲线和反对称曲线.

针对存在于P∈ (P_cr2, P_cr1)区间内的分岔解，对应的特征值问题(18)都有一个负的特征值存在，如图 4中的虚线所示，表明在压强控制下的分岔解同样是不稳定的.

接下来分析b < b_cr的情况，取a=1, b=0.4时，在P-V曲线的下降段只存在一个均匀解，其对应的特征值问题(18)只存在一个负的特征值β=-1.092，对应的特征函数曲线为关于θ=π/2平面的对称曲线.进一步减小b时，由于采用的数值打靶法很难求解，故这里暂不考虑此种情况对应的稳定性问题.

针对P-V曲线的两个上升段均匀解的稳定性分析，分别考虑b>b_cr和b < b_cr的情况进行计算，对应的特征值问题(18)并不存在负的特征值，如取a=1，b=0.8，P=1，两个上升段均匀解对应的最小特征值分别为β=4.757和β=11.644.因而，压强控制下橄榄形薄膜的两个上升段中的均匀解均为稳定解.

2.2 南瓜形薄膜结果

与橄榄形薄膜在内压作用下可以发生两种失稳形式不同的是，南瓜形薄膜不管几何形状是圆是扁，当载荷达到某一临界值时总是存在梨形分岔.

例如，当取a=0.8，b=1，P-V曲线下降段的均匀解和分岔解对应的特征值β随内压P的变化曲线如图 6所示，其中实线为均匀解对应的特征值随内压变化曲线，与存在梨形分岔解的橄榄形薄膜情况类似，P-V曲线下降段的均匀解在P∈ (P_min, P_max)区间存在负的特征值，表明这个解是不稳定的.其中，在P∈ (P_min, P_cr2)和P∈ (P_cr1, P_max)存在一个负特征值，P∈ (P_cr2, P_cr1)区间内存在两个负的特征值，说明在此区间均匀构型有向分岔构型转化的趋势，从而达到一个能量较低的状态. 图 6中的虚线为在P∈(P_cr2, P_cr1)区间中的分岔解对应的特征值随内压变化曲线，针对其中的每一个内压均存在一个负的特征值，说明该分岔解也是不稳定的. 图 7中分别给出了P=1.3和P=1对应的特征函数变化曲线.

当取a=0.1, b=1，均匀解和分岔解的稳定性与a=0.8, b=1一致，图 8中分别给出了P=1.6的特征值-0.828和P=1.4的特征值-1.018，-0.008 5对应的特征函数变化曲线.

同样的，对南瓜形薄膜在P-V曲线两个上升段的均匀解进行稳定性分析，如取a=0.8, b=1, P=1，两个上升段均匀解对应的最小特征值分别为β=6.414和β=13.950，进一步验证其他情况，发现南瓜形薄膜的两个上升段的均匀解均为稳定解.

综上所述，在压强控制下，在P-V曲线下降段不管只有均匀解存在还是均匀解与分岔解共存，所有的解均不稳定.在P-V曲线的两个上升段的均匀解则均为稳定解.

3 质量控制下的稳定性分析

这部分考虑理想气体质量控制下的薄膜椭球稳定性.总能量（除以2πH后）可以写为

$ \mathit{\Pi } = \int_0^{\rm{ \mathsf{ π} }} {} \left[{W({\lambda _1}, {\lambda _2})b\sin \theta \chi (\theta )} \right]{\rm{ d}}\theta + \phi (V, M) $

(31)

其中φ(V, M)为气体势能

$ \phi (V, M) =-kM\ln \frac{V}{{{V_0}}} $

(32)

其中，k为一个大于零的常数，M为气体的总质量，在质量控制情况下为固定值，V和V₀分别是当前构型和初始构型下的气体体积.内压P（为真实的内压除以初始厚度H）可以通过$ HP = - \partial (2\rm{\pi } H\phi ) / \partial V = k(2\rm{\pi } HM) / V $计算得到.

用L表示式(31)，可以写为

$ \mathit{\Pi } = \int_0^{\rm{ \mathsf{ π} }} {} L(\mathit{\boldsymbol{u}}, \mathit{\boldsymbol{u'}}){\rm{d}}\theta + \phi (V, M) + \frac{1}{2}P\int_0^{\rm{ \mathsf{ π} }} {} {r^2}z'{\rm{d}}\theta $

(33)

其中L(u, u')的表达式与式(9)中的形式相同.相应的能量二阶变分可写为

$ {\delta ^2}\mathit{\Pi } = \int_0^{\rm{ \mathsf{ π} }} {} ({\bf{v}} \cdot S{\bf{v}} + 2{\bf{v}} \cdot {\bf{R}}{\rm{ }}{{\bf{v}}^\prime } + {{\bf{v}}^\prime } \cdot {\bf{Q}}{\rm{ }}{{\bf{v}}^\prime }){\rm{ }}\int_0^{\rm{ \mathsf{ π} }} {} {\rm{d}}\theta + {\phi _{VV}}{{\dot V}^2} $

(34)

其中$ \phi _{VV} = \partial ^2\phi / \partial V^2\vert _{V = \bar {V}} = \overline P / (2\rm{\pi } \bar {V}) $，V为平衡状态u=u时的气体体积，另外

$ \dot V = \pi \int_0^{\rm{\pi }} {} ({\rm{2}}\bar r \cdot \bar z'{v_{\rm{1}}}{\rm{ + }}{{\bar r}^{\rm{2}}}{{v'}_{\rm{2}}}){\rm{ d}}\theta $

(35)

将式(31)代入式(34)右端，式(34)可改写为

$ \begin{array}{l} {\delta ^2}\mathit{\Pi }-{\phi _{VV}}{{\dot V}^2} = \int_0^{\rm{ \mathsf{ π} }} {} ({a_1}v_1^2 + 2{a_2}{v_1}{{v'}_2} + \\ {a_3}v_2^{'2} + {a_4}v_1^{'2} + 2{a_5}{{v'}_1}{{v'}_2}){\rm{ d}}\theta \end{array} $

(36)

其中，系数a₁, a₂, a₃, a₄, a₅的表达式与式(15)中定义相同.注意到这里φVV > 0，质量控制下的能量二阶变分δ²Π始终大于压强控制下的能量二阶变分，那么在前一节得到的不稳定的解，在质量控制下，可能变为稳定解.而在压强控制下的稳定解，如在P-V曲线两个上升段的均匀解，在质量控制下仍然稳定.下文中重点分析P-V曲线下降段的解的稳定性.

与前一节类似，为了分析式(34)中δ²Π最小值的正负号，结合归一化条件(13)，可得特征值问题

$ [b_2 {v}'_1]^\prime + b_0 v_1 +b_1 \mathfrak{L}[v_1] = 0, \quad 0\le \theta\le\rm{\pi } $

(37)

$ v_1 (0) = v_1 (\rm{\pi } ) = 0 $

(38)

其中

$ b_1 = \frac{a_2 \bar {r}^{\, 2}}{2(a_3 - \beta )} - \frac{1}{2}\left(\frac{a_5 \bar {r}^{\, 2}}{a_3 - \beta }\right)' - \bar {r} \cdot {\bar {z}}' $

(39)

$ \mathfrak{L}[v_1] = \frac{\bar {P}\dot {V}}{\bar {V}} $

(40)

式中b₀和b₂的表达式与式(20)中的定义相同.把$\dot {V} $中的v'₂消去之后，$ \mathfrak{L}[v_1] $可以写为

$ \mathfrak{L}[{v_1}] = -2{\left( {\frac{{\bar V}}{{{\rm{\pi \bar P}}}} + \frac{1}{2}\int_0^{\rm{\pi }} {} \frac{{{{\bar r}^{{\kern 1pt} 4}}}}{{{a_3} -\beta }}{\rm{d}}\theta } \right)^{ -1}} \cdot \int_0^{\rm{\pi }} {} {b_1}{v_1}{\rm{d}}\theta $

(41)

下面分两种情况来考虑特征值问题(37)式的求解.

情况一： $\mathfrak{L}[v_1] = 0 $

这种情况下，特征值问题(37)与压强控制情况下的问题(18)相同，但是只有同时满足$\mathfrak{L}[v_1] = 0 $的解才是当前质量控制下的解.注意到，u=u对应为均匀解时，b₁关于θ=π/2平面对称，所以在压强控制下的关于θ=π/2的反对称模态将自动满足$\mathfrak{L}[v_1] = 0 $.

在前一节的结果中，如图 5（b）、图 7（b）、图 8（b），橄榄形椭球中b>b_cr时P∈ (P_cr2, P_cr1)区间均匀解和南瓜形椭球中P∈ (P_cr2, P_cr1)区间均匀解所对应的特征函数曲线v₁均为关于θ=π/2平面的反对称曲线，满足$\mathfrak{L}[v_1] = 0 $，所以在质量控制下，这两种情况下的解与压强控制下对应均匀解的稳定性一致，均为不稳定解.这与Fu等^[16]在针对球形薄膜稳定性分析中结果一致，当有梨形解和均匀解同时存在时，一般梨形解的能量比均匀解的能量更低.其他情况下解的稳定性将在情况二中讨论.

情况二： $ \mathfrak{L}[v_1] \ne 0 $

此时需要求解特征值问题(37)的解.注意到该特征值问题为线性问题，当且仅当设$\mathfrak{L}[v_1] = 1 $时存在非平凡解，则满足条件(41)下的原特征值问题存在非平凡解，所以，原特征值问题可以简化为一个边值问题

$ [b_2 {v}'_1]^\prime + b_0 v_1 + b_1 = 0, \quad 0 \le \theta \le \rm{\pi } $

(42)

$ v_1 (0) = v_1 (\rm{\pi } ) = 0 $

(43)

且需要满足$\mathfrak{L}[v_1] = 1 $，即

$ f({v_1}): = 2\int_0^{\rm{\pi }} {} {b_1}{v_1}{\rm{d}}\theta + \left( {\frac{{\bar V}}{{{\rm{\pi }}\bar P}} + \frac{1}{2}\int_0^{\rm{\pi }} {} \frac{{{{\bar r}^{{\kern 1pt} 4}}}}{{{a_3}-\beta }}{\rm{d}}\theta } \right) = 0 $

(44)

假设在压强控制下的特征值问题(18)和(19)的最小特征值为β_min，因为质量控制下的能量二阶变分δ²Π始终大于压强控制下的能量二阶变分，所以满足边值问题(42)和(43)的特征值β一定大于β_min，且当β小于零时，说明u=u为不稳定解；反之，如果没有负特征值或零特征值存在，那么u=u为稳定解.因此，现在只需在β=0到β=β_min区间内以很小的步长逐渐减小β的数值，将其代入边值问题(42)和(43)中进行求解，如果所求v₁满足条件(44)，说明存在负特征值β，则u=u为不稳定解；如果在(β_min, 0)区间内不存在满足要求的β，则说明没有负特征值或零特征值存在，那么u=u为稳定解.另外，注意到边值问题(42)中系数表达式(20)和(39)中均涉及到了$ \dfrac{1}{a_3 - \beta } $项，其中a₃ > 0，所以如果当满足式(44)的特征值β小于零时，不会出现分母为零的情况.

边值问题可采用打靶法实现求解，由于在θ=0, π处的奇异性，求解区间取为θ=δ到θ=π-δ.与前一节过程类似，v₁在θ→ 0和π处的渐进展开为

$ \frac{{v}''_1 (0)}{{v}'_1 (0)} = - \frac{2\beta }{3W_{11} }, \quad {\frac{{v}''_1 (\rm{\pi } )}{{v}'_1 (\rm{\pi } )} = \frac{2\beta }{3W_{11} }} $

(45)

$ \dfrac{v_1^{(3)} (0)}{{v}'_1 (0)} $与$ \dfrac{v_1^{(3)} (\rm{\pi } )}{{v}'_1 (\rm{\pi } )}$的表达式太长，在此省略.其中式(45)右边的表达式分别在u=u解中相应的θ=0和θ=π处取值.$ \dfrac{v_1 (\delta )}{{v}'_1 (\delta )} $和$ \dfrac{v_1 (\rm{\pi } - \delta )}{{v}'_1 (\rm{\pi } - \delta )} $的渐进展开表达式和前一节中式(25)和(26)形式相同.

接下来采用打靶求解，首先猜测一个v'₁(δ)的值，通过计算v₁(δ)=V_L(δ)v'₁(δ)得到v₁(δ)，将这两个值作为初始条件在θ=δ到θ=π-δ区间内进行打靶.通过调整v'₁(δ)的猜测值，在一定容差范围内满足在θ=π-δ处的边界条件

$ v_1 (\rm{\pi } - \delta ) - V_R (\delta ){v}'_1 (\rm{\pi } - \delta ) = 0 $

(46)

即可求解.以上计算过程通过Mathmatica软件实现.

3.1 橄榄形薄膜结果

首先分析b>b_cr的情况，针对下降段的均匀解，在$\mathfrak{L}[v_1] = 0 $情况下已经证明了P∈ (P_cr2, P_cr1)区间内均匀解为不稳定解，现在只需分析P∈ (P_min, P_cr2)和P∈ (P_cr1, P_max)区间内均匀解的稳定性.以a=1，b=0.8，P=1.3为例，在前一节压强控制下β_min=-1.298，图 9给出了在(β_min, 0)区间内式(44)左边的表达式f(v₁)随β的变化曲线，表明在此区间内没有满足式(44)的β存在，该解为稳定解.

针对存在于P∈ (P_cr2, P_cr1)区间内的分岔解，如取P=1时对应的f(v₁)随β的变化曲线如图 10所示，同样不存在负的特征值β使得式(42)和式(43)的解同时满足式(44)，故该分岔解为稳定解.

接下来分析b < b_cr的情况，取a=1, b=0.4时，P-V曲线下降段只有一个均匀解存在，用以上过程验证，此解为稳定解.同样的，此处不考虑进一步减小b时的稳定性问题.

3.2 南瓜形薄膜结果

在南瓜形薄膜构型下，在P-V曲线的下降段始终存在区间P∈ (P_cr2, P_cr1)，其中均匀解和分岔解共存.其中在$\mathfrak{L}[v_1] = 0 $情况下已经验证了P∈ (P_cr2, P_cr1)区间内的均匀解为不稳定解；用$\mathfrak{L}[v_1] \ne 0 $情况下的方法验证P∈ (P_min, P_cr2)和P∈ (P_cr1, P_max)区间内的均匀解和P∈ (P_cr2, P_cr1)区间内的分岔解均为稳定解.

综上所述，在质量控制下，只要有分岔解存在，均为稳定解；而没有分岔解出现的区间，均匀解为稳定解，但是当分岔解与均匀解共存时，均匀解为不稳定解.

4 结论

不同几何尺寸的薄膜椭球在一定内压载荷作用下会出现不同形式的分岔现象，根据P-V曲线发现会有同一个载荷对应多解的情况出现，这些解是否稳定则需要通过稳定性的分析进行判断.能量判据是连续介质力学中广泛用于分析保守系统稳定性的判断准则，通过比较某一平衡状态及施加小扰动之后状态的能量来判断当前状态是否稳定，而这两个状态的能量差可通过判断系统能量二阶变分的正负号实现.

本文分别针对压强控制和质量控制两种加载方式下，不同几何形状的薄膜椭球在P-V曲线下降段对应的均匀解和分岔解进行了稳定性分析.结果表明，在压强控制下，下降段的均匀解和分岔解均为不稳定解；但在质量控制下，均匀解在没有分岔现象发生的区间为稳定解，一旦出现分岔现象，则分岔解为稳定解，均匀解则不稳定.另外，P-V曲线两个上升段的均匀解，则在压强控制或者质量控制下均为稳定解.

致谢 英国基尔大学（Keele University）傅依斌教授对本工作思路上给予了指导，并对文章的审阅提出了重要修改意见，在此致以衷心的感谢！