介绍
考虑下图的情形，中间有一圆孔的无限大平板在x方向受到均匀拉力。
这个问题可分解为以下两种情况的叠加：左是孔边受到板为完整时内部对它的力，此时的应力状态与板为完整时完全相同；右是孔边受到相反的力，而x方向的均匀拉力消失。可以想到，右图中的力会使孔周边的应力不再均匀，一些地方的应力会大于板为完整时，这种情况被称为应力集中。设均匀拉力为\(p\)，前人早已求得，孔边与x轴正方向的夹角为\(\theta\)处的切向应力为\(\sigma_t = p(1-2\cos 2\theta)\)，应力在\(\theta=\pm\frac{\pi}{2}\)处有最大值\(3p\)，在\(\theta=0\)或\(\pi\)处有最小值\(-p\)。远离孔处仍然是均匀应力\(\sigma_x=p\)。应力集中产生的最大应力与未受应力集中影响处的应力的比值称为应力集中系数，所以这里的应力集中系数为3。
在这个例子中，应力集中的系数与孔的大小无关。一个小孔就会使最大应力增至三倍，可能使平板变形甚至破坏。这就促使人们思考，能否改变孔的形状，减少应力集中。文献1在充分汇总前人方法的基础上，使用有限元方法对此问题以及相关的一些其他情况做了分析，得到了目前的最好结果。按此文章，在两个方向的最大尺寸相同的情况下，如果孔为鼓形，y方向的尺寸在中间最大，两端收窄约20%，应力集中系数可降至2.15。
前人求解此类问题，大多是像文献1一样从一初始形状开始，以有限元之类方法计算出应力分布，然后调整形状，以期改善应力分布，多次迭代后以一定条件判定达至最优。但也有文章从逆问题的角度去分析，由最优应力分布逆向构造孔的形状。文献2和文献3就对最优应力分布的两种情形做了分析，但方法的限制较大，很多情况下并不存在满足条件的解，就比如本文的这个看起来是最简单的情况。本文尝试给最优应力分布增加一个变量，使得这种逆向方法适用于更多问题，并将此方法的结果与文献1中的相比较，来理解“最优”究竟意味着什么。
要从应力分布逆向得到形状，用到的是平面弹性问题的复函数法。铁摩辛柯和古地尔的《弹性理论》在第六章介绍了穆斯赫利什维利（Muskhelishvili）提出的这一方法。通过将其中的复函数以有限项级数近似，原本的解析问题可以转而用数值方法求解。1951年的文献4是本人找到的最早使用这种级数方法求解平面弹性问题的文章。穆斯赫利什维利本人的《数学弹性理论的几个基本问题》早有中译，不过部头很大，记号与今日通行的不同，对于本文这种程度问题并无必要。当然，掌握此方法需要复变函数知识，这在理工科是一门必修课，但本人已忘掉大部分，只得找资料再学习一遍。
平面弹性复函数法基本公式
平面弹性复函数法的具体推演过程可见《弹性理论》，这里直接给出基本公式。平面应力可由两个复解析函数\(\phi(z)\)和\(\psi(z)\)表示
\[ \sigma_x+\sigma_y=4\mathrm{Re}\left(\phi'(z)\right)\\ \sigma_y-\sigma_x+2i\tau_{xy}=2\bar z\phi''(z)+2\psi'(z) \]
在边界\(L\)上受到的外力分量分别是\(\bar X\)、\(\bar Y\)，记\(f(z)=i\int_L \left(\bar X+i\bar Y\right)dz\)，则边界条件为
\[ \phi(z)+z\bar\phi'(z)+\bar\psi(z)=f(z) \]
孔的应力求解
如前面的图所示，本问题中的应力可表示为两种状态的叠加，可记为
\[ \phi^1(z) = \phi^0(z)+\phi(z)\\ \psi^1(z) = \psi^0(z)+\psi(z) \]
各应力分量中，本问题只需解得切向正应力\(\sigma^1_t\)。正应力之和在坐标变换下保持不变，且法向正应力\(\sigma^1_n=0\)，有\(\sigma^1_t=\sigma^1_x+\sigma^1_y=4\mathrm{Re}\left(\phi^1{'}(z)\right)\)，所以两个复函数中只需求解一个。
无穷大平板在x方向受到均匀拉力\(p\)，\(\sigma_x^0=p\)，其他应力为0，由前面公式有
\[ \phi^0(z) = \frac14pz\\ \psi^0(z) = -\frac12pz \]
孔边上无外力，\(f^1(z)=0\)，对于\(\phi(z)\)、\(\psi(z)\)，边界条件为
\[ \begin{aligned} \phi(z)+z\bar{\phi}'(z)+\bar\psi(z)&=-\left(\phi^0(z)+z\bar{\phi^0}'(z)+\bar\psi^0(z)\right)\\ &=-\frac{p}2(z-\bar z) \end{aligned} \]
为了考察一般孔形下的应力，设孔的曲线是由解析函数\(z=\omega(\zeta)\)从单位圆映射而来，对于\(\zeta\)前式变为
\[ \phi(\zeta)+\omega(\zeta)\frac{\bar{\phi}'(\zeta)}{\bar\omega'(\zeta)}+\bar\psi(\zeta)=-\frac{p}2(\omega(\zeta)-\bar \omega(\zeta)) \]
令\(\Phi(\zeta)={\phi}'(\zeta)/\omega'(\zeta)\)，有\(\sigma_x+\sigma_y=4\mathrm{Re}\left(\Phi(\zeta)\right)\)，而上式变为
\[ \int\omega'(\zeta)\Phi(\zeta)\mathrm{d}\zeta+\omega(\zeta)\bar\Phi(\zeta)+\bar\psi(\zeta)=-\frac{p}2(\omega(\zeta)-\bar \omega(\zeta)) \]
在本问题中，只需要映射改变孔的附近，远处仍可保持平直，即\(|\zeta|\to\infty\)时，\(\omega(\zeta)\to R\zeta\)，而应力应趋向于0。进一步地，本问题只需考虑孔的形状，与大小无关，\(R\)可取为1。于是可设
\[ \omega(\zeta)=\zeta+\sum_{k=1}^\infty c_k\zeta^{-k}\\ \Phi(\zeta)=\sum_{k=1}^\infty u_k\zeta^{-k}\\ \psi(\zeta)=\sum_{k=1}^\infty v_k\zeta^{-k} \]
再利用当\(|\zeta|=1\)时，\(\bar\zeta=1/\zeta\)，代入前式
\[ \int\left(1+\sum_{k=1}^\infty-kc_k\zeta^{-k-1}\right)\sum_{k=1}^\infty u_k\zeta^{-k}\mathrm{d}\zeta+\left(\zeta+\sum_{k=1}^\infty c_k\zeta^{-k}\right)\sum_{k=1}^\infty \bar u_k\zeta^{k}\\ +\sum_{k=1}^\infty \bar v_k\zeta^{k}=-\frac{p}2\left((\zeta+\sum_{k=1}^\infty c_k\zeta^{-k})-(\frac{1}{\zeta}+\sum_{k=1}^\infty \bar c_k\zeta^{k})\right) \]
只看\(\zeta^{-n}\)的系数，等式为
\[ -\frac1n (u_{n+1}-\sum_{k=1}^{n-1} kc_ku_{n-k})+\sum_{k=n+1}^\infty c_k\bar u_{k-n}=-\frac{p}2c_n \]
\(n=1\)时，右边再加上\(-p/2\)。通过只考虑负幂次项，得到了与\(\psi(\zeta)\)无关、只含有\(\omega(\zeta)\)和\(\Phi(\zeta)\)展开系数的等式。这些展开系数可通过傅立叶级数或施瓦茨–克里斯托费尔映射等方法计算。如果截取前若干项，已知\(c_k\)、\(u_k\)中的一组，便得到关于另一组的线性方程组。如果需要，还可继续利用\(\zeta^{n}\)的系数，得到含有\(v_k\)等式，在已知\(c_k\)、\(u_k\)后，求解出完整的应力函数。从数值求解的角度看，这样得到的矩阵是非稀疏的，不便于求解，更关键的是这方法无法对局部做细化，提高精度比较困难，所以用来处理一般问题当然是不如通行的有限元方法。但是这个方法不仅可以从形状求应力，还可以逆向从应力求形状，这个优势是其他方法没有的。
如果孔有对称性，等式还可进一步简化。如果孔关于x轴对称，这些系数均为实数。如果关于y轴对称，\(c_k\)只有下标为奇数的项，而\(u_k\)只有下标为偶数的项。将前式两边乘上\(-n\)，并应用前述简化得到
\[ u_{2n}-\sum_{k=1}^{n-1} (2k-1)c_{2k-1}u_{2n-2k}-(2n-1)\sum_{k=n+1}^\infty c_{2k-1}u_{2k-2n}=\frac{p}2(2n-1)c_{2n-1} \]
同样地，\(n=1\)时，右边另加上\(-p/2\)。这便是本文的核心等式。
以椭圆孔为例
\(\omega(\zeta)\)的展开仅有一项系数时，\(\omega(\zeta)=\zeta+c_1\zeta^{-1}\)，此时它将单位圆映射为长半轴在x轴上的椭圆，半轴长为\(1+c_1\)和\(1-c_1\)。\(c_1\)的值在0和1之间，取0时形状就是圆，而取1时退化为宽为0的缝。由前面的等式
\[ \begin{aligned} u_2&=\frac p2(c_1-1)\\ u_{k+2}&=c_1u_k \end{aligned} \] 于是 \[ \Phi(\zeta)=\frac p2(c_1-1)\zeta^{-2}(1+c_1\zeta^{-2}+\cdots)=\frac p2(c_1-1)\frac{\zeta^{-2}}{1-c_1\zeta^{-2}} \]
\(\sigma^1_t=p+4\mathrm{Re}\left(\Phi(\zeta)\right)\)，最小值恒为\(-p\)，最大值在\(\zeta=\pm i\)处，等于\(\frac{3-c_1}{1+c_1}p\)。在\(c_1=0\)时，就是前面已知的圆孔时应力集中系数为3，随着孔变椭，应力集中系数降低，到了\(c_1=1\)这一退化情形，则没有应力集中现象。
由前式进一步求解，\(\phi'(\zeta)=\frac p2(c_1-1)\zeta^{-2}\)，\(\phi(\zeta)=\frac p2(1-c_1)\zeta^{-1}\)
\[ \begin{aligned} \phi^1(\zeta)&=\frac p4(\zeta+c_1\zeta^{-1})+\frac p2(1-c_1)\zeta^{-1}\\&=\frac p4 \zeta+\frac p4(2-c_1)\zeta^{-1} \end{aligned} \]
以上结果与《弹性理论》中用其他方法对椭圆孔的计算相同。
逆问题的求解
根据前面的等式，已知\(\Phi(\zeta)\)的情况下便可求得\(\omega(\zeta)\)，但得到图形可能会自重叠，也就是\(\omega(\zeta)\)在单位圆的外部并不是单射。要避免此情形，就要求\(\omega'(\zeta)\)在单位圆外不为0。从\(\omega(\zeta)\)的级数表示，求\(\omega'(\zeta)\)的零点可转为求多项式的根，而所有根的模都必须小于1。通过对一些逆问题的计算，本人发现，对于某个\(\Phi(\zeta)\)，只有一小段区间内的\(s\)可从\(s\Phi(\zeta)\)得到可接受的\(\omega(\zeta)\)。随着级数项次的增加，区间会变得很小。这使人猜想，对于某个\(\Phi(\zeta)\)，只有唯一的\(s\)使得\(s\Phi(\zeta)\)是有物理意义的。可是这里出现了本文的大缺陷，不要谈证明，本人甚至无力用严格的方式表述这一猜测。于是本文只能是用数值方法寻找合适的\(s\)：由\(s\Phi(\zeta)\)得到\(\omega(\zeta)\)，再计算\(\omega'(\zeta)\)的根，得到根的最大模，以数值优化的方法找到令最大模最小的\(s\)。
最优应力分布
应力的优化应减小最大应力，使得应力尽量均匀。最优的情况当然是应力处处相等。在本问题中，\(\Phi(\zeta)\)在闭曲线上的积分为0，如果处处相等，只能是\(\Phi(\zeta)=0\)，此时\(\omega(\zeta)=\zeta+\zeta^{-1}\)，就是前面提过的宽为0的缝，这就是没有其他约束下的最优孔形，此时完全没有应力集中现象。
参考前面对椭圆孔的求解，可以知道，孔在y方向和x方向上的最大尺寸的比值越小，应力集中的程度越低，所以更有意义的优化目标应是，在此比值一定的情况下，寻找应力集中程度的最小的孔形。既然问题是从圆孔引出，本文就先在两个方向上的最大尺度相同的情况下寻找最优。此时应力必然有正有负，要应力尽量均匀，最极端的情况便是，正值都是同一个值，而负值也是如此，而按前述的路径积分，均值应为0。只看第一象限，\(\theta\in[0,\frac\pi2]\)，\(\zeta=e^{i\theta}\)处的应力可由在0和1之间的\(a\)和大于0的\(b\)表示为
\[ \frac{\sigma_x+\sigma_y}{p}=\begin{cases} b &\theta>a\frac\pi2 \\ b\frac{a-1}a &\theta \leq a\frac\pi2 \end{cases} \]
\(\Phi(\zeta)\)展开的系数只有偶数项，可通过傅立叶展开求得
\[ u_{2k}=\frac{p}{4}\frac{4}{\pi}\left(\int^{a\pi/2}_0b\frac{a-1}a\cos 2k\theta \mathrm{d}\theta+\int^{\pi/2}_{a\pi/2}b\cos 2k\theta \mathrm{d}\theta\right)\\ =-\frac{p}{4}\frac{4}{\pi}\frac b{2ka}\sin 2ka\pi \]
最优孔形
有了\(u_{2k}\)的表达式，依照前面逆问题的求解方法，对于某个\(a\)的值，可以找到\(b\)的值，得到有物理意义的孔形。计算孔在y方向和x方向上的最大尺寸的比值，就得到在此比值下的最优孔形。以比值等于1作为目标，用数值优化方法寻找\(a\)，求得\(a=0.387\)，\(b=1.13\)。下图右边就是本文得到的最优孔形，图中的颜色表示的是正应力之和，而左边是圆孔的情形。可以看到，得到的孔形通过使得受拉区域和受压区域各处于均一应力下，明显地减小了最大拉应力，而最大压应力也改善至圆孔的0.79倍。
由\(b\)的值可知此孔形的应力集中系数为2.13，比文献1中的结果2.15提高了0.02。文献1中优化得到的鼓形孔可使受拉区域的应力均等，本文的孔形令左右两边略向内凹，进一步使受压区域的应力均等，因为应力沿边界的积分为定值，故而改进压应力分布亦可减小最大拉应力，但从结果上看改进是很微小的。文献1中还对在y方向和x方向上的最大尺寸的比值更小的几种情况给出了最优结果。可以想到，此时受压区域的占比更小，本文的优化所能带来的效果也会更小，所以本文不再做相关的计算。
感想
最终结果带来的改进是微小的，相较于这点改进，本人更愿意写两点感想。一是今日已极少使用的老方法在某些特定问题上仍然很有魅力。这方法实现起来比基于有限元的优化简单很多，数学意义更为明确，虽然有重要步骤本人未能透彻理解。二是得到的结果是有尖角的，尽管本文的级数近似并不能真正表示尖角。书中讲到应力集中时，总是说要避免尖角，平滑过渡。至少从数学上看，从高应力区跳变到低应力区才有利于减小最大应力，而不是平滑过渡。
结语
统而言之，本文采用一种今日已较为小众的方法对带孔的无限大平板受拉问题进行了分析。这一方法不仅可以由孔的形状计算应力，还可以由应力计算孔的形状，因此本文可以直接由最优的应力分布求得最优孔形。得到的结果确实能超过前人的成绩，但提升是微小的。不过本人自认为此过程仍有一定价值，故撰文以示众人。
参考文献
- Burchill, M., & Heller, M. (2004). Optimal free-form shapes for holes in flat plates under uniaxial and biaxial loading. The Journal of Strain Analysis for Engineering Design, 39(6), 595-614.
- Bjorkman, G. S., & Richards, R. (1976). Harmonic holes—an inverse problem in elasticity. Journal of Applied Mechanics, 43(3), 414-418.
- Vigdergauz, S. B., & Cherkayev, A. V. (1986). A hole in a plate, optimal for its biaxial extension-compression. Journal of Applied Mathematics and Mechanics, 50(3), 401-404.
- Gray, C. A. M. (1951). Polynomial approximations in plane elastic problems. The Quarterly Journal of Mechanics and Applied Mathematics, 4(4), 444-448.
- 铁摩辛柯, 古地尔. 徐芝纶译. 弹性理论.