闫 棣， 苏 祺， 李四平

(上海交通大学 船舶海洋与建筑工程学院， 上海 200240)

固体力学研究领域中，工程结构的非线性屈曲问题一直是非常活跃的话题[1].在使用通用有限元软件进行屈曲模拟的过程中，一些文献[2-3]直接使用软件中提供的线性屈曲分析模块进行屈曲临界载荷的预测与屈曲模态的判断.线性屈曲分析方法相当于弹性屈曲分析，然而材料的非线性、结构形状的初始缺陷以及几何大变形使得实际结构并不在理论弹性范围内屈曲.因此，特征值屈曲分析所得的解为非保守解，一般不能用于实际工程的分析[4-5].除此之外，线性屈曲分析方法对材料非线性等问题无法解决，有一定的局限性.一些文献[6-7]使用线性屈曲分析中所得的一阶模态或前几阶模态进行缺陷的引入，而在引入缺陷时所用的比例因子参数没有定论.这种方法通过人为调控比例因子的大小可能得到不同的结果，计算结果仍需要通过进一步的试验比对来调整参数，而且在线性屈曲方法无法计算的材料非线性等问题上[8]，继续使用此方法进行屈曲的分析，所得结论是不可信的.

本文提出一种更为普适的模拟方法，即通过随机的方法对理想结构模型引入初始缺陷进行屈曲的分析.这种方法不仅能够得到上述文献中所用方法计算得到的多阶屈曲模态和其对应临界力，而且将特征值问题转化为几何非线性静力变形问题进行计算，克服了线性屈曲方法对材料非线性等问题的局限性，获得了更多中间加载过程的信息.此外，以弹性地基上梁的屈曲问题和横向纤维作用下穿透型矩形脱层屈曲问题为例，计算讨论随机缺陷法的普适性和以上文献中所使用方法的局限性.同时，通过对比不同大小缺陷对计算结果的影响，提出合理的随机缺陷引入方法.

1 随机缺陷法

1.1 理论分析

初始缺陷是指与理论计算模型相比，实际构件在尚未受到载荷作用前就已经存在的各种缺陷，其中包含力学初始缺陷以及几何初始缺陷.力学初始缺陷主要是由各种原因引起的初始残余应力；而几何初始缺陷则包含整体初弯曲，初偏心，局部初始缺陷，等等.

考虑如图1所示的两端铰支的细长杆，当它处于微弯曲平衡的临界状态时，若杆内应力不超过材料的比例极限，则压杆任一截面m-n处的弯曲变形与截面弯矩的关系为

(1)

式中：E为杆件弹性模量；I为杆件截面惯性矩；M(x)为杆件在距离下支座x处的由外力F引起的弯矩.

图1 两端铰支细长杆Fig.1 Schematic of hinged-hinged slender beam

考虑模型有初始缺陷，由于初始缺陷是一个随机的无穷小量，忽略高阶小量，做叠加处理得

(2)

式中：M0(x)为初始缺陷引起的初弯矩，是一个随机的非连续小量.由于M(x)=-Fy，所以式(2)可写为

(3)

(4)

解二阶非线性常系数微分方程，其通解及特解分别为

Y=asinkx+bcoskx

(5)

(6)

式中：D为微分算子.代入杆端边界条件(x=0时，y=0，M0(0)=0；x=l时，y=0，M0(l)=0)，可以得到式(4)的解为

(7)

(8)

即当M0(x)趋近于0时，式(7)同式(1)的解等价，表示当加入的随机缺陷足够小时并不会影响模型整体的屈曲结果.

然而理想不含初始缺陷的构件在外力作用下并不会发生分支点屈曲，分支点屈曲发生的原因也正是由于实际构件中包含一定的初始缺陷.在通用有限元软件线性屈曲分析中，将屈曲问题作为特征值问题进行求解.而随机缺陷法，通过随机的方法在理想构件中引入缺陷将屈曲作为几何非线性静力计算问题进行求解，这更符合实际情况.

在对模型施加缺陷时，随机缺陷法强调随机性，充分模拟实际构件中缺陷发生的随机性.缺陷可以是构件的几何缺陷，也可以是力学初始缺陷，引入这些小缺陷的目的就是触发构件屈曲.缺陷大小和随机性是计算结果精度的主要影响因素.

1.2 有限元实现方法的探讨

本文算例中采用随机选取模型网格节点进行初始偏移的方法随机地添加缺陷，根据计算结果可以得知偏移的大小取梁截面径向尺寸的千分之一或更小为宜.这是一种简单方便的几何缺陷引入方法.除此之外，本文也尝试通过程序随机选取结构中某一微小单元并赋予较小的模量来模拟材料的空心几何缺陷.力学缺陷上，可以使用虚拟的触发载荷来启动不稳定性，也可以通过不均匀的温度场致使屈曲的发生.

在保证随机性的前提下，引入初始缺陷的方法有很多，其根本的目的是为了触发屈曲.初始缺陷越小，有限元模型越接近理想模型，计算结果与理论结果也越接近，但此时为触发屈曲，模型的计算量也会相应增加.所以有限元计算中触发屈曲时，初始缺陷的大小、模型的计算量以及计算结果的精度需要有一个平衡的过程.

2 两端铰支细长杆多阶屈曲模拟

设有长l=100 cm，截面面积A=4 cm2，惯性矩Ix=Iy=1.33 cm4的两端铰支细长杆件，长细比λ=173.2，如图2所示.材料的弹性模量E=200 GPa.

图2 两端铰支压杆Fig.2 Schematic of hinged-hinged compression bar

使用通用有限元软件进行建模，采用梁单元进行模拟，并将梁沿轴向划分为100个单元.通过线性屈曲方法得到前5阶屈曲模态及临界载荷，结果以第1阶和第2阶模态为例，如图3所示，其中u2表示杆件径向位移.参照Euler公式

(9)

式中：Fcr为临界力；n为屈曲模态数.线性屈曲方法模拟结果与Euler公式理论结果(Fcr1=26.32 kN，Fcr2=105.28 kN)较为吻合.

图3 线性屈曲方法所得压杆屈曲前两阶模态Fig.3 The first two order modes of compression bar by linear buckling method

图4 随机缺陷法计算压杆屈曲第1阶模态Fig.4 The first mode of compression bar by stochastic defect method

使用随机缺陷法进行模拟.通过编写程序操作计算文件，随机选取杆件上节点施加10-5R0(R0为横截面径向尺寸)的横向偏移，加载位移载荷，使其从0开始缓慢增长，其第1阶屈曲模态如图4所示.

为保证随机性不影响最终结果，进行了多次计算，每次选取的偏移节点均由程序随机生成，所得结果相同.在使用节点偏移法的同时，尝试使用空心几何缺陷，通过实体单元进行建模，随机选择模型中某一位置划分出微小的单元块，并赋予较低的材料模量来模拟空心，载荷施加不变，所得到的结果与前述方法相同.

在第1阶模态基础上，第2次计算时不改变随机缺陷位置，使位移载荷在加载开始时跳过第1阶的临界位移载荷，而后缓慢增长，寻找第2阶模态的临界载荷，计算结果如图5所示.

图5 随机缺陷法计算压杆屈曲第2阶模态Fig.5 The second mode of compression bar by stochastic defect method

通过以上算例可以看出，随机缺陷法不仅可以计算模型屈曲的最低阶模态，而且可以通过加载参数的改变来计算高阶模态，计算结果与理论结果较为吻合.

3 弹性地基上梁的屈曲

3.1 嵌入于地基中的梁

将上文中的杆件嵌入于弹性地基中.地基使用Winkler地基模型进行模拟，弹簧劲度系数β=20 MN/m，模型如图6所示.

图6 Winkler弹性地基梁Fig.6 Schematic of Winkler foundation model

使用通用有限元软件进行建模，采用梁单元进行模拟，并将梁沿轴向划分为100个单元，在除两端外的每个网格节点处添加线性弹簧单元[9-10].

使用线性屈曲方法进行计算，在梁的一端施加单位位移载荷，所得结果如图7所示.计算特征值为7.655 86×10-2，所以临界应变εcr1=7.655 86×10-2.

3.2 搁置于地基表面的梁

在以上模型的基础之上，加入非线性弹簧，即模拟梁未嵌入地基，而是放置在地基表面上的情况.此时，梁在向上弯曲时，与地基表面脱离；而向下弯曲时，则受到来自地基的压力.模型相当于两端铰支的压杆稳定问题.在建模过程中弹簧压缩时的劲度系数β=20 MN/m，而拉伸时劲度系数趋近于0.

使用线性屈曲方法进行计算，发现计算结果与使用线性弹簧时的结果相同.线性屈曲方法会选取弹簧初始刚度或者上一分析步结束时的刚度作为线性刚度进行计算，而忽略后续弹簧参数的变化[11]，即不能计算非线性弹簧或者非线性材料等问题.而在此基础之上通过影响因子的方式将线性屈曲方法的计算结果作为缺陷引入到模型中再进行静力计算，显然也有问题，其结果的有效性有待验证.同时可以认为，对于材料非线性等问题，使用此方法仍然会失效.

使用随机缺陷法进行计算，同样地随机选取梁上节点引入10-5R0的横向偏移，其计算结果如图9与10所示.随机缺陷法计算所得的临界载荷Fcr=26.31 kN，与Euler公式的计算结果吻合， 说明随机缺陷法对非线性弹簧问题同样适用，相比于线性屈曲方法有较强的普适性.

考虑不同大小初始缺陷对最终临界载荷Fcr计算结果的影响.对模型分别引入不同大小初始缺陷，所得计算结果如表1所示.其中:δ=R/R0，R为初始缺陷.

图9 搁置地基表面梁受压端轴向位移与整体支反力关系Fig.9 Relation of axial displacement of foundation surface beam compression end and integral bearing reaction

图10 搁置地基表面梁中点处横向位移与整体支反力关系Fig.10 Relation of transverse displacement at midpoint of foundation surface beam and integral bearing reaction

表1 杆初始缺陷与临界载荷关系Tab.1 Relation table of stochastic defect and critical load

通过对比分析可以发现，当δ=10-3或更小时，所得计算结果的精度较高.

4 横向纤维搭桥下穿透型矩形脱层屈曲

设有长2a=200 mm，宽2b=100 mm的穿透型矩形脱层如图11所示，受到两侧均布载荷Px的作用.脱层的厚度h=1 mm，在横向纤维的影响下，受到均布载荷q(x)的作用抑制屈曲.材料的弹性模量E=2.3 GPa，泊松比μ=0.3.文献[12]中也进行了类似的计算，其所使用的方法是首先使用线性屈曲模块，然后将1阶屈曲模态位移乘以比例因子来进行初始缺陷的引入，最终计算结果与理论结果没有很好的对应关系.

图11 横向纤维搭桥下矩形脱层Fig.11 Schematic of rectangular delamination bridged by fibers

根据文献[13]所述理论进行分析，

q(x)=βω(x),β>0

(10)

(11)

(12)

式中：q(x)为纤维约束力；β为纤维假设为弹簧时的劲度系数；ω(x)为屈曲挠度；D为抗弯刚度；a0为屈曲时特征长度.

根据表2中计算结果与理论值的对比可以发现，使用随机缺陷法所得到的计算结果有较高的可信度.这种方法通过随机地改变初始缺陷，可充分模拟实际构件的缺陷形式，降低人为因素对计算结果的干扰，对工程应用中分支点屈曲问题的模拟有较高的参考价值.

图12 屈曲模态图Fig.12 Schematic of buckling mode

5 结论

通过使用线性屈曲方法和随机缺陷法分别对压杆稳定、弹性地基上梁的屈曲以及横向纤维搭桥下穿透型矩形脱层屈曲问题进行了有限元模拟.得到以下几点结论：

(1) 随机缺陷法将屈曲作为几何非线性静力计算问题进行有限元模拟，改进了线性屈曲计算中将屈曲作为特征值问题进行求解的不足，与实际情况更为相符.

(2) 讨论了随机缺陷法的理论依据，并给出了通用有限元软件使用初始随机缺陷法的几种思路.

(3) 随机缺陷法具有更强的普适性，解决了线性屈曲分析无法计算非线性弹簧等问题.同时，随机缺陷法同样可以计算模型高阶的屈曲模态.

(4) 通过对比不同大小初始随机缺陷的计算结果，发现初始随机缺陷取构件截面径向尺寸千分之一或更小的值时，有限元计算结果与实际情况较为吻合.