（计算数学专业论文）自由液面不可压流体运动的legendre谱元法.pdf

中文摘要 中文摘要 带自由液面的流体流动问题在工程和自然界中大量存在。这类问题的难点在于 如何准确的跟踪模拟自由液面随时间的变化过程。此类问题的控制方程为粘性不可 压n a v i e r - s t o k e s 方程过去人们采用多种数值方法研究带自由液面的流体问题l 州，主要 有差分法，有限元法和边界元法等。本文采用l e g e n d r e 谱元法进行数值模拟。在固定区 域谱元法的基础上作了如下的扩展：用全粘性压力张量格式使得牵引边界条件( 自由液 面) 能很自然的应用；用a l e 方法来精确的跟踪自由液面的位置；用半隐式时间离散格 式和改进的u z a w a 算法一定程度上分离自由液面的演变和区域内部的演变；用s t o k e s 算 子代替弹性算子来控制网格速度。 为了使问题易于分析且更明确，本文主要研究沿平板下落薄膜流动的线性稳定性问 题。沿平板下落薄膜流动问题是一种很典型的带自由液面的流体。在数值试验中本文采 用一个有效的计算格式，在空间方向上采用l e g e n d r e 谱元法，在时间方向上采用半隐式 方法；对移动网格的控制采用s t o k e s 方程。数值结果表明这种计算格式得到的结果不仅 与o r r - s o m m e r f e l d 方程线性稳定性的结果相吻合，两且相比于椭圆弹性算子，s t o k e s 算 子可较好的保证谱元体积守恒，从而证实了a l e 谱元法用于求解此类问题的有效性。 关键词：移动网格；n a v i e r - s t o k e s 方程；l e g e n d r e 谱元法；a l e 方法；线性稳定性。 v a b s t r a c t a b s t r a c t u n s t e a d yi n c o m p r e s s i b l ev i s c o u sf r e e - s u r f a c ef l o w sa r el a r g e l ye n c o u n t e r e di ne n g i - n e e r i n ga n dn a t u r a ls y s t e m s t h ed i f f i c u l t yo ft h i sk i n do fp r o b l e m sr e l i e so nh o wt o a c c u r a t e l yt r a c kt h ep o s i t i o no ft h ef r e es u r f a c ev i at h et i m ee v o l u t i o n i nt h ep a s tt h e p e o p l eu s e dm a n yn u m e r i c a lm e t h o d st oi n v e s t i g a t et h eb e e - s u r f a c ef l o w si n c l u d i n gf d m ， f e m b e m 1 卅8 j 1 ds oo n i nt h i sp a p e rt h el e g e n d r es p e c t r a le l e m e n tm e t h o di sp r o p o s e d t os o l v et h ev i s c o u si n c o m p r e s s i b l ef r e e - s u r f a c ef l o w s t h em a i ni n g r e d i e n t so fo u rm e t h o d c o n s i s t ：u s eo ft h ef u l lv i s c o u ss t r e s st e n s o rf o rn a t u r a li m p o s i t i o n so ft h et r a c t i o n ( s u r f a c e t e n s i o n ) b o u n d a r yc o n d i t i o n s ；u s eo fa l em e t h o d sf o ra c c u r a t er e p r e s e n t a t i o no fm o v i n g b o u n d a r i e s ；u s eo fs e m i i m p l i c i tt i m e - s t e p p i n gp r o c e d u r e sa n di m p r o v e du z a w am e t h o d t op a r t i a l l yd e c o u p l et h ef r e e - s u r f a c ee v o l u t i o na n di n t e r i o rn a v i e r - s t o k e se q u a t i o n s ；a n d u s eo fs t o k e so p e r a t o rt oc o n t r o lt h eg r i d d i n gm o t i o n f o rp u r p o s e so fa n a l y s i sa n dc l a r i t yo fp r e s e n t a t i o n ，t h ea t t e n t i o ni sf o c u s e do nt h e l i n e a rs t a b i l i t yo ft h et h i nf i l mf l o wd o w na ni n c l i n e dp l a n et h a ti sg e o m e t r i c a l l yt h e s i m p l e s tf r e es u r f a c ep r o b l e m i nt h en u m e r i c a lt e s t sw ed e v e l o pa ne f f i c i e n tm e t h o db a s e d o ns e m i - i m p l i c i ts c h e m ei nt i m ea n dl e g e n d r es p e c t r a le l e m e n tm e t h o di ns p a c e ；s t o k e s e q u a t i o ni sa p p l i e dt os o l v et h eg r i d d i n gv e l o c i t y t h en u m e r i c a lr e s u l t so b t a i n e da r e n o to n l yi nc l o s ea g r e e m e n tw i t ht h er e s u l t so ft h el i n e a r i z e ds t a b i l i t yo fo r r - s o 】m e r f e l d e q u a t i o n ，b u ta l s os h o wt h a ts t o k e so p e r a t o rb e t t e re n s u r e st h ec o n s e r v a t i o no ft h e v o l u m eo ft h es p e c t r a le l e m e n t st h a nt h eu s u a le l a s t i co p e r a t o r t h i sp r o v e st h a tt h e a l es e mi sag o o dm e t h o dt ot h ef r e e - s u r f a c ef l o wp r o b l e m s k e yw o r d s ：m o v i n gm e s h ；n a v i e r - s t o k e se q u a t i o n s ；l e g e n d r es p e c t r a le l e m e n tm e t h o d ； a l em e t h o d ；l i n e a rs t a b i l i t y 厦门大学学位论文原创性声明 兹呈交的学位论文，是本人在导师指导下独立完成的研究成果。本人 在论文写作中参考的其他个人或集体的研究成果，均在文中以明确方式标 明。本人依法享有和承担由此论文产生的权利和责任。 声明人( 签名) ： 用寥 砌洚占月7 日 厦门大学学位论文著作权使用声明 本人完全了解厦门大学有关保留、使用学位论文的规定。厦门大学有 权保留并向国家主管部门或其指定机构送交论文的纸质版和电子版，有权 将学位论文用于非赢利目的的少量复制并允许论文进入学校图书馆被查 阅，有权将学位论文的内容编入有关数据库进行检索，有权将学位论文的 标题和摘要汇编出版。保密的学位论文在解密后适用本规定。 本学位论文属于 1 、保密 () ，在年解密后适用本授权书。 2 、不保密( , j 作者签名：闽容 日期：聊万年b , 9 导9 币签名：做匠日期： 孵阳 夕e 1 7 e l 第一章引言 第一章引言 1 1 带自由液面的流体运动 具有自由液面的流体流动问题在海洋工程，石油化工，水利工程和机械工程中大量 存在。自由液面可以是固液交面，液液交面或气液交面，它随着时间在移动变化。这种 移动边界的存在使得流体的运动区域不再是固定的，而是变形的，如何有效地跟踪模拟 自由液面随时间的变化过程正是研究此类问题的困难所在。带自由液面流体运动的基 本方程为n a v i e r - s t o k e s 方程，它的解析解在实际中极难得到，所以实验和数值模拟成为 研究它的主要手段。特别是随着计算技术的发展，数值计算的优势愈发突显，已经成为 研究流体问题的重要手段。过去人们采用了多种数值方法来研究带自由液面的流体问 题 郴】，主要有差分法，有限元法和边界元法等，并且取得很大的成果。本文主要讨论 用a l e ( a r b i t r a r y l a g r a n g e - e u l e r ) 谱元法做数值分析。a l e 方法结合t l a g r a n g e 方法 和e u l e r 法的优点，既可以精确描述自由液面随时间的变化过程，又具有解决大幅晃动和 高雷诺数的特点。谱元法是有限元法和谱方法的结合，它既继承了有限元方法对复杂区 域的适应性i 又保持了谱方法的高精度等优点。 为了分析方便，本文将主要研究沿平板下落薄膜流动问题。薄膜流动就是一层厚度 很小的流体在固体衬体上的流动，它是一种典型的交界面流动。从几何角度看薄膜流动 是一种最简单的自由液面流动，但却有着非常复杂丰富的动力学行为，是一种很典型的 自由液面流体。对它的研究有助于对自由液面流体流动的理解，具有重要的现实意义。 本文将a l e 与谱元法结合数值求解二维状态下沿平板下落薄膜流动线性扰动问题，得到 了与实验结果相吻合的数值结果，从而证明了谱元法求解带自由液面流体问题的有效 性。 1 2a l e 方法 在对自由液面的流体流动问题的数值研究中，最关键的问题就是自由液面的处理。 为了提高自由液面的模拟精度，人们研究了多种自由液面模拟方法。传统的方法主要 有e u l e r 法和l a g r a n g e 法。l a g r a n g e 方法可以精确描述自由液面随时间的变化过程，而且 第一章引言 求解过程不存在对流项，可简化数值求解过程，但是流体流动导致的网格过分扭曲会增 加求解误差或使求解失败。因此l a g r a n g e 方法一般只适合于小幅晃动，低雷诺数和无漩 涡问题。而e u l e r 方法不易模拟自由液面，但可以用于解决大幅晃动和高雷诺数流体问 题。所以就产生了将两种方法结合的a l e ( a r b i t r a r y - l a g r a n g e - e u l e r ) 方法。 a l e 方法最早出现在有限差分中。1 9 6 4 年n o h t l 首先提出了a l e 描述的基本思 想并通过有限差分程序实现，用于求解带移动边界的二维流体动力学问题。后 来d o n e a 和b e l y s c h k o s ，9 1 将a l e 方法和有限元结合，求解了瞬态流固耦合问题，较好的 解决了液面的大幅晃动和边界的大幅运动等问题。文献f 1 0 】利用a l e 有限元法得到了在 任意形状三维容器中非线性自由表面振荡的数值解。文献【1 1 】贝n 用a l e 有限元法和分步 求解结合研究了三维流体晃动和非线性流固耦合问题。至今人们用a l e 方法和有限元结 合求解了很多带自由液面的流体流动问题。 1 3 研究动机 液态薄膜厚度通常只有0 o l i m m ，而它的流向尺度却有i m 的量级，这样一个长度 比的流动不仅对实验测量是一个挑战，还使直接数值模拟变得非常困难。薄膜流动具有 空间上的不规则性和时阃上缓慢变化的特征。有限元是薄膜流动直接数值模拟的主要方 法b a c h l l l 。k h e s h g i l 2 1 ，和m a l a 玎1 a t 蕊【3 l 最先采用拉格朗日有限元法，并利用动边界对物 理变形有限元区域进行重新划分s a l a m o n l l 2 i 的有限元方法将控制方程建立在表面波 速的惯性坐标下，采用了同样的方式处理薄膜的运动边界问题。c a i r n c r o s s t x3 】利用有限 元法求解了三维情况下带自由液面的粘性不可压流体流动问题。 本文旨在证明谱元法对此类问题同样有效，并且可以发挥高精度的优点。谱元法 p a t e r a 1 4 等人发展起来，它是谱方法和有限元方法的结合，其基本思想是将复杂区域 分解成若干个较小的四边形或六面体区域，然后在每个小区域上使用谱方法，得到单元 矩阵，最后利用有限元法的思想，将这些单元矩阵合成为总体矩阵，进而求得问题的近 似解。这样，谱元法既继承了有限元方法对复杂区域的适应性，又保持了谱方法的高精 度等优点。近年来，谱元法获得了很大的发展，用于求解不可压的n a v i e r - s t o k e s 方程方 面已有不少的工作。并已成功地解决了一些流动及对流换热问题。 h o f l 5 曾经采用l e g e n d r e 谱元法，结合任意的拉格朗日欧拉方程形式处理薄膜 流动问题。h o 认为他的数值模拟结果与o r r s o m m e r f e l d 线性稳定性结果吻合，并且 恭i k a p i t z a 1 6 1 在实验中得到的波形和波速进行了比较。本文与h o 的工作涉及相同的问 2 第一章引言 题，但在如下方面做了改进。首先本文采用s t o k e s 算子控制网格运动，能有效的控制网 格变形，而h o 用的二阶椭圆弹性算子没有这个性质。其次是对u z a w a 算法的改进，大大加 快了压力求解速度。最后一节中的数值结果与h o 的结果一致，并且将网格的两种控制方 程对谱元体积的影响作了比较，发现采用s t o k e s 可以较好的保证计算体积守恒。 本文的结构如下：第二章主要介绍沿平板下落薄膜流动问题的控制方程n a v i e r - s t o k e s 方程和它的变分过程，以及移动网格的实现。第三章介绍在空间方向上的谱元离 散及时间方向上的离散方法。第四章先介绍沿平板下滑薄膜流动o r r - s o m m e r f e l d 线性 稳定性理论，然后利用我们提出的方法对薄膜流动进行模拟。我们将模拟得到的结果 与o r r - s o m m e r f e l d 线性稳定性的结果进行比较，验证方法的有效性。最后考察了网格的 两种控制方程对模拟结果的影响。 3 第二章沿平板下落薄膜流动 第二章沿平板下落薄膜流动 2 1 薄膜流动问题 本文考虑二维情况下薄膜流动问题。设平板与平面的夹角为p ，流体密度为p 。假 设在流体流动方向上( 。) 液体薄膜长度是无界的，在交面方向上( ) 薄膜的平均厚度 为h 。薄膜流动的控制方程为n a v i e r - s t o k e s 方程： 卜酉o n + 肚v u = v 7 - + f i n ) ， i v u = 0 ，i na c t ) ， 其中7 - 为应力张量，u = ( u l ，t 2 ) r ，他为戤方向上的流体速度，f 为体积力，q ( t ) 为与时间 有关的流体区域。 对于牛顿流而言，它的应力张量下可表示成 r = 一p i + p ( v u + v u t ) 其中p 为压力，p 为粘性系数。以 为长度尺度，为速度尺度，其中= g h 2 2 p 为重 力加速度) ，h u o 为对n r 度，p 瑶为压强尺度。通过无量纲化过程，可得 i 瓦o n + u v u 却叮+ f 枷， ( 2 - 1 ) 【 v u = 0 ，i nq ( t ) 此时的应力张量r 为： r = - p i + 壶( v u + v u r ) ( 2 - 2 ) 兵中r e = u o h p # 为雷话数 区域边界由以下三部分构成；o f t ( t ) = o a , u o f t q u a q 。( t ) ， 1 在弛上设周期边界条件 t ( o 删= 啦( 等删， p ( o 删叫等删 其中q = 警，q 为波数，入为波长，i = 1 ，2 。 4 ( 2 - 3 a ) ( 2 3 b ) 第二章沿平板下落薄膜流动 2 在a f 乙上设无滑移边界条件 讹= 0 3 在自由液面a q 仃( ) 上设表面张量牵引边界条件 s t 【叫+ 去( v u + v u t ) 】n = o ， n t 一p i + - 去( v u + v u t ) 】n = k ， ( 2 - 4 ) ( 2 - 5 a ) ( 2 5 b ) 其中n = ( 礼l ，n 2 ) t 和s = ( s l ，s 2 ) r 分别为a q 矿( t ) 上的单位法向量和切向量，k 为曲率， 为w e b e r 数，w = , , h u g ，仃为表面张力。在二维情况下，根据7 1 1 8 1 + n 2 5 2 = 0 ，可 将( 2 5 a ) 和( 2 5 b ) 简化成 卜p i + 击( v u + v u t ) i n = w a n ， ( 2 6 ) 总的来说薄膜流动问题的特征由o t ，p ，r e ，w 这四个无单位的参数决定。注意在以下篇幅 中为方便我们将去记为，j ，w i 己为c r 2 。2 变分方程 将动量方程和连续性方程分别乘以地和口，在q ( ) 区域上做积分，并且对动量方程做 分部积分，可得到如下的等价变分形式：找( u ，p ) ，u 础，札【q ( t ) 】2 ，且p 己2 m ( t ) 】，使 得 小象奴+ 小( u v u 肛厶v 刊v u + v u r ) 一，奴 + v p l u ( v u + v u ? ) 】n 打一 f v d x = 0 ， v v - l 。1 ，觚【q ( 亡) 】2 ， ( 2 - 7 ) ，8 瓤t ) j f t ( t ) g v u d x = 0v q l 2 【q ( 亡) 】( 2 - 8 ) d r y ( t ) 其中硪，a 【q ( t ) 】和l 2 【q ( t ) 】为q ( t ) 上的s o b o l e v 空间，l 2 【q ( ) 】为在q ( t ) 上平方可积的 空间，凰，a 【q ( t ) 】表示在q ( 亡) 上为硪【q ( t ) 】空间且在边界上满足方程( 2 4 ) 将( 2 7 ) 具 5 第二章沿平板下落薄膜流动 体展开来如下； 厶m ( 鲁+ - 丽o u l + 砌面o u l , ，, 一p 鲁+ 2 p 瓦o v l 瓦。u l + p 哿( 哿+ 等) - 嘶 叔 一厶卜印，+ 2 p 磐他- 邯c 等+ 尝胁，研打= 。， 厶州警+ u 警+ 抛哿卜p 篱+ 2 p 等等+ 丝o xr 丝o x + 等，一忱f 2 d x 一厶( t ) 【哪z + 2 p 万o u 2 n 2 + p 【瓦o u 2 - o 吲u 1 ) n 忱打扎 q ( t ) 的边界由三部分组成， v 卜p i + p ( v u + v u t ) 】n 如= v 【一p i + p ( v u + v u t ) 】n 打 j a n ( t )，a f b + v 卜一+ 弘( v u + v u r ) n d o r + v 卜- p i + z ( v u + v u t ) 】n 出 _ ，a n i j o f z ，( t ) 现在将上述的边晃条件用于边界积分项 1 在a q 缸上v = 0 ，因此厶v 【一皿+ t z ( v u + v u r ) n d a = 0 2 在a 上根据边界条件( 2 - 3 a ) 和( 2 - 3 b ) 得kv 【一p i + p ( v u + v u ? ) 】n 打= 0 3 在a ( t ) 上根据边界条件( 2 5 a ) 和( 2 5 b ) 得k ( t ) v - p i + 卫( v u + v u t ) 】n 如= k ( t ) o c v n 打 那么可得到： v 卜p i + p ( v u + v u ? ) 】n 如= 矿胛n 如 j a a ( t )j a n ， 在对上述的变分形式进行全离散之前需要对时间项( 方程( 2 7 ) 中的第一项) 做离 散。 警g 筚i l 2 瓦譬矿4 2lz 其中上标+ 1 ) 和n 分别表示当前步和上一个时阔步。如果将( 2 - 7 ) 直接做时阒离散的 话，将会出现上一个时间步的流体速度婶在当前步的网格上做积分，如果网格是固定的 这不会有问题，但对于移动网格既当网格速度不为。的时候，积分中的地基 会引起混 6 第二章沿平板下落薄膜流动 乱，逐柙个一毁任目j 鸵会严笙很大网呙教误爱为1 j 避毙返一臂现，爵蛩掰网。t 网俩寻 移到积分号外面，以下就是这一过程。首先对( 2 - 7 ) 的第一项做分部积分得到( 1 = 1 ，2 ) ： 知甏奴= 幻警地面o ? g i - 胁( 2 - 9 ， 下面我们先介绍一下r e n o l d s 输运定理。 定理2 1 ：对于任意使下列各项积分有意义的妒。成立 旦d t 厂( 0 妒( x 奴一上( 。害叔+ z ( t ) w ，v 妒+ 妒( v ，w ) 奴 ( 2 - 1 。) 证明：这里w = 鲁，记= o 时的流体质点坐标为x o ，ljl 为j a c o b i 行列式 旦d t 厂( 0 肼) 奴= 丢厶脚i 奴0 = 上。蕊d ( 出 t ) i ji ) d x 。 = 厶( 等+ 妒掣 因为髻( x ，t ) = 甓+ w v 妒，掣( x ，t ) = ljl ( v w ) ，将它们代入上式就可得到 丢z 妒( 即) 叔= 毒 篱l 7l + ( w v 纠 ，l + 妒( 阢w ) f 歹| 蕊 = 幻警+ w v 妒州玑w ) ) d x 对( 2 - 7 ) 式右边第一项用上述的r e y n o l d s 输运定理司得 厶警奴；象厶。姚奴一鼬警+ 警溉( 2 - n ， 其中w = ( w l ，毗) t 为网格速度。注意到本文中毗和是独立无关的变量，所以可通过选 择控制挑( i = 1 ，2 ) 的方程避免计算区域内网格产生大的畸变，从而保证计算的可靠性， 这也是a l e g 弓- 法的优点之一。定义仇随流体质点( 记为x ) 运动的时间导数为0 ，即 警l x _ 警+ 叫1 5 挑i + i u 2 筹_ o 、面l x 嚣否f + 叫+ 否歹2 u 将上式代入( 2 1 1 ) 式可得到 z ( 0o u i ，奴= 丢z ( t ) v t u i 叙一f n o u i w l - + 是笋，拿 陪1 2 ， 7 第二章沿平板下落薄膜流动 将( 2 - 1 2 ) 代入( 2 - 7 ) 得到：找( u ，p ) ，u 硪搠。f q ( ) 】2 ，r p l 2 1 2 ( t ) ，使得 丢上v u 叔+ z v f ( u w ) v u - u ( v w ) 】a b c 一上。v ，v p ( v u + v 一) v v ) 叔 一z 仃腓n 打一上f v d x = 。， ， k 书u 奴= 0 ， 在二维情况下曲率代可表示成 摊。8 2 ， s 1 3 2 8 1 ， v v 硪8 n 。【q ( 卯；( 2 - 1 3 ) v 口驴【q ( ) 】 ( 2 - 1 4 ) 其中s = ( s l ，s 2 ) t 为a q 口( t ) 上的单位切向量，f 为自由液面丽上的曲线坐标，口，6 为曲线 的两端点，8 i ，= a s i o f , 。 下面将( 2 - 1 3 ) 中的曲面积分做变换如下 pp a a v l n l d a 。 a ( s 2 ，s l 一8 2 8 1 。) o l n l d o ja o 旧 ，b ( 钟 一 a ( s x ，e s ；一8 1 8 2 8 2 。) v l d a j 笋k ( ”广 。o s t , e v l d o 一0 8 1 ( 8 1 8 1 ，+ 3 2 8 2 ， ) 打 j 讲k ( 1 )j 俄k ( t ) 因为s l s l ， + 8 2 8 2 ，= i 1 硪0 a 6 2 l + s ；) = 壤( 1 ) = 0 ，所以 l 。( t ) o k u l l q l 打= 卜讯加， 2 对( 2 1 5 ) 的右边项做分部积分并应用u l 毒= 秒l ，i 8 i ，得到曲线积分的最终形式 j = l 类似的有： z a , , ( t ) a t c v l n l 如= 一卜z ( 磐s ，+ 万o r l s 2 ) 曲巾郴1 ) l l , , ( oa r g v 2 n 2 出一上b 口蹦, 面d r 2 s + 等s 2 ) 如+ ( 盯忱s z ) l ： 注意( 2 1 7 ) 只适合二维的情况，三维的情况f 1 7 1 要复杂的多。 8 ( 2 - 1 5 ) ( 2 - 1 6 ) ( 2 - 1 7 ) 第二章沿平板下落薄膜流动 2 。3 移动网格 在a l e 谱元法中，自由液面的跟踪实质上是要保证自由液面上的网格点始终要与液 面保持一致，即要保证流体质点不会越过流体与空气的交界面。在二维情况下这就要求 在讹矿( ) 的法向方向上，网格速度戤需满足 t l ，1 l + t 屹，b = ”1 佗】+ u 2 n 2 ( 2 1 8 ) 有了上面的条件我们就可以利用流体的速度来确定自由液面上的网格速度。在自由液面 切向方向上，为了保证网格变形尽量小，可令 i 3 1 3 1 + w 2 3 22 0 解( 2 1 8 ) 和( 2 - 1 9 ) 方程组得到w 在自由液面上的计算公式； ( 2 1 9 ) 铷l = ( u l n l + u 2 n 2 ) n l ，w 2 = ( u l n l + 砌扎2 ) 纷2 ( 2 - 2 0 ) 在其它与时间无关的边界上，我们令 w = ( w 1 ，地) t = 0 ( 2 2 1 ) 一旦网格速度毗( = l ，2 ) 在整个边界上被确定下来，就可设法将它扩展到内部区 域。h o 1 5 】提出一种简单的方法，将内部区域的网格速度定义为满足边界条件( 2 2 0 ) ( 2 - 2 1 ) 的椭圆型方程的解： v ( v w ) ；0 ， i n q ( t ) ( 2 - 2 2 ) 由( 2 - 2 0 ) ，( 2 2 1 ) ，( 2 2 2 ) 确定网格变化的一个主要缺点是无法有效控制网格变形，于 是本文提出用s t o k e s 算子来定义网格速度方程： 这里张量于定义为： i n q ( t ) ， i n q ( 亡) 亍= 一皿+ p ( v w + v w t ) ( 2 2 3 ) ( 2 2 4 ) 其中庐和豇分别为网格压力和网格运动粘度。将q ( t ) 中的流体质点的位置用l a g r a n g i a n 坐 9 o 0 = l i 吁 w v v 、-j l l 第二章沿平板下落薄膜流动 标表不，记为x = ( z l ，x 2 ) 。，我们有以f 的运动条件： 要：w ( 2 - 2 5 )一= w ，一z nl 出 r 7 那么考虑的计算区域随时间的变化过程可通过上式来刻画。 那么最终得到的n a v i e r s t o k e s ：y 的变分形式为：a ( u ，p ) ，u 础，a 吼【q ( t ) 】2 ， 勘l 2 f q ( ) 】，使得 爰上泐v + u d x + 上绀v 【( u w ) v u - u ( v w ) 】d x 一上辩( p v - v - t ( v u + v u t ) v v ) d x + 。水) 仃s ( s v v ) 如一吖sl ( t ) 一上( t ) f v d x = 。， v v 瑶a n 。瞅圳2 ；( 2 - 2 6 ) 厶一u 叔她 v qel 2 a ( 纠( z - 2 7 ) 1 0 第三章数值方法 第三章数值方法 3 1 l e g e n d r e 谱元法 假设q 可以表示成髟个小区域q 知的并，即q = u 鉴l 霓南。把每个小区域称为一个单 元，这些单元需满足以下条件： 1 不同的单元之间不能有重叠，即q in 踢= 仍( i j ) ； 2 二维情况下如果r = 盈n 锄( i 歹) ，则r 要么是一个空集，要么是单元上的一条完 整的边。 本文采用谱元法进行空间离散，第二节得到的变分形式是空间离散的基础。 定义q = a 天为标准单元，其中a = ( 一l ，1 ) 为一维标准区间。假设r 是从q 到m 的映射。标准单元里的坐标用( ns ) 表示，那么有： r ：q ( 7 ，s ) 爿q 南( 。，) ， 也即：z = 硪( r ，s ) ，秽= 露( n s ) 。对于带自由液面的流体流动问题，因为要考虑到网格 速度，所以区域单元随时间演变后就不再是标准的四边形，而是平面上任意的单元。 定义在q 上的分片多项式空间如下： n , k ( f 1 ) = 圣五2 ( q ) ；西。凡( r ，8 ) 尸( q ) ，1 七k 助一2 ( q ) = 圣l 2 ( q ) ；圣or ( r ，s ) p 一2 ( q ) ，1 南) p ( q ) 表示q 上关于变量n s 的阶数均不超过的多项式空间。需要指出的是：虽然 称尸( q ) 和尸一2 ( q ) 为分片多项式空间，但由于r 的作用，它们中的函数限制在q 七上不 一定是多项式。 现在定义速度和压力的谱元近似空间分别如下； = ( 碓a n ，( t ) ( q ) np ，k ( q ) ) 2 ， 知= 舻( q ) n 霸_ 2 k ( q ) 1 1 ( 3 - 1 a ) ( 3 - 1 b ) 第三章数值方法 我们注意到对速度空间用阶多项式，而对压力用的是n 一2 阶多项式，这就是传统 的p p 一2 方法。在高阶的谱元法中压力空间的维数比速度空间低两阶，为了使离散 问题的解存在唯一，这是谱元法鞍点理论中的l b b 条件f 1 8 j 。 变分方程可写为：找u n x n ，p n m n ，使满足 丢( u n , v n ) + c ( u ，w ，v ) 一6 ( v ，p ) + 。( u n ，v n ) = ( f 知，v ) + d ( s n , v n ) ，v v j ，； 其中 b ( u ，q n ) = 0 ，v 如m n ( 3 2 a ) ( 3 - 2 b ) c ( u n ，w n ，v n ) = v t 【( u n w ) v u n u n ( v w n ) d x ，( 3 - 3 a ) b ( v n ，p l y ) = p n v v n d x ，( 3 - 3 b ) a ( u n ，v n ) = p ( v u n + v a t ) v v n d x ， ( 3 - 3 c ) d ( s n ，v n ) = 一 o s n ( 8 n v v n ) d o + o w n s l a n ，( t ) ( 3 - 3 d ) 下面定义g l l ( g a u s s - l o b a t t o - l e g e n d r e ) 和g l ( g a u s w l e g e n d r e ) 积分公式： ( 1 ) 二维标准单元q 上的g l l 积分公式为： 厶认邵) d r d s 丕烈6 舭( 3 - 4 ) 其中& ，分别为一维标准区间a 上的g l l 积分节点和权值。积分公式( 3 4 ) 对恳n - 1 ( 盎) 中 的多项式精确成立记( 砖，谚) = 最( & ，白) ，( 。各，蝣) 被称作单元呸上的局 部g l l 点，i ，歹一般称作) 粥g l l 点的局部编号。q 知上的g l l 积分公式定义为； 小删咖= 五( 加脚 s ) 删s 叠( 懈确， 其中以为映射f k 的j a c o b i 矩阵，定义如下s j = 、l， 如一c；曲一c； 圪 如一加以一务 ，i一一、 ，的逆矩阵为： j _ 1 耋、1 ；高f 一耋一耋 o y o r如 由链式法则在q 七单元上对z ，的导数可以用标准区域中的坐标7 ，s 表示为： 88 r88 s8 瓦。瓦一a r 十蕊一o s ，8 笺8 茁l8 笔 那么对v v 有： 88 i r88 88 万2 万爵十瓦蕊8 馨8 笥8 f j 8 譬8 s ( v v n ) of 七= j - t v t ( v of 七) ， ( v v ) of 丘= j 一1 v ，( v o 严) ， ( 2 ) 二维标准单元佥上的g l 积分公式为： 么地s ) d r d s 一1 慨，白) 瓯奶， ( 3 - 5 ) 其中己，觑分别为一维标准区间人上的g l 积分节点和权值。记( 孟易，踮) = 最( 磊，弓) ，( 磅，姥) 被称作单元q 膏上的局部g 厶点，那么上的g 三积分公式可定义为： z ( z ，暂) 如匆= z ( | 五i ) 。冗( ns ) 办幽 一1 ( 1 五 ) ( 镑，鲢) 觑岛 i j = l 那么( 3 - 2 a ) ，( 3 - 2 b ) 式的数值积分形式为：找t i n ，p n m n ，使满足 忌= l【爰( u ，v ) + c ( u ，w n v n ) f 一6 ( v 舢) ；：+ 。( u ，v ) 期= ( 知，v ) ：+ k f f i l七= l d ( s n ，v n ) ， v v n x n ； 耳 b ( u n ，卵) 赛= 0 ， k = l 1 3 v g 五知 ( 3 - 6 ) ( 孓7 ) 阶一如夙一如 归 第三章数值方法 其中 ( u ，v ) = ( v 。最) ) ( u 。最) ( 勖) i 五渤) l ， i d = o e ( u ，w ，v ) = ( v 。最) ( 如) 【( ( u 一w n ) 。疋) ( 石r ( u n 。最) ) i j = 0 一( u o 最) ( 石1 v ，( w o 最) ) ) ( 岛) i 以( 岛) j 蚴， 一l 6 ( v ，p r ) ： = 。忍) ( 白) ( 石1 v r ( v 。忍) ) ( 毛) i 以( 岛) i ， i d = l n a ( u 。v ，v ) 盏= p ( 石r v r ( t i n 。r ) ) + ( j r r ( v r ( u 。最) ) ) 丁) t j = 0 ( j f r v r ( v n or ) ) ( 岛) l 以( 岛) l ， ( f ，v ) = ( v ，o 最) ( 钿) ( 轴。风) ( 岛) l 厶( 锄) i ， i j = o d ( s - v ，v ) 惫一一o ( s n 。最) ( s 。r ) ( 石t v ，( v 。凡) ) ) ( 6 ) l 砖( & ) i 蛾， i = 0 一l 6 ( u ，q ) 2 = ( 。r ) ( 岛) ( 石1 v r ( u 。疋) ) ( 白) l ( 己) i 为了简单起见，假设自由液面的边映射到标准区域对应s = l 的矩形边，那么这条边上 的g l l 结点为 ，“) ， = 0 ，n 。此外为表面j a u c o b i ，定义如下： 七( “) = 略纠】2 + 【骞( “) 2 标准单元q 里对应于g 己乞点( ，岛) 的l a g r a n g e 基函数可记作：h i i ( r ，s ) = ( r ) ( s ) ，其 中也( r ) 为维标准区间天上的g l l 点6 对应的l a g r a n g e 基函数。面g l 点磊对应的的基函 数记为：( f ，) = 风( 严) ( ) 。令巧1 为凡的逆映射，并记坞( z ，y ) = h v ko 巧1 ( z ，矽) ，类 似可定义 嚣( z ，耖) 。 各( z ，可) 被称作局部基，它是局部g l l 点( $ 乞，掣刍) 对应的l a g r a n g e 插 值多项式；k ( f ，；) 是局部g l 点( 磕，诒) 对应的l a g r m n g e 插值多项式。由于巧1 的作 用， 各( z ，可) 不具有b 那样的变量分离形式，一般也不是q 七上的多项式。与一维情 况类似，对v v x ，有： v ( x ，掣) h = 吨磕( 删) ( 孓8 ) i j = o 1 4 第三章数值方法 定义点集： 乃= ( z ，耖) i ( 。，可) = ( 。，赡) ，1 忌k ，0 冬i ，歹n且；( z ，y ) 罾铀) ， 死= ( z ，) l ( z ，秒) = ( 孟，鲢) ，1s 七k ，1 ，jsn 一1 ) ， 它们分别是所有内部g l l 点和g l 点的的集合。基于这些点可以构造和a 知的一组全 局节点基，对应于点( 姒，y 1 ) a 的基函数幽( z ，耖) 可以表示成： 讥c z ，材，l n 。= ：，z y ) i ，( 统) = ( z 各，药) ， e l s e 对应于点( 磊，巍) 死的基函数磊( z ，毛) 可以表示成： 西c $ ，i 。= 等z 可) l ，( 戤，轨) = ( 诌，筠) ， 对任意的v ，它可以用g l l 或g l 全局基表示： v ( x ，v ) = v ( x l ，y 1 ) 妒l ( x ，可) ( ，认) 死 ( 3 - 9 ) v ( 刚) = v ( x z ，玑) 磊( z ，) ( 3 4 0 ) 仁，m ) e t 2 那么u n ，v n ，轴可以和( 3 - 9 ) 一样用g l l 的全局基展开，p 用g l 点的全局基展开。 现在分别定义g l l 节点基和g l 节点基的导数算子如下【1 9 】； ；掣，d o - 掣，( 3 - 1 1 ) 记勉( 邑) = b ，l 以( 锄) i = lji 各，l ( 己) i l 歹l 各；对v v = ( m ，吨) t x n ，记v 各= v or ( 岛) ，v - 巧k ；v o 最( 玉) ，那么有 ( 町1 v ( v n 。r ) ) ( 岛) 罩【( g - ) 每( u - ) + ( g 2 ) 易m 渤) 1 ，( 3 - 1 2 ) m n ：= 0 或者 n - 1 ( 石1 v ( v 。最) ) ( 己) = f ( 台z ) 各舢( 魂) + ( 岛) n 渤) 袅n 】，( 3 - 1 3 ) 仇。扎= l 1 5 其中： ( g t ) 岛册= ( g 2 ) 各肭= ( 0 t ) 嚣撇= ( 0 2 ) 易m n = n p = 0 p = 0 n l f = i 一l p = l 丽1 ( 。拥绣一d j n 绣) ， ( 3 - 1 4 a ) 赢( 一和易+ d i p d j n 略) ，( 3 - 1 4 b ) 方甫西概锄砖n 站一碥西印岛n 殇) ， ( 3 - 1 4 c ) 瓜1 ( 一d t m 如砖n 矛暴+ k 蕊岛n 矛轨 n i ，i 嚣= ( d p 名易b q 磕一岛口z 岛绉) ， p ，q = o 一l 歹晦= ( 瓯绉岛。皱一岛。童乞西劬踢) p ，q = l ( 3 - 1 4 d ) ( 孓1 4 e ) ( 3 - 1 4 f ) 记u 各= u n o 兄( 岛) ，w 嚣= w