基于不连续Petrov-Galerkin方法的非线性弹性问题：最速下降法

《Computer Methods in Applied Mechanics and Engineering》：Nonlinear elasticity with the discontinuous Petrov-Galerkin method: The steepest descent method

　　乔纳森·张 | 莱谢克·德姆科维茨
美国奥登计算工程与科学研究所

摘要
本报告是对最小残差方法研究的延续，重点探讨不连续佩特罗夫-加莱金方法在非线性弹性问题中的应用。报告的主要目标是研究使用经典最速下降法求解最小残差问题，尤其将其与另一种可行方法——高斯-牛顿法进行比较。这两种方法的区别在于：高斯-牛顿法假定线性化算子满足上确界-下确界条件，而最速下降法则没有此要求。这一假设在处理某些实际重要问题（如结构屈曲）时尤为重要。本报告将这两种方法应用于二维非线性弹性问题的多凸储能函数的一个最成功的变分形式（即文献[1]中所述的公式IV）。该公式为鞍点型公式，将行列式?u+I视为另一个未知量，通过拉格朗日函数建立其与u之间的关系，并通过对积分的分部处理来放宽这一约束。报告还讨论了与迹未知量相关的实现细节。针对若干人工构造的解以及经典的库克膜基准问题，报告给出了数值结果。这些结果证实了最速下降法的线性收敛阶，但也表明其实际收敛速度不及高斯-牛顿法。

引言
预备知识
使用最速下降法（更一般地，基于梯度的方法）求解最小残差问题的思路并非新事物，本报告也绝非首次尝试。J. Perrier领导的法国团队已在[2][3]中发表相关研究成果，将其应用于跨音速流问题。此后，不连续佩特罗夫-加莱金方法作为一种具有诸多优点的强大有限元技术应运而生。德姆科维茨和戈帕拉克里希南[4]对线性问题的相关理论进行了深入研究，而卡尔斯滕森等人[5]则通过分析非线性混合问题，为所谓规则点发展出了更一般的非线性DPG方法理论。文献[5]的显著成果包括证明了混合问题与最小残差问题的等价性，以及给出了混合问题的先验误差估计。坎廷和霍伊尔[6]提出了适用于强单调算子的DPG框架，他们将非线性问题简化为带有额外未知量的线性问题，再将非线性关系视为约束条件，通过惩罚项来实现该非线性约束的强制。兹杜内克等人在其研究[7]中对弹性问题的当前最优布布诺夫-加莱金方法进行了详尽的综述。他们进一步开发了一种用于不可压缩有限弹性的新型四场布布诺夫-加莱金方法，并证明该方法具有压力鲁棒性。

本报告属于探索性质。我们的主要目标是研究最速下降法及其在DPG最小残差问题中的应用，而非与现有文献中的先进公式相竞争。在我们看来，现在就进行此类性能比较和宣称还为时过早。在之前的研究[1]中，我们主要致力于为二维非线性弹性问题的多凸储能函数子类开发和分析各种变分公式。核心思路是将某些量（如变形梯度的行列式）视为额外的独立未知量，并将其与原始解变量（如位移）建立关联。这一思路源于约翰·鲍尔[8]的理论工作，他在证明极小值存在性时正是依靠了这种关联方式。其他采用类似变量分离思想的研究还包括[7][9][10][11][12][13][14][15][16][17][18]。不过，这些公式的实现仅限于高斯-牛顿法。因此，自然会思考最速下降法是否也能带来优势。实际上，它早已在[19]等研究中被用于非线性扩散问题等领域。

本报告的其余部分结构如下：本节剩余部分将介绍最小残差问题、理想DPG方法、高斯-牛顿法以及最速下降法。第2节将介绍待分析的模型问题以及所选的变分公式。第3节将讨论与最速下降法实现相关的若干重要细节，并给出一些数值结果。如需简要了解DPG方法的实际应用方面，建议读者参阅附录B。

最小残差方法
关于最小残差方法的讨论从最一般的层面开始。我们需要解决的问题是求解一般方程B(x)=y，其中B: X → Y是两个实巴拿赫空间X与Y之间的一个可微映射，y ∈ Y。假设x ∈ X是方程B(x)=y的一个规则点，即x满足方程，且x处的弗雷歇导数是X到Y′之间的双射。在无限维度情况下，最小残差问题表现为寻找x∈arg min w∈X∥y?B(w)∥Y。根据假设，最小值必然存在，其值为零，且就在x处取得。对于某个（有限维）子空间Xh ? X，最小残差方法旨在通过寻找xh∈arg min wh∈Xh∥y?B(wh)∥Y，为真实解提供Xh中的近似值。显然，极小值点可能并非唯一。此外，还有一个隐含假设：残差的最小化实际上意味着误差的最小化，即如果x满足B(x)=y，则∥y?B(wh)∥Y→0?∥x?wh∥X→0。

在有限元领域，算子B通常由一个变分表述来确定，例如：{Find u∈U: b(u;v)=?(v) ?v∈V}，其中U表示试验空间，V表示检验空间，?:V→R是线性函数，而算子b:U×V→R可能在第一个参数上为非线性（用分号表示），在第二个参数上为线性。该变分表述直接对应着等价的算子问题：Find u∈U: B(u)=?，此处算子B: U → V′由〈B(u),v〉V′×V:=b(u;v)定义，对所有v ∈ V成立（?·,·?表示对偶配对）。对于这类问题，最小残差方法将试图找到uh∈Uh ? U，使得uh∈arg min wh∈Uh∥??B(wh)∥V′。对偶范数可通过上确界以标准方式定义，即∥φ∥V′:=sup∥v∥V=1|φ(v)|。但实际上，这个上确界往往难以计算。如果检验空间是一个带有内积(·,·)V的希尔伯特空间，那么里斯算子RV: V → V′使得最小残差问题可以用V中的范数来重写。该算子的特性表现为〈RVv,δv〉V′×V=(v,δv)V ?δv∈V，特别地，∥v∥V=∥RVv∥V′。结合范数的非负性，原问题（6）可以等价地替换为uh∈arg min wh∈Uh J(wh)，其中J(wh)=1/2∥RV?1(??B(wh))∥V2，J:U→R，这样的形式更便于处理，同时也能利用希尔伯特空间的结构。与现有的DPG文献[4]一致，我们定义符号ψ:=RV?1(??B(u))，并将其称为“误差表示函数”。

注1：即使检验空间不是希尔伯特空间，上述观点依然适用。如果V是一个类似Lp的空间，1 < p < ∞，那么逆里斯算子会被非线性对偶映射所替代。这种情况已在[20][21]中得到研究，并在[22][23][24][25][26]中应用于线性问题。其推广方式如下：设φ ∈ V′。如果V是希尔伯特空间，那么里斯代表元RV?1φ∈V就是最小化问题RV?1φ=arg min v∈V[1/2∥v∥V2?φ(v)]的解。在类似Lp的巴拿赫检验空间中，里斯代表元的类似物可以通过考虑(9)式的p阶版本来定义，即ψ:=arg min v∈V[1/p∥v∥Vp?φ(v)]。映射V′?φ?ψ ∈ V被称为对偶映射。需要注意的是，与RV的定义不同，对偶映射的定义方向是相反的。更多细节见附录C。

问题（8）的一阶最优性条件为〈duhJ,δu〉U′×U=?(RV?1(??B(uh)),RV?1duhBδu)V=0 ?δu∈U。为简洁起见，后续在方便时将使用符号b′(uh;δu,·):=〈duhBδu,·〉V′×V。上述方程中出现了两个逆里斯算子。通过保留其中一个并改变另一个以用于对偶配对，最小残差问题可以有两种不同的解释。

最优测试
通过定义vδu:=RV?1duhBδu为“与方向δu相关的最优测试函数”，最小残差问题可以被视为具有最优测试函数的非线性佩特罗夫-加莱金问题，即{Find uh∈Uh: b(uh;vh)=?(vh) ?vh∈Vhopt:=RV?1duhBUh}。需要指出的是，除非算子是线性的，否则最优测试空间会依赖于解uh。对于线性问题，“最优”这一术语源于这些测试函数能够使得上确界在上下界常数中达到最大值。在非线性情况下，这些函数则能使线性化算子duB的上确界常数达到最大值。事实上，在当前状态固定时，duB是一个从U → V′的线性算子，且supη∈V b′(u;ξ,η)∥η∥V=∥duBξ∥V′=∥RV?1duBξ∥V=(RV?1duBξ,RV?1duBξ)V∥RV?1duBξ∥V=(RV?1duBξ,vδu)V∥vδu∥V=b′(u;ξ,vδu)∥vδu∥V。

混合问题
另一方面，如果在(11)式中用误差表示函数的表达式替换，最小残差问题就会变成混合问题。此时条件(11)变为0=(ψ,RV?1duBδu)V=(RV?1duBδu,ψ)V=〈duBδu,ψ〉V′×V，而混合问题则为{Find ψ∈V,uh∈Uh(ψ,v)V+b(uh;v)=?(v) ?v∈V b′(u;δu,ψ)=0 ?δu∈Uh}。在实际应用中，DPG方法总是按照这种理解来编码实现。

注2：卡尔斯滕森等人在[5]中通过拉格朗日方法也提出了同样的混合问题。对于线性问题，理想佩特罗夫-加莱金方法的这三种理解是等价的[4]。“理想”一词指的是此时检验空间尚未被离散化。一旦V被替换为Vh ? V，这些方法就被称为“实用”的PG方法。在离散化后，这些理解仍然保持等价。重要的是，当算子B是线性时，duB=B，此时最优测试空间不再依赖于解。此外，最优测试函数的极大化性质意味着原始双线性形式的离散稳定性常数至少与连续形式相当。对于非线性问题，文献[5]已经证明了在规则点条件下，最小残差问题（8）与混合问题（15）是等价的。

无论采用哪种理解，当算子B是非线性时，都需要某种额外的线性化策略来解决该问题。为此，高斯-牛顿法和最速下降法成为两种自然的候选方法。以下段落将介绍这两种方法在最小残差问题中的应用，并分析它们之间的差异。

高斯-牛顿法
如[27][28][29]中所描述，高斯-牛顿法似乎是首选方法。给定一个初始点u0，该方法将算子B近似为B(u)≈B(u0)+du0BΔu，然后求解增量Δuh,n的最小化问题Δuh,n=arg min wh∈Uh[1/2∥??B(uh,n?1)?duh,n?1Bwh∥V′2]。这种最小化过程等同于使用牛顿-拉夫森法求解混合问题（15），即令b(uh,n;v)≈b(uh,n?1;v)+b′(uh,n?1;Δuh,n,v)，然后求解{Find ψ∈V,Δuh,n∈Uh(ψ,v)V+b′(uh,n?1;Δuh,n,v)=?(v)?b(uh,n?1;v) ?v∈V b′(uh,n?1;δu,ψ)=0 ?δu∈Uh}，从而得到下一个猜测值uh,n=uh,n?1+Δuh,n。因此，在这种情况下，这两个术语可以互换使用。需要注意的是，线性化后的混合问题（18）的左右两边都通过线性化算子b′依赖于之前的解猜测值uh,n?1。

最速下降法
应用最速下降法首先需要计算成本函数J的微分。通过初等代数运算以及对里斯算子的操作，可以得出〈duJ,δu〉U′×U=?(RV?1(??B(u)),RV?1duBδu)V=?(RV?1duBδu,RV?1(??B(u)))V=?〈duBδu,RV?1(??B(u))〉V′×V=?〈E[RV?1(??B(u))],duBδu〉V″×V′=?〈(duB)?E[RV?1(??B(u))],δu〉U′×U，其中E: V → V″是将希尔伯特空间V与其对偶空间V″对应起来的赋值映射，(duB)?: V″ → U′是转置算子。由于两者是一一对应的，因此通常会省略E[v]与v之间的区分。因此，微分结果为duJ=?(duB)?ψ∈U′。梯度是微分的Riesz表示，如果试验空间也是Hilbert空间，那么可以通过在试验空间中应用逆Riesz算子来计算它，即U??uJ=RU?1duJ=?RU?1(duB)?RV?1(??B(u))=?RU?1(duB)?ψ。现在，可以应用典型的最速下降更新方法（并配合合适的线搜索算法来确定步长），即un=un?1?αn?1?un?1J，其中αn?1>0。总体而言，对于该方法的每次迭代，都必须进行以下计算：1. 反演Riesz算子RV以计算误差表示函数ψ；2. 对问题进行线性化处理，即计算duB；3. 反演Riesz算子RU以计算梯度。关于线性问题的收敛性讨论，首先需要回顾对偶配对的定义。定义1：定义在两个Banach空间上的双线性（半双线性）形式b(u, v)，如果∥u∥U=sup∥v∥V=1|b(u,v)|且∥v∥V=sup∥u∥U=1|b(u,v)|，则称其为广义对偶配对。以下定理[4]、[19]描述了当试验范数和测试范数处于对偶配对关系时的收敛性特征。定理1：对于线性问题，如果测试范数和试验范数处于对偶配对关系，那么最速下降法（采用单位步长）无需初始猜测，即可一步得到最小残差解。如果算子B是线性的，那么梯度的表达式可简化为?uJ=?RU?1B?RV?1(??Bu)。此时，对于任意初始猜测u0，有u1=u0??u0J=u0+RU?1B?RV?1(??Bu)=RU?1B?RV?1?+(I?RU?1B?RV?1B)u0。如果试验空间采用能量范数，即∥u∥E=∥Bu∥V′，那么U中的Riesz算子可表示为(u1,u2)U=〈RUu1,u2〉U′×U=(Bu1,Bu2)V′=〈Bu2,RV?1Bu1〉V′×V=〈BTRV?1Bu1,u2〉U′×U，即RU=BTRV?1B。利用这个公式，u1=(B?RV?1B)?1B?RV?1?+(u0?(B?RV?1B)?1B?RV?1Bu0)︸=0，从而一步得到DPG（最小残差）解。需要注意的是，通过为试验空间引入能量范数，试验范数和测试范数本质上就处于对偶配对关系。如果试验范数和测试范数不处于对偶配对关系，其收敛性可由以下结果[4]描述，此处不再证明。定理2：设γ和M分别表示算子B的inf-sup常数和连续性常数。令C???M/γ。那么，对于任意q∈[(C2?1)/(C2+1),1)且α∈[(1?q)/γ2,(1+q)/M2]，采用步长α的最速下降迭代会得到一个收缩常数≥q的收缩映射。这两种方法各有优缺点。在比较它们之前，首先值得探讨的一个问题是这两种方法是否有可能重合。高斯-牛顿最小化问题(17)对应的正规方程为(dunB)?RV?1dunBΔu=(dunB)?RV?1(??B(un))。这个表达式与最速下降更新公式(24)（采用单位步长时）极为相似。实际上，当RU=(dunB)?RV?1dunB时，两种更新方法就会重合，而这种情况恰好发生在试验空间采用与步长相关的能量范数∥Δu∥E:=∥duBΔu∥V′时。备注3：上述观察表明，至少从理论上讲，最速下降法有可能与高斯-牛顿法重合，也就是说，此时会出现二次收敛阶。不过，这在实际问题中是否真的会发生，则是一个更为复杂的问题。高斯-牛顿法的核心优势在于其在最优解附近具有二次收敛阶，而最速下降法仅具有线性收敛阶。尽管如此，高斯-牛顿法隐含地假设线性化算子满足inf-sup条件，这一点可以从公式(18)左侧出现的duB看出。虽然这一假设看似只是一个小小的数学调整，但它对物理应用有着重大影响。某些物理现象（如屈曲）会导致解的不唯一性，进而可能使线性化刚度矩阵出现病态。相比之下，最速下降法中线性化算子仅出现在载荷向量的右侧，此时的刚度矩阵则是条件良好的Gram矩阵。尽管存在一定的计算效率问题，但使用梯度方法的优势依然显著，因为它能够直接施加约束，包括那些仅具一次可微性的约束（例如在极端应变条件下可能重要的非自穿透条件）。或许最重要的是，由于存在Riesz算子RU和RV，最速下降法要求同时指定测试范数和试验范数。从数学角度来看，这一要求似乎并不重要，但实际上却有着重大意义，尤其是在试验范数难以计算的情况下。关于实现细节的进一步讨论将放在后面的章节中。模型问题最小残差法、最速下降法以及高斯-牛顿法已被用于解决非线性弹性领域的模型问题[1]。一般来说，非线性弹性的数学研究旨在寻找满足如下形式的积分的最小化函数u:Ω→RN：E(u)=∫ΩW(x,F)+?(x,u)dx?∫Γtguds，其中W是储能函数，?表示体力势，g表示所有施加的拉力，而F:=?u+I则被称为变形梯度。在本节中，我们将详细介绍最速下降法的实现及相关理论和实践细节。处理试验范数最速下降法要求必须指定试验范数，这一要求适用于所有的未知量，尤其是包括迹未知量在内的那些量，这就带来了第一个挑战。这些变量的范数是一种最小能量扩展范数，虽然其定义已知，但计算起来却相当困难。幸运的是，通过利用变量之间的线性依赖关系，可以简化计算过程。结论在本报告中，我们介绍了最速下降法，并探讨了如何将其用于通过DPG方法解决一类二维非线性问题的最小残差问题。我们推导出了微分和梯度的表达式，并且在试验范数和测试范数处于对偶配对关系时，给出了该方法的某些收敛性结果。与通过最小化线性化残差的高斯-牛顿法相比，最速下降法具有诸多优势，比如能带来更好的条件数。致谢作者贡献说明Jonathan Zhang：负责初稿撰写。Leszek Demkowicz：负责审稿与编辑工作。利益冲突声明作者声明存在以下可能被视为潜在利益冲突的财务利益/个人关系：Jonathan Zhang表示获得了美国国家科学基金会的资金支持；Leszek Demkowicz也表示获得了美国国家科学基金会的资金支持。如果有其他作者，他们则声明自己没有已知的可能影响本文所述工作的财务利益或个人关系。