准确高效地处理几何非线性问题对于悬臂梁的大变形等实际工程至关重要．不同于小变形分析中应变位移满足线性关系的假设，大变形问题往往存在大位移或剧烈应变，此时应变位移不再是线性关系，无法满足小变形假设，须要采用几何非线性模型进行计算[1]．严格来讲，自然界所有问题都是非线性的，为了更准确地求解大变形结构，陈嵩涛等[2]研究认为当描述弹性体的非线性增量形式时，以现时构形作为参考构形(更新拉格朗日法)，不断更新构形，会明显降低计算效率，因此参考构形的选择直接影响到分析效率与结果的精确性．丁敏等[3]通过解析法，给出了圆拱结构几何非线性位移计算方法，但仅适用于圆拱结构，而数值解法借助于计算机强大的计算能力，能够轻松求解不同结构的大变形问题．有限元法(finite element method，FEM)作为数值解法，因其通用性与有效性，被广泛应用于固体力学问题求解中．目前在几何非线性分析中，有限元法通常首先结合线性应力-应变关系与非线性格林应变-位移关系得到非线性刚度矩阵，然后将刚度矩阵分为线性和非线性部分，通过网格单元离散求解域，用迭代方法求解待解方程组，最终得到较为精确的解．陈政清[4]基于有限元法建立了梁杆结构几何非线性计算方法，基于虚功方程实现了数值求解．但传统有限元法采用均质单元，计算得到的应力场存在人工不连续现象[5]，同时在处理几何非线性问题中，往往会遇到单元畸变等问题，难以进行力学计算[6]．与有限元法相比，有限体积法(finite volume method，FVM)实施方法简单且在求解域具有强守恒性，并具备明确的物理意义，这对于许多工程应用是十分重要的[7]．根据控制体的构造及变量的存储方式将有限体积法分为格点型有限体积法(cell-vertex finite volume method，CV-FVM)和格心型有限体积法(cell-centred finite volume method，CC-FVM)．有限体积法多用于求解流体问题[8]，1990年代开始应用于求解几何非线性问题．Moosavi等[9]发展了一种等几何无网格有限体积法，该方法的优点是较高的计算效率和精度．许秋怡等[10]结合纽马克数值算法和牛顿-拉夫逊迭代法求解几何非线性梁的动力学平衡方程．Fallah等[11]将格点型有限体积法和有限元法在二维几何非线性上的应用进行了比较，结果表明基于四边形网格，格点型有限体积法能以更高的网格密度达到有限元法的结果．Cardiff等[12]针对正交各向异性材料的大应变和大转动问题，提出并验证了一种格心型有限体积法．刘琦等[13]应用格心型有限体积法研究非线性大旋转问题，提出了一种改进的距离加权插值算法的大应变有限体积方法(LS-FVM)，结果表明LS-FVM能明显提高边界网格移动的数值精度．本研究将适用于线弹性问题的格点型有限体积法进一步推广应用于三维几何非线性问题的求解．结合线性应力-应变关系和非线性应变-位移理论，导出经格点型有限体积法离散的基尔霍夫应力描述的控制方程，以线弹性计算结果为初始近似值，迭代求解非线性方程组，直至达到给定的收敛标准．在整个分析过程中所有变量的表达格式都采用全拉格朗日型，并通过相关算例证明了本文方法的有效性和可靠性．1 控制方程选定图1中初始构形V0为参考构形(全拉格朗日法)，给出用基尔霍夫应力描述的控制方程∂(Fαγσβγ)/∂Xβ+fα=0, (1)式中：F为变形梯度；σ为基尔霍夫应力；X为初始构形坐标；f为单位初始体积的体积力．10.13245/j.hust.250613.F001图1大变形示意图线性基尔霍夫应力-格林应变关系和非线性格林应变-位移关系为σβγ=2Gεβγ+λεααδβγ;εβγ=12∂uβ∂xγ+∂uγ∂xβ+∂uα∂xβ∂uα∂xγ,式中：G和λ为拉梅常数，可通过弹性模量E和泊松比μ计算得到；ε为格林应变；δβγ为克罗尼克尔符号，当β=γ时δβγ=1，当β≠γ时δβγ=0；u为位移；x为现时构形坐标．又有：G=E/[2(1+μ)];λ=μE/[(1+μ)(1-2μ)].边界条件假设如下tα0=tbα0=Fαγσβγnβ    (边界S0t);uα=ubα    (边界S0u),式中：S0t和S0u分别为施加外力tbα0和规定位移ubα的边界；nβ为初始构形中面单位外法线矢量．2 数值离散2.1　控制体的建立如图2所示，采用四节点四面体、六节点棱柱、五节点金字塔或八节点六面体网格划分三维计算域．网格单元由粗实线及虚线围成，其节点用空心圆表示．围绕网格节点(节点1)用虚线依次连接相邻边长中点、面中心及体中心建立控制体，其边长中点、面中心及体中心用实心圆表示．10.13245/j.hust.250613.F002图2CV-FVM三维控制体将待解变量(位移)定义在网格节点上，并假设变量在控制体内均匀分布．已知物理量(如材料参数、体积力等)均定义在网格体中心，并假设在网格单元内均匀分布．2.2　控制方程的离散为便于程序开发，将上述方程的张量形式用矩阵形式重新表示．基于控制体V对式(1)计算体积分，并通过高斯定理将体积分转化为初始构型下面积分，得∫S0FTDεdS-R=0, (2)式中：F为变形梯度矩阵，F=I+H，I为单位矩阵；T为面单位外法线矢量矩阵；D为弹性矩阵；S0为初始构形面积．又有：H=∂(u)/∂(X)∂(u)/∂(Y)∂(u)/∂(Z)∂(v)/∂(X)∂(v)/∂(Y)∂(v)/∂(Z)∂(w)/∂(X)∂(w)/∂(Y)∂(w)/∂(Z);T=nX000nZnY0nY0nZ0nX00nZnYnX0;D=Cλλ000λCλ000λλC0000002G0000002G0000002G;R=-∫V0fdV-∫S0ttb0dS,式中C=2G+λ．结合单元内节点位移的型函数，第一项离散为∫S0FTDεdS=∑i=1nc∑j=1nn∫S0iFiTiDi(Bijuij)dS,式中：nc为当前节点周围分控制体(Vi)总数；nn为Vi所在单元节点总数；Bij为应变矩阵；u=[u，v，w]T为各节点的位移列向量；下标i为Vi所在单元中心的变量值，下标ij为Vi所在单元内第j个节点上的变量值．应变矩阵Bij分为线性和非线性两部分，线性部分应变矩阵为Bij0=∂Nij/∂X000∂Nij/∂Y000∂Nij/∂Z0∂Nij/∂Z∂Nij/∂Y∂Nij/∂Z0∂Nij/∂X∂Nij/∂Y∂Nij/∂X0.非线性部分应变矩阵为BijN=AiGij/2.矩阵A可用应变位移梯度函数θ表示为AT=θXT000θZTθYT0θYT0θZT0θXT00θZTθYTθXT0.以θX为例，有：θX=∂u/∂X,∂v/∂X,∂w/∂XT;Gij=∂Nij∂XI3×3,∂Nij∂YI3×3,∂Nij∂ZI3×3T,式中N为型函数[14]．任意单元i的应变、应力可以表示为：εi=∑j=0niBij0+12AiGijuij=∑j=0niBijz-12AiGijuij;σi=Diεi.第二项基于初始构形离散为R=-∑i=1ncfαiV0i/(ncni)-∑i=1nttpα0S0it,式中：fαi为当前节点所在单元α方向体积力；nt为当前节点周围外力边界网格面的个数；S0it为各外力边界网格面在初始构型中的面积．当给定位移边界ubα时，采用置大数法[15]处理约束位移边界．3 迭代格式及切线刚度矩阵F矩阵包含位移项，导致F矩阵是非线性的，采用牛顿-拉夫逊法求解，将F矩阵表示为待解变量u的泛函f(u)=∑i=1nc∑j=1nn∫S0iFiTiDiεidS-R.设f (u)为具有一阶导数的连续函数，把f (u)在Δu处泰勒展开，保留前两项，忽略第三项及后续高阶项，得f(u+Δu)=f(u)+f(u)'Δu.令f (u)ʹ为切线刚度矩阵Kt，当Δu足够小时，KtΔu=-f(u).当已知一个近似解un时，可以求得Kt和f(u)，进而算出Δu，un+1=un+Δu.基于线弹性结构问题格点型有限体积法求解程序[16]，得到线弹性计算结果uL，将其作为迭代初始值u0，迭代至达到给定收敛标准TO=1×10-6，即f(u)/RTo. (3)具体迭代求解过程如下：步骤1 读取计算域网格及控制信息；步骤2 小变形求解器求解得uL并作为大变形求解器迭代初始值u0；步骤3 根据un建立切线刚度矩阵Ktn并计算不平衡量f (un)；步骤4 计算节点位移增量Δun =f (un)/Ktn并叠加节点位移un+1=un+Δun；步骤5 判断是否达到收敛标准，若达到收敛标准则输出计算结果，否则更新位移变量，即un=un+1，重复步骤3～5．步骤3中Kt可表示为[11]Kt=∑i=1nc∑j=1ni∫s0FiTiDiBijz+TiePiGijdS,式中：Tie=nXiI3×3,nYiI3×3,nZiI3×3;Pi=σXXiI3×3σYXiI3×3σZXiI3×3σXYiI3×3σYYiI3×3σZYiI3×3σXZiI3×3σYZiI3×3σZZiI3×3.F，B，G，P，D均可由初始位移u0或上一迭代步位移un和定义在网格单元中心的已知物理量计算得到，因此Kt=∑i=1nc∑j=1niFiSiDiBijz+SiePiGij,式中S和Se均为控制体边界面的面积矩阵[14]．4 数值方法验证将本研究提出的几何非线性格点型有限体积公式应用于求解三维几何非线性问题，所有计算均采用C++编写的求解程序．为了提高非线性问题求解的收敛性，在进行大量迭代求解之前，须要对初始近似进行一定的评估，在良好的初始近似下，可以实现非常快速的收敛．因此，本研究用线弹性有限体积法的解作为非线性问题的初始近似，并通过牛顿-拉夫逊迭代求解非线性方程组．下面通过3个算例来说明本文方法的性能特点．4.1　点力作用下两端固定梁的大变形分析考虑两端固定梁的大变形问题，梁的长度L=200 mm，高度h=10 mm，厚度b=10 mm，如图3所示．材料参数分别为：弹性模量E=2.1×1011 N/m2，泊松比μ=0.3．采用如图4所示的四面体网格模型，具体网格信息见表1．10.13245/j.hust.250613.F003图3中心点力作用下的两端固定梁10.13245/j.hust.250613.F004图4两端固定梁四面体网格模型10.13245/j.hust.250613.T001表1两端固定梁四面体网格信息网格模型节点数/103单元数/103a62.954335.584b15.48967.516c4.85422.361d0.9723.670根据变形情况选取合适的点力值，最终梁的变形形状如图5所示．本研究格点型有限体积法解和有限元法解获得的中心点C处的位移计算结果如图6所示，从图中可以清楚地观察到，本算法的结果与有限元法的结果一致．对比基于不同网格得到的计算结果，格点型有限体积法具有与有限元法相当的网格收敛性．同时可以看到，在点力值超过一定限度后，线性位移解不再符合实际，表明非线性效应对解的影响较大，必须采用非线性求解理论，本研究格点型有限体积法和有限元法都能得到正确的非线性解．线性解计算公式为w=[F/(48EIe)](4y3-3Ly2),式中：w为z方向的位移；F为点力大小；Ie为惯性矩，Ie=bh3/12．10.13245/j.hust.250613.F005图5中心点力作用下两端固定梁的位移云图10.13245/j.hust.250613.F006图6中心点力作用下两端固定梁的位移曲线基于不同网格模型分析初始位移采用固定值0(方法1)和线弹性结果(方法2)的计算耗时．比较两种方法得到收敛解需要的迭代步数及计算时间，见表2和表3，表中：I，t和Δt分别为收敛步数、计算耗时和两种方法计算耗时之差，即Δt=|t1-t2|．基于相同的网格，采用方法1需要的迭代步数和收敛时间比采用方法2多，且网格单元数越大，这一现象越明显；因此，用线弹性计算结果作为非线性问题的初始近似可以实现较快速的收敛，节省计算时间．与有限元法相比，本文方法整体耗时更多，后续可通过进一步优化计算方法和计算程序缩短计算耗时．10.13245/j.hust.250613.T002表2基于网格模型a的计算耗时F/kN有限元法/s方法1方法2Δt/sI1t1/sI2t2/s10.5002144307.823241.7266.1013.1252164308.543243.4865.0615.7502124308.223246.6561.5621.0002785379.264309.7869.4842.0003266495.135426.5068.6310.13245/j.hust.250613.T003表3基于网格模型b的计算耗时F/kN有限元法/s方法1方法2Δt/sI1t1/sI2t2/s10.500594138.643111.4227.2213.125664139.543109.1430.4015.750644141.653109.1432.5021.000655173.244139.7733.4642.0001276195.515168.5926.924.2　集中载荷作用下悬臂梁的大变形分析如图7所示，对集中载荷作用下的悬臂梁进行大变形分析，计算域被划分为2.299×104个六面体单元．梁的长度L=5 m，横截面b=1 m，h=0.1 m，弹性模量E=3 MN/m2，泊松比取μ=0．施加的集中载荷为F(0 N≤F≤50 N)．10.13245/j.hust.250613.F007图7受集中载荷悬臂梁表4和表5为各载荷水平下自由端中心C点的y向和z向位移结果对比(Qr为误差值)，梁的变形如图8所示．根据误差计算公式，位移最大误差分别为2.79%和1.33%．图8表明：随着集中力的增加，悬臂梁弯曲的程度越来越大，在这个过程中格点型有限体积法计算结果与解析解符合良好，但较低的高宽比(h/b)将导致网格划分量的增加．10.13245/j.hust.250613.T004表4受集中载荷悬臂梁的y向位移vF/Nv/mmQr/%CV-FVMFEM解析解5-78-79-802.7910-276-275-2801.2920-791-809-8001.2930-1 257-1 259-1 2751.3840-1 630-1 649-1 6450.9450-1 923-1 849-1 9400.8510.13245/j.hust.250613.T005表5受集中载荷悬臂梁的z向位移wF/Nw/mmQr/%CV-FVMFEM解析解5-801-810-8101.0610-1 494-1 494-1 5101.0820-2 450-2 477-2 4700.7930-3 002-3 007-3 0150.4440-3 338-3 355-3 3500.3550-3 523-3 523-3 5701.33Qr=Q-QA / QA×100%, (4)式中：Q为计算得到的变量值；QA为解析解值．10.13245/j.hust.250613.F008图8受集中荷载悬臂梁的位移云图4.3　均布载荷作用下壳体的大变形分析考虑图9所示的两端固支θ=7.1°的壳体，其半径R=100 m，宽度a=12.422 m，厚度h=0.5 m，弹性模量E=3×107 N/m2，泊松比取μ=0.3．计算中采用4.882 2×104个四面体网格进行离散．为了验证本研究格点型有限体积法在计算壳体大变形问题中的准确性，同时给出了C点z向位移w解析解和有限元法解作为参照．图10位移曲线表明，计算结果与参照解符合良好．图11给出了该壳体在p=200 Pa时AB两点连线间的x向应力σx，y向应力σy，z向应力σz曲线，进一步证明了本文方法在计算该类型问题中的准确性．10.13245/j.hust.250613.F009图9均布荷载作用下的壳体10.13245/j.hust.250613.F010图10均布荷载作用下的壳体位移曲线10.13245/j.hust.250613.F011图11200 Pa均布荷载作用下的壳体应力曲线5 结论本研究拓展了用于求解弹性体问题的格点型有限体积法，将其拓展应用于求解三维几何非线性问题．基于格点型有限体积法离散三维几何非线性控制方程，以线弹性结果为初始值，迭代求解非线性方程组．通过数值模拟结果与解析解及有限元法结果对比，得到以下结论．a．数值模拟不同点力作用下的悬臂梁，点力超过一定限度后，线性解不再符合实际，因此采用非线性计算理论．另一方面，基于相同网格，采用线弹性计算结果为迭代初始场比0初始位移场需要的迭代步数和收敛时间少，且网格数量越大，节省计算时间越明显，计算效率越高．b．数值模拟不同集中力作用下的两端固支梁，本研究计算的y向和z向位移结果与解析解误差小于2.79%和1.33%，数值结果表明本文方法能够有效地模拟三维几何非线性问题的位移情况．c．数值模拟不同均布载荷作用下的壳体，位移及应力结果均与解析解及有限元法结果符合良好，表明该方法可以用于几何非线性问题的数值预测，并作为有限元法的有效补充方法．