用于控制约束线性弹性方程的原始和混合无锁定非协调虚拟单元方法

《Computers & Mathematics with Applications》：Primal and mixed locking-free non-conforming virtual element methods for control constrained linear elasticity equations

　　

**摘要**

本文提出了针对控制约束线性弹性方程原始形式和混合形式的非协调虚拟单元方法（VEMs）。我们遵循"先优化后离散化"方案，对控制变量采用变分离散化方法，而状态方程和伴随方程则采用VEM离散化进行近似。我们证明了所提方法在准不可压情形下均匀收敛，并推导了状态变量和伴随变量在L2范数和H1范数下的最优收敛估计，以及最优控制在L2范数下的收敛估计。数值实验验证了理论结果，并展示了该方法在一般多边形网格上的无锁定性能。

**引言**

最优控制问题（OCPs）广泛出现在科学和工程应用中，如导航、渠道流工程、晶体生长、热传递以及其他目标导向过程。在这类问题中，物理系统的行为通过适当的控制措施加以影响，以优化预设的性能准则。这种将复杂过程导向预期结果的内在需求为最优控制问题的数学分析和数值逼近提供了强有力的动机，过去几十年间已发展出广泛的数值方法；例如参见[1]、[2]、[3]、[4]、[5]、[6]、[7]、[8]、[9]、[10]、[11]、[12]、[13]、[14]、[15]、[16]及其中的参考文献。

在固体力学的背景下，由线性弹性方程控制的最优控制模型尤为重要，因为它们允许系统地设计力或激励策略，在物理限制条件下将弹性体驱动至预设构型；例如参见[17]、[18]。这类模型出现在弹性及弹塑性结构的变形控制、生物力学中的反向加载和力的重构问题，以及具有严格规格要求和有界激励的制造过程中。在这些场景中，分布体力控制提供了一种自然且灵活的建模框架，而逐点控制约束则考虑了执行器饱和、安全限制或生理边界[18]。

设Ω?R2为边界为?Ω的凸多边形区域。在第一种情况下，对于期望位移yd∈L2(Ω)，考虑如下最小化问题：

$$\min_{(y,q) \in V \times Q_{ad}} J_1(y,q) = \frac{1}{2}\int_\Omega |y - y_d|^2 dx + \frac{\kappa}{2}\int_\Omega |q|^2 dx,$$

受限于以下平面线性弹性问题（以原始位移形式表述）：

$$-\text{div}(2\mu\epsilon(y) + \lambda(\text{div}(y))I) = f + q \quad \text{在} \Omega \text{中},$$
$$y = 0 \quad \text{在} \partial\Omega \text{上},$$

**情形1**

作为第二种情形，对于期望压力pd∈L2(Ω)，考虑如下优化问题：

$$\min_{(y,p,q) \in V \times P \times Q_{ad}} J_2(y,p,q) = \frac{1}{2}\int_\Omega |y - y_d|^2 dx + \frac{1}{2}\int_\Omega |p - p_d|^2 dx + \frac{\kappa}{2}\int_\Omega |q|^2 dx,$$

受限于Herrmann（混合）形式[19]表述的线性弹性问题：

$$-\text{div}(2\mu\epsilon(y)) + \nabla p = f + q \quad \text{在} \Omega \text{中},$$
$$\text{div}(y) + \frac{1}{\lambda}p = 0 \quad \text{在} \Omega \text{中},$$
$$y = 0 \quad \text{在} \partial\Omega \text{上},$$

**情形2**

其中，V和P分别为位移空间和压力空间，将在后文定义，可容许控制空间Qad定义为：

$$Q_{ad} := \{v \in L^2(\Omega) : a \leq v(x) \leq b \text{ 在} \Omega \text{中几乎处处成立}\},$$

其中a, b∈R2，不等式按分量理解，κ > 0为控制代价参数，f∈L2(Ω)为源函数。对于v=(v?, v?)，我们定义以下记号：

$$\epsilon(v) = \frac{1}{2}(\nabla v + (\nabla v)^T),$$
$$\text{div}(v) = \frac{\partial v_1}{\partial x_1} + \frac{\partial v_2}{\partial x_2},$$

其中div(·)表示张量场的散度。Lamé参数μ和λ定义为：

$$\mu = \frac{E}{2(1+\nu)}, \quad \lambda = \frac{E\nu}{(1+\nu)(1-2\nu)},$$

其中E > 0表示杨氏模量，ν∈(0, 0.5)表示泊松比。极限λ→∞对应准不可压情形，这已被证明会带来重大的分析和数值挑战。混合形式已在[20]、[21]、[22]及其参考文献中得到广泛研究。

为该问题设计数值方法的一个主要困难是所谓的"锁定现象"，通常在λ变得非常大时（即ν→1/2时）出现，对应于准不可压材料。在这种情况下，标准离散化方法可能会丧失精度甚至无法收敛。例如，分片线性协调有限元方法（CFEM）在λ→∞时不能均匀收敛。为恢复与λ无关的最优收敛率，CFEM需要三次以上的多项式阶[23]。相比之下，[23]中表明，至少二阶的非协调FEM能够关于λ获得最优的均匀收敛率。然而，线性的非协调元不能直接应用，因为在这种情况下离散Korn不等式不成立。其结果是，相关的离散双线性形式可能不具强制性。为克服这一困难，Hansbo等人提出了一种稳定化的线性非协调FEM，通过引入附加的稳定化项恢复强制性[24]。

虚拟单元方法（VEM）在[25]中被引入，作为经典FEM向一般多边形或多面体元的自然扩展。在每个网格单元上，局部逼近空间包含直到预定阶数的所有多项式，并用额外的非多项式函数进行丰富。VEM的一个关键特性是不需要显式地给出形函数形式。相反，该方法仅依赖于精心选择的自由度，例如网格顶点处的值、边或面上的矩，以及多边形或多面体单元上的矩。这使得方法的构造既灵活又稳健。此外，VEM可以方便地扩展到高阶逼近。凭借这些优势，VEM已被成功应用于广泛的问题；例如参见Adak等[26, 27]、Ayuso de Dios等[28]、Beir?o da Veiga等[29, 30]、Bertoluzza等[31]、Boffi等[32]、Cangiani等[33]、Credali等[34]、da Veiga等[35]、Dassi等[36]、Foligno等[37]、Huang和Wang[38]、Marcon和Mora[39]、Mora等[40]、Mora和Rivera[41]、Mora等[42]、Zhao等[43]及其参考文献。

在[35]中，提出了用于弹性问题的至少二阶的协调VEM，并建立了关于Lamé参数λ的最优均匀收敛率。随后，在[44]中引入了至少二阶的非协调VEM，同样证明了关于λ的均匀且最优的收敛性。然而，最低阶的非协调VEM仅考虑了纯位移形式，因为最低阶情形下离散Korn不等式可能失效。最近，Kwak等人在[45]中发展了两类最低阶VEM，证明了这些方法在准不可压情形下均匀收敛，并在精确解具有适当正则性假设下获得最优收敛率。在[46]、[47]、[48]、[49]、[50]、[51]、[52]、[53]中有很多使用协调和非协调虚拟单元方法求解Stokes和Navier-Stokes方程的工作。[21]中的作者提出了一种新的VEM，用于准不可压线性弹性问题的数值近似。在他们的方案中，位移和乘子均使用等阶或不等阶的虚拟单元空间进行逼近，所得方案被证明是无锁定的。但到目前为止，Herrmann形式表述的线性弹性方程尚未使用非协调VEM进行研究。

近年来，虚拟单元方法已成为求解最优控制问题的灵活而强大的离散化框架，其有效性已在多项研究中得到证实；例如参见Li等[8]、Tushar等[14]、Wang和Zhou[15]、Brenner等[54]、Tushar等[55, 56, 57, 58]、Wang和Zhou[59, 60]、Yang等[61]。尽管VEM在广泛的偏微分方程和优化问题中的应用取得了成功，但其在控制约束线性弹性方面的应用却鲜有关注。现有的关于弹性VEM的结果主要局限于施加预设载荷的正向问题，而最优控制模型通过引入伴随弹性方程和与逐点控制约束相关的变分不等式，从根本上改变了问题的分析结构。在此设置中，必须同时确保状态问题和伴随问题关于Lamé参数的均匀稳定性，且离散控制的精度取决于底层离散化的参数鲁棒性。据我们所知，控制约束线性弹性特别是最低阶非协调情形及准不可压情形下的系统性VEM分析尚不可用。

下文总结本文的主要贡献和创新之处：

- **原始形式和混合形式的VEM离散化**：我们采用Kwak等人[45]提出的最低阶无锁定虚拟单元方法来离散化原始形式和混合形式的状态方程和伴随方程。控制变量采用变分离散化方法进行近似，并结合"先优化后离散化"策略，这允许隐式处理控制约束，而无需引入额外的离散控制空间。

- **参数鲁棒的先验误差分析**：在标准网格正则性假设下，我们推导了关于Lamé参数λ的最优先验误差估计。特别是，我们建立了状态变量和伴随变量在L2范数和离散H1范数下的收敛率，以及最优控制在L2范数下的收敛率。

- **一般多边形网格上的数值验证**：我们在多种凸和非凸多边形网格上进行了数值实验，以验证理论误差估计。数值结果进一步证明了所提方案在大λ值（对应准不可压情形）下的无锁定行为。

本文余下部分结构如下：在第2节中，我们引入最优控制问题分析所需的记号和功能空间，并讨论适定性及最优性条件。在第3节中，我们引入非协调虚拟单元方法来离散化连续问题，并推导原始形式的最优先验误差估计。第4节讨论混合形式的非协调VEM及其先验误差分析。最后，在第5节中，我们给出一些数值实验和基准问题以验证理论估计。

**功能空间与适定性**

设Ω?R2为边界为?Ω的凸多边形区域。对于任意给定开子集D?R2，Hs(D)表示具有范数∥·∥s,D和半范|·|s,D的标准Sobolev空间。当s=0时，我们定义H0(D):=L2(Ω)，配备标准的L2内积(·,·)D和范数∥·∥0,D。粗体记号用于表示相应的向量值空间。若D=Ω，则省略下标。我们将问题中使用的空间V和P定义为：

$$V := \{v \in H^1(\Omega)$$

**情形1的非协调VEM离散化**

本节中，我们提出一种非协调VEM（NCVEM）来离散化前文推导的最优性系统。我们采用隐式的变分离散化方法来近似控制变量。

**情形2的非协调VEM离散化**

本节中，我们在第3.2节引入的多边形网格假设和Kouhia-Stenberg型NCVEM空间下离散化混合最优性系统(2.18)。

除了局部位移虚拟空间V(T)外，对压力变量我们定义离散空间为：

$$P_h := \{\psi_h \in L^2(\Omega) : \psi_h|_T \in P_0(T) \ \forall T \in \mathcal{T}_h\}.$$

对于每个T∈Th，定义局部离散双线性形式bhT: V(T)×P0(T)→R和chT: P0(T)×P0(T)→R为：

$$b_h^T(v_h, \psi_h) := -(\Pi_{0,T}^0(\text{div}_h(v_h)), \psi_h)_T,$$
$$c_h^T(p_h, \psi_h) := \frac{1}{\lambda}(\Pi_{0,T}^0 p_h, \Pi_{0,T}^0 \psi_h)_T,$$

其中divh

**数值实验**

本节中，我们在单位正方形区域Ω=(0,1)2上进行两个数值实验，并在L形区域上进行一个实验，控制代价参数κ=1，以支撑我们在图1所示一般多边形网格上的理论结果。所有数值实验均在配备Apple M2处理器、8GB内存和256GB存储的MacBook上，使用MATLAB R2023b完成。离散最优性系统使用投影梯度算法[7]求解。我们用L2范数误差表示状态

**结论**

本文中，我们研究了一种非协调无锁定虚拟单元方法，用于求解以原始形式和混合形式表述的平面线性弹性方程控制的最优控制问题。使用控制的变分离散化方法，推导了状态变量和伴随状态变量在L2范数和H1范数下以及最优控制在L2范数下的最优先验误差估计。所有估计均通过不同族的数值结果得到了验证。

Harpal Singh | Sumit Mahajan | Arbaz Khan

印度罗尔基印度理工学院数学系，Roorkee, 247667, Uttarakhand, India