摘要

弹性力学中常研究含椭圆边界的区域，以阐明圆不对称性的作用，这类问题还关联着诸多经典的闭合形式解。本文提出一种统一的半解析方法，求解各类椭圆几何变体（如椭圆圆柱、共焦椭圆环、无限平面中的椭圆孔）在任意合理边界条件下的应力与位移场。该方法采用椭圆坐标系下截断为有限项的艾里应力函数广义表达式，并通过配点法推导相关系数以满足边界条件。通过求解多种算例并结合独立的有限元模拟验证，证明了该方法的正确性与有效性。

引言

在弹性力学中，含椭圆边界的二维区域是一种常用的几何模型，因为它便于分析偏离圆对称性的影响，并且存在明确的闭合形式解的可能性 [1-3]。例如，Inglis [4,5] 关于无限平面在远场载荷作用下椭圆孔对其应力场影响的研究，以及二维Eshelby夹杂问题 [6]，就是椭圆几何具有深远影响的两个典型例子。

一种强大且被广泛采用的椭圆线性弹性域解析解构造方法是复变函数法（CVM） [5,7-9]。其求解过程包括：i) 通过保角映射将椭圆边界转换为圆边界；ii) 在圆域上求解目标问题；iii) 最后将解转换回原椭圆几何。通过该方法，椭圆域的复杂性被转化为圆域上修正后的边界条件 [5,7-12]。相关文献众多，难以一一列举，但复变函数法确实已被应用于多种载荷与边界条件，包括利用解析延拓处理无限平面中椭圆孔/夹杂以及实心椭圆圆柱的混合边界条件（MBCs）。在椭圆几何背景下，复变函数法的一个重要应用是推导无限弹性平面中椭圆孔表面受集中力作用的基本解，并用于研究不同尺寸和形状孔洞之间的相互作用 [13]。然而，复变函数法与解析延拓的强大组合无法应用于有限椭圆环域问题 [14]。

另一种寻求解析解的途径是在椭圆坐标系中表达弹性力学问题，并尝试推测合适的应力势，即著名的艾里应力函数。选择艾里应力函数作为候选函数，源于其控制方程（即椭圆坐标系下的双调和方程）的解，该解由Timpe [15] 给出，并记录于文献 [3,16] 中。利用艾里应力函数求解无限平面中椭圆孔问题的闭合形式解示例可见于文献 [1,3,16]。然而，将艾里应力函数方法应用于含椭圆边界的有限域的研究却寥寥无几。本文中，艾里应力函数被表示为Timpe解的线性组合，通过求解系数以满足包含应力和/或位移的边界条件。与圆域的Michell解不同，Timpe的解并非三角函数的线性形式，因此无法利用正交性原理唯一确定系数。取而代之的是，在若干后续称为“配点”的点上选择性地满足边界条件。通常，配点的数量远多于未知系数的数量，从而得到一个超定方程组，并通过最小二乘法 [17]、超定配点法 [17-19,20-23] 求解。Goree [17] 采用最小二乘法处理了含光滑刚性椭圆嵌入物的无限板在均匀远场载荷作用下的接触应力确定问题。Jones和Hozos [21] 利用超定配点技术分析了有限薄矩形板在单向和双向载荷作用下椭圆孔周围的应力分布。Alexandrakis [22] 通过考虑双向载荷作用下薄矩形板中的多个椭圆孔，研究了孔间相互作用对应力集中的影响。Velazquez和Kosmatka [24] 推测了艾里应力函数的形式，并推导了半椭圆曲梁在端部受横向力作用下的应力场，使面力边界条件在积分意义上得到满足。

显然，艾里应力函数方法在含椭圆边界的有限域中的应用，尚未在所有可能的边界条件和几何变体方面得到充分探索。现有文献与当前研究之间的差距，源于计算机技术的巨大进步以及有限元法等稳健数值方法的发展，这也导致了人们对弹性力学边界问题半解析方法的兴趣相对不足。但目前，MATLAB®、Maple®、Mathematica®等复杂且易用的软件及其开源版本的出现，正促使人们重新审视经典弹性力学问题，这也正是本文研究的动机。

论文剩余部分的内容安排如下：下一节将介绍几何模型、符号定义，并建立求解椭圆域线性弹性问题所需的数学框架与数值细节。在第3节中，将通过一系列算例验证本文方法的有效性。最后，在结论部分对当前工作进行总结。

问题表述

图1展示了一个由共焦椭圆面$\epsilon_1$和$\epsilon_2$所界定的线弹性、各向同性、均匀平面域。$x$轴和$y$轴分别为椭圆$\epsilon_1$和$\epsilon_2$的长轴与短轴。为便于说明，图中显示内侧曲面受法向和面力$p$和$q$作用，而外侧曲面受混合边界条件作用。考虑到边界形状以及对半解析提取应力与位移的需求，采用图2所示的曲线椭圆$\epsilon-\eta$坐标系最为合适。

图1：线性弹性域由共焦椭圆曲面1和2分别界定，其中x轴和y轴分别为长轴和短轴。作为示例，上图显示前一个曲面受到牵引边界条件，而后一个曲面施加混合边界条件。

图2：椭圆（$\epsilon-\eta$）坐标系

椭圆坐标$(\epsilon,\eta)$与笛卡尔坐标$\{x,y\}$的关系为 [1,3]：

其中$c$为具有长度量纲的实常数，$\epsilon$为正实数，$\eta$为取值范围在$[0,2\pi]$内的实数。当$\epsilon$为常数时，从式(1)中消去$\eta$可得椭圆方程：

其焦点位于$x=\pm c$处。类似地，当$\eta$为常数时，从式(1)中消去$\epsilon$可得双曲线方程：

该双曲线具有相同的焦点$\pm c$。如图2所示，改变$\epsilon$（或$\eta$）可得到共焦椭圆（或双曲线）。

图1所示的平面椭圆域假设处于准静态平衡状态。假设产生的变形为微小变形，其形式可根据物体的面外尺寸与面内尺寸的相对大小，以及载荷在面外方向的变化，表现为平面应力或平面应变状态。设$\{u_\epsilon, u_\eta\}$分别表示位移在$\epsilon$和$\eta$方向的分量。对于平面应变变形，面外位移$u_z$为零；而对于平面应力情况，面外位移可由面内位移推导得出。设$\{\sigma_\epsilon, \sigma_\eta\}(\{\varepsilon_\epsilon, \varepsilon_\eta\})$和$\sigma_{\epsilon\eta}(\varepsilon_{\epsilon\eta})$分别表示椭圆坐标系下的面内法向和切向应力（应变）分量。在无体力的情况下，适用于平面应力和平面应变变形的、椭圆坐标系中沿$\epsilon$和$\eta$方向的非平凡平衡方程如下 [1,3]：

其中

为比例因子。应力分量可通过一个称为艾里应力函数$\phi(\epsilon, \eta)$的标量势表示，表达式如下 [1,3]：

这些表达式可使平衡方程（4）-（5）自动满足。在小变形假设下，椭圆坐标系中面内应变分量与位移分量的关系如下 [1,3]：

上述面内应变分量并非相互独立，而是受以下协调方程约束 [1,3]：

在线弹性、各向同性的椭圆域（图1）中，面内应力与应变分量服从广义胡克定律 [1]：

其中$E$为弹性模量，$\nu$为泊松比，$\kappa$为Kolosov常数，在平面应力和平面应变变形下分别等于$(3-\nu)/(1+\nu)$和$(3-4\nu)$。联立式(7)、(8)、(9)、(11)、(12)，可推导出椭圆坐标系下控制艾里应力函数$\phi$的著名双调和方程 [3]：

式(13)的广义解由Timpe [15]给出，并可在文献[3,16,21,22]中查阅。以下给出从文献[3]中采用的解的系统表格化形式：

其中$\{a_0, b_0, c_n, d_n, e_n, f_n, g_n, h_n, i_n, j_n\}$为基于边界条件确定的未知系数。应力分量$\{\sigma_\epsilon, \sigma_\eta, \sigma_{\epsilon\eta}\}$可通过将式(7)–(9)应用于式(14)中$\phi$的解来推导。随后，将应力代入式(12)可得面内应变分量，再将其代入式(10)并积分即可得到位移分量。积分常数不影响应变和应力，代表以下刚体项：

其中$\{A_0, B_0, C_0\}$为任意系数，以下称为刚体系数。作为示例，以下给出$\phi = \epsilon$时不含刚体项的应力与位移分量 [3]：

同理，式(14)中其余艾里应力函数对应的应力与位移分量表达式是存在的，但由于结果表达式过长，此处不再列出。感兴趣的读者可参考文献[3]中清晰的表格化内容。

回顾图1，需再次强调，本文旨在建立一个统一的半解析求解框架，以求解椭圆域中的应力与位移场，而无论边界条件的类型（面力、位移或混合边界条件）。对于圆域，Michell解表明，应力与位移的解可表示为正弦和余弦函数的线性组合。结合正弦与余弦函数的正交性，已证明对于圆域上的各类边界条件，均可得到闭合和/或半解析解 [25,26]。椭圆坐标$\epsilon$和$\eta$类似于圆坐标的半径$r$和角度$\theta$。此外，势函数$\phi$及相关的应力与位移分量也包含$\eta$的余弦和正弦函数。但文献[25,26]中针对圆域的方法无法直接应用于椭圆域，原因很简单：$\phi$及相应的应力与位移分量中并不存在线性形式的余弦和正弦函数，这一点可通过式(16)的示例得到验证。下一小节将详细介绍一种数值方法，以弥补椭圆域中余弦和正弦函数正交性无法直接应用的不足。

数值方法

作为第一步，式(14)中$\phi$的广义解中的无穷级数在$n=N$处截断，使得包括刚体系数在内待确定的系数总数为$(8N+5)$。由此，可得到包含待定系数的位移与应力分量的截断表达式。求解这些表达式的目的是在最大程度上满足给定的边界条件。边界条件分为三类：

1. 面力边界：整个表面的面力已知。根据柯西定理 [3]，椭圆面$\epsilon = \{\epsilon_2, \epsilon_1\}$的法向与切向面力分量分别等于$\pm\sigma_\epsilon$和$\pm\sigma_{\epsilon\eta}$。因此，椭圆边界上的面力边界条件可简化为分别指定$\sigma_\epsilon$和$\sigma_{\epsilon\eta}$分量。
2. 位移边界：整个表面的位移$\{u_\epsilon, u_\eta\}$已知。
3. 混合边界：在给定表面的一部分上指定位移$\{u_\epsilon, u_\eta\}$，其余部分指定面力$\{\sigma_\epsilon, \sigma_{\epsilon\eta}\}$。

在曲面$\epsilon = \{\epsilon_1, \epsilon_2\}$上施加上述任意一种边界条件，都会生成四个关于变量$\eta \in [0, 2\pi)$的代数方程，其中包含$8N+5$个未知系数。这些方程在内外表面上的$m_i$（$m_o$）个等间距$\eta$值处进行计算。在引言部分，这些选定的$\eta$值被称为配点，记为$n_c$。对于混合边界条件，将区间$[0, 2\pi)$按面力和位移的指定区域分成两部分并各自等分，使得配点总数根据表面是内表面还是外表面分别达到$m_i$或$m_o$。因此，针对图1的内外边界条件，共生成$2(m_i + m_o)$个方程，所得方程组可写成矩阵形式：

其中$[A]$为大小$2(m_i + m_o) \times (8N + 5)$的矩阵，$\{X\}$为包含未知系数的大小$(8N + 5) \times 1$的向量，$\{b\}$为包含$\eta = n_c$处边界（面力或位移）信息的大小$2(m_i + m_o) \times 1$的向量。过去，式(17)所代表的超定方程组通常采用最小二乘法或多元线性回归（MCR）求解 [17-23]。但研究发现，当边界条件存在不连续性时，MCR方法的解精度不足。相比之下，基于MATLAB® 中mldivide求解器 [28] 得到的式(17)的解精度更高，且在边界条件是否光滑或具有任意特性时均表现良好。因此，本文采用MATLAB® 中的mldivide求解器 [28] 进行计算，所有结果均基于该求解器。

$N$的取值受限于：在计算过程中，矩阵$[A]$中包含$N$的正指数项不应在计算平台上出现无穷大或NaN。另一方面，增大$m_i$和$m_o$的值显然能更精确地捕捉边界条件，但也可能因矩阵$[A]$的病态性和稀疏性等不利特性导致解$\{X\}$出现伪解。在载荷呈现对称性、域包含原点或为无限域的情况下，可对$N$的取值和艾里应力函数的选择进行优化。若载荷关于x轴和/或y轴对称，则无需考虑$\eta$的反对称函数；类似地，若域包含原点（或为无限域），则需剔除艾里应力函数中的负（或正）指数项，以保证应力和位移分量的有界性。

下一节将通过求解椭圆域（图1）上的多种边值问题，并结合闭合形式解析解（若存在）或在COMSOL® 中进行的独立有限元（FE）分析 [27] 验证结果，以展示本文方法的有效性。

结果和讨论

除非另有说明，本节讨论的所有问题均采用平面应力条件，弹性模量$E=1$，泊松比$\nu=0.3$。

已知椭圆坐标系（$\epsilon-\eta$）下的位移与应力场，可通过变换规则推导得到笛卡尔坐标系（$x-y$）下的对应分量，其中（$x,y$）与（$\epsilon,\eta$）通过式(1)关联，变换关系如下 [1,3]：

其中$\alpha=\tan^{-1}(\coth\epsilon\tan\eta)$，$\{u,v\}$和$\{\sigma_x,\sigma_y,\sigma_{xy}\}$分别表示笛卡尔坐标系下的面内位移和应力分量。上述关系可直接逆推，由笛卡尔坐标系下的信息得到椭圆坐标系下的位移与应力分量。

在接下来的小节中，几何模型将从含椭圆孔的无限平面过渡到有限尺寸的共焦椭圆环，最终到实心椭圆盘。对于每种几何模型，所处理的边界条件包括均匀、可变、不连续载荷，以及纯位移边界条件和混合边界条件。

图3：无限平面（$E=1$, $\nu=0.3$）含无面力椭圆孔，在远场等双轴单位面力作用下，椭圆孔表面（$\epsilon_1=1$, $c=1$）切向应力$\sigma_\eta$的变化。闭合形式解为Inglis解 [1,3]，有限元（FE）解通过COMSOL® [27] 获得。“本文研究”的解基于$m_i=m_o=4$, $N=2$。由于问题关于$x$轴和$y$轴对称，解仅展示$0<\eta<\pi/2$的部分。

含椭圆孔的无限平面

基于图1，令$\epsilon_2 \to \infty$即可得到无限平面。为保证$\epsilon \to \infty$时应力有界，式(14)中与系数$\{d_n, f_n\}$和$\{h_n, j_n\}$相关的级数项，其$N$值分别限制为2和1。

椭圆孔无面力且远场受等双轴拉伸

椭圆孔由$\epsilon_1 = 1$、$c = 1$描述（见式(2)），其无面力条件意味着：

远场（$\epsilon \to \infty$）作用单位大小的等双轴载荷，可得：

结合边界条件式(23)和(24)，式(14)中与系数$\{c_n, e_n, g_n, i_n\}$相关的级数项取$N=2$。每个边界的插值点数量相同，即$m_i = m_o = 4$，因此针对式(23)和(24)的边界条件共生成16个方程。对于后者的边界条件，需令$\epsilon_2 \to \infty$以生成关于未知系数的方程。组装后的方程组通过式(17)及前一节讨论的方法求解。

图3将本文计算得到的孔边$\epsilon = \epsilon_1$处切向应力$\sigma_\eta$的变化，与Inglis解析解 [1,3] 及COMSOL有限元计算结果进行了对比，吻合良好。本文结果用圆形标记表示，记为“Present Study”；解析解和有限元结果分别用实线和方形标记表示，记为“Closed Form Solution”和“FE Solution”。由于问题关于$x$轴和$y$轴对称，图3仅展示第一象限$0 < \eta < \pi/2$的结果。有限元计算同样利用了对称性，仅对平面应力条件下尺寸为$100 \times 100$的薄板的四分之一进行了COMSOL® 模拟 [27]。所选尺寸足够大，可反映原问题的无限域特征。网格包含6442个三节点三角形单元、3330个节点，最小单元尺寸为0.0075。四分之一椭圆边界采用7个单元离散。

图4：无限平面（$E=1$, $\nu=0.3$）中，椭圆孔表面（$\epsilon_1=1.3$, $c=1$）施加面力$\sigma_\epsilon=\cos(4\eta)$、$\sigma_{\epsilon\eta}=\sin(4\eta)$时，孔边的a 切向应力$\sigma_\eta$、b 位移分量$u_\epsilon$、c 位移分量$u_\eta$的变化。有限元（FE）解通过COMSOL\textregistered\ [27] 获得。“本文研究”的解基于$m_i=m_o=110$, $N=10$。由于问题关于$x$轴和$y$轴对称，解仅展示$0<\eta<\pi/2$的部分。

椭圆孔边界施加$\eta$的正弦与余弦形式连续面力

由$\epsilon_1 = 1.3$、$c = 1$描述的椭圆孔（见式(2)），其边界施加$\eta$的连续三角函数形式面力，对应应力分量需满足：

随着考察点远离孔表面，应力逐渐衰减，最终在远场满足：

需注意，所选面力载荷为自平衡载荷。式(14)中以$\{c_n, e_n, g_n, i_n\}$为系数的级数项取$N=10$，表面$\epsilon_1$和$\epsilon \to \infty$上的配点数量$m_i$和$m_o$均取110。由于载荷与几何的对称性，与前一案例类似，有限元计算仅考虑尺寸为$100 \times 100$薄板的四分之一。网格包含25,195个三节点三角形单元、12,805个节点，最小单元尺寸为0.002。四分之一椭圆边界采用8个单元离散。图4中本文方法的预测结果（圆形标记）与有限元解（方形标记）的高度吻合，进一步验证了本文方法的有效性。

椭圆孔表面的混合边界条件：孔的部分区域完全约束，其余区域无面力，且远场受x方向均匀面力

椭圆孔由$\epsilon_1 = 1.3$、$c = 1$描述（见式(2)），并施加以下混合边界条件：

无限平面在远场受x方向单位大小的均匀单轴面力作用，当$\epsilon \to \infty$时边界条件为：

式(14)中以$\{c_n, e_n, g_n, i_n\}$为系数的级数项取$N=150$，表面$\epsilon_1$和$\epsilon \to \infty$上的配点数量均取$m_i = m_o = 160$。在表面$\epsilon_1$的160个配点中，40个点用于位移约束区域$\eta \in [2\pi - 0.822, 2\pi] \cup [0, 0.822]$。有限元模拟采用尺寸为$200 \times 200$的完整矩形薄板，含椭圆孔（$\epsilon_1 = 1.3$, $c = 1$）。网格包含28,388个三节点三角形单元、14,411个节点，最小单元尺寸为0.06。椭圆边界采用32个单元离散。

从图5中孔表面（$\epsilon_1 = 1.3$, $c = 1$）的位移分量$u_\epsilon$、$u_\eta$和应力分量$\sigma_\epsilon$、$\sigma_{\epsilon\eta}$的曲线可知，本文解正确反映了混合边界条件式(27)–(28)。约束区域$\eta \in [2\pi - 0.822, 2\pi] \cup [0, 0.822]$内应力解的振荡是三角函数展开的固有缺陷，可通过对解$\{X\}$施加滤波方案消除（见[25,26]）。在边界条件性质发生变化的过渡点（$\eta = \{0.822, 2\pi - 0.822\}$）处，应力解的不连续性和突变是弹性力学解中奇异性的标志。在线弹性中，从位移边界到面力边界的过渡会导致振荡型平方根奇异性[29,30]。无论是有限元中的细网格，还是本文方法中改变$N$和$m$的取值，都无法消除该突变。实际上，突变的不收敛性可用于判断奇异性的存在。但除了检测奇异性，两种方法均无法确定奇异性的精确数学性质。

共焦椭圆环

通过为图1中界定域的椭圆截面指定有限值，使几何范围变为有限。具体而言，本小节处理的算例中，考虑$\{\epsilon_1 = 1.3, c_1 = 1\}$和$\{\epsilon_2 = 1.6, c_2 = 1\}$的情况。

图5：无限平面（$E=1$, $\nu=0.3$）中，椭圆孔表面（$\epsilon_1=1.3$, $c=1$）施加混合边界条件[式(27)-(28)]，并在远场受x方向单位大小的均匀单轴拉伸（式(29)）时，孔边面内分量：a 位移分量$u_\epsilon$、b 位移分量$u_\eta$、c 法向应力$\sigma_\epsilon$、d 切应力$\sigma_{\epsilon\eta}$、e 切向应力$\sigma_\eta$的变化。有限元（FE）解通过COMSOL\textregistered\ [27] 获得。“本文研究”的解基于$m_i=m_o=160$, $N=150$。

内边界施加不连续法向面力、外边界无面力

假设外表面$\epsilon = \epsilon_2$无面力，而内表面$\epsilon = \epsilon_1$（见图1）上切向面力为0，法向面力在部分边界上为均匀非零值1。因此，外表面的边界条件为：

内表面（$\epsilon = \epsilon_1$）的边界条件如下：

式(14)中所有级数项取$N=9$，表面$\epsilon_1$和$\epsilon_2$上的配点数量分别取$m_i=80$和$m_o=40$。有限元模拟中，考虑到关于$x$轴和$y$轴的对称性，构建并采用四分之一模型。网格包含7712个四节点四边形单元、8005个节点，最小单元尺寸为0.00155。

从图6可以明显看出，本文方法的结果与有限元计算结果高度吻合。法向面力的不连续性并未反映在图6所示内表面$\epsilon_1$的切向应力$\sigma_\eta$和位移分量$\{u_\epsilon, u_\eta\}$中。

外表面$\epsilon=\epsilon_2$的混合边界条件：部分边界完全约束，其余部分无面力，且内表面（$\epsilon=\epsilon_1$）施加面力

内表面$\epsilon=\epsilon_1$，$\eta \in [0, 2\pi]$施加面力载荷，对应边界条件为：

而外表面$\epsilon=\epsilon_2$施加混合边界条件：

式(14)中所有级数项统一取$N=11$，内表面$\epsilon_1$和外表面$\epsilon_2$的配点总数分别为160和120，后者在位移约束区域和面力自由区域中均匀分布。考虑到关于$x$轴和$y$轴的对称性，有限元模拟采用四分之一模型，网格参数与前一案例相同。图7对比了本文方法与有限元方法得到的内表面切向应力$\sigma_\eta$和位移分量$u_\epsilon$、$u_\eta$，由于对称性，结果仅展示$0 < \eta < \pi/2$的部分。

内表面$\epsilon=\epsilon_1$完全约束且外表面$\epsilon=\epsilon_2$施加反对称面力

内表面$\epsilon=\epsilon_1$完全约束，边界条件为：

外表面$\epsilon=\epsilon_2$施加反对称载荷，满足：

图6：共焦椭圆环（$E=1$, $\nu=0.3$）中，外表面（$\epsilon_2=1.6$, $c=1$）无面力（式(30)）、内表面（$\epsilon_1=1.3$, $c=1$）部分区域施加大小为1的均匀法向面力（式(31)-(32)）时，内表面的面内分量：a 位移分量$u_\epsilon$、b 位移分量$u_\eta$、c 法向应力$\sigma_\eta$的变化。有限元（FE）解通过COMSOL\textregistered\ [27] 获得。“本文研究”的解基于$m_i=80$, $m_o=40$, $N=9$。

式(14)中所有级数项统一取$N=9$，内表面$\epsilon_1$和外表面$\epsilon_2$的配点总数相同（$m_i = m_o$），均为60。有限元模拟采用完整模型，网格包含56330个三节点三角形单元、28845个节点，最小单元尺寸为0.00155。图8对比了本文方法与有限元方法得到的完全约束内表面的面内应力分量，结果仅展示$0 < \eta < \pi/2$的部分，利用反对称性可扩展至全范围。

图7：共焦椭圆环（$E=1$, $\nu=0.3$）中，内表面（$\epsilon_1=1.3$, $c=1$）施加式(33)的面力载荷、外表面（$\epsilon_2=1.6$, $c=1$）施加式(34)-(35)的混合边界条件时，内表面的a 切向应力$\sigma_\eta$、b 位移分量$u_\epsilon$、c 位移分量$u_\eta$的变化。有限元（FE）解通过COMSOL\textregistered\ [27] 获得。“本文研究”的解基于$m_i=160$, $m_o=120$, $N=11$。

实心椭圆盘

实心椭圆盘可通过在图1中令$\epsilon_1 = 0$得到（见式(2)）。

针对该几何模型，仅研究一个混合边值问题：外边界为椭圆$\epsilon_2$（$\epsilon_2 = 1.3$, $c_2 = 1$）的圆盘，其部分边界完全约束，其余部分施加对称面力，混合边界条件表示为：

在本文数值方法中，式(14)所有级数项取$N=20$，外边界$\epsilon_2$的配点数量取$m_o = 160$。由于问题的对称性，有限元计算仅采用实心盘的四分之一模型，网格参数与3.2.2小节的问题类似。外边界上$\eta = 0.822$的奇异点可从图9a中有限元解的$\sigma_\eta$突变直接识别，而在本文方法的解中，由于$N$取值较小，该突变并不明显，这一点也适用于图5c-e的结果。尽管如此，所选的$N$值足以准确预测奇异点附近的解，这可从面力加载区域$0 < \eta < 0.822$内，本文方法与有限元方法得到的切向应力$\sigma_\eta$、位移分量$\{u_\epsilon, u_\eta\}$的高度吻合得到验证。如前所述，本文方法得到的图9中$0.822 < \eta < \pi/2$范围内的振荡现象，是使用三角函数正弦和余弦的必然结果。

数值参数$N, m_i, m_o$的影响

所有算例中报告的$N, m_i, m_o$值，均确保了与解析解或有限元数值解的最佳吻合。为阐明这些数值参数对应力与位移场的影响，重新研究3.2.2节中处理的共焦椭圆环混合边值问题，其载荷与边界条件由式(33)–(35)给出。选择该算例的原因在于，它是一个综合性问题，覆盖了本文方法的所有方面。

以图7中采用的$N=11$、$m_i=160$、$m_o=120$作为基线参数，通过扰动这些参数来检查应力与位移场的变化。具体而言，与图7对比，评估内表面切向应力$\sigma_\eta$的变化。图10a展示了保持$m_i=160$、$m_o=120$不变时，改变$N$的影响：显然，随着$N$增大，结果收敛于期望的有限元解。图10b、c展示了保持$\{N=11, m_o=120\}$和$\{N=11, m_i=160\}$不变时，改变$m_i$和$m_o$对$\sigma_\eta$的影响：尽管$m_i$或$m_o$发生大幅变化，$\sigma_\eta$的变化却很小。其他场量也观察到类似趋势。

一般而言，针对特定问题确定$\{N, m_i, m_o\}$值的准则如下：

1. 对于给定的$\{N, m_i, m_o\}$，将式(17)得到的系数$\{a_0, b_0, c_n, d_n, e_n, f_n, g_n, h_n, i_n, j_n\}$代入位移与应力场表达式，检查边界条件是否被正确捕捉。若满足，则进行步骤3；否则，执行步骤2。
2. 调整$\{N, m_i, m_o\}$的取值，重复上述过程，直到边界条件被满意地捕捉。
3. 在边界条件被恢复后，计算特定位置（如应力集中位置）的应力和/或位移，或全局量（如单位面外长度的应变能），以确保结果收敛。若感兴趣的量存在显著偏差，则调整$\{N, m_i, m_o\}$的值。

图8：共焦椭圆环（$E=1$, $\nu=0.3$）中，内表面（$\epsilon_1=1.3$, $c=1$）完全约束、外表面（$\epsilon_2=1.6$, $c=1$）施加式(37)的反对称面力时，内表面的a 法向应力$\sigma_\epsilon$、b 切应力$\sigma_{\epsilon\eta}$、c 切向应力$\sigma_\eta$的变化。有限元（FE）解通过COMSOL\textregistered\ [27] 获得。“本文研究”的解基于$m_i=m_o=60$, $N=9$。

图9：椭圆盘（$E=1$, $\nu=0.3$）外表面（$\epsilon=1.3$, $c=1$）的$0 < \eta < \pi/2$部分，在式(38)–(39)所示混合边界条件下，a 切向应力$\sigma_\eta$、b 位移分量$u_\epsilon$、c 位移分量$u_\eta$的变化。有限元（FE）解通过COMSOL\textregistered\ [27] 获得。“本文研究”的解基于$m_o=160$, $N=20$。

本文采用等间距配点的选择，确保了算法的统一性，而与边界条件和/或载荷的性质无关。但缺点是，矩阵规模的增大会导致数值成本增加，最终延长计算时间。一种替代方案是采用沿边界非均匀间距的配点，这通常由边界条件发生剧烈变化的区域决定（例如混合边界条件下的过渡点、不连续载荷或载荷随位置呈复杂函数变化等情况）。重新回顾3.2.2节中讨论的混合边值问题算例，通过内表面切向应力$\sigma_\eta$分布的变化，评估外边界总配点数$m_o$的非均匀分配的影响。在外边界的位移[面力]指定区域，配点数分别取为$m_o(u)$和$m_o(\sigma)$，确保所有情况下均满足$m_o = m_o(u) + m_o(\sigma)$。计算采用相同的数值参数：$N = 11$，$m_i = 160$，$m_o = 120$。针对$(m_o(u), m_o(\sigma))$考虑三种情况：$(20, 100)$、$(60, 60)$和$(100, 20)$。中间情况代表配点沿整个边界均匀分布，其余两种情况分别对应面力指定边界和位移指定边界的偏态分布。需注意，在边界的给定部分（即位移或面力指定区域）内，配点仍为等间距排列。

图10：共焦椭圆环（$E=1$, $\nu=0.3$）中，内表面（$\epsilon=1.3$, $c=1$）施加式(33)的面力载荷、外表面（$\epsilon=1.6$, $c=1$）施加式(34)-(35)的混合边界条件时，数值参数a $N$、b $m_i$、c $m_o$的变化对其内表面切向应力$\sigma_\eta$的影响。数值参数基线为$N=11$, $m_i=160$, $m_o=120$，每张图中仅改变一个参数，其余两个参数取基线值。有限元（FE）解通过COMSOL\textregistered\ [27] 获得。

图11：共焦椭圆环（$E=1$, $\nu=0.3$）中，内表面（$\epsilon=1.3$, $c=1$）施加式(33)的面力载荷、外表面（$\epsilon=1.6$, $c=1$）施加式(34)-(35)的混合边界条件时，外边界位移指定区域和面力指定区域配点数$m_o(u)$和$m_o(\sigma)$的相对分配对其内表面a 切向应力$\sigma_\eta$、b 位移分量$u_\epsilon$、c 位移分量$u_\eta$的影响。数值参数为$N=11$, $m_i=160$, $m_o = m_o(u) + m_o(\sigma) = 120$。有限元（FE）解通过COMSOL\textregistered\ [27] 获得。

基于图11a可推断，应力分量$\sigma_\eta$对$m_o$配点的取值及分布并不敏感，但图11b、c中的位移曲线清楚表明，需要为面力控制的边界部分分配足够数量的配点$m_o(\sigma)$。因此，重新分配配点可有效获得精确且数值高效的解。沿边界给定部分采用非均匀间距配点是另一个值得探索的方向，留待未来研究。

总结

本文提出了一种统一的策略，用于求解由椭圆界定的线弹性、各向同性、均匀平面域在面力、位移或混合边界条件下的应力与位移场。研究涵盖了三种几何情况：含椭圆孔的无限平面、共焦椭圆环和实心椭圆盘。所有二维（平面应变、平面应力）分析均在椭圆坐标系中进行，以精确满足边界条件及其伴随的约束。

二维弹性问题采用艾里应力函数方法求解，以双调和方程为控制方程。椭圆坐标系下的广义解作为出发点，通过施加边界条件求解相关系数作为未知量。边界条件在沿边界等间距分布的配点上精确满足。除了每个边界的配点数量外，该方法还需要指定双调和方程级数解中考虑的项数$N$。根据边界条件的性质，调用艾里应力函数的广义解得到应力和/或位移分量的表达式，并生成线性方程组，通过MATLAB® 中的mldivide求解器 [28] 求解。该方法在所有椭圆几何和边界条件变体上进行了严格测试，结果通过COMSOL® [27] 中的有限元计算以及可用的闭合形式解析解进行了充分验证。

传统上，复变函数法结合保角映射被用于处理椭圆域问题 [9]，但该方法无法处理椭圆环的混合边值问题 [25]。此外，当边界条件变得复杂时，复变函数法与保角映射的组合代数运算会愈发繁琐，抵消了解析解形式带来的优势。有限元和边界元法等方法虽然通用，但属于纯数值方法，因此本文的半解析性质仍具有重要价值。

在近期的未来工作中，本研究可扩展至断裂和接触问题。借鉴文献[26]的思路，还可将该方法应用于非均匀介质。另一个研究方向是探索配点位置的优化，以获得数值精确且高效的解。