圆轴是极为常见的零件，一般情况下，它的受力分析是三维问题，但对于拉伸、扭转、纯弯三种常见的受力情形，问题可以简化为二维。虽然以今日计算机的能力，直接处理三维问题也不是太难，但化简仍很有意义。关于弹性力学和有限元的资料已有很多，本文中先简要介绍，再将其应用到圆轴的分析中。
圆柱坐标系下的弹性力学方程
圆轴几何上为轴对称，使用圆柱坐标系较为方便。圆柱坐标下，应变与位移的表达式如下 \[ \begin{pmatrix} \varepsilon_r \\ \varepsilon_\theta \\ \varepsilon_z \\ \gamma_{r\theta} \\ \gamma_{\theta z} \\ \gamma_{zr} \end{pmatrix} = \begin{pmatrix} \frac{\partial}{\partial r} & 0 & 0 \\ \frac{1}{r} & \frac{1}{r}\frac{\partial}{\partial \theta} & 0 \\ 0 & 0 & \frac{\partial}{\partial z} \\ \frac{1}{r}\frac{\partial}{\partial \theta} & \frac{\partial}{\partial r}-\frac{1}{r} & 0 \\ 0 & \frac{\partial}{\partial z} & \frac{1}{r}\frac{\partial}{\partial \theta} \\ \frac{\partial}{\partial z} & 0 & \frac{\partial}{\partial r} \end{pmatrix} \begin{pmatrix} u_r\\u_\theta\\u_z \end{pmatrix} \] 记为\(\boldsymbol{\varepsilon}=\mathbf{Su}\)。应力与应变的关系和直角坐标类似 \[ \begin{pmatrix} \sigma_r \\ \sigma_\theta \\ \sigma_z \\ \tau_{r\theta} \\ \tau_{\theta z} \\ \tau_{zr} \end{pmatrix} = \frac{E}{(1+\nu)(1-2\nu)} \begin{pmatrix} 1-\nu & \nu & \nu & 0 & 0 & 0 \\ \nu & 1-\nu & \nu & 0 & 0 & 0 \\ \nu & \nu & 1-\nu & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} \end{pmatrix} \begin{pmatrix} \varepsilon_r \\ \varepsilon_\theta \\ \varepsilon_z \\ \gamma_{r\theta} \\ \gamma_{\theta z} \\ \gamma_{zr} \end{pmatrix} \] 记为\(\boldsymbol{\sigma}=\mathbf{D}\boldsymbol{\varepsilon}\)，其中\(E\)为弹性模量，\(\nu\)为泊松系数。将体积力的分量记为\(F_r\)、\(F_\theta\)、\(F_z\)，则圆柱坐标下的平衡方程为 \[ \begin{aligned} & \frac{\partial \sigma_{r}}{\partial r} + \frac{1}{r}\frac{\partial \tau_{r\theta}}{\partial \theta} + \frac{\partial \tau_{zr}}{\partial z} + \frac{1}{r}(\sigma_{r}-\sigma_{\theta}) + F_r= 0 \\ & \frac{\partial \tau_{r\theta}}{\partial r} + \frac{1}{r}\frac{\partial \sigma_{\theta}}{\partial \theta} + \frac{\partial \tau_{\theta z}}{\partial z} + \frac{2}{r}\tau_{r\theta} + F_\theta = 0 \\ & \frac{\partial \tau_{zr}}{\partial r} + \frac{1}{r}\frac{\partial \tau_{\theta z}}{\partial \theta} + \frac{\partial \sigma_{z}}{\partial z} + \frac{1}{r}\tau_{zr} + F_z = 0 \end{aligned} \]
有限元
由以上关系，应力的各个分量用位移表达，代入平衡方程中，得到三个关于位移的方程，未知的位移分量也是三个，便可求得相应的解。但是，前述方程的表达式非常复杂，除极简单的情形外，无法得到解析解，所以应用中通常使用数值方法。在数值解法中，有限元方法是一种主要方法，在弹性力学领域更是占有绝对地位。实际上，有限元方法求解的并不是偏微分方程，而是使用加权残差法将偏微分方程转化为另一种形式的方程。这种转化可对应于力学中的虚功原理：设有虚位移\(\delta\mathbf{u}\)，相应的虚应变\(\delta\boldsymbol{\varepsilon}=\mathbf{S \delta u}\)，产生的虚应变能为\(\int_\Omega \delta\boldsymbol{\varepsilon}^T\boldsymbol{\sigma}dV=\int_\Omega \delta\boldsymbol{\varepsilon}^T\mathbf{D}\boldsymbol{\varepsilon}dV\)，应等于外力在虚位移上作的功，即体积力的功\(\int_\Omega \delta\mathbf{u^T F} dV\)加上边界上力的功\(\int_{\partial\Omega} \delta\mathbf{u^T f} dS\)。
有限元法中，需将求解域划分为若干单元，再在单元的边界上取一些结点，将域上的位移近似为结点位移\(\mathbf{u}_i\)的插值 \[ \mathbf{u} = \sum_i N_i\mathbf{u}_i \] \(N_i\)是插值函数，只在结点相邻的单元内不为零。将此近似表示代入到虚功原理的公式中，便可得到关于\(\mathbf{u}_i\)的线性方程组。在建立方程的过程中，需要处理形如\(N_iN_j\)，或是它们微分乘积的积分，而这只有当结点\(i\)、\(j\)都属于同一单元时才不为零，所以这些积分可以一个单元一个单元地计算，而得到的方程中，大多数系数都是零，所以，可以高效地构建和求解方程。这一优势是有限元方法得到广泛应用的重要原因。从方程解得\(\mathbf{u}_i\)后，按插值公式又可进一步算出应变、应力等。
接下来转到本文的主题，讲述三种不同受力情形的求解。其中的关键是确定相应的方程，也就是找到相应的矩阵\(\mathbf{S}\)和\(\mathbf{D}\)。
扭转
在此情形中，位移与\(\theta\)无关，且只有\(u_\theta\)不为零，所以应变分量中只有\(\gamma_{r\theta}\)、\(\gamma_{\theta z}\)不为零 \[ \begin{pmatrix} \gamma_{r\theta} \\ \gamma_{\theta z} \end{pmatrix} = \begin{pmatrix} \frac{\partial}{\partial r}-\frac{1}{r}\\ \frac{\partial}{\partial z} \end{pmatrix} u_\theta \] 应力分量相应地只有两项 \[ \begin{pmatrix} \tau_{r\theta} \\ \tau_{\theta z} \end{pmatrix} = \frac{E}{2(1+\nu)} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} \gamma_{r\theta} \\ \gamma_{\theta z} \end{pmatrix} \] 这里的系数即剪切模量\(G=\frac{E}{2(1+\nu)}\)。原问题转为只有一个未知量\(u_\theta\)的二维问题。端部受力的边界条件可按照材料力学中的设定，正比于\(r\)，比值为扭矩除以截面极惯性矩。
对这一问题还可作一补充。注意到 \[ \gamma_{r\theta} = r\frac{\partial (u_\theta/r)}{\partial r} ,\, \gamma_{\theta z} = r\frac{\partial (u_\theta/r)}{\partial z} \] 如果以\(\psi=u_\theta/r\)为未知量，就可以将方程写得更为对称。虚应变能可表示为 \[ 2\pi G\iint_\Omega \left(\frac{\partial \delta\psi}{\partial r}\frac{\partial \psi}{\partial r} + \frac{\partial \delta\psi}{\partial z}\frac{\partial \psi}{\partial z}\right)r^3drdz \] 这一形式实则对应于传导方程。如果对于这种对应不熟悉，也可以从平衡方程推导。由于只有\(\tau_{r\theta}\)、\(\tau_{\theta z}\)不为零，三个平衡方程只需考虑第二个，将应力用\(\psi\)表达，代入方程 \[ \frac{\partial}{\partial r}\left(Gr\frac{\partial \psi}{\partial r}\right) + \frac{\partial}{\partial z}\left(Gr\frac{\partial \psi}{\partial z}\right) + \frac{2}{r}\left(Gr\frac{\partial \psi}{\partial r}\right) = 0 \] 两边除以\(G\)，再乘上\(r^2\) \[ r^3\frac{\partial^2 \psi}{\partial r^2} + 3r^2\frac{\partial \psi}{\partial r} + r^3\frac{\partial^2 \psi}{\partial z^2} = 0 \] 即 \[ \frac{\partial}{\partial r}\left(r^3\frac{\partial \psi}{\partial r}\right) + \frac{\partial}{\partial z}\left(r^3\frac{\partial \psi}{\partial z}\right) = 0 \] 便是二维传导方程，只是传导率不恒定，而等于\(r^3\)。由于这一关系，将板的厚度取为\(r^3\)，通过电传导模拟，在数值方法成熟以前，便可对轴的扭转作精确的分析。
拉伸
拉伸的情形正好与前面扭转的情形与互补。位移同样与\(\theta\)无关，不过\(u_\theta\)为零，而\(u_r\)、\(u_z\)不为零。此时，切应变中只有\(\gamma_{zr}\)不为零 \[ \begin{pmatrix} \varepsilon_r \\ \varepsilon_\theta \\ \varepsilon_z \\ \gamma_{zr} \end{pmatrix} = \begin{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} \begin{pmatrix} u_r \\ u_z \end{pmatrix} \] 应力亦相应简化 \[ \begin{pmatrix} \sigma_r \\ \sigma_\theta \\ \sigma_z \\ \tau_{zr} \end{pmatrix} = \frac{E}{(1+\nu)(1-2\nu)} \begin{pmatrix} 1-\nu & \nu & \nu & 0 \\ \nu & 1-\nu & \nu & 0 \\ \nu & \nu & 1-\nu & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2} \end{pmatrix} \begin{pmatrix} \varepsilon_r \\ \varepsilon_\theta \\ \varepsilon_z \\ \gamma_{zr} \end{pmatrix} \] 原问题转为关于\(u_r\)、\(u_z\)的二维问题。端部受力的边界条件为均匀力，等于总力除以截面积。商用有限元软件都提供专门处理这种几何和受力均为轴对称情形的方法。
弯曲
此情况比上面两种复杂不少，位移的分量都不为零，且随\(\theta\)改变。为寻找思路，先看一种已得到解析解的简单情形。圆柱杆的截面极惯性矩\(I_p=\frac{\pi d^4}{64}\)，在力矩\(M\)下的弯曲曲率\(\kappa=\frac{M}{EI_p}\)。设弯曲在\(xz\)平面，弯向\(x\)轴的负方向，则位移的解为\(u_x = -\frac12\kappa(z^2+\nu(x^2-y^2))\)，\(u_y = -\kappa\nu xy\)，\(u_z = \kappa xz\)。转为柱坐标则是\(u_r = -\frac12\kappa(z^2+\nu r^2)\cos\theta\)，\(u_\theta = \frac12\kappa(z^2-\nu r^2)\sin\theta\)，\(u_z=\kappa zr\cos\theta\)，而应力\(\sigma_z=\frac{M}{I_p}x=\frac{M}{I_p}r\cos\theta\)，其他的应力都为零。
对于一般的轴受弯问题，前面的解启发我们假定：\(u_r = \bar u_r\cos\theta\)，\(u_\theta=\bar u_\theta\sin\theta\)，\(u_z = \bar u_z\cos\theta\)，带上划线的项与\(\theta\)无关，可得应变 \[ \begin{aligned} \varepsilon_{r} & = \frac{\partial \bar u_r}{\partial r}\cos\theta \\ \varepsilon_{\theta} &= \frac{1}{r}\left(\bar u_\theta\cos\theta + \bar u_r\cos\theta\right) \\ \varepsilon_{z} &= \frac{\partial \bar u_z}{\partial z}\cos\theta \\ \gamma_{r\theta} & = -\frac{\bar u_r}{r}\sin\theta + \frac{\partial \bar u_\theta}{\partial r}\sin\theta - \frac{\bar u_\theta}{r}\sin\theta \\ \gamma_{\theta z} & = \frac{\partial \bar u_\theta}{\partial z}\sin\theta - \frac{1}{r}\bar u_z\sin\theta \\ \gamma_{zr} &= \frac{\partial \bar u_r}{\partial z}\cos\theta + \frac{\partial \bar u_z}{\partial r}\cos\theta \end{aligned} \] 可以看出，每一项只与\(\cos\theta\)或只与\(\sin\theta\)相关，而进一步得到的\(\sigma_z\)只与\(\cos\theta\)相关，也和边界条件一致，所以这样的假定是可行的。在计算虚功时，对\(\theta\)积分会出现\(\int^{2\pi}_{0}\cos^2\theta d\theta\)或\(\int^{2\pi}_{0}\sin^2\theta d\theta\)，它们的值都是\(\pi/2\)，从而可以从等式中消去，最终得到的方程与\(\theta\)无关。原问题转化为关于\(\bar u_r\)、\(\bar u_\theta\)、\(\bar 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\)对应于柱坐标下的微体积\(2\pi ydxdy\)，所以方程的积分中要乘上\(y\)，\(2\pi\)这个常数可以从方程两边消去。
另外，前述应变公式中，有的分量中出现了\(1/r\)，在\(r=0\)时会出现奇异。所以，在作几何定义时，最好避开\(r=0\)处，就是说，即使是实心轴，也可以在中心留一个微小的孔。
