（计算数学专业论文）一类四阶发展型变分不等式的数值方法.pdf

一类四阶发展型变分不等式的数值方法中文摘要 一类四阶发展型变分不等式的数值方法 中文摘要 障碍问题及动态障碍问题在物理、力学及工程中有着广泛的应用。其数学描述分 别为具有障碍约束的四阶椭圆变分不等式和发展型变分不等式。目前四阶椭圆变分不 等式的数值方法多集中在有限差分法，有限元方法、罚方法、对偶方法、区域分解方 法等等。而发展型变分不等式多采用时间差分格式，如向前差分格式或n e w m a r k 方 法等处理时间项，化为椭圆变分不等式讨论。本文介绍了一种求解动态障碍问题的数 值方法。通过改进二重网格方法，给出了四阶椭圆变分不等式问题及发展型变分不等 式问题的两套数值方法。本文的主要工作如下： 第二章通过二重网格方法，构造了求解四阶椭圆变分不等式的两种数值方法。采 用罚方法将四阶椭圆变分不等式转化为等价的非线性罚方程。方法一用 m a r c h u k y a n e n k o 格式将罚问题分解为两个子问题求解。方法二直接采用n e w t o n 法 将罚方程转化为线性问题求解。章末通过二维数值算例说明了两种方法是有效的，并 分析了障碍区域、罚因子等参数对解的影响。 第三章讨论了具有障碍约束的四阶发展型变分不等式的两种数值方法。给出了两 套方法的算法流程图，通过数值算例模拟了动态障碍问题的变化过程，章末同样讨论 了障碍区域等参数的影响。比较了两种方法的计算效率。 关键词：四阶椭圆变分不等式；四阶发展型变分不等式；u z a w a 方法；罚方法；有限 元方法 作者：孙海波 指导老师：丁睿( 教授) 一类四阶发展型变分不等式的数值方法 英文摘要 t h en u m e r i c a lm e t h o df o rak i n do f e v o l u t i o n a r y v a r i a t i o n a li n e q u a l i t yo fo r d e r4 a b s t r a c t o b s t a c l ep r o b l e ma n dt h ed y n a m i co b s t a c l ep r o b l e ma r ea p p l i e dw i d e l yi np h y s i c s ， m e c h a n i c sa n de n g i n e r i n g t h em a t h e m a t i c a lm o d e l so ft h e s et w o p r o b l e m sa r ef o r m u l a t e d a sa l le l l i p t i cv a r i a t i o n a li n e q u a l i t yo fo r d e r4a n da ne v o l u t i o n a r yv a r i a t i o n a li n e q u a l i t yo f o r d e r4w i t ho b s t a c l ec o n s t r a i n tr e s p e c t i v e l y t h en u m e r i c a lm e t h o d sf o rt h i sk i n do f e l l i p t i cv a r i a t i o n a li n e q u a l i t ym a i n l yf o c u so nf e m ，d u a l i t ym e t h o d ，u z a w am e t h o da n d d o m a i nd e c o m p o s i t i o nm e t h o d o r d i n a r yw eo f t e nu s et i m es p l i t t i n gs c h e m et oc h a n g e t h ee v o l u t i o n a r yv a r i a t i o n a li n e q u a l i t yi n t oa ne l l i p t i cv a r i a t i o n a li n e q u a l i t y , f o ri n s t a n c e f o r w o r dt i m ed i f f e r e n c es c h e m eo rn e w m a r km e t h o de t c i nt h i sp a p e r , t w ok i n d so f n u m e r i c a lm e t h o d sf o rv a r i a t i o n a li n e q u a l i t ya n de v o l u t i o n a r yv a r i a t i o n a li n e q u a l i t ya r e d i s c u s s e d ，t h i st h e s i si sc o m p o s e do ft h ef o l l o w i n gs e c t i o n s ： i nt h ec h a p t e rt w o ，u s i n gp r e s e n t st w o g r i d sp r o j e c t i o na l g o r i t h m ，t w ok i n d so fm e t h o d f o rt h ee l l i p t i cp r o b l e mo fv a r i a t i o n a li n e q u a l i t yo fo r d e r4 t h ep r o b l e mi sr e d u c e di n t oa n e q u i v a l e n tn o n - l i n e a re q u a t i o nb yt h ep e n a l t ym e t h o d m a r c h u k y a n e n k os c h e m ei s a p p l i e di nt h ef i r s tm e t h o d i nt h es e c o n dm e t h o dw ed i r e c t l yu s en e w t o nm e t h o dt os o l v e t h en o n - l i n e a rp e n a l t ye q u a t i o n t w on u m e r i c a le x a m p l e sa r ei m p l e m e n t e d ，t h er e s u l t s s h o wt h ee f f i c i e n c yo ft h em e t h o d i ta n a l y z e st h ee f f e c to ft h eo b s t a c l er e g i o n ，p e n a l t y f a c t o rt ot h es o l u t i o n s i nt h ec h a p t e rt h r e e ，t w ok i n d so fm e t h o d sf o r t h ee v o l u t i o n a r yv a r i a t i o n a li n e q u a l i t yo f o r d e r4a r eo b t a i n e d t h ef l o wc h a r t sa r ep r e s e n t e di nt h i sc h a p t e r t w od i m e n s i o n a l n u m e r i c a le x a m p l e sp r o v et h ee f f i c i e n c yo ft h i sm e t h o d w ea l s od i s c u s st h ee f f e c to f o b s t a c l er e g i o n ，p a r a m e t e r st ot h es o l u t i o no fe v o l u t i o n a r yv a r i a t i o n a li n e q u a l i t y i n a d d i t i o n a l ，w ec o m p a r e dt h ee f f i c i e n c yo ft h e s et w om e t h o d s k e yw o r d s ：v a r i a t i o n a li n e q u a l i t yo fo r d e r4 ：e v o l u t i o n a r yv a r i a t i o n a li n e q u a l i t yo fo r d e r 4 ；u z a w am e t h o d ；p e n a l t ym e t h o d ；f i n i t ee l e m e n tm e t h o d ( f e m l i i w r i t t e n b ys u nh a i b o s u p e r v i s e db yp r o f d i n gr u i 苏州大学学位论文独创性声明及使用授权的声明 学位论文独创性声明 本人郑重声明：所提交的学位论文是本人在导师的指导下，独立进 行研究工作所取得的成果。除文中已经注明引用的内容外，本论文不含 其他个人或集体己经发表或撰写过的研究成果，也不含为获得苏州大学 或其它教育机构的学位证书而使用过的材料。对本文的研究作出重要贡 献的个人和集体，均已在文中以明确方式标明。本人承担本声明的法律 责任。 研究生签名： 墨尘! i 盈i 垒一日 期： 蝉! ! ： 学位论文使用授权声明 苏州大学、中国科学技术信息研究所、国家图书馆、清华大学论文 合作部、中国社科院文献信息情报中心有权保留本人所送交学位论文的 复印件和电子文档，可以采用影印、缩印或其他复制手段保存论文。本 人电子文档的内容和纸质论文的内容相一致。除在保密期内的保密论文 外，允许论文被查阅和借阅，可以公布( 包括刊登) 论文的全部或部分 内容。论文的公布( 包括刊登) 授权苏州大学学位办办理。 研究生签名： 型! 垡遮 日期： 型：! 三 导师签名： 要盎 日 期：立学。址 一类四阶发展型变分不等式的数值方法第一章前言 1 1 综述 第一章前言弟一早月i j 苗 障碍问题1 1 , 2 , 5 及动态障碍问题1 1 , 2 , 5 1 在物理、力学、工程中有大量的应用。其数学 描述分别为具有障碍约束的四阶椭圆变分不等式和发展型变分不等式。目前四阶椭圆 变分不等式的数值方法多集中在有限差分法，有限元方法，罚方法，对偶方法，区域 分解方法【1 0 1 等。而发展型变分不等式多采用时间差分格式，如向前差分格式或 n e w m a r k 方法 i h 等处理时间项，将其化为椭圆变分不等式讨论。文【4 】在处理由二阶 发展型变分不等式描述的障碍问题时，引入了二重网格投影方法。该方法具有明显的 高计算效率。文 6 】6 中将这一方法推广到动态弹塑性扭转问题。 这一方法是否可以推广到四阶椭圆变分不等式及四阶发展型椭圆变分不等式? 我 们在实际计算中发现了一些困难( 见第二章节2 1 ) 。为此对方法做了两个变形。方法 一中首先采用罚方法将原来含有约束项的四阶椭圆变分不等式化为相应的罚方程，然 后采用m a r c h u k y a n e n k o 格式将问题分解成两个子问题分别求解，考虑到投影误差较 大，这里均采用同一网格，避免了二重网格之间的投影，对非线性方程采用n e w t o n 方法处理。针对是否有必要采用m a r c h u k y a n e n k o 格式的疑问，构造了利用n e w t o n 方法直接求解罚方程的方法二。最后将两种方法推广到四阶发展型变分不等式中。分 别实现了数值算例，比较了两套方法的优劣。 1 2 论文的主要研究内容 第二章介绍了障碍问题及其描述形式，给出了具障碍约束的不等式问题的相关理 论结果。通过二重网格方法，引入了求解四阶椭圆变分不等式的两种数值方法。采用 罚方法将四阶椭圆变分不等式转化为等价的非线性罚方程。方法一用 m a r c h u k y a n e n k o 格式将罚问题分解为两个子问题求解。方法二直接采用n e w t o n 法 将罚方程转化为线性问题求解。分别给出了两种方法的实现步骤和程序流程图。通过 二维数值算例说明了方法是有效的并且分析了障碍区域、罚因子等参数对解的影响， 比较了两种方法的优劣。 第三章给出了动态障碍问题的数学描述，讨论了具有障碍约束的四阶发展型变分 一类四阶发展型变分不等式的数值方法第一章前言 不等式的两种数值方法。给出了两套方法的实现步骤，章末通过数值算例模拟了动态 障碍问题的变化过程，同样讨论了障碍区域等参数的影响并且比较了两种方法的计算 效率。 最后给出了结论和未来工作的展望。 本文的内容得到了国家自然科学基金项目( n o 1 0 2 0 1 0 2 6 ，n o 1 0 6 7 2 1 1 1 ) 的资助，作 者在此表示衷心的感谢1 2 一类四阶发展型变分不等式的数值方法第二章一类p n q 阶椭圆变分不等式的数值方法 第二章一类四阶椭圆变分不等式的数值方法 本章中，首先给出了四阶椭圆变分不等式问题的背景介绍和现有的常用算法，其 次给出了该问题的两种数值方法和步骤，最后给出了数值算例比较了本章所用的两种 不同算法的结果。 2 1 问题介绍 设i c ，p ( q ) ，f r ( q ) ，k 为硪( q ) 中的闭凸子空间，q 是r ”( n = l ，2 ，3 ) 中一个 有界闭区域或多边形区域，e 是q 的一个多边形子区域，四阶椭圆变分不等问题p l 描 述如下： 扣( 卜“) d x 互a f ( v - u ) d x v v k( 2 1 ) i “k = 1 ，iv 1 - 1 2 0 ( q ) ，v ( x ) 沙( x ) ，a e v x e ) 以二维情况为例，不等式( 2 1 ) 描述了一个水平放置的弹性薄板在垂直方向小变形。其 边界r 固定，整个薄板位于高度是y 的障碍上方，在密度为厂的垂直力场作用下，薄 板受力后始终位于障碍上方，问题( 2 1 ) 模拟了该物理过程。 在( 2 1 ) 中令 口( ) = ，a u a v d x ，( 厂，v ) = s i v , & 贝j j ( 2 1 ) 可写为：求u k ，使得 口( 甜，v 一甜) ( 厂，v = - 甜) ， v v k ( 2 1 ) 【“k = viv h ；( q ) ，v ( x ) ( x ) ，口e v x e ) 直接核验椭圆变分不等式解的存在唯一性定理1 2 1 ，可知 定理2 1 问题( 2 1 ) 的解存在且唯一。 定理2 2 问题( 2 1 ) 等价于以下极小化问题： j 山( “) 厶( 1 ，) ，v 1 ，k( 2 2 ) fu k 。 其中能量泛函为山( v ) = 丢l a v2 d x l 声出= j 1 口( v ，v ) 一( 厂，v ) 证明： ( j ) 由双线性形式a ( u ，v ) 的定义知： 一类四阶发展型变分不等式的数值方法第二章一类阴阶椭圆变分不等式的数值方法 a ( v ，v ) 一a ( u ，材) = 2 a ( u ，一“) + a ( 1 ，- u ，v - u ) 则： 山( v ) 一厶( “) = 争1 口( v ，) 一日( “，“) 一( ，v 一甜) = 圭【2 口( u , v - u ) + 口( v 一“v 一“) 卜( 厂，v 一材) = 口( z f ，v 一“) 一( 厂，一“) + 三1 口( v - u , v - u ) 由( 2 1 ) j 3 乏a 的强制性可知： 厶( v ) 一山( 甜) 0 ，v v k ( 乍) 任取，k ，u k 是( 2 2 ) 的解，由k 的定义知： u + 五( 1 ，一u ) k ，v 2 0 ，1 y , n y os o ( z f ) 是泛函s o 的极小值，所以有： 山( “) s o ( “+ 旯( v 一“) ) 令f ( 五) = s o ( “+ 旯( v 一“) ) 则上式为f ( 0 ) f ( 2 ) ，v a , 【o ，1 】，即0 是f ( 2 ) 的极小值点，所以，( 0 ) 0 而 f ( o ) = 烛业业掣 + o 五 = 驻n 1 【口( “，v 一“) 一( 厂，v 一甜) + 2 a ( v u ，v - u ) 】 a - - + u = a ( u ，1 ，一 ) 一( f ，一“) 由此可得( 2 1 ) 证毕 ( 2 1 ) 这类四阶椭圆变分不等式的数值方法主要集中在有限元方法，对偶方法u z a w a 方法0 , 2 1 ，多重网格法【9 1 等。 4 一类四阶发展型变分不等式的数值方法 第二章一类网阶椭圆变分不等式的数值方法 2 2 一类四阶椭圆变分不等式的数值方法 文献 4 】中给出了一类二阶发展型变分不等式的数值方法一二重网格投影方法。我 们希望引入该方法来处理问题( 2 1 ) 。为此下面先介绍我们的思路。 2 2 1 计算思路 二重网格算法的基本思想【4 1 是：先利用罚方法将原椭圆变分不等式问题转换为 一个非线性罚形式的变分方程；再由m a r c h u k y a n e n k o 格式川将罚方程转化为两个嵌 套求解的子问题，即一个经典的椭圆方程和一个非线性方程。对于经典的椭圆方程采 用分段常函数逼近的有限差分法【7 】或有限体积法【7 1 等这类非连续近似来求解，对非线 性方程采用分段线性有限元近似，再用n e w t o n 方法求解非线性代数方程。两个问题 的网格不同且在各自单元上近似函数一个非连续，另一个连续，为此引入投影方法， 利用投影方法建立非连续网格近似函数与另一种连续网格近似函数之间的联系。 本章希望采用这套思路求解( 2 1 ) 这种四阶椭圆不等式，但是在求解中主要有以下 几个困难： 1 四阶非线性不等式采用分段线性有限元法是不合适的，至少要采用二次元上 逼近。 2 若四阶线性椭圆方程再采用常单元可能会造成投影时的积累误差加大。 3 由于问题( 2 1 ) 中含有约束条件，因而还需要先采用相应方法处理约束项。 为了克服以上困难，首先通过罚方法将含约束条件的椭圆变分不等式问题转化为罚形 式的变分问题；然后采用m a r c h u k y a n e n k o 格式将罚方程转化为两个嵌套的子问题， 对于线性椭圆方程采用有限差分法求解，对于非线性方程采用二次元离散后再用 n e w t o n 方法求解。采用投影方法建立两个解之间的联系。 在实际计算中我们发现高次元与常单元之间的投影误差较大。此外在采用罚形式 后是否有必要将罚方程化为两个方程求解。为此我们对原二重网格法做了两个变形， 用以求解四阶椭圆变分不等式。 方法一：由于投影过程必然会在整体计算中对计算结果和算法性能产生一定的影 响，因此当利用m a r c h u k y a n e n k o 格式将罚方程转化为两个嵌套求解的子问题时，若 对两个子问题都采用有限元离散求解，这样也就不存在粗细网格近似函数之间连续性 一类四阶发展型变分不等式的数值方法第二章一类四阶椭圆变分不等式的数值方法 1 空间r ( q ) ，h 2 ( q ) ，磁( q ) 的近似 设q 是r 2 中的多边形区域，l 是q 上的三角形剖分，且l 满足正则剖分的条件， 昱是所有次数不大于2 的两变元多项式的集合。则r ( q ) ，h 2 ( q ) 的近似为： 圪= i 屹c 。( 西) ，i r 只，v 丁r 。 瑶( q ) 的近似为： v o 。= i 圪，屹= o ，d 甩r 2 ( 2 1 ) 式的罚问题 对( 2 1 ) 使用罚方法【1 1 ，将( 2 1 ) 转化为相关的罚形式 i 求“。日2 ( q ) ，s j 。僦一* 。计诎钔肌职q ) q 3 其中占是罚因子，这里取g = 1 0 - 1 ，( 口) - 定义为( 口) 一= m a x ( 一口，0 ) 在( 2 3 ) 中令罚函数( ) = ( 一缈) ：，口( ，v ) = a v d x ，则( 2 3 ) n - f 写成： f 求“。h 2 ( q ) ，s j 1 砸) 一( 触a v ) = ( 加) ，执鲋2 ( q ) 直接验证罚函数的性质，由【1 】1 中第一章节3 2 可知 定理2 3 罚形式( 2 3 ) - 与n 问题( 2 1 ) 等价。 一类四阶发展型变分不等式的数值方法第二章一类四阶椭圆变分不等式的数值方法 定理2 - 4 收敛性定理，即娥“。= “，其中u e 是罚形式的解，材是原问题( 2 1 ) 的解。 3 问题( 2 3 ) 的m a r c h u k y a n e n k o 方案 利用m a r c h u k - y a n e n k o 格式处r ( 2 3 ) ，则可以等到与( 2 3 ) 等价的两个嵌套的子问 题： 畿篡警吐机鳓 亿4 ， 1 “：诎+ 眺= ( 厂，v ) ，v v 日2 ( q ) 毕 f 求“”+ 1 日2 ( q ) ，s ， h 弦抛。训：眺= o ，v v 谢( 哟 q 巧 在用m a r c h u k y a n e n k o 分离方案对罚问题进行处理后，得到了两个不同特性的问 题。问题( 2 4 ) 是一个经典的椭圆方程，因而采用分段常函数逼近的有限差分法或有限 体积法等这类非连续近似都可以很好地用来求解此问题，且该近似问题还是适定的。 另一方面，非线性问题( 2 5 ) 对于非连续逼近格式是不适定的，必须采用连续近似， 如对求解区域进行三角剖分，并且在每一个三角形单元上采用分段线性有限元近似， 再用n e w t o n 方法求解。由第三步可知问题( 2 4 ) ，( 2 5 ) 是按顺序求解的两个相互耦合 的问题，如果求解问题( 2 4 ) 采用一种非连续的分段常函数离散格式，而( 2 5 ) 采用另一 种连续分段线性函数近似，两个问题的网格不同，并且在各自单元上近似函数一个非 连续，另一个连续，导致问题( 2 4 ) 的解无法直接带入问题( 2 5 ) ，同样问题( 2 5 ) 的解也 无法直接带入问题( 2 4 ) 。：3 1 a 投影方法，即利用投影方法将问题( 2 4 ) 求出的一种非连 续网格近似函数与问题( 2 。5 ) 求出的另一种连续网格近似函数之间建立联系。具体的 投影规则可见文献 4 】。由于实际计算发现，对于四阶问题双重网格之间投影误差很 大，因而这里没有选择投影方法，而是采用节2 2 1 中的变形，即方法一，二。 4 n e w t o n 方法求解问题( 2 5 ) 由于( 2 5 ) 是非线性问题，而非线性问题的求解通常较困难，一般采用非线性数值 方法将非线性问题转化为线性问题求解，这里采用n e w t o n 法。 问题( 2 5 ) 的左端可以写成如下形式： c ( “) = 0 一类四阶发展型变分不等式的数值方法第二章一类四阶椭圆变分不等式的数值方法 其中c ( 材) = 一甜：眺一去l ( “一) ：眺，v v 磁( q ) ， 则应用n e w t o n 方法可以得到： 已知u o h o ( q ) 对于后0 ，若u 女已知，则通过下式计算 + 。= 一c ( ) - 1 e ( ) ， 令u = 一+ ，通过计算( ) ( 一+ ，) = 0 可以得到关于+ ，的化线性方程，即： f 求小满足 t 一詈“刊一融= p 诎一“训v v 圪 g 届 以下是方法一的程序设计流程图： 原问题( 2 1 ) 罚方法 1 罚形式( 2 3 ) n k 。格式 椭圆方程( 2 4 ) 非线性方程( 2 5 ) n e w t o n 法解非线 有限元法解线性代数方程组有限元法 性代数方程组 r 代入 1 r ( 2 4 ) t 懈u ： 线性方程( 2 6 ) l lk = k + 1 l 输出结果 y 旷，i i ：多n 一类四阶发展型变分不等式的数值方法第二章一类四阶椭圆变分不等式的数值方法 方法二： 对( 2 3 ) 应用n e w t o n 法，令“( t ) 是n e w t o n 迭代的第k 次的迭代数值，令 u k = 材( i ) 一甜( 川) ，则问题转化为如下迭代形式：求甜( ) ，满足 对线性化方程( 2 7 ) 采用有限元法，然后对有限元离散后的矩阵方程应用共轭梯度法 求解。 以下是方法二的程序设计流程图： i 原问题( 2 1 ) l 罚力 1 l 罚形式c 2 - 3 ) 法 9 诎 虹 一” 眈 一、， l 力 v 1 一 ， 乞 0 一 以 撇 一共四阶发展型变分不等式韵数值方法第二章一共四阶椭圈变分不等式的数值方法 2 3 数值算倒 四阶椭圆变分不等式( 2 1 ) 中，l 玟f = 1 0 ，妒= o ，q = 1 - 1 ，1 】x 卜i ，1 】，罚因子= l o 一。 为了体现障碍区域对解的影响，分别对。的闭凸子区域e 的不同情况进行计算。采用 m a t l a b 编程计算，计算时步长 2 亩，n e 州加选代的终止条件为 恢 + l 】一叶”k 1 矿，迭代采用s o r 方法，松弛困子= 仉o 8 , 1 4 5 采用本文的两种 方法分别计算，其中m 表示第t 步n e w t o n 选代的数值解 第一种情况：e = f 2 ，计算迭代所取的松弛因子= 0 , 0 8 , i 4 5 的情况。 v 圈1 m = 0 时的三维图像图2 m 印8 时的三维图像 图3 m = 14 5 时的三维图像 第二种情况：e = ( l y ) l e + ，2e o 5 2 ，松弛因子= 0 , 0 ，8 ，14 5 的情况 一共四阶发展型变分不等式的数值方法 第= 章一共四阶喃囡变分不等式的数值方法 9 图4 m = 0 时的三维图像图5 m = 08 时的三维图像 图6 = 1 4 5 时的三维图像 第三种情况：= 嘛 l ，+ y 2 - s o 2 5 2 ，松弛困子m = 0 , 0s , 14 5 的情况 “0 渗 图7 o = o 时的三维图像图8 = 0 8 时的三维圈像 纛 一娄四骱发展型变分不等式的数值方法第= 章一类四阶椭圆变分t 等式的数值方法 图9 街= 14 5 时的三维图像 通过观察圈1 - 9 可以发现，第一种情况时( 图1 - 3 ) ，此时障碍区域是整个q 表面， 可以看到q 表面几乎都有突起；第二种情况时( 图4 - 6 ) ，障碍区域是n 内部以05 为 半径的圆域，可以看到突起部分集中与q 内部的圆域上，经过仔细观察突起部分圆域 半径约为0 5 ，而此外的平面几乎看不出有突起；第三种情况时( 图7 9 ) ，障碍区域 是n 内部以02 5 为半径的圆域，而突起部分集中于半径约为02 5 的圆域内。这是由于 在外力，恒定的情况下受力区域越小，受力的网格节点就越少，因此产生变化的网 格节点就越少，与该问题描述的物理问题是吻合的。 表1 总结了当障碍区域e 取为各种不同情况时的计算误差，解的最大值等计算结 果。 表1 不同障碍区域e 时解的最大值及误差 m 缸 e r r o r 】 00 1 0 0 44 1 e 一5 n 080 1 0 0 73 5 一5 l4 5 0 1 0 0 62 9 一5 000 6 6 1 32 一5 0 5 为半径的圆盘 0 800 6 1 822 一5 】4 500 6 1 52 1 e 一5 00 0 2 4 516 e 一5 o2 5 为半径的盟盘 0 80 0 2 4 61 _ 4 e 一5 i 4 500 2 4 51 l e 一5 一类四阶发展型变分不等式的数值方法第二章一类【r q 阶椭圆变分不等式的数值方法 其中误差e r r o r l = 愀“一u k u k 表示n e w t o n 迭代的第七步的数值解，u k + l 表示终止 步的数值解。 由表1 可以看出在障碍区域e 从q 逐渐减小为0 2 5 为半径的圆盘的过程中，由于 受力区域越来越小，所以数值解的最大值是依次下降的。 表1 中还列出了计算不同情况的障碍区域时，松弛因子的不同取法所导致数值结 果相对误差的变化情况。 表2 给出了两种不同方法在计算同一问题时误差以及时间的比较： 表2 两种方法在计算e 的第一种情况时的比较 计算方法方法1方法2 收敛准则 0 “一一甜”1 i i 1 0 一 1 1 2 e r r o r l 4 1 e 一57 o e 一5 e r r o r 2 0 7 9 计算时间( s ) 1 2 0 l1 4 4 6 其中e 仃( ) r 2 = 些矗牡，表示方法1 的得到的数值解，材：表示方法2 得到的数值解。 慨0 由表2 可以看出在相同的算法收敛条件下，第二种方法所需要的计算时间更长， 多出约4 分钟。尽管从编程角度方法- - l l 较麻烦，因为需要求解两个耦合问题，但是 由于第一个问题实现简单，且其中第二个子问题的非线性方程的形式也较方法二中要 简洁，因而实现的时间代价要小于第二种方法。这说明在求解四阶椭圆变分不等式时 采用m a r c h u k y a n e n k o 格式将罚方程分成两个子问题求解，其计算效率要高于直接求 解罚方程。 表3 给出了计算同一问题时罚因子s 的不同对求解的影响。 一类四阶发展型变分不等式的数值方法第二章一类四阶椭圆变分不等式的数值方法 表3 计算e 的第一种情况时罚因子对迭代次数、误差的影响 占 1 0 11 0 一2l o 一3 迭代次数 2 53 75 2 e r r o r l 5 6 2 e 一5 4 9 3 e 一54 1 l e 一5 由表3 可以看出随着罚因子s 的减小，伴随着迭代的次数的增多，求解的精度也 会增加，这一数据符合本章定理2 4 所陈述的收敛性。 1 4 一类四阶发展型变分不等式的数值方法第三章类阴阶发展型变分不等式的数值方法 第三章一类四阶发展型变分不等式的数值方法 本章中，首先给出四阶发展型变分问题的意义，其次给出了一类四阶发展型变分 不等式的数值方法，最后给出了数值算例和结果讨论。 3 1 问题介绍 动态障碍问题f l ，2 j 】是力学中的一个重要问题，在第二章叙述的障碍问题上如果考 虑时间效应，则可描述为一类四阶发展型变分不等式。 设y r ( q ) ，f r ( q ) ，k 为m o ( q ) 中的闭凸子空间，q 是r ”( ，? = 1 ，2 ，3 ) 中一个 有界闭区域或多边形区域，是q 的一个多边形子区域，则四阶发展型变分不等式1 1 1 如下： l 昙“( 吼v 一“( 嘞+ l “( 分( v 一“( ，) ) 出厂( 以v 一“( 呦出，v v k “k = vv 乒瑶( q ) ，v ( x ) 之v ( x ) ，a e v x e ) “( o ) = u o ( k ) ，f 【o ，t 】 令 施) = 警，嘶，v ) 2 p 掀，( 加) 2 驴 则( 3 1 ) 可写为：求u ( t ) k ，使得 ( 3 1 ) ( “o ) ，v 一“) + 口( “，1 ，一“) ( ，v 一“) ，v v k z f k = vv ，露( q ) ，v ( x ) y ( x ) ，a e v x e )( 3 1 ) u ( o ) = u o ( k ) ，f 【0 ，t 】 直接核验双线性形式的性质，由 2 】中发展型变分不等式解的存在唯一性定理可知： 定理3 1 问题( 3 1 ) 的解存在且唯一。 3 2 四阶发展型变分不等式的数值方法 ( 3 1 ) 这类四阶发展型变分不等式的数值方法主要集中在时间离散格式下的有限元 一类四阶发展型变分不等式的数值方法第三章一类四阶发展型变分不等式的数值方法 方法，u z a w a 方法，对偶方法等【l 2 1 ，这里我们也希望将第二章所述的两种方法用来 解决该问题。 方法一： 1 问题( 3 1 ) 离散时间的格式 设a t ( 0 ) 是时间步长，则对( 3 1 ) 离散时间，后的格式如下： u o = u o ； ( 3 2 ) 对于刀0 ，若u ”已知，则可通过以下问题求解u 肘1 f “”“k ，v v k t 牟旷1 胁似州胆 州，o 门 由变分不等式解的存在唯一性定理1 2 直接可得： 定理3 2 问题( 3 3 ) 的解存在且唯一。 对于v n 0 ，( 3 3 ) 等价于以下问题： j “( v 一“) d x + a t i a “( v 一“) d x _ a ti f ( v - u ) d x ，v v k ；u k ( 3 4 ) 2 问题( 3 3 ) 的罚形式 问题( 3 2 ) ，( 3 3 ) 是具有唯一解的，因此假设占0 ，并且记m a x ( o ，一v ) 为v - ，( 3 2 ) ，( 3 4 ) 的近似罚问题形式为： j 严v 出+ p an + l 心一一p 1 州v 出= 诎 5 ， 【v v 日；( q ) 3 问题( 3 4 ) 的m a r c h u k - y a n e n k o 方案 利用m a r c h u k y a n e n k o 方案处n ( 3 4 ) ，这表明在每一个时间步长，( 3 4 ) 等价于求 解以下两个嵌套的子问题： ( 1 ) 求材2 h 2 ( q ) ，s t 严写等吡+ 眺= ，川眺, w e h 2 c c 3 固 1 6 一类p q 阶发展型变分不等式的数值方法第三章一类四阶发展型变分不等式的数值方法 ( 2 ) 求“：“h 2 ( q ) ，s j 严眈巧1 1 卅诎= 。 7 ， 4 有限元方法求解问题( 3 6 ) 5 n e w t o n 方法求解l - j 题( 3 7 ) 由于( 3 7 ) 是非线性问题，而非线性问题的求解通常较困难，一般采用非线性数值 方法将非线性问题转化为线性问题求解，这里采用n e w t o n 法。 i 口- j 题( 3 7 ) 可以写成如下形式： c ( ”) = 0 其中胁) = 严竽眺巧1 p n 2 眺，v v 硪( q ) ， 则应用n e w t o n 方法可以得到： 已知u o 明( q ) 对于m 0 ，若u 。己知，则通过下式计算u 。+ 。 u 。+ l = u 。一c ( ) q e ( ) ， 令u = u 。一“川，通过计算( 甜。) ( “。- - u r n + 。) = 0 可以得到关于“，+ 。的化线性方程，即： f mv 出一半c 一吐五m 池2p 。诎一等量一妙吡一，川陀诎c 3 固n 。 q0。nq 以下是方法一的程序设计流程图： 1 7 一类四阶发展型变分不等式的数值方法 第三章一类四阶发展型变分不等式的数值方法 1 8 一类四阶发展型变分不等式的数值方法第三章一类四阶发展型銮坌至箜茎堕鏊! 量! 壁 方法二： 和第二章一样，对罚方程( 3 5 ) 不使用m a r c h u k y a n e n k o 方案，而是直接用n e w t o n 法将非线性问题转化为线性问题。 ( 3 4 ) 的罚形式为： 廿s 咖址卜以诎q 一肛叫) v d x = q j f v d x( 3 9 ) q qn n o j j i v v 瑶( q ) 对( 3 9 ) 用n e w t o n 法转化为线性问题： f 眺+ 址二诚一- 薹，u m - y ) 一五眺 nn 。q ( 3 1 0 ) = p 。慨+ 出，a u m a v 出一去且一) ：v 出一j 声出 其中二：一“。+ 。，”。是n e 吼o n 迭代的第所次迭代数值解，“川是第m + 1 次数值解。 对( 3 1 0 ) 采用有限元离散，对离散后的矩阵方程用共轭梯度法求解。 以下是方法二的程序设计流程图： 1 9 一类四阶发展型变分不等式的数值方法 第三章一类四阶发展型变分不等式的数值方法 3 3 数值算例 四阶发展型变分不等式( 3 1 ) 中，取厂= 1 0 s i n ：r t ，沙= o ，o f 1 ，q = 卜1 ，l 】 一1 ，1 】， 二差婴坠墨墨型茎竺尘兰墨堕塾堡查堡 墨三兰二差璺堕垄墨里壅坌至塑苎箜墨堕查些 t = l 。 采用m a t l a b 编程，空间步长 = 0 2 5 ，时间步长f = 0 0 1 ，算法外迭代的收敛准则 为p h h l 2 1 0 。，算法内迭代的收敛准则为帆) 一h ”k 1 0 。下面的图示依次 给出了e 取为不同区域时各时间点的计算结果： 第一种情况：e = q 图1 0 2 时刻的三维图像 图2 0 4 时刻的三维图像 i 气 图3 0 6 时刻的三维图像 图4 0 8 时刻的三维图像 图5 0 9 时刻的三维图像 2 l 一娄四阶发展型变分不辞式的数值方法第! 章一娄四阶发展型变丹不等式的数值方法 第二种情况：= 幢，) i ，2 + y 2 o 5 2 。i 问 图6 仉2 时刻的三维图像图7 0 4 时刻的三维图像 i 9 图8 0 6 时刻的三维图像图9 0 8 时刻的三维图像 图1 0 0 9 时刻的三维图像 第三种情况：占= o y ) i x 2 + y 2 s o 2 5 2 一娄四聃发展型变分不等式的数值方法 第三章一娄四阶发展型变分不等式的教值方法 图1 10 2 时刻的三维图像图1 2 0 4 时刻的三维图像 ： 图1 3 0 6 时刻的三维图像 图1 4 0 8 时刻的三维图像 图1 5 09 时刻的三维图像 一娄四阶发展型变分不等式的数值方法第= 章一类口阶发展型变分不等式的数值方法 表1 第一种障碍区域时的计算结果 时刻t 04 080 9 节点数8 l 三角形单元数1 2 8 外选代次数1 6，81 6 n e w t o n 选代次数 8 31 0 7 表1 统计了计算第一种障碍区域时的所计算的节点数，三角形单元数，外迭代次数 以及n e w t o n 迭代次数。 表2 不同障碍区域时各时间点数值解的最大值 时刻f 02040 6n 80 9 n0 0 5 2 40 0 6 7 0o0 8 1 5 0 5 为半径的圆盘0 0 2 6 80 0 3 5 70 0 4 4 70 0 5 3 6 02 5 为半径的圆盘0 0 1 0 70 0 1 2 4 表2 比较了不同障碍区域时每个时间点的解的最大值，可以看出当障碍区域相同 时随着时间的增长，该问题的数值解的最大值是逐渐增大的。而在时问点相同而障碍 区域不同时，障碍区域越大，数值解的最大值就越大，这是由于在相同时刻，外力， 是相同的，所以受力区域越小，受力的网格节点就越少，因此产生变化的网格节点就 越少。 图1 6 第一种障碍区域时对角线y = z 上各时问点数值解的比较 反 一类日阶发展型变井不等式的盐值方法第! 章一共四阶发展型变分不等式的数值芹珐 图1 7 第二种障碍区域时对角线y = x 上各时间点数值解的比较 图1 8 第三种障碍区域时对角线y = x 上各时间点数值解的比较 图1 6 - 1 8 分别计算了三种不同障碍区域下不同时刻对角线上数值解的比较，从三 张图中可以看出，当时刻逐渐增大时该问题的解也是逐渐增大的。 表3 两种方法在计算第一种障碍区域时的比较 计算方法方法1方法2 收敛准则 i i “。一“”4 | 【- 1 0 。 e r r o r 21 2 计算时n ( s 1 3 2 7 2 一类四阶发展型变分不等式的数值方法第三章一类四阶发展型变分不等式的数值方法 其中e m 咖甘，其札表示方法- 的删的数值解心表示方法2 矧的 数值解。和第二章两种方法的比较结果一样，方法二在计算该问题时的时间消耗 要多于方法一，多出约6 分钟。 2 6 一类四阶发展型交分不等式的数值方法 第p l j 章结论 第四章结论 动态障碍问题在物理、力学及工程应用中是一个很重要的问题。目前求解此类发 展型变分不等式的常用数值方法主要集中在时间离散下的有限元法、对偶法、u z a w a 方法等。本文介绍了一种求解动态障碍问题的数值方法。通过改进二重网格方法，给 出了四阶椭圆变分不等式问题及发展型变分不等式问题的两套数值方法。通过数值算 例分别实现了两套方法。比较了方法的精度、计算效率，说明了障碍区域大小对解的 影响。分析了两套方案的优劣。 本文的主要工作如下： 1 介绍了障碍问题及其描述形式，给出了具障碍约束的不等式问题的相关理论 结果。通过二重网格方法，引入了求解四阶椭圆变分不等式的两种数值方法。 采用罚方法将四阶椭圆变分不等式转化为等价的非线性罚方程。方法一用 m a r c h u k y a n e n k o 格式将罚问题分解为两个子问题求解。方法二直接采用 n e w t o n 法将罚方程转化为线性问题求解。给出了两种方法的实现步骤和流程 图。通过二维数值算例说明了方法是有效的，分析了障碍区域、罚因子等参 数对解的影响，比较了两种方法的优劣。 2 给出了动态障