基于Grobner基的圆薄板非线性大挠度分析

由于重量轻，易加工，又可承受较大的载荷，卞结构在工程中使用广泛。为了准确地评估其工作性能，常常需要考察卞在载荷作用下的力学行为，圆卞非线性大挠度弯曲就是其中的-个重要问题 。

圆卞大挠度弯曲问题有多种数学模型，其中卡门 (Von Karman)方程组由于形式相对简单，又合理地考虑了非线性效应，因而得到广泛采纳。该方程组包括-对耦合的非线性偏微分方程，各国学者对其求解方法进行了大量研究[1-5]0Vincent 以载荷为摄动参数，首先求解了均布载荷作用下的圆板大挠度弯曲问题，结果的误差较大。钱伟长 I3 以中心位移为摄动参数求解了同样的问题，计算结果非秤近实验值，并且求解过程比Vincent的方法收敛更快。叶开沅等改进了摄动法的计算程序，提出了修正迭代法4，使得各阶摄动方程不再耦合，计算过程得到较大简化。上述方法均需要多次迭代求解微分方程组，以逐步提高近似解的精度，这要求工程设计人员有较高的数学分析和运算能力，不便于工程应用。

变分法不直接求解卡门方程组，而是将大挠度弯曲问题转化为等价的弹性势能极值问题，通过选取适 当的挠曲面函数，得到-个代数方程组，进而求出近似解，求解过程的数学分析较为简便。但是，当挠曲面函数待定系数增多时，将形成多元高次非线性代数方程组。常用的数值方法 ，/lNewton.raphson迭代法，需要给出好的初始值才能保证非线性方程组的收敛。而解析法中的GrObner基法 ，通过将-组非线性多项式方程化简为-个等价的三角化方程组，可求得封闭形式的解析解，不会产生增根或少根，能够最大限度地满足方程组的求解精度，克服了数值方法的不足。

本文根据能量原理，用里茨法推出等厚度轴对称圆卞非线性弯曲问题的变分列式，结合符号计算软件Mathematica，应用GrObner基法求解了非线性方程组，得到了圆卞受均布载荷时的最大挠度♂果表明，采用Gr6bner基法求解大挠度问题，计算效率较高，且能求得高精度的近似解。

1 里茨法变分列式在轴对称分布载荷作用下，圆卞大挠度弯曲时，板的形变势能 包括弯曲形变势能 和中面应变势能 ：[8： 1 (1)。J1 p( 吉(等] dw dw]却 (2) f等 ) 等 c3，其 中， p为极坐标；W为板的中面挠度方收稿日期：2012-12-22基金项目：中央高校基本科研业务费专项资金项目(2009B30014)作者简介：张迅炜 (1972-)，男，江苏常州人，讲师，从事机械结构设计与分析研究。

82 第35卷 第3期 2013-03(下) I 匐 似程： 为中面内点的径向位移； 为泊松比；弯曲刚度DE8 /[12(1- )]，式中：6为板厚，E为弹性模量。

周边固支时，板的位移边界条件及轴对称条件分另 为(Up：： ：W)pa-0 0( )。； ：。 ) JMathematica等)中该算法及其改进版本均已实现，可直接以函数方式调用。

本文用改进的Buchberger算法对式(8)、(9)组成的非线性方程组进行求解，步骤为：1)假设n，得到关于 与Cm的非线性方程组(4) G；)p 。㈡ 。 (5)应用里茨法求解时，要求假设的位移函数同时满足式 (2)、 (3)，因此将中面内各点的径向与轴向位移分别取为：( - (w Cmf - a2 )mO (6)(7)其中， 、 为互不依赖的待定系数；n为非负整数。

根据最小势能原理，圆卞轴对称弯曲问题的解应满足下列二组方程：盟 ：0 (8)：2n Iqf1- r (9) 0L- a其中，q为载荷集度假设式(6)、(7)中的参数值n，并代入式(8)、(9)，得到-个关于 与Cm的非线性方程组，它包括2n1个方程～该方程组的解代入式(6)、(7)即得到板中面各点的径向与轴向位移的近似解，据此可进-步求出应力、内力等其他未知量。随着参数n取值的增大，待定系数的数量相应增加，各级近似解将不断趋近于精确解。

2 Gr0bner基法简介与求解步骤Gr6bner基法 是-种建立非线性多项式系统标准基的理论和方法。与求解线性方程组的高斯消去法思路相类似，Gr6bner基法将非线性系统转化为等价的三角形方程组，从而可用回代的方法求解 。该方法在原非线性多项式系统所构成的多项式环 内，通过对多项式变元的适 当排序，求多项式系统中多项式对的S-多项式，并进行彼此的消元 ，最终生成-个与原 系统完全等价的标准基 。Buchberger算法是构造多项式系统Gr6bner基的-种有效算法，在计算代数系统 (如2)将G中所有多项式赋于H；3) 确 定H中的 变 元 的字 典 法 排 序 ，如c0 P P P， P P Alp 。。

4)用GrObner基函数将H化为最简Gr6bner标准基，进而得到与原多项式方程组等价的三角化方程组 ：g。(c0)0g：(Co，G)0g (Co，C1，， ， ，4 。)05)求解其中-元多次方程得到C。，回代方程逐个求出其它变元。

3 算例与程序设等厚度轴对称圆板，半径为 a2m，厚度为60.03m，受均布载荷q40kN/m 作用，周边固 支 ，材 料 的弹 性 模 量 E200GPa，泊松 比 0.3，如图1所示。求圆板的最大挠度。

g 上、 ， 上D ， r口1图1周边固支轴对称圆卞由式(7)可知，P0时圆板的中心挠度也就是最大挠度 ∑ (10)设 0，式 (6)、 (7)分别表示为： 1-口P/l口O- - (Ao4导) c ，wcof -将式(11)、(12)代入式 (1) (2) (3)，得到形变势能 表达式，再代入式 (8)、 (9)得到第35卷 第3期 2013-o3(下) I83] l 匐 似非线性方程组：1.03569x10 6.21414x109Al-7.03612x10 204 10Sc-40000n/3 8 28552 10 co 1 40722 10 CO。。 (13)- . × - . ×1.13104x10 4c06.31278x10 0应用Gr6bner基法，按前述步骤求解并转换得到的等价三角化方程组为：I-0.00003764280.00186146c01. 301. -0.604365c0 (14)1.4o894o48 0由式(14)的第-个方程式可求得CO的三个根。

其中唯-的实根，据式 (10)，就是 取0时板的最大挠度：WmaxCo≈0.01739m再设 1，则：1-旦a]j旦a k 4导 ]wC0f - clf - 2)类似地，求得最大挠度：WmaxCoC ≈0.01867-0.001552≈0.01712m取不同值时，求解GrObner基和最大挠度的Mathematica程序如下：( 设置参数n的值 )n1：( 例题中的已知参数值)[Mu3/10.；a2；[Delta]3/100；q40000；( 求弯曲刚度，建立近似位移方程 )FR2 10 11 [Delta]3/(12 (1-[ 2))；W(1-[Rho] 2/a 2) 2 Sum[c[i (1-[Rho] 2/a2)i，i，0，n)]；u[Rho](1-[Rho]/a) [Rho]/a Sum[A[i] ([Rho]/a) i，i，0，n1)]；( 求弯曲应变能 )V[Epsilon]1Pi FR Integrate[[Rho] (D[w，[Rho]，2)2lA[Rho] (D[w，[Rho]) 22 [Mu] D[w，[Rho]，21 Dw，[Rho]，[Rh0]，0，a；V[Epsilon]2Pi FR 12A[Delta]2 Integrate[((D[u[Rho]，[Rho]]l/2 (D[w[Rho]) 2) 2u[Rho] 2A[Rho] 22[Mu u[Rho]A[Rho](D[u[Rho]，[Rho](1/2 (D[w，[Rho])2)) Rho)[Rho，0，a)]；( 生成非线性方程组)[841 第35卷 第3期 2013-03(下)equlSimplify[Table[D[V[Epsilon]lV[Epsilon]2，A[i]]，i，0，rl1]；equ2Lef-Simplify[Table[D[V[EpsilonlV[Epsilon]2，c[i]，i，0，n)]]；equ2RigSimplify[Table[Integrate[2 Pi q (1-[Rho] 2/a2) (i2) Rho]'[Rho]，0，a)]，i，0，n]；equ2equ2Lef-equ2Rig；equSysJoin[equl，equ2]( 设定字典序 )tabReverse[Join[Table[c[i]，i，0，n]，Table[A[i]，i，0，n1]( 生成Gr6bner基 )equGroGroebnerBasis[equSys，tab( 逐个回代求解三角形方程组，-c)solSolve[equGro[1]0，c[0]；solc[0->Select[c[0]/.sol，#[Element]Reals&][[1]]；I f[n1，S O l 1S O 1 V e(e q u G r o[[2]]/。

so1)0，c[1]]，sollc[1]->0)]；( 求最大挠度 )wmax(c[0]/.so1)(c[1]/.sol1)/N为了检验GrObner基法以及Mathematica程序的求解精度，对算例中的q取不同值，求得最大挠度wmax，结果汇总于表1。

表1 不同载荷集度q下的最大挠度比较计算结果表明：1)应用里兹法和GrObner法求解等厚度圆卞最大挠度，计算过程收敛速度快 ， ：0时，与最大挠度摄动解的误差不超过2%。2)随着参数n的增大，中面各点的位移表达式待定系数相应增加，n1时，最大误差减小为0.97%，最大挠度更加趋近于摄动解，且结果已能达到工程应用的精度要求。

4 结束语根据能量原理，本文结合里茨法和Gr6bner基法，应用符号计算软件Mathematica求解了等厚度圆卞的大挠度问题。该方法原理清晰，求解步骤易于编程实现，无需应用复杂的数学分析方法和技巧，就可得到结果的解析符号表达式，或者代入参数求出高精度的数值结果，便于工程应用。本文的方法和程序可推广到其它载荷与约束条件下圆卞的轴对称大挠度弯曲问题。