圆轴的有限元分析

Posted 2020-11-24 wildabc

圆轴是极为常见的零件，一般情况下，它的受力分析是三维问题，但对于拉伸、扭转、纯弯三种常见的受力情形，问题可以简化为二维。虽然以今日计算机的能力，直接处理三维问题也不是太难，但化简仍很有意义。关于弹性力学和有限元的资料已有很多，本文中先简要介绍，再将其应用到圆轴的分析中。

圆柱坐标系下的弹性力学方程

圆轴几何上为轴对称，使用圆柱坐标系较为方便。圆柱坐标下，应变与位移的表达式如下
[ egin{pmatrix} ? ? varepsilon_r \ varepsilon_ heta \ varepsilon_z ? ? gamma_{r heta} \ gamma_{ heta z} \ gamma_{zr} end{pmatrix} = egin{pmatrix} ? ? frac{partial}{partial r} & 0 & 0 ? ? frac{1}{r} & frac{1}{r}frac{partial}{partial heta} & 0 ? ? 0 & 0 & frac{partial}{partial z} ? ? frac{1}{r}frac{partial}{partial heta} & frac{partial}{partial r}-frac{1}{r} & 0 ? ? 0 & frac{partial}{partial z} & frac{1}{r}frac{partial}{partial heta} ? ? frac{partial}{partial z} & 0 & frac{partial}{partial r} end{pmatrix} egin{pmatrix} ? ? u_r\u_ heta\u_z end{pmatrix} ]
记为(oldsymbol{varepsilon}=mathbf{Su})。应力与应变的关系和直角坐标类似
[ egin{pmatrix} ? ? sigma_r \ sigma_ heta \ sigma_z ? ? au_{r heta} \ au_{ heta z} \ au_{zr} end{pmatrix} = frac{E}{(1+ u)(1-2 u)} egin{pmatrix} ? ? 1- u & u & u & 0 & 0 & 0 ? ? u & 1- u & u & 0 & 0 & 0 ? ? u & u & 1- u & 0 & 0 & 0 ? ? 0 & 0 & 0 & frac{1-2 u}{2} & 0 & 0 ? ? 0 & 0 & 0 & 0 & frac{1-2 u}{2} & 0 ? ? 0 & 0 & 0 & 0 & 0 & frac{1-2 u}{2} end{pmatrix} egin{pmatrix} ? ? varepsilon_r \ varepsilon_ heta \ varepsilon_z ? ? gamma_{r heta} \ gamma_{ heta z} \ gamma_{zr} end{pmatrix} ]
记为(oldsymbol{sigma}=mathbf{D}oldsymbol{varepsilon})，其中(E)为弹性模量，( u)为泊松系数。将体积力的分量记为(F_r)、(F_ heta)、(F_z)，则圆柱坐标下的平衡方程为
[ egin{aligned} ? ? & frac{partial sigma_{r}}{partial r} + frac{1}{r}frac{partial au_{r heta}}{partial heta} + frac{partial au_{zr}}{partial z} + frac{1}{r}(sigma_{r}-sigma_{ heta}) + F_r= 0 ? ? & frac{partial au_{r heta}}{partial r} + frac{1}{r}frac{partial sigma_{ heta}}{partial heta} + frac{partial au_{ heta z}}{partial z} + frac{2}{r} au_{r heta} + F_ heta = 0 ? ? & frac{partial au_{zr}}{partial r} + frac{1}{r}frac{partial au_{ heta z}}{partial heta} + frac{partial sigma_{z}}{partial z} + frac{1}{r} au_{zr} + F_z = 0 ? end{aligned} ]

有限元

由以上关系，应力的各个分量用位移表达，代入平衡方程中，得到三个关于位移的方程，未知的位移分量也是三个，便可求得相应的解。但是，前述方程的表达式非常复杂，除极简单的情形外，无法得到解析解，所以应用中通常使用数值方法。在数值解法中，有限元方法是一种主要方法，在弹性力学领域更是占有绝对地位。实际上，有限元方法求解的并不是偏微分方程，而是使用加权残差法将偏微分方程转化为另一种形式的方程。这种转化可对应于力学中的虚功原理：设有虚位移(deltamathbf{u})，相应的虚应变(deltaoldsymbol{varepsilon}=mathbf{S delta u})，产生的虚应变能为(int_Omega deltaoldsymbol{varepsilon}^Toldsymbol{sigma}dV=int_Omega deltaoldsymbol{varepsilon}^Tmathbf{D}oldsymbol{varepsilon}dV)，应等于外力在虚位移上作的功，即体积力的功(int_Omega deltamathbf{u^T F} dV)加上边界上力的功(int_{partialOmega} deltamathbf{u^T f} dS)。

有限元法中，需将求解域划分为若干单元，再在单元的边界上取一些结点，将域上的位移近似为结点位移(mathbf{u}_i)的插值
[ mathbf{u} = sum_i N_imathbf{u}_i ]
(N_i)是插值函数，只在结点相邻的单元内不为零。将此近似表示代入到虚功原理的公式中，便可得到关于(mathbf{u}_i)的线性方程组。在建立方程的过程中，需要处理形如(N_iN_j)，或是它们微分乘积的积分，而这只有当(i)、(j)都属于同一单元时才不为零，所以这些积分可以一个单元一个单元地计算，而得到的方程中，大多数系数都是零，所以，可以高效地构建和求解方程。这一优势是有限元方法得到广泛应用的重要原因。从方程解得(mathbf{u}_i)后，按插值公式又可进一步算出应变、应力等。

接下来转到本文的主题，讲述三种不同受力情形的求解。其中的关键是确定相应的方程，也就是找到相应的矩阵(mathbf{S})和(mathbf{D})。

扭转

在此情形中，位移与( heta)无关，且只有(u_ heta)不为零，所以应变分量中只有(gamma_{r heta})、(gamma_{ heta z})不为零
[ egin{pmatrix} ? ? gamma_{r heta} \ gamma_{ heta z} end{pmatrix} = egin{pmatrix} ? ? frac{partial}{partial r}-frac{1}{r}?? frac{partial}{partial z} end{pmatrix} u_ heta ]
应力分量相应地只有两项
[ egin{pmatrix} ? ? au_{r heta} \ au_{ heta z} end{pmatrix} = frac{E}{2(1+ u)} egin{pmatrix} ? ? 1 & 0 \ 0 & 1 end{pmatrix} egin{pmatrix} ? ? gamma_{r heta} \ gamma_{ heta z} end{pmatrix} ]
这里的系数即剪切模量(G=frac{E}{2(1+ u)})。原问题转为只有一个未知量(u_ heta)的二维问题。端部受力的边界条件可按照材料力学中的设定，正比于(r)，比值为扭矩除以截面极惯性矩。

对这一问题还可作一补充。注意到
[ gamma_{r heta} = rfrac{partial (u_ heta/r)}{partial r} ,, gamma_{ heta z} = rfrac{partial (u_ heta/r)}{partial z} ]
如果以(psi=u_ heta/r)为未知量，就可以将方程写得更为对称。虚应变能可表示为
[ 2pi Giint_Omega left(frac{partial deltapsi}{partial r}frac{partial psi}{partial r} +? frac{partial deltapsi}{partial z}frac{partial psi}{partial z} ight)r^3drdz ]
这一形式实则对应于传导方程。如果对于这种对应不熟悉，也可以从平衡方程推导。由于只有( au_{r heta})、( au_{ heta z})不为零，三个平衡方程只需考虑第二个，将应力用(psi)表达，代入方程
[ frac{partial}{partial r}left(Grfrac{partial psi}{partial r} ight) + frac{partial}{partial z}left(Grfrac{partial psi}{partial z} ight) + frac{2}{r}left(Grfrac{partial psi}{partial r} ight) = 0 ]
两边除以(G)，再乘上(r^2)
[ r^3frac{partial^2 psi}{partial r^2} + 3r^2frac{partial psi}{partial r} + r^3frac{partial^2 psi}{partial z^2} = 0 ]
即
[ frac{partial}{partial r}left(r^3frac{partial psi}{partial r} ight) + frac{partial}{partial z}left(r^3frac{partial psi}{partial z} ight) = 0 ]
便是二维传导方程，只是传导率不恒定，而等于(r^3)。由于这一关系，将板的厚度取为(r^3)，通过电传导模拟，在数值方法成熟以前，便可对轴的扭转作精确的分析。

拉伸

拉伸的情形正好与前面扭转的情形与互补。位移同样与( heta)无关，不过(u_ heta)为零，而(u_r)、(u_z)不为零。此时，切应变中只有(gamma_{zr})不为零
[ egin{pmatrix} ? ? varepsilon_r \ varepsilon_ heta \ varepsilon_z ? ? gamma_{zr} end{pmatrix} = egin{pmatrix} ? ? frac{partial}{partial r} & 0 ? ? frac{1}{r} & 0 ? ? 0 & frac{partial}{partial z} ? ? frac{partial}{partial z} & frac{partial}{partial r} end{pmatrix} egin{pmatrix} ? ? u_r \ u_z end{pmatrix} ]
应力亦相应简化
[ egin{pmatrix} ? ? sigma_r \ sigma_ heta \ sigma_z ? ? au_{zr} end{pmatrix} = frac{E}{(1+ u)(1-2 u)} egin{pmatrix} ? ? 1- u & u & u & 0 ? ? u & 1- u & u & 0 ? ? u & u & 1- u & 0 ? ? 0 & 0 & 0 & frac{1-2 u}{2} end{pmatrix} egin{pmatrix} ? ? varepsilon_r \ varepsilon_ heta \ varepsilon_z ? ? gamma_{zr} end{pmatrix} ]
原问题转为关于(u_r)、(u_z)的二维问题。端部受力的边界条件为均匀力，等于总力除以截面积。商用有限元软件都提供专门处理这种几何和受力均为轴对称情形的方法。

弯曲

此情况比上面两种复杂不少，位移的分量都不为零，且随( heta)改变。为寻找思路，先看一种已得到解析解的简单情形。圆柱杆的截面极惯性矩(I_p=frac{pi d^4}{64})，在力矩(M)下的弯曲曲率(kappa=frac{M}{EI_p})。设弯曲在(xz)平面，弯向(x)轴的负方向，则位移的解为(u_x = -frac12kappa(z^2+ u(x^2-y^2)))，(u_y = -kappa u xy)，(u_z = kappa xz)。转为柱坐标则是(u_r = -frac12kappa(z^2+ u r^2)cos heta)，(u_ heta = frac12kappa(z^2- u r^2)sin heta)，(u_z=kappa zrcos heta)，而应力(sigma_z=frac{M}{I_p}x=frac{M}{I_p}rcos heta)，其他的应力都为零。

对于一般的轴受弯问题，前面的解启发我们假定：(u_r = ar u_rcos heta)，(u_ heta=ar u_ hetasin heta)，(u_z = ar u_zcos heta)，带上划线的项与( heta)无关，可得应变
[ egin{aligned} ? ? varepsilon_{r} & = frac{partial ar u_r}{partial r}cos heta ? ? varepsilon_{ heta}? &= frac{1}{r}left(ar u_ hetacos heta + ar u_rcos heta ight) ? ? varepsilon_{z}? &= frac{partial ar u_z}{partial z}cos heta ? ? gamma_{r heta} & = -frac{ar u_r}{r}sin heta + frac{partial ar u_ heta}{partial r}sin heta - frac{ar u_ heta}{r}sin heta ? ? gamma_{ heta z}? & = frac{partial ar u_ heta}{partial z}sin heta - frac{1}{r}ar u_zsin heta ? ? gamma_{zr} &= frac{partial ar u_r}{partial z}cos heta + frac{partial ar u_z}{partial r}cos heta end{aligned} ]
可以看出，每一项只与(cos heta)或只与(sin heta)相关，而进一步得到的(sigma_z)只与(cos heta)相关，也和边界条件一致，所以这样的假定是可行的。在计算虚功时，对( heta)积分会出现(int^{2pi}_{0}cos^2 heta d heta)或(int^{2pi}_{0}sin^2 heta d heta)，它们的值都是(pi/2)，从而可以从等式中消去，最终得到的方程与( heta)无关。原问题转化为关于(ar u_r)、(ar u_ heta)、(ar u_z)的二维问题。端部受力的边界条件同样可按照材料力学中的设定，正比于(r)，比值为弯矩除以截面惯性矩。

具体实现

能够自定义方程的有限元软件也有不少，本文中使用的是开源的FreeFem++。它的文档比较详细，老版本的文档有中文版，基本的用法与现版本是相同的。在手册中也给出了弹性力学的例子，不过没有使用矩阵来定义问题。个人觉得使用矩阵会更易懂，而且在FreeFem++中实现也很简单。关键点是，要用func而不是matrix来声明方程定义中需要的矩阵。以本文中最复杂的弯曲问题为例

func D = E/(1+nu)/(1-2*nu)*
            [[1-nu, nu, nu, 0, 0, 0],
            [nu, 1-nu, nu, 0, 0, 0],
            [nu, nu, 1-nu, 0, 0, 0],
            [0, 0, 0, (1-2*nu)/2, 0, 0],
            [0, 0, 0, 0, (1-2*nu)/2, 0],
            [0, 0, 0, 0, 0, (1-2*nu)/2]];

而后使用宏来定义应变

macro epsilon(uz, ur, ut) [dx(uz), dy(ur), (ut+ur)/y, dy(ut)-(ur+ut)/y, dx(ut)-uz/y, dx(ur)+dy(uz)] //

在定义方程时域上的虚功积分便可写为

int2d(Th) (epsilon(uz,ur,ut)'*D*epsilon(vz,vr,vt)*y)

这里是以(x)轴为柱坐标中的(z)轴，(y)轴为柱坐标中的(r)轴。平面上微元(dxdy)对应于柱坐标下的微体积(2pi ydxdy)，所以方程的积分中积分中要乘上(y)，(2pi)这个常数可以从方程两边消去。

另外，前述应变公式中，有的分量中出现了(1/r)，在(r=0)时会出现奇异。所以，在作几何定义时，最好避开(r=0)处，就是说，即使是实心轴，也可以在中心留一个微小的孔。