（流体力学专业论文）基于高阶边界元的三维数值波浪港池.pdf

上海交通大学 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定， 同意学校保留并向国家有关部门或机构送交论文的复印件和电子 版，允许论文被查阅和借阅。本人授权上海交通大学可以将本学位 论文的全部或部分内容编入有关数据库进行检索，可以采用影印、 缩印或扫描等复制手段保存和汇编本学位论文。 保密口，在一年解密后适用本授权书。 本学位论文属于 不保密口。 ( 请在以上方框内打“4 ”) 学位论文作者签名： 指导教师签名： 日期：年月e t 日期：年月日 上海交通大学 学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指导下， 独立进行研究工作所取得的成果。除文中已经注明引用的内容外，本 论文不包含任何其他个人或集体已经发表或撰写过的作品成果。对本 文的研究做出重要贡献的个人和集体，均已在文中以踞确方式标明。 本人完全意识到本声明的法律结果由本人承担。 学位论文作者签名： 日期：年月日 基于高阶边界元的三维数值波浪港池 摘要 本文初步建立了一个基于高阶边界元的三维数值波浪港池，它可 以用来模拟完全非线性波浪的运动。本文的数值模型使用高阶三维边 界元方法和可调节时间步长的基于二阶显式泰勒展开的混合欧拉。拉格 郎日时间步进来求解带自由表面的完全非线性势流方程。在港池端 可以造出非线性的周期性波浪，在另一端为反射边界或通过人工粘性 消除反射波的吸收边界。f 在所使用的边界元方法中，边界几何形状和 流场变量采用g r i l l i ( 2 0 0 1 ) 所提出的1 6 节点三次滑移四边形单元来表 示，这种单元在单元内具有高阶的精度同时在单元之间具有良好的连 续性。对于这些单元的数值积分进行了有效的处理以保证具有较高的 精度。使用速度势连续性条件妥善处理了对数值解的稳定性有重要影 响的自由表面和边壁交界处的边界条件。在局部曲线坐标系中采用2 5 节点划移四次( g r i l l i ( 2 0 0 1 ) ) 单元计算县赢高阶精度的空间导数以用于 计算自由表面上节点的时间步进。j 模型取得了对于体积和能量守恒的 ， 相当高的计算精度，在计算过程中无需数值光滑和网格重划技术而没 有出现锯齿状不稳定性。本文通过计算二维波浪问题对该数值波浪港 池进行了充分的校核。f 首先给定初始的孤立波波形和其上的速度势， 计算了沿常水深的孤立波的传播。然后在港池的一端模拟推板造波板 造波，分别计算了另一端为反射边界时驻波的形成和另一端加上消波 段时均匀波列的形成。计算结果表明，该数值波浪港池可以有效地模 拟孤立波的传播和周期性非线性波浪的运动，具有良好的数值精度和 第1 页 摘要 随着时间空间离散细化的收敛性。港池的消波段取得了比前人结果更 好的消波性能，说明该港池具有实际应用前景。最后本文计算了一个 三维问题，初步验证了本港池可以用于模拟三维问题。，、1 关键词边界元，数值波浪港池，孤立波，造波，消波 一一_ _ _ _ h _ 一 一一 一 第2 页 摘要上海交通大学硕士论文 a3 dn u m e r l c a lw a v et a n kb a s e do n hig h e r 一0 r d e rb e m a b s t r a c t at h r e ed i m e n s i o n a ln u m e r i c a lw a v et a n k f n w t ) b a s e d o nah i g h e r o r d e r b o u n d a r ye l e m e n tm e t h o d ( b e m ) i s d e v e l o p e d ，w h i c hc a nb eu s e d t os i m u l a t e f u l l y n o n l i n e a r w a v e s t h em o d e ls o l v e s f u l l y n o n l i n e a r p o t e n t i a lf l o we q u a t i o n sw i t haf r e es u r f a c eu s i n gah o b e ma n dam i x e d e u l e r i a n l a g r a n g i a n t i m e u p d a t i n g ，w h i c h i sb a s e do ns e c o n d o r d e r e x p l i c i tt a y l o rs e r i e se x p a n s i o n sw i t h a d a p t i v et i m es t e p s n o n l i n e a r p e r i o d i cw a v e sc a nb eg e n e r a t e da to n ee n do fn w t , a n dr e f l e c t i v eo r a b s o r b i n gb o u n d a r yc o n d i t i o n ss p e c i f i e do nt h eo t h e re n da n dl a t e r a l b o u n d a r i e s i nt h e b e m ，b o u n d a r yg e o m e t r ya n df i e l dv a r i a b l e sa r e r e p r e s e n t e db y1 6 。n o d ec u b i c s l i d i n g q u a d r i l a t e r a le l e m e n t sd e v e l o p e d b yg r i l l i ( 2 0 0 1 ) t h e s ee l e m e n t s p r o v i d eh i g h e r o r d e r a c c u r a c y i n t h e m s e l v e sa n dg o o di n t e r - e l e m e n t c o n t i n u i t y a c c u r a t e a n de f f i c i e n t n u m e r i c a l i n t e g r a t i o n s a r e d e v e l o p e df o rt h e m d i s c r e t i z e db o u n d a r y c o n d i t i o n sa ti n t e r s e c t i o n sb e t w e e nt h ef r e es u r f a c eo rt h e b o r o ma n d 第3 页 捅璺 j f _ j m 、丁o l a t e r a lb o u n d a r i e sa r ew e l l p o s e d i na l l c a s e so fm i x e d b o u n d a r y c o n d i t i o n s h i g h e r o r d e rt a n g e n t i a l d e r i v a t i v e s ，r e q u i r e d f o rt h et i m e u p d a t i n g ，a r ec a l c u l a t e di n al o c a lc u r v i l i n e a rc o - o r d i n a t es y s t e m ，u s i n g 2 5 一n o d e s l i d i n g f o u r t h - o r d e rq u a d r i l a t e r a le l e m e n t s ( g r i l l i ( 2 0 0 1 ) ) v e r y h i g ha c c u r a c yi sa c h i e v e di nt h em o d e lf o rm a s sa n de n e r g yc o n s e r v a t i o n n os m o o t h i n go ft h es o l u t i o ni s r e q u i r e d i nt h ec a l c u l a t i o nt oa v o i d s a w t o o t hi n s t a b i l i t y t h en w t i st e s t e dt h r o u g hc a l c u l a t i n g2 dw a v e s b ys p e c i f y i n gs o l i t a r yw a v ep r o f i l ea n di t s v e l o c i t yp o t e n t i a l ，s o l i t a r y w a v e p r o p a g a t i o no v e r c o n s t a n td e p t hi s s i m u l a t e d b yg e n e r a t i n gw a v e s a to n ee n do f n w t s t a n d i n gw a v e sa n ds t e a d yw a v ep r o p a g a t i o na r e s i m u l a t e dw i t h r e f l e c t i n g a n d a b s o r b i n gb o u n d a r y o nt h eo t h e re n d r e s p e c t i v e l y t h er e s u l t ss h o w t h a tt h e p r e s e n t n 耳c a r ls i m u l a t es o l i t a r y w a v ea n dp e r i o d i cn o n l i n e a rw a v ep r o p a g a t i o n e f f i c i e n t l y w i t h g o o d n u m e r i c a l a c c u r a c y a n d c o n v e r g e n c e w i t har e f i n e d s p a f i o - t e m p o r a l d i s c r e t i z a t i o n t h e s p o n g el a y e r i nt h en w ta c h i e v e sb e r e rw a v e a b s o r b i n gp r o p e r t yt h a ni np r e v i o u sa t t e m p t sr e p o r t e di nt h el i t e r a t u r e ， w h i c hm e a n st h en w t h a sa p r o m i s i n gp r o s p e c ti np r a c t i c a la p p l i c a t i o n f i n a l l y , at h r e ed i m e n s i o n a lc a s ei sc a l c u l a t e di nt h e3 d n w t , w h i c h s h o w st h ec u r r e n tn tc a nb eu s e d t os i m u l a t e3 dc a s e s k e yw o r d sb o u n d a r ye l e m e n tm e t h o d , n u m e r i c a l w a v et a n k ，s o l i t a r y $ 膏 第4 页 摘要 上海交通大学硕士论文 w a v e s ，w f l v eg e n e r a t i n g ，w a v ea b s o r b i n g 第5 页 第一章论文背景上海交通大学硕士论文 1 1 研究概况 第一章论文背景 在过去的二十多年里，发展更为精确、高效的数值方法来模拟强非线性自由 表面波浪直是海洋、海岸工程与科学里的一项挑战和任务。确实如此，因为波 浪动力学控制着这些领域里发生的大多数物理现象( 例如大气、海洋的相互作用， 波浪靠岸和破碎) 和所使用的工程方法( 如破碎波对海岸结构物的冲击力的计 算) 。 在处理破碎前的波浪时，不管是理论上还是数值上，到目前为止最为成功的 方法都是基于势流理论的，即忽略波浪流的粘性和有旋性。在这种情况下，控制 方程一连续性方程，是流场速度势函数的l a p l a c e 方程。这个线性偏微分方程可 以用边界积分方程进行有效求解。但是非线性项出现在运动学和动力学边界条件 上，不同的方法因处理这些非线性边界的不同而各有不同的精度，适用范围和数 值效率。 基于上述的大多数波浪理论和数值模型的一个传统做法就是定义所谓的“小 参数”( 例如波陡、水深与波长比) ，然后将自由表面的边界条件和几何形状的截 断序列展开式表达成这些小参数的表达式。这些方法在处理海岸地区的浅水和中 等水深问题时取得了相当好的结果。 但是当处理强非线性波浪时( 例如那些接近破碎或开始破碎的波浪) ，这些 基于小参数的方法便不再有效了。在这种情况下，我们必须直接求解控制方程的 原始形式，并且使用混合欧拉拉格郎日法来跟踪自由表面的运动。这个方法最 早是由l o n g u e t - h i g g i n s 和c o k e l e t ( 1 9 7 6 ) 7 i 入用来模拟二维波浪。以后发展了很 多基于此方法的数值解，它们大多数是二维情况下的。例如，d o m m e r m u t h 等 ( 1 9 8 8 ) 和s k y n e r ( 1 9 9 6 ) 对深水和中等水深下波浪翻卷地模拟，g r i l l i 等( 1 9 9 4 ， 1 9 9 7 ) 和l i 和r a i e h l e n ( 1 9 9 8 ) 对斜坡下波浪地爬坡与破碎地模拟。这些计算结 果取得了相当好的精度。值得一提的是，近年来随着计算机计算能力的日益强大， 流体体积法f v o f ) 得到了越来越多的使用。这些方法直接求解自由表面流动的 纂6 页 m n s 方程，可以用来模拟破碎后的波浪运动。但是它们需要耗费大量的计算时间 ( 特别是在三维的情况下) ，而且由于数值耗散导致波浪长距离传播时产生非真 实的波浪能量损失。因此，除非要模拟破碎后的波浪，非线性势流理论是比较好 的模型。 因而到目前为止大多数模拟波浪的研究都是基于非线性势流理论的。而且， 在大多数的模型里，l a p l a c e 方程都是通过高阶边界元方法( b e m ，见b r e b b i a 6 ) 来求解，这些方法或者基于g r e e n 函数，或者基于c a u c h y 积分。自由表面的跟 踪( 即混合拉格郎日欧拉法) 则使用时间步进方法( r u n g e k u t t a 法) 或预测一 矫正方法( a d a m s b a s h f o r t h - m o u l t o n ) ，或泰勒展开方法( d o l d 和p e r e ：g r i n e ( 1 9 8 6 ) ) 。早期的计算都是二维的，局限于常水深、空间周期性的波浪( 例如 l o n g u e t h i g g i n s ( 1 9 7 6 ) ，d o l d 和p e r i g r i n e ( 1 9 8 6 ) ) ；更进一步的二维模型可以处 理随机入射波和复杂的底部地形变化( 例如，g r i l l i ( 1 9 8 9 ) ，c o i n t e ( 1 9 9 0 ) ， o h y a m a 和n a d a o k a ( 1 9 9 1 ) ) 。大部分现在的模型也能处理实际的空间域，即数 值波浪港池( n w t ) 一入射波由造波端产生，然后在其它边界被反射、吸收或 辐射( 例如g r i l l i 和s u b r a m a n y a ( 1 9 9 6 ) ) 。 尽管在这个领域已经做了很多工作，但是我们仍需要发展出更为高效和稳定 的数值模型来满足科学研究和工程应用的需要。而且到目前为止，对于三维问题 的数值模拟进行得还比较少。这其中只有少数研究模拟一般的三维数值波浪水池 里的非线性波浪( 例如c e l e b i 等( 1 9 9 8 ) ，c r d l l i 等( 2 0 0 1 ) ) 。 要建立一个一般性的数值波浪港池，主要有几个重要问题需要解决： 1 如何有效地跟踪位置事先未知的自由表面： 2 如何提高计算精度，减少计算量，控制各种数值误差，使计算有良好的 稳定性和收敛性： 3 如何有效地模拟造波和消波。 第一个问题前面已经述及，即使用拉格郎日欧拉混合方法来跟踪自由表面。 在每个时间步里，首先在欧拉场下用边界元求解给定边界条件的l a p l a c e 方程， 然后通过非线性运动学和动力学边界条件对时间积分来更新自由表面上质点的 位置和边界值的大小。这个方法被广泛采用，被证明是行之有效的。 第二个问题是最为重要的，也比较难处理。数值上误差是由很多地方共同产 第7 页 第一章论文背景 上海交通大学硕士论文 生的。对于自由表面的数值离散不可避免要产生误差。边界元插值函数的选取十 分重要，要提高精度和稳定性，高阶的插值函数是必需的，而且单元之间要有好 的连续性。在最新的文献里，x u e 等( 2 0 0 1 ) 采用高阶二次边界元，g r i n 等( 2 0 0 1 ) 采用三次划移四边形单元，都取得了非常好的计算结果。自由表面上物理量的数 值误差在一定时间后会逐渐放大，导致自由表面出现锯齿状不稳定性，所以通常 每隔一定时间就采用网格重划技术( 例如d o m m e r m u t h 和y u e ( 1 9 8 7 ) ) 或数据 平滑技术( 例如x u e 等( 2 0 0 1 ) ) 来控制误差。 数值误差和不稳定出现的另一个重要原因来自于对几个边界交界处的处理， 即所谓的角点问题。对于这个问题，g r i l l i ( 1 9 9 0 ) 做了比较详细地讨论。对于 数值模拟来说，主要问题是如何处理好角点处所谓双节点问题，即两个具有相同 坐标，不同边界条件的点。比较有效的方法是利用速度势函数和速度的连续性条 件来处理这些交界点。另一个问题是角点附近积分精度的问题，特别是对于波浪 爬坡问题。交界处的两个边界距离十分接近，被积函数的变化十分大，用少数的 高斯点来积分难以反映实际的变化情况，造成比较大的误差。g r i l l i ( 1 9 9 0 ) 给 出了一个解决办法，能够有效地解决这类问题。 有效地造波和消波也是大家一直在研究的课题。这方面的文献很多，也提出 了许多的方法。用造波板造波的主要问题是开始时造波板的突然启动会产生解的 奇异性。为了消除奇异性，通常的做法是给造波板的运动乘上一个函数，保证在 计算的最初几个时间步内板的加速度保持比较小。另一个造波的方法是事先给定 初始的边界条件。对于消波，要建立一个具有一般性的数值水池，消波段应该能 对各种不同的来波有效地消波。在目前已有的消波方法中( 可参见r o m a t e ( 1 9 9 1 ) 的总结) ，使用人工粘性( 例如c e l e b i 等( 1 9 9 8 ) ) 可以比较有效地消去各种波长、 各种类型的来波。但是人工粘性对于消波段的长度有要求，要求其长度至少要大 于一个来波的波长。这对于波长较长的来波来说比较浪费计算时间。考虑到活塞 式消波边界( c l e m e n t ( 1 9 9 6 ) ) 可以方便有效地消去低频的来波，c l e m e n t ( 1 9 9 6 ) 将二者耦合起来使用以达到优化的效果。g r i u i 等( 1 9 9 7 ) 对其进行了进一步的 改善。o h y a m a ( 1 9 9 1 ) 则将消除长波入射波的s o m m e r f i e l d 辐射条件的近似表达 式和人工粘性偶合起来使用。他同时指出波浪遇到障碍物后的反射和绕射传播到 造波端时会受到造波的影响，为了可以分析非线性和随机波浪，他在造波端也加 第8 页 第一章论文背景上海交通大学硕士论文 上了消波条件。 1 2 本文工作和意义 本文初步建立了一个有效的、普遍的三维数值波浪港池并验证了该港池可以 用于来模拟和研究真实情况下波浪的运动变化规律。 本文采用混合拉格郎日欧拉法来跟踪自由表面，时间步进采用二阶显式泰勒 展开法。自由表面的边界单元离散本文采取g r i l l i ( 2 0 0 1 ) 的方法一1 6 点三次划 移四边形单元。这种边界单元的插值函数阶数高( 三阶) ，单元之间各变量的连 续性好，是十分理想的插值函数。用这个方法，在大多数情况下不需要进行数据 光滑或网格重划，数值稳定性十分好。本文对角点问题做了仔细的处理。 本文用二维情况下的算例对三维数值波浪港池进行了校核。首先给定自由面 上孤立波的波形和速度势函数的初值，计算了孤立波的传播。然后模拟造波板造 波产生非线性波浪，计算了在港池另一端为全反射边界时驻波的形成和在另一端 使用人工粘性消波来模拟周期性波浪在开边界的传播。 本文最后给出了一个三维算例，初步分析了本数值波浪港池用于模拟三维问 题的适用性。 第9 页 第二章数学描述与数值方法上海交通大学硕士论文 第二章数学描述与数值方法 2 1 控制方程和边界条件 带自e h 表面的完全非线性势流的控制方程如下所示。在笛卡尔坐标系 ( x ，y ，z ) 下用速度势 ，) 来表达不可压无粘无旋三维流动( 见图2 1 ) ，z = o 为 未受扰动的水面位置。 f i g u r e2 1s k e t c ho fc o m p u t a t i o n a ld o m a i nf o r3 d b e mo f f n p f e q u a t i o n s t a n g e m i a lv e c t o r s a tp o i n tr ( oo ft h ef r e es u r f a c e 0 ( f ) a r ed e f i n e d 哺m ) a n d o u t w a r dn o r m a lv e c t o r n d e p t h i ss e tt oh 咄( ) y ) 图2 1 非线性势流方程组3 d b e m 求解的计算域示意图在自由表面r ，( f ) 上任意一点r ( t ) 处的切向矢量定义为( $ ) t r l ) ，与它们垂直指向外的矢量定义为 水深的变化为h - b ( x ，y ) 速度定义为 “= v 妒= ，v ，w ) ( 2 1 ) 荣 0 页 第二章数学描述与数值方法 上海交通大学硕士论文 边界为r ( ，) 的流场区域q ( ，) 的连续性方程是速度势函数的l a p l a c e 方程 v 2 矿= 0 ，在q ( ，) 内 对应的三维空间里的g r e e n 函数定义为 g ( x , x i ) = 去并( 圳= 一石17 r f ! 这里 ，= 工一x - ( 2 2 ) ( 2 3 ) ( 2 4 ) r 。为边界面上的点x s z ) 到边界面上的参考点工，s ( x ，y ，毛) 的距离，弹 代表x 点处边界面指向外的单位法向矢量。 g r e e n 第二函数将方程( 2 2 ) 变换为边界积分方程 a ( x 脚( x 沪r 凳( x ) g ( ，) - 妒( 工) 罢( ，) ) d r ( 2 5 ) r f x o no n 这里口( 上，) = 舅( 4 n - ) ，岛为x l 点处边界的外实体角( 对于一个光滑边界是 2 n ，参见b r e b b i a 6 ) 。 边界r 按照不同的边界条件划分成不同的部分。在自由表面r ，( f ) 上，妒满 足非线性运动学和动力学边界条件 面d r = 牲= v 妒 在r ，( f ) 上 ( 2 6 ) 警= 一g z + 三v 印妒一鲁在r ，( f ) 上 ( 2 7 ) 这里胃是自由表面上流体质点的位置矢量，g 是重力加速度，p l 是大气压，p 是 流体密度，随体导数定义为 尝；昙v ( 2 8 ) d fa f 有不同的方法可以用来造波，这到后面再详述。这里说明一下造波壁的边 界条件的给法。比如，当波由模拟计算域的边界r ，上的造波板运动产生时，造 波板的运动和速度为 x 2 x ， 凳：h ，以 在r ，( f ) 上 荔刈， 征i ，1 ( f ) 上 ( 2 9 ) 纂i l 页 第二蕈数学描述与数值方法上海交通大学硕士论文 上划线表示给定的值。 沿着底部l 和其它的静止边界r ，：，边界条件为 娑：o ，在l 和r ，：上 d 所有的边界为 r s o u l ur r 2 u l 。 2 2 时间步进 ( 2 1 0 ) 自由表面边界条件( 2 6 ) 和( 2 7 ) ：i m i 立在时间r 积分建立下一个时刻o + a t ) ( a t 是可变时间步长) 自由表面的新位置和边界条件。 根据g r i l l i 等人( 1 9 8 9 ) 在二维模型里所采用的方法，本文采用混合欧拉 拉格郎日方法的二阶显示泰勒展开来表达自由表面的新位置置“+ a t ) 和 ( 胄( ，+ ，” r ( t + a t 埘等+ 竽等+ 。】 眨 勰+ r ) ) - ) + ，警+ 竽等+ o 【( f ) 3 】 ( 2 1 2 ) 这些泰勒展开式的系数被表达成速度势及其导数和它们沿着自由面的法向和切 向导数的函数( 计算的细节见后面的部分) 。 更具体的讲，一阶的系数由方程( 2 6 ) 和( 2 7 ) 给出，需要计算自由面上 的西和0 0 o t 。这可由求解时刻t 的边界积分方程( 2 5 ) 得到，边界条件为( 2 9 ) ， ( 2 1 0 ) 和( 2 1 2 ) ( 详见下部分) 。二阶的系数由方程( 2 6 ) 和( 2 ，7 ) 的随体 导数得到，需要计算时刻t 的a 妒a t 和a 2 o t o n 。这可由求解类似于方程( 2 5 ) 的边界积分方程得到。这第二个边界积分方程的自由表面边界条件可由求解对 于西的边界积分方程后从b c m o u u i 方程( 2 7 ) 中得到 詈= 一髟一j 1v 妒v 一鲁，在。( 。上 ( 2 1 3 ) 簟1 2 页 第二章数学描述与数值方法上簿变通大学坝士论文 对于由造波板产生的波浪，方程( 2 9 ) 给出 a 2 a ( u p 押) 瓦万一瓦一 对于底部和其它静止边界 旦卫：o 在l 和：上 ( 2 1 4 ) ( 2 1 5 ) 注意到对于毋和a 妒a t 的边界积分方程都是在时刻r 求解，因而对应于相同 的边界几何形状，具有相同的离散形式( 见下一部分) 。所以第二个边界积分方 程的求解只需要占用求解第一个边界积分方程所需时间的一小部分。这使得这 种时间步进方法非常有效，特别是相对于在其它文献中所采用的高阶 r u n g k u t t a 法或a d a m s - b a s h f o r t h - m o u l t o n 法，这些方法要求在每个时间步长的 几个中间时间里对边界积分方程进行多次运算。这个方法的另一个好处是它是 显式的，并且在计算，+ f 时刻自由表面的量时用到变量的空间导数。这给数值 解提供了更高的稳定性，使得有可能用较大的时间步长得到差不多的精度，使 整个的求解更加有效。 2 3 边界的离散 2 3 1 边界兀概述 速度势矽的边界积分方程( 2 5 ) 和其对应的a # a t 的边界积分方程通过边 界元( 见b r e b b i a ( 1 9 7 8 ) ) 来求解。边界被离散成r 个节点，用m ，个高阶边 界单元在这些节点中的m 个中插值。这样，在第_ | 个边界单元r 。里，边界的 几何位置和物理量( 为了简单起见，令”= 妒，q = 却a 拧) 都通过形函数来离散。 这里，形函数由一个解析的表达式给出，为在参考单元0 ，上的高阶多项式。这 个参考单元是这m ，个任意形状的边界单元通过一个变量转换得到的。参考元素 上的局部坐标表示成( 孝，玎) 【1 ，1 】。每个单元七上的几何量和物理量的变化由单 纂1 3 页 第二量数掌描述与数值方法上海交通大学硕士论文 元的节点值工，“，q ，和形函数( 善，7 7 ) 来描述 x ( 善，7 7 ) 。n ，( 孝，v ) x ( 2 1 6 ) “( 善，7 ) 2n ，( 吉，1 7 ) u ，g ( 孝，7 7 ) = ，( 毒，r ) q j ( 2 1 7 ) 这里产，m ，m 是每个单元r 。的节点个数，重复的下标表示求和。 形函数为( 善，叩) 的多项式，其系数由要求方程中的“( 孝，叩) 在节点x ，上取节 点的值“气即在方程( 2 1 7 ) 中 “( 善( t ) ，叩( 墨) ) = n ( 孝，仇) “，= “f 这样，对于一个有m 个节点的参考单元的第i 个节点，形函数必须满足 g 。，仉) = j f ，u 2 1 ，m 在参考单元0 。上 ( 2 1 8 ) 这里屯是k r o n e c k e r 符号。对于给定的单元元素形状和阶数，方程组( 2 1 8 ) 的解给出形函数的近似表达式( 见下一部分) 。 232 三维鬲阶边界元 等参单元能在单元内部提供高阶的近似，但是在单元之间的公共节点上只 具有几何量和物理量的c 。阶连续性。所以等参单元单元间的连续性不好，不是 十分理想的三维边界单元。在这里，本文采用g r i l l i ( 2 0 0 1 ) 所使用的方法，定义 一种新的单元。这种单元既有高阶的精度，又有很好的连续性。 边界单元是4 节点的四边形，用三次形函数插值，插值函数用到本单元的 节点和在各个方向的相邻节点，共m = 1 6 个点。这样，在计算方程( 2 5 ) 的边 界单元积分时只用到三次形函数的一部分( 一般是中间的部分，参见图2 2 ) 。 羹l 页 第二章数学描述与数值方法上海交通大学硕士论文 z x 卜 p v ，i o ， + 1 1 01 1 + 1 l 1 f r 6 7 l 孝，彘) 6 f i g u r e2 2s k e t c ho f1 6 - n o d ec u b i cb o u n d a r ye l e m e n t f a n d c o r r e s p o n d i n gr e f e r e n c ee l e m e n t 0 t h e b o u n d a r ye l e m e n ti sl o c a t e di nt h em i d d l ea n dt h es u r r o u n d i n gr e f e r e n c ee l e m e n t sa n d n o d e sa r ef o rc a l c u l a t i o no f v a r i a b l e so f t h eb o u n d a r ye l e m e n t t h ec u r v i l i n e a rc o o r d i n a t es y s t e m ( 品埘，一) h a sb e e nm a r k e d a tp o i n t ，o f t h ee l e m e n t 图2 2 1 6 节点三次插值的边界单元r 的示意图和其对应的参考单元0 ，口中问的单元对 应于一个边界单元，四周的8 个单元和相应的节点为相邻的用于插值的单元和节点单元 上的点r 的曲线坐标( s ，挣b 捏) 标示见图。 蔫砖页 第二童数学描述与数值方法 上海交通大学硕士论文 二维的双三次形函数定义于参考单元上，为两个一维的三次形函数。( ) ( c = 1 4 ) 的乘积，即 n ，( 孝，r ) = n 6 r ，( ( f ，彘) ) n d r ，( ( 玑) ) ( 2 1 9 ) 这里6 d = - i ，4 ；j = 4 ( d - 1 ) + b 。由( 2 1 8 ) ，有 n ：( u j ) = 6 t ，“i = ( 2 i 一5 ) 3 ( 2 2 0 ) 这里i = 1 ，4 。求解方程( 2 2 0 ) ，得到 ，( ) = ( 1 - a ) ( 9 a 2 _ 1 ) ， ， 0 2 ( ) 2 素( 1 一2 ) ( 1 3 a ) ， ， 0 3 ( ) 2 素( 1 一卢2 ) ( 1 + 3 p ) ， 。( ) = ( 1 + a ) ( 9 a 2 - 1 ) ， ( 2 2 1 ) 从转化到参考单元上的局部坐标( 善，7 7 ) 的表达式为 ( z ，z 。) = z 。+ i 1 ( 1 + z ) ( 2 2 2 ) 这里z = 善或r ；屁= 氏( 或r o ) = - 1 ，- 1 3 ，或1 3 ，按照选取由肌2 1 6 个节点 定义的9 个四边形中( 见图2 2 ) 的哪一个而定。 2 3 3 边界积分方程的离散 边界积分方程( 2 5 ) 的积分转化成沿边界单元的积分之和，在每个单元上 积分都转化到参考单元。上计算。这个转化关系为 妻= 掣= 半 秘例c ) ) ) ，a a f。i 。a la f “ 坠o r l o r l = 牮胁例洲h ， 这里户l ，肌，在单元r 。上( 七= 1 ，m r ) 。由方程( 2 2 2 ) ，得 第1 6 页 第二章数学描述与数值方法 上海交通大学硕士论文 塑：丝：三 a a q 3 相应的切向单位矢量定义为 ，( ：笔， 力，d 亡 啊= 斟 埘( 酬：士罢 n 2d 叩 铲例 第三个矢量定义为同一点的垂直方向上的矢量 堡：堡。鱼 8 a 8 q 对应的单位矢量定义为 蟛肋= 壶毒= ( 2 2 4 ) ( 2 2 5 ) ( 2 2 6 ) ( 2 2 7 ) 这个矢量指向计算域的外部，只要定义所考虑的单元的( s ，珊) 的叉积指向外。 单元转化的j a c o b i a n 矩阵定义为 髂新哪7 在边界积分方程中第k 个单元积分中的要用到的j a c o b i a n 矩阵的行列式的大小 为 a r k ( 亭，可) l = ，h ：，后= 1 ，m ， ( 2 2 8 ) 田万程( 2 2 1 ) 利【2 2 3 ) 一( 2 2 5 ) ，e j 以征早兀内的仕俐一个点处解析地 计算出来。 通过转化，方程( 2 5 ) 中的积分表达成以下离散的形式 j l 附，a - 鲁n g i d r = 鍪誊，( 善川) g ( ( 工( 参破_ ) 旷( 最，7 ) p 乒班筹( x ，) = 塾，等= 争，。等 簟1 7 页 第二章数学描述与数值方法 上海交通大学硕士论文 h 妒鲁一篱 n r ，r = 。，( 善，叩) 里鱼兰! j ；孑堕型1 ，( 孝，7 ) i d 考d 叩) 妒( x ，) e # k 丸= 置f ”办 ( 2 - 3 0 ) 这里1 = i ，n ，d 和k 分别代表当地( 即对于单元k 来说) 和整体( 即组装 后的) 的d i r i c h l e t 矩阵，e 和置4 为n e u m a n 矩阵。注意这里下标，为边界上 的整体节点编号，表示单元七上的节点值。在方程( 2 2 9 ) 和( 2 3 0 ) 中的g r e e n 函数、形函数和j a c o b i a n 行列式分别由方程( 2 3 ) ，( 2 1 9 ) ，( 2 2 1 ) 和( 2 2 8 ) 给出。 利用方程( 2 2 9 ) ，( 2 3 0 ) ，边界积分方程( 2 5 ) 最终表达成 ； 口= 艺。4 9 ，一k 。”j ( 2 3 1 ) ，= l 这里f 1 ，。 方程( 2 3 1 ) 对应的边界条件是：( 1 ) 甜= 谚或a 妒a t 的d i r i c h l e t 条件( 方 程( 2 1 2 ) 或( 2 1 3 ) ) ；( 2 ) q = 却锄或a 2 ( k o t o n 的n e u m a n n 条件( 方程( 2 9 ) ， ( 2 1 0 ) 或( 2 1 4 ) ，( 2 1 5 ) ) 。通过将未知的节点变量移到方程的左边，将己知 的移到右边，最终得到的线性代数方程组为 c p 一置，“) “p d 置叫4q ，= 置一。一 + 置一n 扣g ( 2 3 2 ) 这里i = 1 ，n ，表示节点的总体编号；g = l ，n 。，代表自由表面f 上的具 有d i r i c h l e t 条件的点；p = 1 ，n ，代表边界r ，1u r r 2u l 上具有n e u m a n n 条 件的点。c 是一个由系数组成的对角矩阵。 23 4 刚体模型方法( r i g i dm o d em e t h o d ) 方程( 2 3 2 ) 的系数g 可以由纯几何的方法得到，即计算节点处的实体 角q 。这些系数可以间接地由一种叫做刚体模型的方法( 参见b r e b b i a 6 ) 的方 纂l i 页 第二章数学描述与数值方法上海交通大学硕士论文 法更为精确和有效地得到。 考虑一个均匀的d i r i c h l e t 问题，这里在整个边界上都己知“，而且 百= c o n s t 0 。这时法向的量。必须在每个节点为零。方程( 2 3 2 ) 简化为 c ，+ ”阳= o ( 2 3 3 ) 要求对所有的，花括号内的值均为零。这样，通过分离方程左边的对角元素， 我们得到 n ， + 玩“) = 一b “，z = 1 ，r ( 2 3 4 ) ，( t ，j e l 由这个方程可以解出方程( 2 3 3 ) 的每一行方程的对角元素项，其值即它的非 对角系数的和的相反数。在离散形式的方程组( 2 3 2 ) 中，这些对角项被直接 替换掉。 这个方法可以提高数值解的精度；而且从物理意义上讲，对于势流，这也 对应于使离散的问题在特定条件下精确满足零流量条件。 23 5 离散后边界条件在角点的处理 边界条件和法向方向在边界之间的交界处( 例如自由表面或底部和边壁间 的交界处，见图2 1 ) 通常是不同的。这些相交的地方被称做边缘，对应的离散 的节点称做角点。要在数值模型中反映这些不同，本文采用双节点方法( g r i l l i ( 1 9 9 0 ) ) 来表达这些点，它们具有相同的坐标，但是不同的法向。这样，在每 个节点上，有两个不同的离散的边界积分方程来表达。 对于d i r i c h l e t - n e u m a r m 条件，我们有方程，例如( i ) l = p ，在造波边界r r ， 上；( i i ) l - - f , 在自由表面r ，上。因为速度势必须是单值的，即在给定的点上 唯一，在最后的方程组中，这两个方程中的一个必须改为满足条件九2 ，( 即 速度势连续) ，方程的形式变为 蟊二= o ，v j p ，k ”= 肘，l ，= m 妒， m = 慨“r ( 2 3 5 ) 这里l p 表示方程的右端的值，吖为权数，等于方程左端系数绝对值的最大值， 纂1 9 页 第二毒数学描述与数值方法 上海交通大学硕士论文 使用这个参数可以保证方程组求解的性能。 对于n e u m a n n - n e u m a r m 条件，我们有，例如( i ) l = p ，在造波边界r ，l 上 ( i i ) l = b ，在底部l 上。这是速度势连续条件写成九一九= 0 ，这两个量都是 未知的。在最终的方程组里，将两个方程中的个换成表达这个条件，另一个 方程形式不变。改变后的方程为 置= 0 ，j p ，b ，置盏。一m ，k 月4 2 m ，l ，2 0 ， m = i k u n j 一( 2 3 6 ) 类似的连续性条件同样使用于相应的边界积分方程在角点处的a a t 。 在三个边界面的交界处，用三角点来表示这些点。在最终的方程组中，对 应于这些点的三个边界积分方程中的两个替换成表示速度势及其对时间取偏导 的连续性条件。 23 6 网格的生成 网格的生成采用十分简化的办法。初始的网格等距划分。在自由面上，网 格节点固结在流体质点上，随着时间的步进，网格节点随着流体质点的移动而 移动。在其它的面上，网格节点在z 方向上均匀分布，随着自由面高度的起伏 变化，每个时间步后重新划分。由于所采用的高阶边界元十分有效，在计算中 自由表面无需进行数据光滑或网格的重划。注意到切向矢量的求法( 见后) ，在 每个方向上单元的数量不得少于4 个。 2 4 数值积分 由于解析形式十分复杂，方程( 2 2 9 ) ，( 2 3 0 ) 中的积分相对于每个节点由 数值积分得出。对于每个单元，这些积分表达成当地矩阵d 口，e 口。 当参考的节点不属于要积分的单元时，使用一般的g a u s s 积分来积分即可。 当参考节点属于要积分的单元时( 图2 3 ，= ，) ，g r e e