第二节　槽身局部稳定性分析

一、柱壳有限条元的位移模式

采用柱壳有限条元法分析U形薄壳渡槽结构的稳定性问题，首先要划分单元，用一系列过圆柱轴线的平面把圆柱壳结构切成有限个柱壳条元，条元之间由节点与线相连。如图2-1所示，其单元e的曲率半径为R，壳体厚度为t，条元宽度弧长为B。结构失稳时，设单元中曲面内任一点的挠度为W。为了描述结构的失稳位形，需要由单元节线位移函数插值构成，于是设单元e的节线为i、j，而单元节线i、j与柱壳横截面的交点i、j作为单元e的节点。其节点位移用Wi、αi、Wj、和αj来表示，写成向量形式为

式中　Wi、Wj——单元节点i、j的挠度；

θi、θj——单元节点i、j绕x坐标轴的转角。

由弹性力学中极坐标系下的几何方程，可知

式中　β——圆心角，其变化范围为β∈［-α，α］。

位移函数选择为半解析型，沿y轴方向（周向）的径向位移选用Y（η）函数，不仅要求径向位移连续，而且要求y方向的一阶导数连续，这样才能保证失稳屈曲后的管壳与yz平面的交线为光滑曲线，也就是转角位移连续。因此Y（η）可以用三阶的插值多项式函数向量［Φ（η）］与节点位移向量{δ}的乘积来表示，其具体表达式如下

沿x轴方向（轴向）的径向位移选用解析函数X（ξ），其表达式为

其中　ξ=x/L

式中　ξ——无因次坐标，ξ∈［-1/2，1/2］；

L——待定的失稳半波长。

由此可见，计算模型选取了一系列失稳波形中的一典型失稳波段，并考虑到计算方便，纵向坐标（即x坐标或ξ坐标）原点取至波峰处。这种半解析位移函数的选取有利于减少结构自由度，提高计算精度。

柱壳条元任意一点的径向位移可用如下的位移函数来描述：

式中　ξ——x轴方向无因次坐标，ξ∈［-1/2，1/2］；

η——y轴方向无因次坐标，η∈［-1，1］。

二、柱壳有限条元的刚度矩阵

1.几何方程

根据壳体分析的诺瓦日洛夫理论，考虑到薄壳的大挠度应变-位移关系，略去高阶小量，可得柱壳有限条元的几何方程

式中　{εl}——柱壳弯曲的变形分量矩阵，它与z坐标有关，中曲面上其值为零，在壳体的上下表面其值为最大；

{εN}——由大挠度所引起的非线性应变，它沿壳体厚度方向是均匀分布的。

对于大挠度问题，单元内任意一点的应变应当包括这两部分的叠加，因此，式（2-8）可简写为

在后面的公式推导中，为了使运算清晰简洁，在此约定了下面一些符号：

因此，式（2-8）就可以表示成

式中　“*”号——两矩阵对应元素相乘形成一个同阶矩阵的运算。

2.物理方程

在壳体的稳定性分析中，认为材料是线弹性的，并考虑壳体在横向荷载作用下所引起的正应力σz与平行于中曲面上的正应力相比是可以忽略的，因此它具有和平面应力问题一样的应力-应变关系，可写成下面的矩阵形式

式中　［D］——弹性矩阵；

E——弹性模量；

μ——泊松比。

3.柱壳有限条元弹性变形能公式

根据能量原理，可知单元的应变能公式为

式中　U——单元体的应变能；

V——单元体的体积。

现在把式（2-9）、式（2-14a）代入式（2-15），可得

式（2-16）等号右侧的第三项积分是非线性应变的二次项，它是高阶微量，可以忽略。则应变能公式主要由两部分组成，即

式中　Ul——单元的线应变能；

UN——非线性应变引起的应变能。

下面分别对两个积分进行计算。

4.柱壳有限条元弹性刚度矩阵的推导

在式（2-17）中，首先考虑第一个积分Ul，把式（2-12a）代入式（2-18a），并整理得

则柱壳有限条元的弹性刚度矩阵［kE］为

把式（2-14b）代入式（2-20）中，可得

然后把式（2-11a）、式（2-11b）和式（2-11c）代入式（2-21），并沿周向、法向积分，可得

在式（2-24）中右端的四项分别记为［kE1］、［kE2］、［kE3］和［kE4］，则此式可简写为

要计算［kE1］、［kE2］、［kE3］和［kE4］，分别把多项式向量［Φ］、［Φ′］和［Φ″］的表达式（2-4）、式（2-10a）和式（2-10b）代入式（2-25），并进行矩阵乘法计算，然后再对每个元素积分，其计算过程如下：

把式（2-26a）～式（2-26d）四式代入式（2-24），可得到［kE1］、［kE2］、［kE3］和［kE4］如下式所示：

把式（2-27a）～式（2-27d）四式代入式（2-25），就得到了柱壳有限条元的弹性刚度矩阵［kE］。

5.柱壳有限条元几何刚度矩阵的推导

再考虑单元应变能表达式中的第二项，即式（2-18b），并做如下变换：

将式（2-12b）代入式（2-28），可得

如果把列阵σ｛｝改写成对角阵，如

因此，式（2-29）就可以写成如下形式

将式（2-31）展开，可得

由此可知，柱壳有限条元的几何刚度矩阵［kG］可分为［kGx］、［kGy］和［kGxy］三部分，下面将分别进行计算。

将式（2-11d）和式（2-11e）代入式（2-32），并沿柱壳轴向、径向和周向积分可得

由此可以看出［kGxy］为一个零矩阵。［kG］称为单元的几何刚度矩阵，其中［kGx］是与轴向压力有关的几何刚度矩阵；而［kGy］是与周向压力有关的几何刚度矩阵，把式（2-34b）与式（2-34c）对η积分可得

由此可知，柱壳有限条元的几何刚度矩阵［kG］=［kGx］+［kGy］，对于渡槽结构只考虑轴向应力σξ，因此，矩阵［kGy］为零矩阵，渡槽结构有限条元的几何刚度矩阵［kG］=［kGx］。

至此，就推导出了U形薄壳渡槽局部稳定性分析的有限条元的弹性刚度矩阵和几何刚度矩阵。

三、U形槽局部稳定性的计算

1.渡槽截面应力分布特征

要想求得渡槽局部稳定性临界应力，先要清楚渡槽截面的应力分布特征，渡槽上任一点的应力可表示为

2.有限条单元刚度方程

设单元的等效节点荷载向量为｛P｝，由卡氏第一定理可得

把式（2-19）、式（2-32）代入上式，并考虑矩阵［kE］、［kG］的对称性，可得

这就是柱壳条元的刚度方程，由于单元的几何刚度矩阵［kG］中含有应力向量，可知柱壳元刚度方程是一个非线性方程。

3.渡槽结构的总体刚度矩阵

利用刚度集成法形成结构的总体弹性刚度矩阵［KE］和总体几何刚度矩阵［KG］，其叠加公式分别为

式中　［KE］e、［KG］e——与结构总体刚度矩阵同阶扩大的单元弹性刚度矩阵和单元几何刚度矩阵。e表示单元编号，结构刚度矩阵中的任意一个元素为

式中　脚标r、s——刚度矩阵的第r行和第s列。

4.渡槽结构局部稳定性问题的特征方程

形成结构的总体刚度矩阵以后，就可以写出结构的刚度方程

式中　｛δ｝——结构的节点位移向量；

｛P｝——结构的节点荷载向量。

在结构刚度方程中，引入约束条件，在结构刚度矩阵中划掉与约束位移相关的行和列，就得到一个新的方程。为简单记，此处仍然用式（2-42）表示修改后的刚度方程。

由于式（2-42）中的几何刚度矩阵中含有应力项，因此刚度方程是一个非线性方程。如果把应力项从几何刚度矩阵中提出，写成应力项的显式，则式（2-42）就可改写成

这样矩阵［KG］*就不再具有“刚度”的量纲，在此称它为当量刚度矩阵，应力σξ就是这个当量刚度矩阵的比例因子。

柱壳结构在外压作用下，在初始的平衡位置壳体处于无矩状态，但是当柱壳结构一旦丧失稳定性就会突然发生巨大的弯曲变形。此时，描述柱壳结构平衡的方程式（2-43）是不稳定的，即矩阵（［KE］+σξ［KG］*）奇异，也就是

这就是描述柱壳结构局部稳定性问题的特征方程，式（2-44）可改写成

式中　［I］——单位矩阵。

这样就可以把求解柱壳结构局部稳定性问题转化为求实矩阵（-［KE］-1［KG］*）的最大特征根（1/σξ）max问题。

在此，设［A］=-［KE］-1［KG］*，λ=1/σξ。则式（2-45）化成了标准的特征值问题

然后，就可以利用幂法求解特征方程的最大特征根λmax。其倒数（1/λmax）=σξ，min，就是所求的柱壳结构最小临界荷载σcr。

四、U形槽局部屈曲规律分析

1.上缘宽度对渡槽局部屈曲的影响

设U形薄壳渡槽上缘厚度、竖直段厚度和半圆段厚度都为t1=t2=0.2m，竖直段高度h=3.0m，半圆段平均半径R=1.5m。此时渡槽的高宽比H/B为1.5，现分析渡槽上缘宽度在0.2～1.2m范围内变化下，其局部稳定性分析结果参见表2-1，渡槽上缘宽度与局部屈曲临界应力关系曲线如图2-2所示。

2.高宽比对渡槽侧局部屈曲的影响

设U形薄壳渡槽上缘厚度、竖直段厚度和半圆段厚度都为t1=t2=0.2m，上缘宽度b=1.0m，半圆段平均半径R=1.5m。现分析渡槽高宽比在1.0～2.0范围内变化下，其侧局部定性分析结果参见表2-2，渡槽高宽比与局部屈曲临界应力关系曲线如图2-3所示。

从表2-1和图2-2可以看出，随着渡槽上缘宽度的增大，局部屈曲临界应力减小，其中渡槽上缘宽度对运行工况的临界应力影响较大，并且渡槽上缘宽度对其局部屈曲的半波长影响很小。

从表2-2和图2-3可以看出，随着渡槽高宽比的增大，局部屈曲临界应力减小，并且渡槽运营工况和施工工况的局部屈曲临界应力变化趋势大致相同。

3.壁厚对渡槽局部屈曲的影响

设U形薄壳渡槽上缘宽度b=1.0m，竖直段高度h=3.0m，半圆段平均半径R=1.5m。此时渡槽的高宽比H/B为1.5，现分析渡槽上缘厚度t1、竖直段厚度和半圆段厚度t2在0.10～0.30m范围内变化时，其局部稳定性分析结果参见表2-3，渡槽壁厚与局部屈曲临界应力关系曲线如图2-4所示。

4.不同壁厚、高厚比时渡槽局部屈曲分析

设U形薄壳渡槽上缘宽度与槽宽之比b/B=0.3，高宽比H/B为1.5。现分析渡槽壁厚t1=t2=t和高厚比H/t2变化时，其局部稳定性分析结果参见表2-4，渡槽壁厚、高厚比与局部屈曲临界应力关系曲线如图2-5所示。

从表2-3和图2-4可以看出，随着渡槽上缘厚度的增大，其局部屈曲临界应力逐渐减小，并且变化得很平缓。随着渡槽竖直段和半圆段厚度的增大，其局部屈曲临界应力增大，并且变化得比较显著。但是渡槽壁厚对其局部屈曲的半波长影响很小。

图2-5是根据表2-4中的壁厚t1=t2=t=0.2m，在不同高厚比下渡槽局部屈曲临界应力做出的曲线。从图2-5可以看出，随着渡槽高厚比的增大，其局部屈曲临界应力减小，并且高厚比在10～30范围内，局部屈曲临界应力变化得较显著，高厚比在30～50范围内，局部屈曲临界应力变化的趋于平缓。从表2-4可以看出，渡槽壁厚取0.10m、0.20m和0.30m的不同值时，只要渡槽高厚比相同，其局部屈曲临界应力非常接近。

综上所述，可得到以下一些结论：

（1）在渡槽的截面尺寸中，渡槽上缘宽度对运营工况的临界应力影响较大，并且随着渡槽上缘宽度的增大，局部屈曲临界应力减小。但上缘宽度对其局部屈曲的半波长影响很小。

（2）在渡槽的壁厚中，渡槽侧墙厚度对其局部屈曲临界应力影响最显著，而渡槽上缘厚度对其局部屈曲临界应力影响最小，并且随着渡槽侧墙厚度和底板厚度的增大，其局部屈曲临界应力增大。但是渡槽壁厚对其局部屈曲的半波长影响很小。

（3）随着渡槽高宽比、高厚比的增大，其局部屈曲临界应力减小，高厚比在10～30范围内，局部屈曲临界应力变化得较显著。并且在高厚比确定的情况下，即使渡槽壁厚变化，其局部屈曲临界应力非常接近。

（4）影响渡槽局部屈曲临界应力的主要因素有渡槽上缘宽度、壁厚、高宽比和高厚比等，其中对局部屈曲临界应力影响最大的因素是渡槽高厚比。

五、沙河渡槽局部稳定性分析

根据上述半解析壳有限条元法的理论和方法编制了薄壳结构稳定性分析的专用计算机软件（WDF2），沙河梁式渡槽的局部稳定性问题采用WDF2软件进行分析计算。

1.槽顶缘板、侧壁板局部稳定性分析

采用WDF2软件计算出沙河梁式渡槽的失稳半波长L=5.30m，失稳波形如图2-6所示。由图中可以明显看出，失稳首先发生在侧壁板处。求得的临界失稳应力σcr=93.86MPa，该值为混凝土抗压强度标准值（27.0MPa）的3.5倍。

根据有限元分析结果，静力荷载作用下槽顶侧壁板的最大压应力为3.664MPa，动力荷载作用下槽顶侧壁板的最大压应力为5.037MPa，可见结构失稳临界应力远远超出了在使用荷载作用下该处的最大压应力，U形薄壳渡槽上缘板及侧壁板处不会发生局部失稳屈曲破坏。

2.槽底局部稳定性分析

由WDF2软件计算出的槽底失稳半波长L=17.5m，失稳波形如图2-7所示。由图中可以看出，结构局部失稳的部位发生在槽底，这时的失稳临界应力值σcr=111.0MPa，该值为混凝土抗压强度标准值的4.1倍。根据有限元分析结果，静力荷载作用下槽底的最大压应力为1.929MPa，可见该工况下，槽身底部有较大的稳定性安全贮备，槽底在预应力施加过程中不会发生局部失稳，结构是安全的。