郭孟武
, 钟宏志
清华大学 土木工程系, 北京 100084
收稿日期：2015-12-03
基金项目：国家自然科学基金资助项目（51378294）
作者简介：郭孟武 (1991—), 男, 博士研究生
通信作者：钟宏志, 教授, E-mail：
hzz@tsinghua.edu.cn 摘要：在模型检验研究中，针对离散误差的后验误差估计扮演着重要角色。在工程有限元分析中，有关问题解答场的标量性质的目标输出量是后验误差分析中人们所关注的。在已有的面向目标误差估计技术中，有2种方法能够提供目标量误差的严格上下界，即本构关系误差法与凸目标函数优化法。该文简要总结了这2种方法，并从计算列式的一致性和基本原理的等价性2个层面论证了2种方法的等价性，给出了2种严格界方法的本质均为余能原理的论断。对2种方法等价性的探讨有助于结合2种表达格式的优势，并拓展到更复杂的问题中，形成简明有效的计算列式。
关键词：微分方程的数值解法后验误差估计面向目标误差估计本构关系误差凸目标函数约束优化严格界余能原理
Equivalence of two strictly bounding approaches for goal-oriented error estimation
GUO Mengwu
, ZHONG Hongzhi
Department of Civil Engineering, Tsinghua University, Beijing 100084, China
Abstract: For model verification which is mainly focused on the control of discretization errors of numerical results, a posteriori error estimation plays an important role in various numerical tools such as the finite element method. The outputs of interest are usually converted into integral functionals over the problem domain for a posteriori error analysis. Among the available techniques for goal-oriented error estimation, only two approaches, the constitutive relation error estimation and the constrained optimization method with convex objective functions, have been claimed to be able to offer guaranteed strict upper and lower bounds for the errors in quantities of interest. The two approaches are briefly reviewed and the equivalence of their formulation and principle is given. Both approaches are shown to be essentially based on the complementary energy theorem. The equivalence of the two approaches is instrumental in simplification of error estimation and extension of applications into more complex problems.
Key words: numerical methods for partial differential equations (PDEs)a posteriori error estimationgoal-oriented error estimationconstitutive relation errorconstrained optimization with convex objective functionsstrict boundscomplementary energy theorem
有限元法作为一种十分有效的数值方法，在工程设计与计算分析中得以广泛应用。为了控制误差以保证计算模拟的质量，模型检验 (model verification) 作为数值方法中的重要步骤已经得到广泛研究。而在各种计算误差的来源中，离散误差具有主导地位。为了评价有限元的离散误差，一系列后验误差 (a posteriori error) 估计
[1-3]技术得以发展，例如：显式型 (残值型) 误差估计方法
[4]、隐式型误差估计方法
[5-6]、恢复型误差估计方法
[7]、阶谱误差估计方法
[8]、本构关系误差 (constitutive relation error) 方法
[9]等；这些方法均可对有限元分析带来的整体离散误差进行评估。
然而在实际的有限元分析中，人们关注的往往是一些有关问题解的目标量，例如平均位移、局部应力等。因此，在后验误差分析时，对目标量的离散误差评价是非常必要的，这种针对目标输出量进行的误差估计称为面向目标误差估计 (goal-oriented error estimation)。面向目标误差估计基于伴随方法
[10-11]，自20世纪90年代中期开始发展
[12-13]以来，各种估计技术相继被提出，主要包括伴随加权残值法
[10, 14-15]、能量模估计方法
[16]、Green函数分解法
[17]、本构关系误差法、凸目标函数优化 (constrained optimization with convex objective functions) 法等。与此同时，面向目标误差估计也在Poisson方程、固体力学问题、特征值问题、时变问题、流体力学问题等诸多问题的数值求解中得到了应用
[15, 18]。
在已有的面向目标误差估计技术中，有2种方法可以提供目标量误差的严格上下界，分别为Ladevèze教授提出的本构关系误差方法
[19-24]及Patera等提出的凸目标函数约束优化法
[25-27]。这2种方法为目标量提供可计算的严格误差界，这种严格界性质及在各种问题中的适用性使得这2种方法在各种面向目标误差估计技术中表现突出。
本文旨在理论上论证这2种方法的等价性。在对2种方法的基本思路进行综述的基础上，从计算列式上说明了2种方法的一致性，并进一步揭示二者的基本原理均为线弹性系统的余能最小值原理，故具有原理上的等价性。
1 模型问题考虑定义在域Ω?
${{\mathbb{R}}^{3}}$上的线弹性体，其Dirichlet边界为
?1Ω，Neumann边界为
?2Ω=
?Ω\
?1Ω。该结构上作用有体力
fp∈[
L2(
Ω)]
3，另在
?1Ω上给定位移
up，在
?2Ω上给定面力
Tp∈[
L2(
?2Ω)]
3。
在本文中，采用如下集合记号：
| $\left\{ \begin{array}{l}U = \left\{ {\mathit{\boldsymbol{u}} \in {{\left[ {{H^1}\left( \mathit{\Omega } \right)} \right]}^3}:\mathit{\boldsymbol{u}}\left| {_{{\partial _1}\mathit{\Omega }}} \right. = {\mathit{\boldsymbol{u}}_p}} \right\},\\V = \left\{ {\mathit{\boldsymbol{v}} \in {{\left[ {{H^1}\left( \mathit{\Omega } \right)} \right]}^3}:\mathit{\boldsymbol{v}}\left| {_{{\partial _1}\mathit{\Omega }}} \right. = {\bf{0}}} \right\},\\S = \left\{ {\mathit{\boldsymbol{\sigma }} \in {{\left[ {{L_2}\left( \mathit{\Omega } \right)} \right]}^{3 \times 3}}:\mathit{\boldsymbol{\sigma = }}{\mathit{\boldsymbol{\sigma }}^{\rm{T}}}} \right\}.\end{array} \right.$ | (1) |
并定义退化能量内积
au(·, ·):[
H1(
Ω)]
3×[H
1(Ω)]
3→
$\mathbb{R}$及相应的半模‖·‖
u，能量内积
aσ(·, ·):
S×
S→
$\mathbb{R}$及相应的模‖·‖
σ以及有界线性泛函
l(·):[
H1(
Ω)]
3→
$\mathbb{R}$如下:
| $\left\{ \begin{array}{l}\begin{array}{*{20}{c}}{{a_u}\left( {\mathit{\boldsymbol{u}},\mathit{\boldsymbol{v}}} \right) = \int_\mathit{\Omega } {\mathit{\boldsymbol{\varepsilon }}\left( \mathit{\boldsymbol{u}} \right):\mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( \mathit{\boldsymbol{v}} \right)} {\rm{d}}\mathit{\Omega },}\\{{{\left\| \mathit{\boldsymbol{u}} \right\|}_u} = \sqrt {{a_u}\left( {\mathit{\boldsymbol{u}},\mathit{\boldsymbol{u}}} \right)} ,}\end{array}\\{a_\sigma }\left( {\mathit{\boldsymbol{\sigma }},\mathit{\boldsymbol{\tau }}} \right) = \int_\mathit{\Omega } {\mathit{\boldsymbol{\sigma }}:{\mathit{\boldsymbol{H}}^{ - 1}}:\mathit{\boldsymbol{\tau }}} {\rm{d}}\mathit{\Omega ,}\\\begin{array}{*{20}{c}}{{{\left\| \mathit{\boldsymbol{\sigma }} \right\|}_\sigma } = \sqrt {{a_\sigma }\left( {\mathit{\boldsymbol{\sigma }},\mathit{\boldsymbol{\sigma }}} \right)} ,}\\{l\left( \mathit{\boldsymbol{v}} \right) = \int_\mathit{\Omega } {{f_p} \cdot \mathit{\boldsymbol{v}}{\rm{d}}\mathit{\Omega }} + \int_{{\partial _2}\mathit{\Omega }} {{\mathit{\boldsymbol{T}}_p} \cdot \mathit{\boldsymbol{v}}{\rm{d}}\partial \mathit{\Omega }} .}\end{array}\end{array} \right.$ | (2) |
其中：
H为材料的Hooke张量；
ε(
u) 表示Cauchy应变张量，
ε(
u)=(▽
u+
u▽)/2。
那么，该模型问题 (原问题) 可如下表述：求解定义在
Ω上的位移场
u和应力场
σ，使得满足以下3组方程。
协调条件：
| $\mathit{\boldsymbol{u}} \in U;$ | (3) |
平衡条件：
| $\mathit{\boldsymbol{\sigma }} \in S,\;\;\;\;{a_\sigma }\left( {\mathit{\boldsymbol{\sigma }},\mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( \mathit{\boldsymbol{v}} \right)} \right) = l\left( \mathit{\boldsymbol{v}} \right),\;\;\;\;\forall \mathit{\boldsymbol{v}} \in V;$ | (4) |
本构方程：
| $\mathit{\boldsymbol{\sigma }} = \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( \mathit{\boldsymbol{u}} \right).$ | (5) |
进一步地，综合式 (3)—(5)，该边值问题的弱形式描述为求解
u∈
U使得
| ${a_u} = \left( {\mathit{\boldsymbol{u}},\mathit{\boldsymbol{v}}} \right) = l\left( \mathit{\boldsymbol{v}} \right),\;\;\;\forall \mathit{\boldsymbol{v}} \in V.$ | (6) |
若采用Galerkin有限元对该问题进行离散，并采用插值空间
Ph，则该问题的有限元格式为求解
uh∈
Uh使得
| ${a_u} = \left( {{\mathit{\boldsymbol{u}}_h},\mathit{\boldsymbol{v}}} \right) = l\left( v \right),\;\;\;\forall \mathit{\boldsymbol{v}} \in {V_h}.$ | (7) |
其中：
Uh=
U∩
Ph，
Vh=
V∩
Ph。而相应的有限元应力解由本构关系式 (5) 得到，即
σh=
H:
ε(
uh)。
基于该模型问题进行面向目标分析，不失一般性地，将目标量
I定义为关于位移场
u的有界线性泛函
lO:[
H1(Ω)]
3→
$\mathbb{R}$，即
| $I = {l^O}\left( \mathit{\boldsymbol{u}} \right).$ | (8) |
而其计算值
Ih=
lO(
uh) 的误差 (
I－
Ih) 是本文关注的重点。下述2种方法均可给出 (
I－
Ih) 的严格上下界，以作为面向该目标量的离散误差估计。
2 本构关系误差方法本构关系误差建立在容许场的概念上，通过本构关系定义线弹性问题的误差；通过有限元解选定容许场后，本构关系误差可为位移误差场的能量模提供严格上界；基于此原理，在面向目标误差估计中，由本构关系误差可得到目标量误差的严格上下界。
2.1 本构关系误差的概念本构关系误差的概念基于模型问题的容许场
$\left( {\mathit{\boldsymbol{\tilde u}},\mathit{\boldsymbol{\tilde \sigma }}} \right)$，其中
$\mathit{\boldsymbol{\tilde u}}$为满足协调条件式 (3) 的位移场，称为协调位移场，
${\mathit{\boldsymbol{\tilde \sigma }}}$为满足平衡方程式 (4) 的应力场，称为平衡应力场。而本构关系误差则通过本构关系式 (5) 进行度量：
| ${e_{{\rm{CRE}}}}\left( {\mathit{\boldsymbol{\tilde u}},\mathit{\boldsymbol{\tilde \sigma }}} \right) = {\left\| {\mathit{\boldsymbol{\tilde \sigma }} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {\mathit{\boldsymbol{\tilde u}}} \right)} \right\|_\sigma }.$ | (9) |
可以证明本构关系误差具有如下性质：
| ${e_{{\rm{CRE}}}}\left( {\mathit{\boldsymbol{\tilde u}},\mathit{\boldsymbol{\tilde \sigma }}} \right) = 0 \Leftrightarrow \left( {\mathit{\boldsymbol{\tilde u}},\mathit{\boldsymbol{\tilde \sigma }}} \right) = \left( {\mathit{\boldsymbol{u}},\mathit{\boldsymbol{\sigma }}} \right)\;\;\;\;{\rm{a}}{\rm{.e}}{\rm{.,}}$ | (10a) |
| $e_{{\rm{CRE}}}^2\left( {\mathit{\boldsymbol{\tilde u}},\mathit{\boldsymbol{\tilde \sigma }}} \right) = \left\| {\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}}} \right\|_u^2 + \left\| {\mathit{\boldsymbol{\sigma }} - \mathit{\boldsymbol{\tilde \sigma }}} \right\|_\sigma ^2,$ | (10b) |
| $\begin{array}{l}{e_{{\rm{CRE}}}}\left( {\mathit{\boldsymbol{\tilde u}},\mathit{\boldsymbol{\tilde \sigma }}} \right) = 2{\left\| {\mathit{\boldsymbol{\sigma }} - {{\mathit{\boldsymbol{\tilde \sigma }}}^ * }} \right\|_\sigma },\\\;\;\;\;\;\;\;\;{{\mathit{\boldsymbol{\tilde \sigma }}}^ * }: = \left( {\mathit{\boldsymbol{\tilde \sigma }} + \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {\mathit{\boldsymbol{\tilde u}}} \right)} \right)/2.\end{array}$ | (10c) |
在有限元分析中，由于
uh∈
Uh?
U，往往取
${\mathit{\boldsymbol{\tilde \mu }}}$=
uh；然而，由于
σh不满足平衡方程式 (4)，故平衡应力场
$\mathit{\boldsymbol{\tilde \sigma }} = {{\mathit{\boldsymbol{\tilde \sigma }}}_h}$需要额外构造。该应力场
${{\mathit{\boldsymbol{\tilde \sigma }}}_h}$可基于有限元应力解
σh通过一种称为延长条件的能量关系
[1, 28]进行构造，即
| $a_\sigma ^E\left( {{{\mathit{\boldsymbol{\tilde \sigma }}}_h} - {\mathit{\boldsymbol{\sigma }}_h},\mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( \mathit{\boldsymbol{v}} \right)} \right) = 0,\;\;\;\forall \mathit{\boldsymbol{v}} \in {\left[ {{H^1}\left( E \right)} \right]^3} \cap {P_h}.$ | (11) |
其中：
E为有限元分析时划分出的每一个单元，而
aσE表示
aσ在单元
E上的限制。
由式 (10b) 易见位移误差场
e:=
u－
uh能量模的严格上界，即
| ${\left\| \mathit{\boldsymbol{e}} \right\|_u} \le {e_{{\rm{CRE}}}}\left( {{\mathit{\boldsymbol{u}}_h},{{\mathit{\boldsymbol{\tilde \sigma }}}_h}} \right).$ | (12) |
2.2 针对目标量的伴随问题如式 (8) 所示，模型问题解的标量目标量
I作为[
H1(
Ω)]
3上的有界线性泛函，常具有如下的形式：
| ${l^O}\left( \mathit{\boldsymbol{u}} \right) = \int_\mathit{\Omega } {{f_d} \cdot \mathit{\boldsymbol{u}}{\rm{d}}\mathit{\Omega }} + \int_{{\partial _2}\mathit{\Omega }} {{\mathit{\boldsymbol{T}}_d} \cdot \mathit{\boldsymbol{u}}{\rm{d}}\partial \mathit{\Omega }} .$ | (13) |
其中：
fd∈[
L2(
Ω)]
3，
Td∈[
L2(
?2Ω)]
3为提取因子。
针对该目标量引入伴随问题：求解定义在
Ω上的位移场
w和应力场
τ，使其满足
协调条件：
| $\mathit{\boldsymbol{w}} \in V;$ | (14) |
平衡条件：
| $\mathit{\boldsymbol{\tau }} \in S,\;\;\;{a_\sigma }\left( {\mathit{\boldsymbol{\tau }},\mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( \mathit{\boldsymbol{v}} \right)} \right) = {l^O}\left( \mathit{\boldsymbol{v}} \right)\;\;\;\forall \mathit{\boldsymbol{v}} \in V;$ | (15) |
本构方程：
| $\mathit{\boldsymbol{\tau }} = \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( \mathit{\boldsymbol{w}} \right).$ | (16) |
该边值问题的弱形式描述为求解
w∈
V使得:
| ${a_u}\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{v}}} \right) = {l^O}\left( \mathit{\boldsymbol{v}} \right)\;\;\;\;\forall \mathit{\boldsymbol{v}} \in V.$ | (17) |
该问题的有限元格式为求解
wh∈
Vh使得:
| ${a_u}\left( {{\mathit{\boldsymbol{w}}_h},\mathit{\boldsymbol{v}}} \right) = {l^O}\left( \mathit{\boldsymbol{v}} \right)\;\;\;\;\forall \mathit{\boldsymbol{v}} \in {V_h}.$ | (18) |
与模型问题 (原问题) 类似，进一步可得有限元应力解
τh，并通过延长条件构造平衡应力场
${{\mathit{\boldsymbol{\tilde \tau }}}_h}$。
2.3 目标量误差的表示与估计由式 (17)，考虑到
e∈
V，并结合Galerkin正交性，即对?
v∈
Vh，
a(
e,
v)=0，可得目标量误差
I－
Ih=
lO(
e) 的表达式为
| $I - {I_h} = {a_u }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{\varepsilon }}} \right).$ | (19) |
其中,
ε:=
w－
wh为伴随问题的位移误差场, 由式 (12) 并借助Cauchy不等式，该误差的严格上下界可由下式给出：
| $\left| {I - {I_h}} \right| \le {e_{{\rm{CRE}}}}\left( {{\mathit{\boldsymbol{u}}_h},{{\mathit{\boldsymbol{\tilde \sigma }}}_h}} \right) \cdot {e_{{\rm{CRE}}}}\left( {{\mathit{\boldsymbol{w}}_h},{{\mathit{\boldsymbol{\tilde \tau }}}_h}} \right).$ | (20) |
除此之外，借由本构关系误差的性质式 (10c) 可进一步推证：
| $\begin{array}{*{20}{c}}{\left| {I - {I_h} - \frac{1}{2}{a_\sigma }\left( {{{\mathit{\boldsymbol{\tilde \sigma }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{u}}_h}} \right),{{\mathit{\boldsymbol{\tilde \tau }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{w}}_h}} \right)} \right)} \right| \le }\\{\frac{1}{2}{e_{{\rm{CRE}}}}\left( {{\mathit{\boldsymbol{u}}_h},{{\mathit{\boldsymbol{\tilde \sigma }}}_h}} \right) \cdot {e_{{\rm{CRE}}}}\left( {{\mathit{\boldsymbol{w}}_h},{{\mathit{\boldsymbol{\tilde \tau }}}_h}} \right).}\end{array}$ | (21) |
则目标量
I的误差的严格、可计算上下界可表示为
| $\left\{ \begin{array}{l}\begin{array}{*{20}{c}}{{I^{{\rm{upper}}}} - {I_h} = \frac{1}{2}{a_\sigma }\left( {{{\mathit{\boldsymbol{\tilde \sigma }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{u}}_h}} \right),{{\mathit{\boldsymbol{\tilde \tau }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{w}}_h}} \right)} \right) + }\\{\frac{1}{2}{e_{{\rm{CRE}}}}\left( {{\mathit{\boldsymbol{u}}_h},{{\mathit{\boldsymbol{\tilde \sigma }}}_h}} \right) \cdot {e_{{\rm{CRE}}}}\left( {{\mathit{\boldsymbol{w}}_h},{{\mathit{\boldsymbol{\tilde \tau }}}_h}} \right) \cdot }\end{array}\\\begin{array}{*{20}{c}}{{I^{{\rm{lower}}}} - {I_h} = \frac{1}{2}{a_\sigma }\left( {{{\mathit{\boldsymbol{\tilde \sigma }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{u}}_h}} \right),{{\mathit{\boldsymbol{\tilde \tau }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{w}}_h}} \right)} \right) - }\\{\frac{1}{2}{e_{{\rm{CRE}}}}\left( {{\mathit{\boldsymbol{u}}_h},{{\mathit{\boldsymbol{\tilde \sigma }}}_h}} \right) \cdot {e_{{\rm{CRE}}}}\left( {{\mathit{\boldsymbol{w}}_h},{{\mathit{\boldsymbol{\tilde \tau }}}_h}} \right) \cdot }\end{array}\end{array} \right.$ | (22) |
3 凸目标函数约束优化方法凸目标函数约束优化方法将所求误差 (
I－
Ih) 表示为凸目标函数，并转化为约束优化问题；再通过Lagrange形式将其转化为无约束优化问题，由此得到误差的严格界。
3.1 Lagrange形式与目标量误差上下界由式 (6) 和 (7)，模型问题精确解
u的未知使得目标量的误差
I－
Ih=
lO(
u－
uh) 可表示为以下约束优化形式：
| $\begin{array}{*{20}{c}}{I - {I_h} = \mathop {\min }\limits_{e \in V} \left[ {{l^O}\left( \mathit{\boldsymbol{e}} \right) + \kappa \left( {{a_u}\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{e}}} \right) - R\left( \mathit{\boldsymbol{e}} \right)} \right)} \right],}\\{{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;{a_u}\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}} \right) = R\left( \mathit{\boldsymbol{w}} \right)\;\;\;\;\;\forall \mathit{\boldsymbol{w}} \in V.}\end{array}$ | (23) |
其中:
κ∈
${{\mathbb{R}}^{+}}$;
R(·):=
l(·)－
au(
uh, ·) 为定义在[
H1(
Ω)]
3上的有界线性泛函，称为残值泛函。且值得注意的是，目标函数
lO(
e)+
κ(
au(
e,
e)－
l(
e)) 具有凸性。
当采用插值空间
Ph进行有限元分析时，将单元划分的具体形式记作
Γh，而将所有内部的单元交界及单元与边界的交界组成的集合记作
Ξ(
Γh)，则函数空间V亦可表示为
| $\begin{array}{*{20}{c}}{V = \left\{ {\mathit{\boldsymbol{v}} \in \hat V:0 = b\left( {\mathit{\boldsymbol{v}},\mathit{\boldsymbol{q}}} \right): = } \right.}\\{\left. {\sum\limits_{{\gamma _h} \in \Xi \left( {{\mathit{\Gamma }_h}} \right)} {\int_{{\gamma _h}} {{{\left[ \mathit{\boldsymbol{v}} \right]}_{{\gamma _h}}} \cdot \mathit{\boldsymbol{q}}\left| {_{{\gamma _h}}{\rm{d}}s} \right.\;\;\;\forall \mathit{\boldsymbol{q}} \in Q} } } \right\}.}\end{array}$ | (24) |
其中：
${\hat V}$为破坏
V在
Ξ(
Γh) 上所有连续性 (包括Dirichlet边界条件) 形成的函数空间；[
v]
γh当
γh为内部交界时为
v的跳跃，而当
γh在边界上时为
v的迹；
Q的定义为
Q={
q∈
L2(
Ξ(
Γh)):
q|
?2Ω=0}。
因此，(
I－
Ih) 可进一步地表示为如下的约束优化问题：
| $\begin{array}{l}I - {I_h} = \mathop {\min }\limits_{e \in T} \left[ {{l^O}\left( \mathit{\boldsymbol{e}} \right) + \kappa \left( {{a_u}\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{e}}} \right) - R\left( \mathit{\boldsymbol{e}} \right)} \right)} \right],\\T = \left\{ \begin{array}{l}\mathit{\boldsymbol{e}} \in \hat V:{a_u}\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}} \right) = R\left( \mathit{\boldsymbol{w}} \right)\;\;\;\;\;\forall \mathit{\boldsymbol{w}} \in V;\\b\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{q}}} \right) = 0\;\;\;\;\forall \mathit{\boldsymbol{q}} \in Q\end{array} \right\}.\end{array}$ | (25) |
在此基础上，通过引入Lagrange乘子将上述约束优化问题转化为无约束优化问题，相应的Lagrange形式为
| $\begin{array}{*{20}{c}}{{L^ \pm }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}},\mathit{\boldsymbol{q}}} \right): = }\\{ \pm {l^O}\left( \mathit{\boldsymbol{e}} \right) + \kappa \left( {{a_u}\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{e}}} \right) - R\left( \mathit{\boldsymbol{e}} \right)} \right) + R\left( \mathit{\boldsymbol{w}} \right) - }\\{{a_u}\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}} \right) + b\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{q}}} \right).}\end{array}$ | (26) |
则误差 (
I－
Ih) 为
| $ \pm \left( {I - {I_h}} \right) = \mathop {\min }\limits_{\mathit{\boldsymbol{e}} \in \hat V} \mathop {\max }\limits_{\mathit{\boldsymbol{w}} \in V,\mathit{\boldsymbol{q}} \in Q} {L^ \pm }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}},\mathit{\boldsymbol{q}}} \right).$ | (27) |
进一步地，对于取定的
wh±∈
V与
qh±∈
Q，
| $\begin{array}{*{20}{c}}{ \pm \left( {I - {I_h}} \right) \ge \mathop {\max }\limits_{\mathit{\boldsymbol{w}} \in V,\mathit{\boldsymbol{q}} \in Q} \mathop {\min }\limits_{\mathit{\boldsymbol{e}} \in \hat V} {L^ \pm }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}},\mathit{\boldsymbol{q}}} \right) \ge }\\{\mathop {\min }\limits_{\mathit{\boldsymbol{e}} \in \hat V} {L^ \pm }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}_h^ \pm ,\mathit{\boldsymbol{q}}_h^ \pm } \right).}\end{array}$ | (28) |
即
| $\mathop {\min }\limits_{\mathit{\boldsymbol{e}} \in \hat V} {L^ + }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}_h^ + ,\mathit{\boldsymbol{q}}_h^ + } \right) \le I - {I_h} \le - \mathop {\min }\limits_{\mathit{\boldsymbol{e}} \in \hat V} {L^ - }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}_h^ - ,\mathit{\boldsymbol{q}}_h^ - } \right).$ | (29) |
给出了 (
I-
Ih) 的严格上下界。
3.2 误差上下界计算方法式 (28) 中的
wh±与
qh±取在Lagrange形式 (26) 在子空间
${{\hat V}_h} \times {V_h} \times {Q_h} \subset \hat V \times V \times Q$中的鞍点处，
| ${L^ \pm }\left( {\mathit{\boldsymbol{e}}_h^ \pm ,\mathit{\boldsymbol{w}}_h^ \pm ,\mathit{\boldsymbol{q}}_h^ \pm } \right) = \mathop {\min }\limits_{\mathit{\boldsymbol{e}} \in {{\hat V}_h}} \mathop {\max }\limits_{\mathit{\boldsymbol{w}} \in {V_h},\mathit{\boldsymbol{q}} \in {Q_h}} {L^ \pm }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}},\mathit{\boldsymbol{q}}} \right).$ | (30) |
其中:
${{\hat V}_h} = \hat V \cap {P_h}$;
Qh={
q∈Π
h(
γh), ?
γh∈
Ξ(
Γh):
q|
?2Ω=0}; Π
h(·) 表示与
Ph对应的多项式空间。该鞍点问题的驻值条件为
| $\begin{array}{*{20}{c}}{b\left( {\mathit{\boldsymbol{v}},\mathit{\boldsymbol{q}}_h^ \pm } \right) \pm {l^O}\left( \mathit{\boldsymbol{v}} \right) + \kappa \left( {2{a_u}\left( {\mathit{\boldsymbol{e}}_h^ \pm ,\mathit{\boldsymbol{v}}} \right) - R\left( \mathit{\boldsymbol{v}} \right)} \right) - }\\{{a_u}\left( {\mathit{\boldsymbol{v}},\mathit{\boldsymbol{w}}_h^ \pm } \right) = 0\;\;\;\;\forall \mathit{\boldsymbol{v}} \in {{\hat V}_h},}\end{array}$ | (31a) |
| $R\left( \mathit{\boldsymbol{v}} \right) - {a_u}\left( {\mathit{\boldsymbol{e}}_h^ \pm ,\mathit{\boldsymbol{v}}} \right) = 0\;\;\;\;\forall \mathit{\boldsymbol{v}} \in {V_h},$ | (31b) |
| $b\left( {\mathit{\boldsymbol{e}}_h^ \pm ,\mathit{\boldsymbol{p}}} \right) = 0\;\;\;\;\;\forall \mathit{\boldsymbol{p}} \in {Q_h}.$ | (31c) |
式 (31c) 表明
eh±∈
Vh，进一步可由式 (31b) 得到
eh±=
0，则式 (31a) 简化为
| ${a_u}\left( {\mathit{\boldsymbol{v}},\mathit{\boldsymbol{w}}_h^ \pm } \right) = \pm {l^O}\left( \mathit{\boldsymbol{v}} \right) - \kappa R\left( \mathit{\boldsymbol{v}} \right) + b\left( {\mathit{\boldsymbol{v}},\mathit{\boldsymbol{q}}_h^ \pm } \right)\;\;\;\;\;\forall \mathit{\boldsymbol{v}} \in {{\hat V}_h}.$ | (32) |
当取?
v∈
Vh?
${{\hat V}_h}$时，该式退化为
au(
v,
wh)=
lO(
v)，而
wh±=±
wh。将该
wh±带入式 (31a) 并令
qh±=
$ \mp $q0h+
κq1h，可将该式分解为
| $b\left( {\mathit{\boldsymbol{v}},{\mathit{\boldsymbol{q}}_{0h}}} \right) = {l^O}\left( \mathit{\boldsymbol{v}} \right) - {a_u}\left( {\mathit{\boldsymbol{v}},{\mathit{\boldsymbol{w}}_h}} \right)\;\;\;\;\;\forall \mathit{\boldsymbol{v}} \in {{\hat V}_h},$ | (33a) |
| $b\left( {\mathit{\boldsymbol{v}},{\mathit{\boldsymbol{q}}_{1h}}} \right) = l\left( \mathit{\boldsymbol{v}} \right) - {a_u}\left( {{\mathit{\boldsymbol{u}}_h},\mathit{\boldsymbol{v}}} \right)\;\;\;\;\;\forall \mathit{\boldsymbol{v}} \in {{\hat V}_h}.$ | (33b) |
由此，考虑式 (29) 中的优化问题，即求解
$\mathit{\boldsymbol{e}} = \mathit{\boldsymbol{\tilde e}}_h^ \pm \in \hat V$使得
L±(
e,
wh±,
qh±) 取最小值，相应的驻值条件为
| $\begin{array}{*{20}{c}}{ \pm {l^O}\left( \mathit{\boldsymbol{v}} \right) + \kappa \left( {2{a_u}\left( {\mathit{\boldsymbol{\tilde e}}_h^ \pm ,\mathit{\boldsymbol{v}}} \right) - R\left( \mathit{\boldsymbol{v}} \right)} \right) - }\\{{a_u}\left( {\mathit{\boldsymbol{v}},\mathit{\boldsymbol{w}}_h^ \pm } \right) + b\left( {\mathit{\boldsymbol{v}},\mathit{\boldsymbol{q}}_h^ \pm } \right) = 0\;\;\;\;\forall \mathit{\boldsymbol{v}} \in \hat V.}\end{array}$ | (34) |
令
$\mathit{\boldsymbol{\tilde e}}_h^ \pm = \mp {{\mathit{\boldsymbol{\tilde e}}}_{0h}}/\kappa + {{\mathit{\boldsymbol{\tilde e}}}_{1h}}$，则上式可分解为
| $\begin{array}{*{20}{c}}{2{a_u}\left( {{{\mathit{\boldsymbol{\tilde e}}}_{0h}},\mathit{\boldsymbol{v}}} \right) = {l^O}\left( \mathit{\boldsymbol{v}} \right) - b\left( {\mathit{\boldsymbol{v}},{\mathit{\boldsymbol{q}}_{0h}}} \right) - }\\{{a_u}\left( {\mathit{\boldsymbol{v}},{\mathit{\boldsymbol{w}}_h}} \right)\;\;\;\;\;\forall \mathit{\boldsymbol{v}} \in \hat V,}\end{array}$ | (35a) |
| $\begin{array}{*{20}{c}}{2{a_u}\left( {{{\mathit{\boldsymbol{\tilde e}}}_{1h}},\mathit{\boldsymbol{v}}} \right) = l\left( \mathit{\boldsymbol{v}} \right) - b\left( {\mathit{\boldsymbol{v}},{\mathit{\boldsymbol{q}}_{1h}}} \right) - }\\{{a_u}\left( {{\mathit{\boldsymbol{u}}_h},\mathit{\boldsymbol{v}}} \right)\;\;\;\;\;\forall \mathit{\boldsymbol{v}} \in \hat V.}\end{array}$ | (35b) |
该式为确定
${{\mathit{\boldsymbol{\tilde e}}}_{0h}}$与
${{\mathit{\boldsymbol{\tilde e}}}_{1h}}$的控制方程。
3.3 可计算的误差上下界表达考虑式 (34) 与
R(
wh±)=0，可得到
L±(
e,
wh±,
qh±) 的最小值为
| $\begin{array}{*{20}{c}}{\mathop {\min }\limits_{\mathit{\boldsymbol{e}} \in \hat V} {L^ + }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}_h^ + ,\mathit{\boldsymbol{q}}_h^ + } \right) = {L^ \pm }\left( {\mathit{\boldsymbol{\tilde e}}_h^ \pm ,\mathit{\boldsymbol{w}}_h^ \pm ,\mathit{\boldsymbol{q}}_h^ \pm } \right) = }\\{ - \kappa \left\| {\mathit{\boldsymbol{\tilde e}}_h^ \pm } \right\|_u^2.}\end{array}$ | (36) |
进一步代入
$\mathit{\boldsymbol{\tilde e}}_h^ \pm = m{{\mathit{\boldsymbol{\tilde e}}}_{0h}}/\kappa + {{\mathit{\boldsymbol{\tilde e}}}_{1h}}$可将式 (29) 重新表达为
| $\begin{array}{*{20}{c}}{\left| {I - {I_h} - 2{a_u}\left( {{{\mathit{\boldsymbol{\tilde e}}}_{0h}},{{\mathit{\boldsymbol{\tilde e}}}_{1h}}} \right)} \right| \le }\\{\kappa \left\| {{{\mathit{\boldsymbol{\tilde e}}}_{1h}}} \right\|_u^2 + \frac{1}{\kappa }\left\| {{{\mathit{\boldsymbol{\tilde e}}}_{0h}}} \right\|_u^2.}\end{array}$ | (37) |
而当
$\kappa = {\left\| {{{\mathit{\boldsymbol{\tilde e}}}_{0h}}} \right\|_u}/{\left\| {{{\mathit{\boldsymbol{\tilde e}}}_{1h}}} \right\|_u}$时该不等式右侧取最小值，故可得到 (I－
Ih) 的严格、可计算上下界如下：
| $\left\{ \begin{array}{l}{I^{{\rm{upper}}}} - {I_h} = 2{a_u}\left( {{{\mathit{\boldsymbol{\tilde e}}}_{0h}},{{\mathit{\boldsymbol{\tilde e}}}_{1h}}} \right) + 2{\left\| {{{\mathit{\boldsymbol{\tilde e}}}_{1h}}} \right\|_u}{\left\| {{{\mathit{\boldsymbol{\tilde e}}}_{0h}}} \right\|_u},\\{I^{{\rm{lower}}}} - {I_h} = 2{a_u}\left( {{{\mathit{\boldsymbol{\tilde e}}}_{0h}},{{\mathit{\boldsymbol{\tilde e}}}_{1h}}} \right) - 2{\left\| {{{\mathit{\tilde e}}_{1h}}} \right\|_u}{\left\| {{{\mathit{\boldsymbol{\tilde e}}}_{0h}}} \right\|_u}.\end{array} \right.$ | (38) |
4 两种方法的等价性如前文所述，在本构误差关系法 (后文简称为本构法) 与凸目标函数优化法 (后文简称为优化法) 的求解步骤中，均采用有限元法求解原问题 (见式 (7)) 和完全一致的伴随问题 (见式 (18) 与式 (32))，得到原问题的近似解
uh与伴随问题的近似解
wh。不难看出，通过进一步论证本构法中对平衡应力场
${{\mathit{\boldsymbol{\tilde \sigma }}}_h}$与
${{\mathit{\boldsymbol{\tilde \tau }}}_h}$的构造同优化法中
${{\mathit{\boldsymbol{\tilde e}}}_{0h}}$与
${{\mathit{\boldsymbol{\tilde e}}}_{1h}}$的求解具有一致性，即可说明2种方法计算目标量严格误差界的列式的一致性。此外，这种列式上的一致性的背后是基本原理上的等价性。本节通过讨论这些内容来说明2种方法的等价性。
4.1 计算列式的一致性在采用本构法进行面向目标误差估计时，对任意单元E，平衡应力场
${{\mathit{\boldsymbol{\tilde \sigma }}}_h}$的构造所遵循的原则包括 (以处理原问题为例)：
1) 单元内部满足平衡方程，且单元边界与Neumann边界重合时满足Neumann边界条件，即
| $\left\{ \begin{array}{l}\nabla \cdot {{\mathit{\boldsymbol{\tilde \sigma }}}_h} + {\mathit{\boldsymbol{f}}_p} = {\bf{0}}\;\;\;{\rm{in}}\;E,\\\mathit{\boldsymbol{n}} \cdot {{\mathit{\boldsymbol{\tilde \sigma }}}_h} = {\mathit{\boldsymbol{T}}_p}\;\;\;\;{\rm{on}}\;\;\;\partial E \cap {\partial _2}\mathit{\Omega }\mathit{.}\end{array} \right.$ | (39) |
或者表示为弱形式，即
| $\begin{array}{*{20}{c}}{a_\sigma ^E\left( {{{\mathit{\boldsymbol{\tilde \sigma }}}_h},\mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( \mathit{\boldsymbol{v}} \right)} \right) = }\\{{l^E}\left( \mathit{\boldsymbol{v}} \right) - \int_{\partial E} {{\mathit{\boldsymbol{p}}_h} \cdot \mathit{\boldsymbol{v}}{\rm{d}}\mathit{\Omega }} \;\;\;\;\forall \mathit{\boldsymbol{v}} \in {{\left[ {{H^1}\left( E \right)} \right]}^3}.}\end{array}$ | (40a) |
其中
l(
v) 在单元
E上的限制
lE(
v) 与单元内部边界面力
ph定义如下：
| $\begin{array}{*{20}{c}}{{l^E}\left( \mathit{\boldsymbol{v}} \right) = \int_E {{\mathit{\boldsymbol{f}}_p} \cdot \mathit{\boldsymbol{v}}{\rm{d}}\mathit{\Omega }} + \int_{\partial E \cap {\partial _2}\mathit{\Omega }} {{\mathit{\boldsymbol{T}}_p} \cdot \mathit{\boldsymbol{v}}{\rm{d}}\partial \mathit{\Omega }} ,}\\{{\mathit{\boldsymbol{p}}_h} = \left\{ \begin{array}{l} - \mathit{\boldsymbol{n}} \cdot {{\mathit{\boldsymbol{\tilde \sigma }}}_h},\;\;{\rm{on}}\;\;\;\partial E\backslash \left( {\partial E \cap {\partial _2}\mathit{\Omega }} \right),\\0,\;\;\;\;\;\;\;\;\;\;{\rm{on}}\;\;\partial E \cap {\partial _2}\mathit{\Omega }\mathit{.}\end{array} \right.}\end{array}$ | (40b) |
2) 单元交界处面力平衡，即对?
γh∈
Ξ(
Γh)\(
Ξ(
Γh)∩
?2Ω)，
γh=
?E∩
?E′，满足
| ${\mathit{\boldsymbol{p}}_h}\left| {_{{\gamma _h} \cap \partial E}} \right. + {\mathit{\boldsymbol{p}}_h}\left| {_{{\gamma _h} \cap \partial E'}} \right. = 0.$ | (41) |
3) 单元内满足延长条件式 (11)，若结合式 (40a)，延长条件还可以表达为如下形式：
| $\begin{array}{*{20}{c}}{\int_{\partial E} {{\mathit{\boldsymbol{p}}_h} \cdot {\mathit{\boldsymbol{\varphi }}_i}{\rm{d}}\mathit{\Omega }} = {l^E}\left( \mathit{\boldsymbol{v}} \right) - a_{^\sigma }^E\left( {{\mathit{\boldsymbol{\sigma }}_h},\mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{\varphi }}_i}} \right)} \right),}\\{i = 1,2, \cdots .}\end{array}$ | (42) |
其中
φi为描述
H1(
E)∩
Ph的每一个形函数；此外，为了保证该平衡应力场构造的可解性，还要求
ph|
γh∩?E∈Π
h(
γh)。
在本构方法的具体实施中，满足以上3个条件即可构造出平衡应力场
${{\mathit{\boldsymbol{\tilde \sigma }}}_h}$。不难看出，式 (41)、(42) 与式 (33b) 完全一致；进一步对比可见，式 (40a) 与式 (35b) 一致。考虑到伴随问题中平衡应力场的构造与原问题同理，2种方法中相应量的等价关系如下：
| $\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{q}}_{1h}}\left| {_{{\gamma _h}}} \right. = {\mathit{\boldsymbol{p}}_h}\left| {_{{\gamma _h} \cap \partial E}} \right.,\;\;\;\;\;{\mathit{\boldsymbol{q}}_{0h}}\left| {_{{\gamma _h}}} \right. = {\mathit{\boldsymbol{s}}_h}\left| {_{{\gamma _h} \cap \partial E}} \right.,}\\{\mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{{\mathit{\boldsymbol{\tilde e}}}_{1h}}} \right) = \frac{{{{\mathit{\boldsymbol{\tilde \sigma }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{u}}_h}} \right)}}{2},}\end{array}$ | (43a) |
| $\mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{{\mathit{\boldsymbol{\tilde e}}}_{0h}}} \right) = \frac{{{{\mathit{\boldsymbol{\tilde \tau }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{w}}_h}} \right)}}{2}.$ | (43b) |
其中
sh为伴随问题中与
ph定义一致的单元边界面力。
由式 (43) 所示的等价关系可见，由本构法得到的误差上下界式 (22) 与由优化法得到的误差上下界式 (38) 是对应相等的。
值得说明的是，在本构法构造平衡应力场时，延长条件是“空降”的，但其具有很优秀的“平衡性质”。1) 由式 (42) 构造的任意单元E的边界力
ph与该单元上给定外力是平衡的：若此单元发生任意形式的刚体位移
r，该刚体位移可通过完备的形函数线性表示，即
r=
∑riφi；考虑到
ε(
r)=
0，由式 (42) 可得:
| ${l^E}\left( \mathit{\boldsymbol{v}} \right) - \int_{\partial E} {{\mathit{\boldsymbol{p}}_h} \cdot \mathit{\boldsymbol{r}}{\rm{d}}\mathit{\Omega }} = 0.$ | (44) |
即意味着边界力与给定外力平衡；这种性质是通过延长条件式 (42) 从平衡方程式 (40a) 处继承来的。2) 如式 (42) 所示，对于所有公用节点的单元，等号左边相加即公用节点边界上
ph在该点处的等效荷载之和，由于共用边界的单元在该边界处的面力等大反向如式 (41)，此加和应为0，而等号右边相加之和为0则是由有限元列式直接决定的。
然而在优化法中，延长条件式 (33) 是由驻值条件式 (32) 得到，并非“空降”，从优化法的基本思路也可进一步认识上述2种“平衡性质”的根源。对于刚体位移
r，
au(
r,
v)=0，?
v∈[
H1(
E)]
3∩
Ph，则式 (33) 可保证单元边界力与外力的平衡；当考察公用节点边界在该点处的等效荷载之和时，这些公用节点单元之间的连续性 (协调性) 没有被破坏，当式 (33) 等号右端限定在这些单元时，有限元列式保证其为零。
4.2 基本原理的等价性在本节中，从另一视角分别看待2种方法，进一步考察2种方法的一致性，并得到二者在原理上的等价性。
对于本构法，目标量误差式 (19) 可以进一步表示为
| $\begin{array}{*{20}{c}}{I - {I_h} = \frac{1}{4}\left\| {\sqrt \kappa \mathit{\boldsymbol{e + }}\frac{1}{{\sqrt \kappa }}\mathit{\boldsymbol{\varepsilon }}} \right\|_u^2 - }\\{\frac{1}{4}\left\| {\sqrt \kappa \mathit{\boldsymbol{e}} - \frac{1}{{\sqrt \kappa }}\mathit{\boldsymbol{\varepsilon }}} \right\|_u^2.}\end{array}$ | (45) |
其中
$\kappa ={{\mathbb{R}}^{+}}$，则 (
I－
Ih) 的上下界可表示为
| $\begin{array}{*{20}{c}}{ - \frac{1}{4}\left\| {\sqrt \kappa \mathit{\boldsymbol{e}} - \frac{1}{{\sqrt \kappa }}\mathit{\boldsymbol{\varepsilon }}} \right\|_u^2 \le I - {I_h} \le }\\{\frac{1}{4}\left\| {\sqrt \kappa \mathit{\boldsymbol{e + }}\frac{1}{{\sqrt \kappa }}\mathit{\boldsymbol{\varepsilon }}} \right\|_u^2.}\end{array}$ | (46) |
由式 (10b) 可得
| $\begin{array}{*{20}{c}}{\frac{1}{4}\left\| {\sqrt \kappa \mathit{\boldsymbol{e}} \pm \frac{1}{{\sqrt \kappa }}\mathit{\boldsymbol{\varepsilon }}} \right\|_u^2 \le }\\{\frac{1}{4}\left\| {\sqrt \kappa \left( {{{\mathit{\boldsymbol{\tilde \sigma }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{u}}_h}} \right)} \right) \pm \frac{1}{{\sqrt \kappa }}\left( {{{\mathit{\boldsymbol{\tilde \tau }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{w}}_h}} \right)} \right)} \right\|_u^2 = }\\{ \pm \frac{1}{2}{a_\sigma }\left( {{{\mathit{\boldsymbol{\tilde \sigma }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{u}}_h}} \right),{{\mathit{\boldsymbol{\tilde \tau }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{w}}_h}} \right)} \right) + }\\{\frac{\kappa }{4}e_{{\rm{CRE}}}^2\left( {{\mathit{\boldsymbol{u}}_h},{{\mathit{\boldsymbol{\tilde \sigma }}}_h}} \right) + \frac{1}{{4\kappa }}e_{{\rm{CRE}}}^2\left( {{\mathit{\boldsymbol{w}}_h},{{\mathit{\boldsymbol{\tilde \tau }}}_h}} \right).}\end{array}$ | (47) |
可见这种形式已与优化法的式 (37) 一致。当
$\kappa = {e_{{\rm{CRE}}}}\left( {{\mathit{\boldsymbol{w}}_h},{{\mathit{\boldsymbol{\tilde \tau }}}_h}} \right)/{e_{{\rm{CRE}}}}\left( {{\mathit{\boldsymbol{u}}_h},{{\mathit{\boldsymbol{\tilde \sigma }}}_h}} \right)$时，得到如式 (22) 所示的可计算的目标误差上下界。
对于优化法，目标量误差式 (23) 亦可表示为
| $ \pm \left( {I - {I_h}} \right) = \mathop {\min }\limits_{\mathit{\boldsymbol{e}} \in V} \mathop {\max }\limits_{\mathit{\boldsymbol{w}} \in V} L_0^ \pm \left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}} \right).$ | (48) |
其中
| $\begin{array}{*{20}{c}}{L_0^ \pm \left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}} \right): = }\\{ \pm {l^O}\left( \mathit{\boldsymbol{e}} \right) + \kappa \left( {{a_u}\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{e}}} \right) - R\left( \mathit{\boldsymbol{e}} \right)} \right) + R\left( w \right) - {a_u}\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}} \right) = }\\{{L^ \pm }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w,q}}} \right) - b\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{q}}} \right).}\end{array}$ | (49) |
而对于伴随问题的解
wh±=±
wh，
| $ \pm \left( {I - {I_h}} \right) \ge \mathop {\min }\limits_{\mathit{\boldsymbol{e}} \in V} L_0^ \pm \left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}_h^ \pm } \right).$ | (50) |
该式右侧对应的驻值条件为
e0±∈
V，
| $\begin{array}{*{20}{c}}{ \pm {l^O}\left( \mathit{\boldsymbol{v}} \right) + \kappa \left( {2{a_u}\left( {\mathit{\boldsymbol{e}}_0^ \pm ,\mathit{\boldsymbol{v}}} \right) - R\left( \mathit{\boldsymbol{v}} \right)} \right) - {a_u}\left( {\mathit{\boldsymbol{v}},\mathit{\boldsymbol{w}}_h^ \pm } \right) = 0}\\{\forall \mathit{\boldsymbol{v}} \in V.}\end{array}$ | (51) |
该驻值点
e0±可用原问题与伴随问题的位移误差场
e与
ε表示如下：
| $\mathit{\boldsymbol{e}}_0^ \pm = \frac{1}{2}\left( {\mathit{\boldsymbol{e}} \mp \frac{1}{\kappa }\mathit{\boldsymbol{\varepsilon }}} \right),$ | (52) |
代入式 (49) 与 (50) 可得
| $ \pm \left( {I - {I_h}} \right) \ge L_0^ \pm \left( {\mathit{\boldsymbol{e}}_0^ \pm ,\mathit{\boldsymbol{w}}_0^ \pm } \right) = - \frac{1}{4}\left\| {\sqrt \kappa \mathit{\boldsymbol{e}} \mp \frac{1}{{\sqrt \kappa }}\mathit{\boldsymbol{\varepsilon }}} \right\|_u^2.$ | (53) |
可见与本构法的式 (46) 完全相同。然而，由于
e与
ε是不可计算的，故需要将上式所示的界进一步放宽。优化法对此作如下考虑：由于
| $\begin{array}{*{20}{c}}{\mathop {\min }\limits_{e \in V} L_0^ \pm \left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}_h^ \pm } \right) = \mathop {\min }\limits_{e \in V} {L^ \pm }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}_h^ \pm ,\mathit{\boldsymbol{q}}_h^ \pm } \right) \ge }\\{\mathop {\min }\limits_{e \in \hat V} {L^ \pm }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}_h^ \pm ,\mathit{\boldsymbol{q}}_h^ \pm } \right),}\end{array}$ | (54) |
故目标量的严格上下界由
${\mathop {\min }\limits_{e \in \hat V} {L^ \pm }\left( {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{w}}_h^ \pm ,\mathit{\boldsymbol{q}}_h^ \pm } \right)} $给出，如式 (36)—(38) 所示。
从上述讨论可见，本构法与优化法都得到相同的误差界式 (46) 或 (53)，并各自采用不同的思路对该误差界进行放松，得到可计算的形式。其中，本构法引入了平衡应力场与本构关系误差，如式 (47) 所示；优化法则继续放松约束，如式 (54) 所示；这2种思路看似相异，但其原理完全一致，即线弹性体的余能最小值原理。
2种方法对于误差界式 (46) 或 (53) 的放松，都归结于求
${\left\| {\sqrt {\kappa \mathit{\boldsymbol{e}}} \mp \mathit{\boldsymbol{\varepsilon /}}\sqrt \kappa } \right\|_u}$的可计算上界。对于
$\mathit{\boldsymbol{\xi }}_\kappa ^ \mp = \sqrt {\kappa \mathit{\boldsymbol{e}}} \mp \mathit{\boldsymbol{\varepsilon /}}\sqrt \kappa $，可视为如下问题的解：求解
$\mathit{\boldsymbol{\xi }}_\kappa ^ \mp \in V$使得：
| ${a_u}\left( {\mathit{\boldsymbol{\xi }}_\kappa ^ \mp ,\mathit{\boldsymbol{v}}} \right) = \sqrt \kappa R\left( \mathit{\boldsymbol{v}} \right) \mp {R^O}\left( \mathit{\boldsymbol{v}} \right)/\sqrt \kappa .$ | (55) |
其中
RO(·):=
lO(·)－
au(
wh, ·)。若将此问题视为一个线弹性问题，则其对应的平衡条件为
| $\begin{array}{*{20}{c}}{E_\kappa ^ \mp : = \left\{ {\mathit{\boldsymbol{\zeta }}_\kappa ^ \mp \in \zeta :{a_\sigma }\left( {\mathit{\boldsymbol{\zeta }}_\kappa ^ \mp ,\mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( \mathit{\boldsymbol{v}} \right)} \right) = } \right.}\\{\left. {\sqrt \kappa R\left( \mathit{\boldsymbol{v}} \right) \mp {R^O}\left( \mathit{\boldsymbol{v}} \right)/\sqrt \kappa \;\;\;\;\forall \mathit{\boldsymbol{v}} \in V} \right\}.}\end{array}$ | (56) |
考虑到该线弹性体系的内力势能与余能值相等，由余能最小值原理可得：
| $\begin{array}{*{20}{c}}{\frac{1}{2}\left\| {\sqrt \kappa \mathit{\boldsymbol{e}} \mp \mathit{\boldsymbol{\varepsilon }}/\sqrt \kappa } \right\|_u^2 = }\\{\mathop {\min }\limits_{\zeta _\kappa ^ \mp \in E_\kappa ^ \mp } \frac{1}{2}\left\| {\mathit{\boldsymbol{\zeta }}_\kappa ^ \mp } \right\|_\sigma ^2 \le \frac{1}{2}\left\| {\mathit{\boldsymbol{\zeta }}_\kappa ^{ \mp *}} \right\|_\sigma ^2.}\end{array}$ | (57) |
其中
$\mathit{\boldsymbol{\xi }}_\kappa ^{ \mp *} \in E_\kappa ^ \mp $。依据前文对本构法与优化法的讨论，2种方法对
$\mathit{\boldsymbol{\xi }}_\kappa ^{ \mp *}$的取法为
| $\mathit{\boldsymbol{\zeta }}_\kappa ^{ \mp * } = \sqrt \kappa \left( {{{\mathit{\boldsymbol{\tilde \sigma }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{u}}_h}} \right)} \right) \mp \frac{1}{{\sqrt \kappa }}\left( {{{\mathit{\boldsymbol{\tilde \tau }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{w}}_h}} \right)} \right),$ | (58a) |
| $\begin{array}{*{20}{c}}{\mathit{\boldsymbol{\zeta }}_\kappa ^{ \mp * } = 2\mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {\sqrt \kappa {{\mathit{\boldsymbol{\tilde e}}}_{1h}} \mp \frac{1}{{\sqrt \kappa }}{{\mathit{\boldsymbol{\tilde e}}}_{0h}}} \right) = }\\{2\sqrt \kappa \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {\mathit{\boldsymbol{\tilde e}}_h^ \pm } \right).}\end{array}$ | (58b) |
前文已述这2种取法完全相等。回到目标误差估计问题，
| $\begin{array}{*{20}{c}}{ \pm \left( {I - {I_h}} \right) \ge - \frac{1}{4}\left\| {\sqrt \kappa \mathit{\boldsymbol{e}} \mp \frac{1}{{\sqrt \kappa }}\mathit{\boldsymbol{\varepsilon }}} \right\|_u^2 \ge }\\{ - \frac{1}{4}\left\| {\mathit{\boldsymbol{\zeta }}_\kappa ^{ \mp * }} \right\|_\sigma ^2.}\end{array}$ | (59) |
在2种方法中，
| $\begin{array}{*{20}{c}}{ - \frac{1}{4}\left\| {\mathit{\boldsymbol{\zeta }}_\kappa ^{ \mp * }} \right\|_\sigma ^2 = - \frac{1}{4}\left\| {\sqrt \kappa \left( {{{\mathit{\boldsymbol{\tilde \sigma }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{u}}_h}} \right)} \right) \mp } \right.}\\{\left. {\frac{1}{{\sqrt \kappa }}\left( {{{\mathit{\boldsymbol{\tilde \tau }}}_h} - \mathit{\boldsymbol{H}}:\mathit{\boldsymbol{\varepsilon }}\left( {{\mathit{\boldsymbol{w}}_h}} \right)} \right)} \right\|_u^2 = - \kappa \left\| {\mathit{\boldsymbol{\tilde e}}_h^ \pm } \right\|_u^2.}\end{array}$ | (60) |
即式 (47) 和式 (36)。
5 结论本文论证了本构关系误差法与凸目标函数约束优化法在对有界线性泛函目标量进行面向目标误差估计时的计算列式的一致性和基本原理的等价性；经过对比讨论，2种方法基本思想的实质均为余能最小值原理。本构关系误差形式简明、力学意义明确，但表达格式受到本构形式的限制，在拓展到其他问题时处理方案不甚直观；凸目标函数约束优化法表述方式比较复杂，但数学本质更为明确，数学上的一般性更为明显，便于向更复杂的问题进行拓展。对2种方法等价性的讨论有助于2种表达格式的互补，以在更复杂的问题中形成实用有效、简明清晰的表达形式。
参考文献 | [1] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Ladevèze P, Pelle J P. Mastering Calculations in Linear and Nonlinear Mechanics[M]. New York: Springer, 2004. |
| [2] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Ainsworth M, Oden J T. A Posteriori Error Estimation in Finite Element Analysis[M]. New York: John Wiley & Sons, 2000. |
| [3] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Gr?tsch T, Bathe K J. A posteriori error estimation techniques in practical finite element analysis[J]. Computers & Structures, 2005, 83(4): 235–265. |
| [4] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Babu?ka I, Rheinboldt W C. Error estimates for adaptive finite element computations[J]. SIAM Journal on Numerical Analysis, 1983, 20(4): 736–754. |
| [5] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Babu?ka I, Rheinboldt W C. A posterior estimates for the finite element method[J]. International Journal for Numerical Methods in Engineering, 1978, 12(10): 1597–1615. DOI:10.1002/(ISSN)1097-0207 |
| [6] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Bank R E, Weiser A. Some a posteriori error estimators for elliptic partial differential equations[J]. Mathematics of Computation, 1985, 44(170): 283–301. DOI:10.1090/S0025-5718-1985-0777265-X |
| [7] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Zienkiewicz O C, Zhu J Z. A simple error estimator and adaptive procedure for practical engineering analysis[J]. International Journal for Numerical Methods in Engineering, 1987, 24(2): 337–357. DOI:10.1002/(ISSN)1097-0207 |
| [8] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Deufhard P, Leinen P, Yserentant H. Concepts of an adaptive hierarchical finite element code[J]. Impact of Computing in Science and Engineering, 1989, 1(1): 3–35. DOI:10.1016/0899-8248(89)90018-9 |
| [9] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Ladevèze P, Leguillon D. Error estimate procedure in the finite element method and application[J]. SIAM Journal on Numerical Analysis, 1983, 20(3): 485–509. DOI:10.1137/0720033 |
| [10] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Becker R, Rannacher R. An optimal control approach to a posteriori error estimation in finite element method[J]. Acta Numerica, 2001, 10: 1–102. |
| [11] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Giles M B, Süli E. Adjoint methods for PDEs: A posteriori error analysis and postprocessing by duality[J]. Acta Numerica, 2002, 11: 145–236. |
| [12] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Estep D. A posteriori error bounds and global error control for approximations of ordinary differential equations[J]. SIAM Journal on Numerical Analysis, 1995, 32(1): 1–48. DOI:10.1137/0732001 |
| [13] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Rannacher R, Suttmeier F T. A feed-back approach to error control in finite element methods: Application to linear elasticity[J]. Computational Mechanics, 1997, 19(5): 434–446. DOI:10.1007/s004660050191 |
| [14] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Fidkowski K J, Darmofal D L. Review of output-based error estimation and mesh adaptation in computational fluid dynamics[J]. AIAA Journal, 2011, 49(4): 673–694. DOI:10.2514/1.J050073 |
| [15] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Bangerth W, Rannacher R. Adaptive Finite Element Methods for Differential Equations[M]. Basel: Birkh?user, 2003. |
| [16] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Oden J T, Prudhomme S. Goal-oriented error estimation and adaptivity for the finite element method[J]. Computers and Mathematics with Applications, 2001, 41(5): 735–756. |
| [17] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Gr?tsch T, Hartmann F. Pointwise error estimation and adaptivity for the finite element method using fundamental solutions[J]. Computational Mechanics, 2006, 37(5): 394–407. DOI:10.1007/s00466-005-0711-4 |
| [18] | Journal of Central South University(Science and Technology), 41(2):649-654.--> Stein E, Rüter M. Finite element methods for elasticity with error-controlled discretization and model adaptivity[C]//Stein E, de Borst R, Hughes T J R. Encyclopedia of Computational Mechanics Ⅱ Solids and Structures. New York: John Wiley & Sons, 2004. |
| [19] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Ladevèze P, Rougeot P, Blanchard P, et al. Local error estimators for finite element linear analysis[J]. Computer Methods in Applied Mechanics and Engineering, 1999, 176(1-4): 231–246. DOI:10.1016/S0045-7825(98)00339-9 |
| [20] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Chamoin L, Ladevèze P. Strict and practical bounds through a non-intrusive and goal-oriented error estimation method for linear viscoelasticity problems[J]. Finite Elements in Analysis and Design, 2009, 45(4): 251–262. DOI:10.1016/j.finel.2008.10.003 |
| [21] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Panetier J, Ladevèze P, Louf F. Strict bounds for computed stress intensity factors[J]. Computers & Structures, 2009, 87(15): 1015–1021. |
| [22] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Ladevèze P. Strict upper error bounds on computed outputs of interest in computational structural mechanics[J]. Computational Mechanics, 2008, 42(2): 271–286. DOI:10.1007/s00466-007-0201-y |
| [23] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Florentin E, Gallimard L, Pelle J P. Evaluation of the local quality of stresses in 3D finite element analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2002, 191(39): 4441–4457. |
| [24] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Ladevèze P, Pled F, Chamoin L. New bounding techniques for goal-oriented error estimation applied to linear problems[J]. International Journal for Numerical Methods in Engineering, 2013, 93(13): 1345–1380. DOI:10.1002/nme.v93.13 |
| [25] | Journal of Central South University(Science and Technology), 41(2):649-654.--> Patera A T, Peraire J. A general Lagrangian formulation for the computation of a posteriori finite element bounds[C]//Barth T J, Deconinck H. Error Estimation and Adaptive Discretization Methods in Computational Fluid Dynamics. Berlin Heidelberg: Springer-Verlag, 2003: 159-206. |
| [26] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Pares N, Bonet J, Huerta A, et al. The computation of bounds for linear-functional outputs of weak solutions to the two-dimensional elasticity equations[J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195(4-6): 406–429. DOI:10.1016/j.cma.2004.10.013 |
| [27] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Sauer-Budge A M, Bonet J, Huerta A, et al. Computing bounds for linear functionals of exact solutions to Poisson's equation[J]. SIAM Journal on Numerical Analysis, 2004, 42(4): 1610–1630. DOI:10.1137/S0036142903425045 |
| [28] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Pled F, Chamoin L, Ladevèze P. On the techniques for constructing admissible stress fields in model verification: Performances on engineering examples[J]. International Journal for Numerical Methods in Engineering, 2011, 88(5): 409–441. DOI:10.1002/nme.v88.5 |
