ISSN１０００Ｇ００５４
CN１１Ｇ２２２３/N JTsinghuaUniv(Sci& Technol),清华大学学报 (自然科学版)２０１７年 第５７卷 第４期２０１７,Vol．５７, No．４ ３６２Ｇ３６８５/１９
两种严格界面向目标误差估计方法的等价性
郭孟武, 钟宏志 (清华大学 土木工程系,北京 １０００８４) 收稿日期:２０１５Ｇ１２Ｇ０３ 基金项目:国家自然科学基金资助项目(５１３７８２９４) 作者简介:郭孟武(１９９１—),男,博士研究生. 通信作者:钟宏志,教授,EＧmail:hzz＠tsinghua．edu．cn 摘 要:在模型检验研究中,针对离散误差的后验误差估计 扮演着重要角色.在工程有限元分析中,有关问题解答场的 标量性质的 目 标 输 出 量 是 后 验 误 差 分 析 中 人 们 所 关 注 的. 在已有的面向目标 误 差 估 计 技 术 中,有 ２种方法能够提供 目标量误差的严 格 上 下 界,即 本 构 关 系 误 差 法 与 凸 目 标 函 数优化法.该文简要总结了这２种方法,并从计算列式的一 致性和基本原理的等 价 性２个层面论证了２种方法的等价 性,给出了２种严格界方法的本质均为余能原理的论断.对 ２种方法等价性的探讨有助于结合２种表达格式的优势,并 拓展到更复杂的问题中,形成简明有效的计算列式. 关键词:微分方程的数值解法;后验误差估计;面向目标误 差估计;本构关系误差;凸目标函数约束优化;严 格界;余能原理 中图分类号:O２４１．８２ 文献标志码:A 文章编号:１０００Ｇ００５４(２０１７)０４Ｇ０３６２Ｇ０７ DOI:１０．１６５１１/j．cnki．qhdxxb．２０１７．２５．００５Equivalenceoftwostrictlybounding
approachesforgoalＧorientederrorestimation
GUOMengwu,ZHONGHongzhi (DepartmentofCivilEngineering,TsinghuaUniversity, Beijing１０００８４,China)Abstract:For modelverification whichis mainlyfocused onthe controlofdiscretizationerrorsofnumericalresults,aposteriorierror estimationplaysanimportantroleinvariousnumericaltoolssuchas thefinite element method．The outputs ofinterestare usually convertedintointegralfunctionalsovertheproblem domainfora posteriori error analysis． Among the available techniques for goalＧorientederrorestimation,onlytwoapproaches,theconstitutive relationerrorestimationandtheconstrainedoptimization method withconvexobjectivefunctions,havebeenclaimedtobeableto offerguaranteedstrictupperandlowerboundsfortheerrorsin quantitiesofinterest．Thetwoapproachesarebrieflyreviewedand theequivalenceoftheirformulationand principleisgiven．Both approachesareshowntobeessentiallybasedonthecomplementary energy theorem． The equivalence of the two approaches is instrumentalinsimplificationoferrorestimationandextensionof applicationsintomorecomplexproblems．
Keywords:numerical methods for partial differential equations
(PDEs);a posteriorierror estimation;goalＧoriented errorestimation;constitutiverelationerror;constrained optimization with convex objective functions;strict bounds;complementaryenergytheorem 有限元法作为一种十分有效的数值方法,在工 程设计与计算分析中得以广泛应用.为了控制误差 以保证计算模拟的质量,模型检验(modelverificaＧ tion)作为数值方法中的重要步骤已经得到广泛研 究.而在各种计算误差的来源中,离散误差具有主 导地位.为了评价有限元的离散误差,一系列后验 误差(aposteriorierror)估计[１３]技术得以发展,例 如:显式型(残值型)误差估计方法[４]、隐式型误差 估计方法[５６]、恢复 型 误 差 估 计 方 法[７]、阶 谱 误 差 估 计 方 法[８]、本 构 关 系 误 差 (constitutiverelation error)方法[９]等;这些方法均可对有限元分析带来 的整体离散误差进行评估. 然而在实际的有限元分析中,人们关注的往往 是一些有关问题解的目标量,例如平均位移、局部 应力等.因此,在后验误差分析时,对目标量的离 散误差评价是非常必要的,这种针对目标输出量进 行的误差估计称为面向目标误差 估 计(goalＧoriented errorestimation).面向目标 误 差 估 计 基 于 伴 随 方 法[１０１１],自２０ 世纪 ９０ 年代中期开始发展[１２１３]以 来,各种估计技术相继被提出,主要包括伴随加权 残值法[１０,１４１５]、能量模估计方法[１６]、Green函数分 解法[１７]、本构关系误差法、凸目标函数优化(conＧ strainedoptimization withconvexobjectivefuncＧ tions)法 等.与 此 同 时,面 向 目 标 误 差 估 计 也 在 Poisson方程、固体力学问题、特征值问题、时变问
郭孟武,等: 两种严格界面向目标误差估计方法的等价性 ３６３ 题、流体力学问题等诸多问题的数值求解中得到了 应用[１５,１８]. 在已有的面向目标误差估计技术中,有２种方 法 可 以 提 供 目 标 量 误 差 的 严 格 上 下 界,分 别 为 Ladevèze教授 提 出 的 本 构 关 系 误 差 方 法[１９２４]及 Patera等提出的凸目标函数约束优化法[２５２７].这２ 种方法为目标量提供可计算的严格误差界,这种严 格界性质及在各种问题中的适用性使得这２种方法 在各种面向目标误差估计技术中表现突出. 本文旨在理论上论证这２种方法的等价性.在 对２种方法的基本思路进行综述的基础上,从计算 列式上说明了２种方法的一致性,并进一步揭示二 者的基本原理均为 线 弹 性 系 统 的 余 能 最 小 值 原 理, 故具有原理上的等价性.
１ 模型问题
考 虑 定 义 在 域 Ω ⊂ R３上 的 线 弹 性 体, 其 Dirichlet边界为Ə１Ω,Neumann边界为Ə２Ω＝ƏΩ\ Ə１Ω.该 结 构 上 作 用 有 体 力 fp∈ [L２(Ω)]３,另 在 Ə１Ω 上 给 定 位 移 up,在 Ə２Ω 上 给 定 面 力 Tp∈ [L２(Ə２Ω)]３. 在本文中,采用如下集合记号: U ＝ {u ∈ [H１(Ω)]３:u|Ə １Ω＝up}, V ＝ {v∈ [H１(Ω)]３:v|Ə １Ω ＝０}, S ＝ {σ ∈ [L２(Ω)]３×３:σ ＝σT}． ì î í ï ï ï ï (１) 并 定 义 退 化 能 量 内 积au(, ):[H１(Ω)]３× [H１(Ω)]３→ R及 相 应 的 半 模 ‖‖ u,能 量 内 积 aσ(,):S×S→R及相应的模 ‖‖σ 以 及 有 界 线性泛函l():[H１(Ω)]３→R如下: au(u,v)＝∫
Ωε(u):H:ε(v)dΩ, ‖u‖u＝ au(u,u), aσ(σ,τ)＝∫
Ωσ:H －１:τdΩ, ‖σ‖σ＝ aσ(σ,σ), l(v)＝∫
ΩfpvdΩ＋∫
Ə２ΩTpvdƏΩ． ì î í ï ï ï ï ï ï ï ï ï ï ï (２) 其中:H 为材料的 Hooke张量;ε(u)表示 Cauchy 应变张量,ε(u)＝(Ñu＋uÑ)/２. 那么,该模型问题(原问题)可如下表述:求解 定义在Ω 上的位移场u 和应力场σ,使得满足以下 ３组方程. 协调条件: u ∈U; (３) 平衡条件: σ ∈S, aσ(σ,H:ε(v))＝l(v), ∀v∈V; (４) 本构方程: σ ＝ H:ε(u)． (５) 进一步地,综合式(３)—(５),该边值问题的弱形式 描述为求解u∈U 使得 au(u,v)＝l(v), ∀v∈V． (６) 若采用Galerkin有限元对该问题进行离散,并采用 插值空 间Ph,则该问题的有限元格式为求解uh∈ Uh使得 au(uh,v)＝l(v), ∀v∈Vh． (７) 其中:Uh＝U∩Ph,Vh＝V∩Ph.而相应的有限元 应力解由本构关系式(５)得到,即σh＝H:ε(uh). 基于该模型问题进行面向目标分析,不失一般 性地,将目标量I定义为关于位移场u 的有界线性 泛函lO:[H１(Ω)]３→R,即 I＝lO(u)． (８) 而其计算值Ih＝lO(uh)的误差(I－Ih)是本文关注 的重点.下述２种方法均可给出(I－Ih)的严格上 下界,以作为面向该目标量的离散误差估计.２ 本构关系误差方法
本构关系误差建立在容许场的概念上,通过本 构关系定义线弹性问题的误差;通过有限元解选定 容许场后,本构关系误差可为位移误差场的能量模 提供严格上界;基于此原理,在面向目标误差估计 中,由本构 关 系 误 差 可 得 到 目 标 量 误 差 的 严 格 上 下界. ２．１ 本构关系误差的概念 本 构 关 系 误 差 的 概 念 基 于 模 型 问 题 的 容 许 场 (u~,σ~),其中 u~ 为 满 足 协 调 条 件 式(３)的位移场, 称为协调位移场,σ~ 为 满 足 平 衡 方 程 式(４)的应力 场,称为平衡应力场.而本构关系误差则通过本构 关系式(５)进行度量: eCRE(u~,σ~)＝ ‖σ~－H:ε(u~)‖σ． (９) 可以证明本构关系误差具有如下性质:eCRE(u~,σ~)＝０⇔(u~,σ~)＝ (u,σ) a．e．,
(１０a) e２ CRE(u~,σ~)＝ ‖u－u~‖２u＋‖σ－σ~‖σ２, (１０b) eCRE(u~,σ~)＝２‖σ－σ~∗‖σ, σ~∗∶＝ (σ~＋H:ε(u~))/２． (１０c)
３６４ 清 华 大 学 学 报(自 然 科 学 版) ２０１７,５７(４) u~＝uh;然而,由于σh不满足平衡方程式(４),故平 衡应力场σ~＝σ~h需要额外构造.该应力场σ~h可基 于有限元应力 解σh通 过 一 种 称 为 延 长 条 件 的 能 量 关系[１,２８]进行构造,即 aE σ(σ~h－σh,H:ε(v))＝０, ∀v∈ [H１(E)]３∩Ph． (１１) 其中:E 为有限元分析时划分出的每一个单元,而 aE σ 表示aσ在单元E 上的限制. 由式(１０b)易见位移误差场e∶＝u－uh能量模 的严格上界,即
‖e‖u≤eCRE(uh,σ~h)． (１２)
２．２ 针对目标量的伴随问题 如式(８)所示,模型问题解的标量目标量I 作 为[H１(Ω)]３ 上 的 有 界 线 性 泛 函,常 具 有 如 下 的 形式: lO(u)＝
∫
ΩfdudΩ＋∫
Ə２ΩTdudƏΩ． (１３) 其中:fd∈ [L２(Ω)]３,Td∈ [L２(Ə２Ω)]３ 为 提 取 因子. 针对该目标量引入伴随问题:求解定义在Ω 上 的位移场w 和应力场τ,使其满足 协调条件: w ∈V; (１４) 平衡条件: τ∈S, aσ(τ,H:ε(v))＝lO(v) ∀v∈V; (１５) 本构方程: τ＝ H:ε(w)． (１６) 该边值问题的弱形式描述为求解w∈V 使得: au(w,v)＝lO(v) ∀v∈V． (１７) 该问题的有限元格式为求解wh∈Vh使得: au(wh,v)＝lO(v) ∀v∈Vh． (１８) 与模型问题(原问题)类似,进一步可得有限元应力 解τh,并通过延长条件构造平衡应力场τ~h. ２．３ 目标量误差的表示与估计 由式(１７),考虑到e∈V,并结合 Galerkin正交 性,即对∀v∈Vh,a(e,v)＝０,可得目标量误差I－ Ih＝lO(e)的表达式为 I－Ih＝au(e,ε)． (１９) 其中,ε∶＝w－wh为伴随问题的位移误差场,由式 (１２)并借助 Cauchy不等式,该误差的严格上下界 可由下式给出:|I－Ih|≤eCRE(uh,σ~h)eCRE(wh,τ~h)．(２０)
除此之外,借由本构关系误差的性质式(１０c)
可进一步推证:
I－Ih－ １
２aσ(σ~h－H:ε(uh),τ~h－H:ε(wh)) ≤
１
２eCRE(uh,σ~h)eCRE(wh,τ~h)． (２１)
则目标量I的误差的严格、可计算上下界可表示为
Iupper－Ih＝ １
２aσ(σ~h－H:ε(uh),τ~h－H:ε(wh))＋
１２eCRE(uh,σ~h)eCRE(wh,τ~h),
Ilower－I
h＝ １２aσ(σ~h－H:ε(uh),τ~h－H:ε(wh))－
１２eCRE(uh,σ~h)eCRE(wh,τ~h)．
ì î í ï ï ï ï ï ï ï ï ï ï (２２)
３ 凸目标函数约束优化方法
凸目标函数约束 优 化 方 法 将 所 求 误 差(I－Ih) 表示为凸目标函数,并转化为约束优化问题;再通 过Lagrange形式将其转化为无约束优化问题,由 此得到误差的严格界. ３．１ Lagrange形式与目标量误差上下界 由式(６)和(７),模型问题精确解u 的未知使得 目标量的误差I－Ih＝lO(u－uh)可表示为以下约束 优化形式: I－Ih＝min e∈V[l O(e)＋κ(a u(e,e)－R(e))], s．t． au(e,w)＝R(w) ∀w ∈V． (２３) 其中:κ∈R＋;R()∶＝l()－a u(uh,)为定义 在[H１(Ω)]３ 上 的 有 界 线 性 泛 函,称 为 残 值 泛 函. 且值得注意 的 是,目 标 函 数lO(e)＋κ(a u(e,e)－ l(e))具有凸性. 当采用插值空 间Ph 进 行 有 限 元 分 析 时,将单 元划分的具体形式记作Γh,而将所有内部的单元交 界及单元与边界的交界组成的集合记作Ξ(Γh),则 函数空间V 亦可表示为 V ＝{
v∈V^:０＝b(v,q)∶＝∑
γh∈Ξ(Γh∫
)γh [v]γhq|γhds ∀q∈Q}
． (２４) 其中:V^ 为 破 坏V 在 Ξ(Γh)上 所 有 连 续 性 (包 括 Dirichlet边界条件)形成的函数空间;[v]γh当γh为 内部交界时为v 的跳跃,而当γh在边界上时为v 的 迹;Q 的定义为Q＝{q∈L２(Ξ(Γh)):q|Ə２Ω＝０}. 因此,(I－Ih)可进一步地表示为如下的约束 优化问题:郭孟武,等: 两种严格界面向目标误差估计方法的等价性 ３６５ I－Ih＝min e∈T[l O(e)＋κ(a u(e,e)－R(e))], T ＝ e∈V^:au(e,w)＝R(w) ∀w ∈V; b(e,q)＝０ ∀q∈Q
{
}
． (２５) 在此基础上,通过引入 Lagrange乘子将上述约束 优化问题转化为无约束优化问题,相应的 Lagrange 形式为 L±(e,w,q)∶＝ ±lO(e)＋κ(a u(e,e)－R(e))＋R(w)－ au(e,w)＋b(e,q)． (２６) 则误差(I－Ih)为 ± (I－Ih)＝ min e∈V^w∈V,maxq∈Q L±(e,w,q)． (２７) 进一步地,对于取定的w± h∈V 与qh±∈Q, ± (I－Ih)≥ max w∈V,q∈Qmine∈V^L ±(e,w,q)≥ min e∈V^ L± (e,w± h,qh±)． (２８) 即 min e∈V^ L＋(e,w＋ h,qh＋)≤I－Ih≤－min e∈V^ L－ (e,w－ h,qh－)． (２９) 给出了(I－Ih)的严格上下界. ３．２ 误差上下界计算方法 式(２８)中的w± h 与qh± 取 在Lagrange形式(２６) 在子空间V^h×Vh×Qh⊂V^×V×Q 中的鞍点处, L± (e± h,wh±,qh±)＝ min e∈V^h max w∈Vh,q∈QhL ± (e,w,q)． (３０) 其中:V^h＝V^∩Ph;Qh＝ {q∈Πh(γh),∀γh∈Ξ (Γh):q|Ə２Ω＝０};Πh()表示与 Ph 对 应 的 多 项 式 空间.该鞍点问题的驻值条件为 b(v,q± h)±lO(v)＋κ(２au(eh±,v)－R(v))－ au(v,wh±)＝０ ∀v∈V^h, (３１a) R(v)－au(eh±,v)＝０ ∀v∈Vh, (３１b) b(e± h,p)＝０ ∀p ∈Qh． (３１c) 式(３１c)表 明e± h ∈Vh,进 一 步 可 由 式 (３１b)得 到 e± h ＝０,则式(３１a)简化为 au(v,wh±)＝±lO(v)－κR(v)＋b(v,qh±) ∀v∈V^h． (３２) 当 取∀v∈Vh⊂V^h 时,该 式 退 化 为 au(v,wh)＝ lO(v),而wh±＝ ±wh.将 该wh± 带 入 式(３１a)并 令 q± h＝∓q０h＋κq１h,可将该式分解为 b(v,q０h)＝lO(v)－au(v,wh) ∀v∈V^h, (３３a) b(v,q１h)＝l(v)－au(uh,v) ∀v∈V^h． (３３b) 由此,考虑式(２９)中的优化问题,即求解e＝e~h±∈ V^ 使 得L±(e,w± h,qh±)取 最 小 值,相 应 的 驻 值 条 件为 ±lO(v)＋κ(２a u(e~h±,v)－R(v))－ au(v,wh±)＋b(v,qh±)＝０ ∀v∈V^． (３４) 令e~± h＝∓e~０h/κ＋e~１h,则上式可分解为 ２au(e~０h,v)＝lO(v)－b(v,q０h)－ au(v,wh) ∀v∈V^, (３５a) ２au(e~１h,v)＝l(v)－b(v,q１h)－ au(uh,v) ∀v∈V^． (３５b) 该式为确定e~０h与e~１h的控制方程. ３．３ 可计算的误差上下界表达 考虑式(３４)与 R(w± h )＝０,可得到 L±(e,wh±, q± h)的最小值为 min e∈V^ L±(e,w± h,qh±)＝L± (e~h±,wh±,qh±)＝ －κ‖ e~± h‖２u． (３６) 进一步 代 入e~h±＝me~ ０h/κ＋e~１h可 将 式(２９)重新表 达为 |I－Ih－２au(e~０h,e~１h)|≤ κ‖ e~１h‖２u＋ １κ ‖ e~０h‖２u． (３７) 而当κ＝ ‖e~０h‖u/‖e~１h‖u 时 该 不 等 式 右 侧 取 最 小值,故可得到(I－Ih)的严格、可计算上下界如下:Iupper－Ih＝２au(e~
０h,e~１h)＋２‖e~１h‖u‖e~０h‖u,
Ilower－I
h＝２au(e~０h,e~１h)－２‖e~１h‖u‖e~０h‖u．
{
(３８)４ 两种方法的等价性
如前 文 所 述,在本构误差关系法(后文简称为 本构法)与凸目标函数优化法(后文简称为优化法) 的求解步骤中,均采用有限元法求解原问题(见式 (７))和完全一致的伴随问题(见式(１８)与式(３２)), 得到原问题的近似解uh与伴随问题的近似解wh.不 难看出,通过进一步论证本构法中对平衡应力场σ~h 与τ~h的 构 造 同 优 化 法 中e~０h与e~１h的 求 解 具 有 一 致 性,即可说明２种方法计算目标量严格误差界的列 式的一致性.此外,这种列式上的一致性的背后是 基本原理上的等价性.本节通过讨论这些内容来说 明２种方法的等价性.３６６ 清 华 大 学 学 报(自 然 科 学 版) ２０１７,５７(４) ４．１ 计算列式的一致性 在采用本构法进行面向目标误差估计时,对任 意单元E,平衡应力场σ~h 的 构 造 所 遵 循 的 原 则 包 括(以处理原问题为例): １)单 元 内 部 满 足 平 衡 方 程,且 单 元 边 界 与 Neumann边界重合时满足 Neumann边界条件,即 Ñσ~h＋fp＝０ inE, nσ~h＝T p onƏE ∩Ə２Ω．
{
(３９) 或者表示为弱形式,即 aE σ(σ~h,H:ε(v))＝ lE(v)－∫
ƏEphvdΩ ∀v∈ [H １(E)]３．(４０a) 其中l(v)在单元E 上的限制lE(v)与单元内部边界 面力ph定义如下: lE(v)＝∫
EfpvdΩ＋∫
ƏE∩Ə２ΩT pvdƏΩ, ph＝ －nσ ~ h, onƏE\(ƏE ∩Ə２Ω), ０, onƏE ∩Ə２Ω．{
(４０b) ２)单元交界处面力平衡,即对∀γh∈Ξ(Γh)\ (Ξ(Γh)∩Ə２Ω),γh＝ƏE∩ƏE′,满足 ph|γh∩ƏE＋ph|γh∩ƏE′＝０． (４１) ３)单 元 内 满 足 延 长 条 件 式 (１１),若 结 合 式 (４０a),延长条件还可以表达为如下形式:∫
ƏEphφidΩ ＝l E(v)－aE σ(σh,H:ε(φi)), i＝１,２,． (４２) 其中φi为描述H１(E)∩Ph的每一个形函数;此外, 为 了 保 证 该 平 衡 应 力 场 构 造 的 可 解 性,还 要 求 ph|γh∩ƏE∈Πh(γh). 在本构方法的具体实施中,满足以上３个条件 即 可 构 造 出 平 衡 应 力 场σ~h.不 难 看 出,式 (４１)、 (４２)与 式 (３３b)完 全 一 致;进 一 步 对 比 可 见,式 (４０a)与式(３５b)一致.考虑到伴随问题中平衡应力 场的构造与原问题同理,２种方法中相应量的等价 关系如下: q１h|γh ＝ph|γh∩ƏE, q０h|γh ＝sh|γh∩ƏE, H:ε(e~１h)＝ σ ~ h－H:ε(uh) ２ , (４３a) H:ε(e~０h)＝ τ ~ h－H:ε(wh) ２ ． (４３b) 其中sh为 伴 随 问 题 中 与ph 定 义 一 致 的 单 元 边 界 面力. 由式(４３)所示的等价关系可见,由本构法得到 的误差上下界式(２２)与由优化法得到的误差上下界 式(３８)是对应相等的. 值得说明的 是,在本构法构造平衡应力场时, 延长条件是“空降”的,但其具有很优秀的“平衡性 质”.１)由式(４２)构造的任意单元E 的边界力ph与 该单元上给定外力是平衡的:若此单元发生任意形 式的刚体位移r,该刚体位移可通过完备的形函数 线性表示,即r＝Σriφi;考 虑 到ε(r)＝０,由 式 (４２)可得: lE(v)－∫
ƏEphrdΩ ＝０． (４４) 即意味着边界力与给定外力平衡;这种性质是通过 延长条件式(４２)从平衡方程式(４０a)处继承来的. ２)如式(４２)所示,对于所有公用节点的单元,等号 左边相加即公用节 点 边 界 上 ph 在 该 点 处 的 等 效 荷 载之和,由于共用边界的单元在该边界处的面力等 大反向如式(４１),此加和应为０,而等号右边相加 之和为０则是由有限元列式直接决定的. 然而在优化法中,延长条件式(３３)是由驻值条 件式(３２)得到,并非“空降”,从优化法的基本思路 也可进一步认识上述２种“平衡性质”的根源.对于 刚体位移r,au(r,v)＝０,∀v∈[H１(E)]３∩Ph,则 式(３３)可保证单元边界力与外力的平衡;当考察公 用节点边界在该点处的等效荷载之和时,这些公用 节点单元之间的连续性(协调性)没有被破坏,当式 (３３)等号右端限定在这些单元时,有限元列式保证 其为零. ４．２ 基本原理的等价性 在本节中,从另一视角分别看待２种方法,进 一步考察２种方法的一致性,并得到二者在原理上 的等价性. 对于本构法,目标量误差式(１９)可以进一步表 示为 I－Ih＝ １４‖
κe＋ １ κε‖
２ u－ １ ４‖
κe－ １κε‖
２ u． (４５) 其中κ∈R＋,则(I－I h)的上下界可表示为 － １４‖
κe－ １ κε‖
２ u≤I－Ih≤ １ ４‖
κe＋ １κε‖
２ u． (４６) 由式(１０b)可得 １ ４‖
κe± １κε‖
２u ≤郭孟武,等: 两种严格界面向目标误差估计方法的等价性 ３６７ １ ４
‖
κ(σ~h－H:ε(uh))± １ κ(τ ~ h－H:ε(wh))‖
２u＝ ± １２aσ(σ~h－H:ε(uh),τ~h－H:ε(wh))＋ κ ４e２CRE(uh,σ~h)＋ １ ４κe２CRE(wh,τ~h)． (４７) 可见这种 形 式 已 与 优 化 法 的 式(３７)一 致.当κ＝ eCRE(wh,τ~h)/eCRE(uh,σ~h)时,得到如式(２２)所示的可计算的目标误差上下界. 对于优化法,目标量误差式(２３)亦可表示为 ± (I－Ih)＝ min e∈V maxw∈VL ± ０(e,w)． (４８) 其中 L± ０(e,w)∶＝ ±lO(e)＋κ(a
u(e,e)－R(e))＋R(w)－au(e,w)＝
L±(e,w,q)－b(e,q)． (４９) 而对于伴随问题的解w± h＝±wh, ± (I－Ih)≥ min e∈VL ± ０(e,wh±)． (５０) 该式右侧对应的驻值条件为e± ０∈V, ±lO(v)＋κ(２a u(e０±,v)－R(v))－au(v,wh±)＝０ ∀v∈V． (５１) 该驻值点e± ０ 可用原问题与伴随问题的位移误差场e 与ε表示如下: e± ０ ＝ １２æèçe∓ １κεö ø ÷, (５２) 代入式(４９)与(５０)可得 ± (I－Ih)≥L０±(e０±,wh±)＝－ １４
‖
κe∓ １ κε‖
２ u． (５３) 可见与本构法的式(４６)完全相同.然而,由于e与 ε 是不可计算的,故需要将上式所示的界进一步放 宽.优化法对此作如下考虑:由于 mine∈VL± ０(e,wh±)＝ min e∈VL ± (e,w± h,qh±)≥ min e∈V^ L± (e,w± h,qh±), (５４) 故目标量的严格上下界由min e∈V^ L±(e,w± h ,qh±)给出, 如式(３６)—(３８)所示. 从上述讨论可见,本构法与优化法都得到相同 的误差界式(４６)或(５３),并各自采用不同的思路对 该误差界进行放松,得到可计算的形式.其中,本 构法引入了平衡应力场与本构关系误差,如式(４７) 所示;优化法则继续放松约束,如式(５４)所示;这 ２种思路看似相异,但其原理完全一致,即线弹性 体的余能最小值原理. ２种方法对于误差界式(４６)或(５３)的放松,都 归结于求‖ κe∓ε/κ‖u的可计算上界.对于ξκ∓＝ κe∓ε/κ,可 视 为 如 下 问 题 的 解:求 解ξκ∓ ∈V 使得: au(ξκ∓,v)＝ κR(v)∓RO(v)/ κ． (５５) 其中RO()∶＝lO()－a u(wh,).若将此问题 视为一个线弹性问题,则其对应的平衡条件为 E∓ κ∶＝{ζκ∓ ∈ζ:aσ(ζκ∓,H:ε(v))＝ κR(v)∓RO(v)/ κ ∀v∈V}． (５６) 考虑到该线弹性体系的内力势能与余能值相等,由 余能最小值原理可得: １ ２ ‖ κe∓ε/ κ‖ ２ u＝ min ζκ∓∈Eκ∓ １ ２ ‖ζκ∓‖σ２≤ １ ２ ‖ζκ∓∗‖２σ． (５７) 其中ζ∓∗ κ ∈Eκ∓.依 据 前 文 对 本 构 法 与 优 化 法 的 讨 论,２种方法对ζ∓∗ κ 的取法为 ζκ∓∗ ＝ κ(σ~h－H:ε(uh))∓ １ κ(τ ~h－H:ε(wh)), (５８a) ζκ∓∗ ＝２H:ε κe~１h∓ １ κe ~ ０h æ è ç ö ø ÷＝ ２ κH:ε(e~± h)． (５８b) 前文已述这２种取法完全相等.回到目标误差估计 问题, ± (I－Ih)≥－ １ ４‖
κe∓ １κε‖
u２ ≥ － １４ ‖ζκ∓∗‖２ σ． (５９) 在２种方法中, － １４ ‖ζ∓∗ κ ‖σ２＝－ １４‖ κ(σ~h－H:ε(uh))∓ １ κ(τ ~ h－H:ε(wh))‖２u＝－κ‖ e~h±‖２u．(６０) 即式(４７)和式(３６).５ 结 论
本文论证了本构关系误差法与凸目标函数约束 优化法在对有界线性泛函目标量进行面向目标误差 估计时的计算列式 的 一 致 性 和 基 本 原 理 的 等 价 性; 经过对比讨论,２种方法基本思想的实质均为余能 最小值原理.本构关系误差形式简明、力学意义明 确,但表达格式受到本构形式的限制,在拓展到其 他问题时处理方案不甚直观;凸目标函数约束优化 法表述方式比较复杂,但数学本质更为明确,数学 上的一般性更为明显,便于向更复杂的问题进行拓３６８ 清 华 大 学 学 报(自 然 科 学 版) ２０１７,５７(４) 展.对２种方法等价性的讨论有助于２种表达格式 的互补,以在更复杂的问题中形成实用有效、简明 清晰的表达形式. 参考文献 (References) [１] LadevèzeP,PelleJP．MasteringCalculationsinLinearand NonlinearMechanics[M]．New York:Springer,２００４． [２] AinsworthM,OdenJT．A PosterioriErrorEstimationin
FiniteElement Analysis [M]．New York:John Wiley & Sons,２０００．
[３] Grätsch T, Bathe K J． A posteriori error estimation techniquesinpracticalfiniteelementanalysis[J]．Computers
&Structures,２００５,８３(４):２３５ ２６５．
[４] BabuškaI,Rheinboldt W C．Errorestimatesforadaptive finite element computations [J]．SIAM Journal on
NumericalAnalysis,１９８３,２０(４):７３６ ７５４．
[５] BabuškaI,Rheinboldt W C．A posteriorestimatesforthe finite element method [J]．International Journal for Numerical Methods in Engineering, １９７８, １２(１０): １５９７ １６１５．
[６] BankRE,WeiserA．Someaposteriorierrorestimatorsfor ellipticpartialdifferentialequations [J]．Mathematics of Computation,１９８５,４４(１７０):２８３ ３０１．
[７] ZienkiewiczO C,ZhuJZ．A simpleerrorestimatorand adaptiveprocedurefor practicalengineering analysis [J]．
International Journal for Numerical Methods in Engineering,１９８７,２４(２):３３７ ３５７．
[８] Deufhard P,Leinen P, Yserentant H． Concepts of an adaptivehierarchicalfinite elementcode [J]．Impactof ComputinginScienceandEngineering,１９８９,１(１):３ ３５． [９] LadevèzeP,Leguillon D．Errorestimateprocedureinthe finiteelementmethodandapplication[J]．SIAM Journalon
NumericalAnalysis,１９８３,２０(３):４８５ ５０９．
[１０]BeckerR,RannacherR．Anoptimalcontrolapproachtoa posteriorierrorestimationinfiniteelementmethod[J]．Acta Numerica,２００１,１０:１ １０２．
[１１]GilesM B,SüliE．AdjointmethodsforPDEs:Aposteriori error analysis and postprocessing by duality [J]．Acta Numerica,２００２,１１:１４５ ２３６．
[１２]EstepD．Aposteriorierrorboundsandglobalerrorcontrol forapproximations ofordinary differentialequations [J]．
SIAM JournalonNumericalAnalysis,１９９５,３２(１):１ ４８． [１３]RannacherR,SuttmeierFT．AfeedＧbackapproachtoerror controlin finite element methods: Application to linear elasticity[J]．ComputationalMechanics,１９９７,１９(５):４３４ ４４６．
[１４]FidkowskiKJ,DarmofalDL．ReviewofoutputＧbasederror estimation and mesh adaptation in computational fluid dynamics[J]．AIAAJournal,２０１１,４９(４):６７３ ６９４． [１５]Bangerth W,RannacherR．AdaptiveFiniteElementMethods forDifferentialEquations[M]．Basel:Birkhäuser,２００３． [１６]OdenJT,PrudhommeS．GoalＧorientederrorestimationand adaptivityforthefiniteelementmethod [J]．Computersand MathematicswithApplications,２００１,４１(５):７３５ ７５６．
[１７]Grätsch T,Hartmann F．Pointwiseerrorestimation and adaptivityforthefiniteelement methodusingfundamental solutions[J]．ComputationalMechanics,２００６,３７(５):３９４ ４０７．
[１８]SteinE,RüterM．Finiteelementmethodsforelasticitywith errorＧcontrolled discretization and modeladaptivity [C]// Stein E,de Borst R, Hughes T J R．Encyclopedia of Computational Mechanics Ⅱ Solids and Structures． New York:John Wiley& Sons,２００４．
[１９]LadevèzeP,RougeotP,Blanchard P,etal．Localerror estimatorsforfiniteelementlinearanalysis [J]．Computer Methodsin Applied Mechanics and Engineering,１９９９, １７６(１ ４):２３１ ２４６．
[２０]ChamoinL,LadevèzeP．Strictandpracticalboundsthrough anonＧintrusiveandgoalＧorientederrorestimationmethodfor linear viscoelasticity problems [J]．Finite Elements in AnalysisandDesign,２００９,４５(４):２５１ ２６２．
[２１]PanetierJ,LadevèzeP,LoufF．Strictboundsforcomputed stressintensityfactors[J]．Computers& Structures,２００９, ８７(１５):１０１５ １０２１．
[２２]LadevèzeP．Strictuppererrorboundsoncomputedoutputs of interest in computational structural mechanics [J]．
ComputationalMechanics,２００８,４２(２):２７１ ２８６． [２３]FlorentinE,GallimardL,PelleJP．Evaluationofthelocal
quality of stresses in ３D finite element analysis [J]．
ComputerMethodsinApplied MechanicsandEngineering, ２００２,１９１(３９):４４４１ ４４５７．
[２４]LadevèzeP,PledF,ChamoinL．Newboundingtechniques forgoalＧorientederrorestimationappliedtolinearproblems [J]．International Journal for Numerical Methods in Engineering,２０１３,９３(１３):１３４５ １３８０．
[２５]PateraAT,PeraireJ．AgeneralLagrangianformulationfor thecomputationofaposteriorifiniteelementbounds [C]// Barth T J,Deconinck H．Error Estimation and Adaptive Discretization Methodsin Computational Fluid Dynamics． BerlinHeidelberg:SpringerＧVerlag,２００３:１５９ ２０６． [２６]ParesN,BonetJ,Huerta A,etal．Thecomputationof
boundsforlinearＧfunctionaloutputsofweaksolutionstothe twoＧdimensionalelasticityequations[J]．Computer Methods inApplied MechanicsandEngineering,２００６,１９５(４ ６):
４０６ ４２９．
[２７]SauerＧBudgeA M,BonetJ,Huerta A,etal．Computing boundsforlinearfunctionalsofexactsolutionstoPoissons equation[J]．SIAM JournalonNumericalAnalysis,２００４, ４２(４):１６１０ １６３０．
[２８]PledF,Chamoin L,Ladevèze P．Onthetechniquesfor constructingadmissiblestressfieldsin modelverification: Performanceson engineering examples [J]．International Journal for Numerical Methodsin Engineering,２０１１, ８８(５):４０９ ４４１．
