（流体力学专业论文）基于高阶边界元的无反射数值波浪水池研究.pdf

哈尔滨工程大学硕七生毕业论文 l itr tli ih i i i i i i i i i i i 置 摘要 在开边界的海洋结构物水动力计算中，远方辐射条件一直是人们所关注 的问题，亦是难以解决的一个问题。有效并准确的建立一个无反射边界条件 就成为水动力计算中亟待解决的一个问题。 本文采用高阶边界元的方法，自由面边界条件采用混合欧拉拉格朗日 ( m e l ) 观点，利用计算机高级语言f o r t r a n 9 0 分别建立了二维线性和全非线 性两个数值波浪水池。对于远方辐射条件，本文将基于多次透射公式的人工 边界条件应用于线性及全非线性数值波浪中，此人工边界条件能够有效的吸 收线性及非线性数值波浪，但此人工边界条件有其局限性，仅能够有效吸收 低频波浪。而人工阻尼区方法能够有效消除高频波浪，却不能有效吸收低频 长波，增加阻尼区长度能够解决这一问题，但计算时间却大幅增加，这在三 维计算中尤为明显。因此本文针对包含了高频和低频的不规则波，提出耦合 人工边界条件与人工阻尼区的方法，从而得到了一个新的有效而精确的无反 射边界条件。 基于多次透射公式的人工边界条件在二维全非线性数值波浪中的应用， 以及耦合阻尼区方法在不规则波中的应用，在国际上均属首创。本文对几个 消波方法进行了详尽的分析和讨论，最终得到了几种消波方法各自的计算参 数和适用条件。从而解决了海洋结构物水动力计算中开边界的难题。最后本 文以哈尔兵工程大学拖曳水池为对象，建立了一个二维摇板数值波浪水池， 并对两个典型的强非线性波浪进行了模拟研究。 关键词：无反射边界条件；人工阻尼；多次透射公式：边晁元；数值波浪 哈尔滨丁程大学硕士生毕业论文 i i i i i i i i i i i 葺置i i i 置i i i 宣i i i l li i e i i i 萱i 宣i i i 宣i i i i i a b s t r a c t p e o p l ep a ya t t e n t i o nt ot h ep r o b l e mo fn o n - r e f l e c t i n gb o u n d a r yc o n d i t i o n s i i lt h eu n b o u n d e dr e g i o na l lt h et i m e ，a n di ti sa l s oad i f f i c u l tp r o b l e mo nt h e c o m p u t a t i o no f o c e a ns t r u c t u r e s h y d r o d y n a m i c s a ne f f i c i e n t a n de x a c t n o n - r e f l e c t i o nb o u n d a r yc o n d i t i o ns h o u l d b es e tu pn e c e s s a r i l yo nd y n a m i c s c o m p u t a t i o n i nt h i sp a p e r , w es i m u l a t el i n e a ra n dn o n l i n e a r2 dg r a v i t yw a v e sp r o p a g a t i n g b yl l i 曲一o r d e rb e ma n dt h em i x e de u l e r - l a g r a n g e ( m e l ) i nt h ef r e es u r f a c e t h r o u g ht h ec o m p u t e rp r o g r a ml a n g u a g ef o r t r a n 9 0 f o rt h ed i s t a n c er a d i a t i o n c o n d i t i o n , w ea p p l ya r t i f i c i a lb o u n d a r yc o n d i t i o n ( a b c ) d e r i v e df r o m m u l t i - t r a n s m i t t i n gf o r m u l a ( m t f ) t oa b s o r bf u l l yn o n l i n e a ra n dl i n e a rn u m e r i c a l w a t e rw a v e s b u tt h i sa b ci sl o c a l ，i ti se f f i c i e n ti nt h el o wf r e q u e n c yo ral i m i t e d f r e q u e n c yr a n g e 。d a m p i n gz o n e ( d z ) m e t h o dc o u l d a b s o r bh i g hf r e q u e n c i e s w a v e se f f i c i e n t l y , b u tn o tv e r ye f f i c i e n tf o rl o wf r e q u e n c i e sw a v e s ，i n c r e a s i n gt h e d a m p i n gz o n el e n g t hc o u l ds o l v et h i sp r o b l e m ，b u ti tw i l lc o s tm u c hm o r et i m e t o c o m p u t a t i o ne s p e c i a l l yi n3 ds i m u l a t i o n t h e r e f o r e ，w ep r e s e n tc o u p l i n gt h ea b c a n dd a m p i n g7 _ z n em e t h o d ，w h i c hg i v eae f f i c i e n ta n de x a c tn o n - r e f l e c t i o n b o u n d a r yc o n d i t i o n u s i n gt h ea r t i f i c i a lb o u n d a r yc o n d i t i o n ( a b c ) d e r i v e df r o mm u l t i - t r a n s m i t t i n g f o r m u l a ( m t f ) t oa b s o r b2 df u l l yn o n l i n e a rn u m e r i c a lw a t e rw a v e si s t h ef i r s t t i m ea th o m ea n da b r o a d ，a n dt h ec o u p l i n gm e t h o di sa l s op r e s e n t e df i r s tt i m e i n t h ep a p e r , w ea n a l y z ea n dd i s c u s st h ea b s o r b i n gm e t h o d se x h a u s t i v e l y , a n dg e tt h e c o m p u t a t i o np a r a m e t e ra n da p p l y i n gc o n d i t i o nf o rt h em e t h o d sr e s p e c t i v e l y , w h i c hs o l v et h ep r o b l e mo fo c e a ns t r u c t u r e sh y d r o d y n a m i c sc o m p u t a t i o ni nt h e u n b o u n d e dr e g i o n 哈尔滨工程大学硕士生毕业论文 k e yw o r d s ：n o n - r e f l e c t i n gb o u n d a r yc o n d i t i o n ；d a m p i n gz o n e ；m u l t i - t r a n s m i t t i n g f o r m u l a ；b e m ；n u m e r i c a lw a v e 哈尔滨工程大学 学位论文原创性声明 本人郑重声明：本论文的所有工作，是在导师的指导下，由 作者本人独立完成的。有关观点、方法、数据和文献的引用已在 文中指出，并与参考文献相对应。除文中已注明引用的内容外， 本论文不包含任何其他个人或集体己经公开发表的作品成果。对 本文的研究做出重要贡献的个人和集体，均己在文中以明确方式 标明。本人完全意识到本声明的法律结果由本人承担。， 作者( 签字 日期： r e t 哈尔滨工程大学 学位论文授权使用声明 本人完全了解学校保护知识产权的有关规定，即研究生在校 攻读学位期间论文工作的知识产权属于哈尔滨工程大学。哈尔滨 工程大学有权保留并向国家有关部门或机构送交论文的复印件。 本人允许哈尔滨工程大学将论文的部分或全部内容编入有关数据 库进行检索，可采用影印、缩印或扫描等复制手段保存和汇编本 学位论文，可以公布论文的全部内容。同时本人保证毕业后结合 学位论文研究课题再撰写的论文一律注明作者第一署名单位为哈 尔滨工程大学。涉密学位论文待解密后适用本声明。 本论文( 囱在授予学位后即可口在授予学位1 2 个月后 口 解密后) 由哈尔滨工程大学送交有关部门进行保存、汇编等。 作者( 签字) 易纪爨守导师( 签字) ：险吱i7 p 日期：万卯哆年乡月，日立夕哆年；月j i 日 哈尔滨工程大学硕士生毕业论文 1 1 研究背景及意义 第1 章绪论 水波是自然界最普遍的现象之一，研究水波问题有许多现实的意义。在 海洋、河流、港口、湖泊中都存在大量与水波有关的闽题。在深海工程方面， 浮式或固定式钻井平台的设计、海底建筑的建造、海底油管的铺设，都必须 考虑海浪压力的传播和作用。不然，有可能造成重大事故。海浪研究在造船 工程方面的应用更是不言而喻的。在河流、湖泊中，泥沙的输运与环境和水 利建设都密切相关。在港口建设中，怎样提高防波的效果，也与水波的研究 有关。 人类对水波的研究有着悠久的历史，发展到今日，己经取得了许多的成 就。但水波问题并没有得到彻底充分的解决，比如，对于湍流的研究，有关 波浪破碎的机理等等。在自然科学的发展之中，理论研究与试验研究是相辅 相成、齐头并进的。在工程应用中，试验手段也发挥着巨大的作用。 与此同时，随着计算机技术的飞跃发展，计算机数值模拟技术在各研究 领域得到了广泛的应用，大到原子弹爆炸、全球气候的模拟，小到物质老化、 微电子运动的计算都离不开计算机的数值模拟。在空气动力学实验领域，科 学家早在8 0 年代就已成功的研制出精确的数字风洞模式替代了传统的实验风 洞装置，并将其应用于飞机制造、汽车设计等重要的工业领域。数值化的波 浪水池便在这种情况下应运面生了。 数值波浪水池是一种计算机的仿真模拟程序，其目的是用来尽可能的模 拟真实的波浪水池的各种功能，并最终能够替代实验波浪水池完成相应的科 学研究与工程应用。i s o p e1 9 9 7 年会议上明确提出了数值波浪水池的必备条 件： 1 必须具有实时模拟能力： 哈尔滨工程大学硕士生毕业论文 2 可在有界区域内进行模拟； 3 能够模拟重力作用下的自由表面波动； 4 要求自由面及固体边界面都满足非线性条件； 5 包含物理造波机制( 运动的边界、变化的压力) 从理论上讲实验水池和数值水池所实现的功能应当是完全一样的。但显 然与试验手段相比，数值模拟有它的许多优点。数值模拟可以节约人力、物 力、财力和时间。数值模拟的可重复性好、条件易于控制，比试验研究更自 由、更灵活，便于优化设计。有些试验上难以完成的问题，可以通过数值模 拟来实现，比如数值模拟可以比较容易地获得试验难以侧量的物理量。另外， 数值模拟的经济风险和人员伤亡与试验手段相比要小得多。人们在对波浪 进行研究时，往往要在水池中进行试验。如果能够通过数值方法来模拟水池 中的这些试验，将能更加方便人们对波浪的研究。数值波浪水池可以帮助人 们实现这一点。正因为有如此巨大的优势，世界各国的水波动力学研究人员 都非常重视波浪数值水池的开发和应用，可见其巨大的潜在发展前景。 1 2 数值造波及数值消波概述 1 2 1 数值造波 与实验水池对应，数值造波机也是数值波浪水池产生波浪的源动力。它 按照原理大致可分为两种类型。一种是仿物理造波技术，即仿照实验水池的 造波机原理，模拟固体振动边界作为扰动源，将其数值化并作为相应的边界 条件嵌入数值水池中。活塞式造波机就是一种最典型的仿物理造波机，它的 应用也最广泛，通过活塞的规律性往复运动向水池的内区激发出波浪，它的 原理简单，易于嵌入主体程序中，造波功能完善，最大的优点的是由于存在 现实模型，非常容易得到检验。类似的还有摇板式造波机、柱塞式造波机等( 参 考b u c h m a n n l 9 9 6 ) ；另- - 种是所谓的纯数值造波技术，可以想象这种造波技术 在现实中是无法实现的，它是利用数值模拟的优势创造出的理想化造波方式。 2 哈尔滨工程大学硕士生毕业论文 其种类包括设置造波边界条件和安放脉动奇点。与物理造波机不同，设置造 波边界条件是指直接控制边界面的一部分作为扰动源，这种边界通常是垂直 的并且固定的，而非运动边界( 实际上前述的区域分解方法的本质也是将相邻 的水池作为造波边界条件) ，般根据行波的解析解或数值解来给定造波边界 上的值。例如对于小振幅波，可使用线性理论的解析解，对于有限振幅规则 波，则可使用s t o k e s 波解、椭圆余弦波解等。若想模拟全非线性的规则波则 需要使用较复杂的波解形式。 另一种数值造波方法是安放脉动奇点。b r o r s e na n dl a r s e n ( 19 8 7 ) 首先使用 了一种水下脉动点作为扰动源，著名数值水池专家c l e m e n t ( 1 9 9 9 ) 使用了一种 独具匠心的水下旋转偶极子作为扰动源，成功的造出了一种向一面传播的不 对称波场。其物理模型相当于在自由面下存在一个旋转的偏心圆柱。由于偶 极子位于水中，因此在模拟时，需要将边界积分方程改变形式。 1 2 2 数值消波 数值消波亦称为无反射边界条件( n o n - r e f l e c t i n gb o u n d a r yc o n d i t i o n ， n r a c ) 在数值波浪水池中是一个十分重要的问题，亦是难以解决的问题。如 果在数值水池远方不能有效的消除波浪，波浪将会反射回来而影响到整个流 场，因此建立一个理想的无反射边界条件就显得至关重要。常用的消波方法 有以下几种： 1 辐射边界条件 在波浪线性频域数值模拟中，使用s o m m e r f e l d 辐射边界条件，是非常有 效的。通过自由面上满足辐射边界条件的g r e e n 函数，可以处理无界区域的 波场。s u g i h a r a & i s s k h ( 1 9 8 0 ) 通过实验分析研究了辐射边界条件在封闭区域 内的线性频域数值模拟的适用性，结果非常理想，与无限水域的理论解相当 吻合。而对于进行时域内的非线性模拟的数值水池来说，使用辐射边界条件 时则需要一些技巧。o r l a n s k i ( 1 9 7 6 ) 首先将其应用于时域模拟的数值水池中， 因此它通常被称为s o m m e r f e l d o f l a n s k i 辐射边界条件。从7 0 年代末到8 0 年 哈尔滨工程大学硕七生毕业论文 代中期，一些学者基于经典的一阶s o m m e r f e l d - o r l a n s k i 辐射边界条件提出了 更高阶的无反射边界条件，其中较为著名的就是e n g q u i s t m a j d a 无反射边界 条件和b a y l i s s - t u r k e l 无反射边界条件。 2 主动消波 对于非线性海浪，仅使用上面介绍的辐射边界条件消波是远远不够 的。于是人们尝试模仿实验波槽中的一些主动消波装置，运用波动信号的反 馈原理，通过边界的运动将传来的波动抵消掉。j a g a n n a t h a n ( 1 9 9 2 ) 、m a i s o n d i e u 和c 1 6 m e n t ( 1 9 9 3 ) 分别以消波器界面所受波浪力作为反馈信号，制造了活塞式 消波器并应用于二维波浪水池中。目前流行的主动消波器是一种被称为“主动 消波控制系统”( a c t i v ew a v e a b s o r p t i o nc o n t r o ls y s t e m ) 的技术，简称a w a c s 。 它使用“i n f i n i t ei m p u l s er e s p o n s e ”( i t r ) 过滤器以消波板的水位作为反馈信号 进行消波。 3 人工阻尼区 人工阻尼区又称为数值沙滩或海绵吸收层，其原理就是将自由面条件中 加入阻尼项，当波浪通过阻尼区时其能量将会被吸收，以达到消波的目的。 b a k e re ta 1 ( 1 9 8 1 ) 首次在其二维边界元模式中引入了阻尼消波区( d a m p i n g z o n e ) 的概念，此后这种方法即被广泛采用。阻尼项可以加在动力学自由面条 件和运动学自由面条件上，亦可以只加在动力学边界条件上或运动学自由面 条件上。 阻尼消波对高频波非常有效，但对于波长超过阻尼区长度的长波吸收则 不理想，增大阻尼区长度可以解决这一问题，但是以计算量激增为代价。 c 1 6 m e n t ( 1 9 9 6 ) 将其与主动消波器结合起来使用，取得了良好的效果，因为主 动消波器对低频波较为有效，耦合了两者后恰好取长补短。 4 人工边界条件 用透射理论直接在人工边界( 此处即远方辐射条件) 上模拟其应该满足 的边界条件，称这种边界条件为人工边界条件( a r t i f i c i a lb o u n d a r yc o n d i t i o n , 简记为a b c ) 。此方法即为本文采用的方法，将在下文中做详细的介绍。 4 哈尔滨工程大学硕士生毕业论文 1 3 本文的主要研究内容 目前世界上许多学者都致力于数值波浪水池的研究，并且提出了各种数 值消波的方法，但是仍然没有一种理想的方法既能够精确的消除不规则波， 又能够不额外耗费大量计算时间。本文通过计算机高级语言f o r t r a n 9 0 编制了 二维线性和全非线性数值波浪水池，并提出了一种新的消波方法，即将基于 多次透射公式的人工边界条件与人工阻尼区相结合的方法。通常来讲，人工 边界条件方法能够有效消除低频波，而人工阻尼则对高频比较有效，将两者 结合起来将是一个较为理想的消波方法。 第一章，绪论。主要阐述了本文的研究背景及目的意义。 第二章，数学模型和理论基础。详细叙述了本文中研究内容的数学模型 和基于高阶边界元法的线性及全非线性数值水池的数学实现 过程。 第三章，无反射边界条件。详述了基于多次透射公式的人工边界条件以 及人工阻尼区消波方法，并对他们的数值模拟结果进行了详尽 的讨论。最后提出了将两者相结合的消波方法。 第四章，数值水池及应用。以哈尔滨工程大学拖曳水池为模型，建立了 一个摇板数值水池，并模拟了两个典型的强非线性波浪，对理 论和工程应用都有一定的参考价值。 5 哈尔滨工程大学硕士生毕业论文 第2 章数学模型和理论基础 本章详述了以边界元( b e m ) 为基础的二维数值波浪水池的建立过程。 自由面采用混合欧拉拉格朗日( m e l ) 观点，重点介绍了b e m 方法的数值 实现过程，以及求解过程中的数值要点。 2 1 控制方程及边界条件 这里的二维造波水池问题参考图2 1 ，笛卡尔坐标系原点建立在静止水面 与推板的交点处，x 轴与静水面一致，y 轴正方向向上。自由面、推板和底面 分别用品、品和& 表示，远方辐射面( 人工边界) 用表品示。 sw 一 d a l m p i n g z 。n e q 。 s 图2 1 数值波浪水池模型 我们假定流体为不可压缩无粘性并且无旋，整个流场的运动可以用一个 速度势矽来表示，并且在计算域q 内满足拉普拉斯方程： v 2 矽= 0 i nq ( 2 1 ) 在推板和底面上，边界条件可以写成 鼍= 痧( f ) 0 瓯( 2 - 2 ) 6 哈尔滨工程大学硕士生毕业论文 i tr i 暑宣薯；i 宣i ；暑i i 鲁宣萱宣；每暑暑暑罱i 育暑宣昌宣 娑：o。以& 锄 。 ( 2 3 ) 式中，u ( f ) 为造波板运动的法向速度，法向相量指向域外。 在人工边界& 上速度势所满足的边界条件由多次透射公式得到，其解法 将在后面的章节介绍。瞬时非线性自由面条件可以写成： o n 品 o n 品 ( 2 - 4 ) ( 2 5 ) 线性自由面条件将忽略运动学条件与动力学条件中一级以上的无穷小量，并 且在静水面上满足 a a 叩 一= 一 加a 砂 饼 等+ = o 0 1 1j ，= 0 o n y = 0 ( 2 6 ) ( 2 7 ) 式中弭g 分别表示时间和重力加速度。非线性自由面条件写成拉格朗日形式 d xa 西 d t苏 d y a 西 d t砂 罢= 一g y + 押2 相应的线性自由面条件写成拉格朗日形式则成为： 脚8 垂 d t 卸 d 西 茜2 一g y 7 ( 2 8 ) ( 2 9 ) ( 2 一1 0 ) ( 2 1 1 ) ( 2 1 2 ) l l 却瓦 卯 型跏扣 铆百 驯 却一钞 缈一a 哈尔滨工程大学硕士生毕业论文 式( 2 8 ) 一( 2 1 2 ) 即为本文所采用的混合欧拉一拉格朗h 型( m i x e de u l e r - l a g r a n g e ) 自由表面条件。 2 2 数值过程 2 2 1 边界元法求解的基本原理 想要对自由面的非线性运动进行连续的模拟，就必须通过已知的初、边 界条件得到自由面在下一时刻的运动速度场。本文采用二维简单格林函数的 高阶奇点分布边界元法。对于已知的初边值问题，一般可分为两个步骤来求 解： 1 在每一固定时刻，给定物面法向速度九，则未知量为速度势痧；给定 自由面及辐射面上的速度势，则未知量为法向速度杰。通过求解边界积分 方程即可得到未知量。 2 根据求得的自由面上质点的法向速度丸，以及已知量速度势矽，利用 m e l 自由表面条件可得到下一时刻的自由面位置和速度势分布，这样就可反 复步进求解。 。 在求解过程中最为关键的两点就是求解边界积分方程的准确性和稳定、 精确的步进模式。下面先介绍边界积分方程的求解过程。 对于本文的数学模型，对每一个时刻，在空间上，采用简单格林函数l n r 的边界元求解。用g r e e n 公式将问题转化为积分方程 叫咖且尝彤薏产= o p ，3 ， 其中qy o u r 算域边界上的一点，为q 点上的法向量，指向域外。c p ) 为p 点处边界的两条半切线之间的夹角在域内的部分。硷为q 的整个边界。 g = 1 1 1 ，= p 一剑，a q 为计算域周界，包括、s w 、& 、但除去 p 点的邻域。 8 哈尔滨工程大学硕士生毕业论文 在品、& 上，矽已知，丸未知；在品、& 上，唬已知，矽未知；这样 积分方程( 2 1 3 ) 即可离散成为线性代数方程来求解。 积分方程求解的精度是保证后面时间步进准确性的关键。由于自由表面 的非线性，对于自由面上速度势的求解，采用二阶以下奇点分布不能保证切 向速度的连续性，因而本文在处理这种非线性自由面的问题时采用三次奇点 分布。方法如下： 在流场边界上取一系列节点，其上的矽和九用分段三次插值表达，插值 点为临近的节点，即以边界弧长为变数的三次差值多项式系数中含有了临近 四个节点的矽和丸，把每个节点当作定点，积分方程就成为以诸节点上矽、么 为未知数的线性代数方程。 积分方程离散过程中，一个基本的工作是求每个高次分布线元的相应积 分，本文采用了贺五洲和戴遗山( 1 9 9 1 ) 一给出的该积分的精确表达式。 初始条件我们给定自由面为静水面，速度势分布为o ，底部边界固定， 推板按照给定的信号运动。 在自由面和物面交点p 处采用二重节点技术，此节点既可以看作是自由 面的点也可以看作是物面的点。p 点作为自由面的点时，速度势( p ) i 。已知， 可知( p ) l 蹄2 矽( p ) k ；所以在p 点只有 一个未知数，人工边界与底 面交点处用同样的方法。在自由面与人工边界交点o 处，只有一个己知量 矽( q ) ，而 均为未知量。由于自由面及人工边界上的速度 势均为已知量，故由自由面和人工边界的切向速度所决定的速度在x 轴的分 量即为 ，这样q 点就贿掣卜糊飙 9 哈尔滨工程大学硕士生毕业论文 对于线性数值水池，其所有满足边界积分方程的边界皆为固定不动。即 自由面条件在未扰动的静水面上满足( y ：o ) ，推板边界则在其简谐振动的中 心满足，人工边界和底部边界亦固定。而对于非线性水池，自由面条件在自 由面上满足( y = 巧) ，并且物面湿表面和人工边界上随时间变化，物面上的 网格亦随湿表面变化而变化。推板边界满足其位移函数的变化规律，人工边 界始终垂直x 轴，并且其位置随自由面与人工边界交点处的质点运动而运动。 人工边界与自由面相交的节点，其位置坐标及速度势值按m e l 方法求 得，与自由面上的其他质点求法相同，而人工边界上其他节点则由多次透射 公式得到。 2 2 2 时间步进模式 在某一时间步求解边界积分方程后，便得到了自由面上所有节点的法相 速度杰，自由面切向速度则可由自由面速度势沿自由面切向的导数求得。 这样就能得到自由面上每个节点的速度，按照拉格朗日型自由面条件便可得 到节点下一时刻的位置和速度势痧值。时间计算若要计算几十个周期，对每 个质点来说，步进过程的稳定和准确性是关键。为保证计算的稳定和准确， 必须使用高阶的步进模式。当自由面需要重新划分节点时，若采用多步法的 步进模式需要前面几个时刻自由面的信息，将不能保证下一时刻的准确性， 因而本文采用具有单步法四阶r t m g e k u t t a 法。 令p 表示三个变量( x ，y ，痧) ，假设p 是f = n a t 时刻的值，则f = ( n + 1 ) a t 时刻 的值p 肿1 可以用下式表示： 1 0 哈尔滨工程大学硕士生毕业论文 宣宣i 暑暑宣；高暑暑宣苦；i i 皇宣；昌i 宣昌i 葺- - - ti i - - = ： 2 一m i m z t t i 毒暑i 宣i 高置皇i i 暑暑i 宣i i i 暑置宣 = 只+ 詈阢+ 2 k 2 + 2 k s + k t ) k 。= s ( t 。，只) 毕水+ 等卅等蜀) 毕水+ 等”等k ：) 在时间上每步进一次，就需要计算一次边界积分方程，也就是说每步迸一个 时间步长要计算四次边界积分方程。但是r u n g e 。k u t t a 法具有四阶精度，而 且稳定性也非常好，d o m m e r m u t h 和y u e ( 1 9 8 7 ) 对使用线性自由边界条件的 r k 时域积分方法作了v o n n e w m a n n 稳定性分析，得到了它的稳定性c o u r a n t 条件为：a t ：旦垒，其中出为时间步长，a x 为相邻节点的空间步长。此稳 充g 定性条件很容易得到满足，本文中线性与非线性模拟均满足此稳定性条件。 2 2 3 自由面的光顺处理 在模拟强非线性自由面时，会产生节点的拉伸和压缩变形，为了避免线 元过大的变形，每隔一定时间步长后需要对自由面进行光滑处理，即自由面 上节点的重新划分。本文采用c z w a n g ，g x w u ( 2 0 0 6 ) 介绍的方法： 假设在一条曲线上有n + 1 个节点霉伉，y 。= 0 , 1 刀) ，在节点霉和置+ ，之 间我们用三次b 样条进行插值。用b 样条统一表示为： 声( 甜) - zq ，g 矿q ( f - - o ，1 以) ( 2 1 5 ) 式中“从霉点材= o 到丘“点z = 1 ；矿町为控制点，岛， ) 为三次b 样条函数 可以表达为： ( 咖吉( 1 训3 ，( 炉l ( 3 u s - 6 u 2 + 4 ) ， 哈尔滨工程大学硕士生毕业论文 垦“垆i 1 ( 一3 u 3 + 3 u 2 + 3 州) ，b 3 , 3 ( 咖吉以 ( 2 1 6 ) 整条曲线需要通过所有的点置，用公式表示为 鼻( o ) = 考( e + 4 k + l + 巧+ 2 ) ( i = 0 , 1 刀) ( 2 - 1 7 ) 可以写成矩阵的形式： 阻妒 ：园( 2 - i s ) 式中l 】为一个三对角线矩阵，陟】和p 】分别为控制点矿和节点置的组成向 县 里 m = 6 0 14 l l 41 0 6 ，纠： 虼 一 k ： - 吒+ 2 f 己 崩：6 i ： b 令曲线的边界点瓦和露+ ：分别设置为豆，和丘卅，而巨。和丘+ 。可从下式中得到： 巨。三2 露一蜃，丘+ 。= 2 丘一露一。 节点磊到节点霉的曲线长度可由下式计算得到： 耻j v k 引d u j 。出( 圳，l ，规) 式中 f ，d s 、 l 瓦j ，2，并且( 妄，亳) = 掣= 赛挈+ ， 我们把整个曲线分成m 份。每一个线元的长度可以根据流场的特性来决定。 如果我们希望接近节点磊的线元小一些，而接近节点磊的线元大一点，则可 以应用下面的过程得到。假定新的节点为只。，( f _ o ，1 m ) ，从节点昂，到节点 只。，的弧长s ，可以由下式获得： 1 2 = 最黼辫忙叫，州。， 式中为常数g f l l 。 除了节点的拉伸和压缩变形，在水波模拟的时域模拟过程中经常会看见 锯齿形的现象。本文用能p _ a 里t 9 的方法来解决这一问题。假设用位置向量 q ( f = o ，1 ，z ) 来定义一系列离散的节点，节点磊在光顺后变为霉。我们定 义通过节点霉( 江0 , 1 ，n ) 的曲线的能量为： 疋= 寺( q + l - q ) 2 ( 2 2 1 ) 式中= 1 1 0 , - 0 , 一。0 ，为点q 与点q 一。之间的距离；巳= ( a 一仍一。) ，其中 只( f = o ，1 ，疗) 为点丘的x 坐标和y 坐标。光顺的过程要确保丘与 q ( f = o ，1 ，功之间的差别尽可能地小，为了达到这个目的，我们定义了一 个目标函数： e = 口疋+ 屏( 乃- q j ) 2 ( 2 - 2 2 ) 式中口和乃是常数，g f ( i = 0 ，1 m ) 为点q 的x 坐标和y 坐标。将式( 2 - 2 1 ) 代 入( 2 2 2 ) 得 枷否, , - t 而1c 半一学2 + 嘉讹训2 上式的第一项反应了曲线的光顺度，第二项反映了光顺前与光顺后的差 别。令c 对b 0 = 0 ，1 ，刀) 的导数为0 ，则可得到这两项的最小值。将其写成 n e q = o 式中彳是一个带宽为5 的矩阵， 1 3 ( 2 2 4 ) 哈尔滨工程大学硕士生毕业论文 彳1 = c o磊 岛q + 呸 。 吃一2 。矾一。 口nb nc n 在本文的模拟中岛被设定为一个统一的值，参数口是j 个关于 ( 扛o ，1 ，, ) 的数或者可以有数值试验得到。我们采用下面的方法得到，假 定( 江o ，1 ，z ) 的最小值为k ，令口= c 3 一，式中c 为一个系数，选为5 到1 0 之间的数可取得较好的结果。令磊二= q = 吃= o 使边界点磊和磊分 别等于露和丘。数值模拟结果见图2 2 ，2 3 图2 2 曲线光顺示例，4 0 个节点 1 4 哈尔滨工程大学硕士生毕业论文 2 2 4 计算精度校核 图2 3 曲线光顺示例，1 0 0 个节点 对于线性问题，计算结果可与线性解析解直接进行比较来检验计算精度。 而对于完全非线性的问题计算，则需要额外的处理。本文采用两种方法来校 核全非线性数值计算的正确性。 1 交点处计算校核 自由面与物面的交点是时域完全非线性步进解计算中特别值得注意的， 保证该点计算的合理性对整个步进过程是很重要的。由于该点既在物面上有 在自由面上，因此它应该满足自由面动力学边界条件p = 0 ，考察步进过程中 交点处压力的变化，便是对交点处物面条件与自由面条件相容与否的一个重 要检验。 2 边界积分方程求解精度校核 根据质量守恒，容易得到一个检验边界积分方程求解合理与否的一个简 单判断，l 。唬出= o ，既校核了边界积分方程求解精度同时也检验了流域的质 1 5 哈尔滨工程大学硕士生毕业论文 量守恒。 2 3 全非线性受力公式 对于水动力的计算，多数研究者采用b e r n o u l l i 方程算的压力，然后沿物 面积分的传统方法。但是在全非线性的时域水动力计算中，由于物面的变动， b e r n o u l l i 方程中第一项速度势关于时间的偏导数的计算一般不太准确也不是 非常方便，特别是在计算波浪和水面浮体的运动耦合问题时，往往由于水动 力计算的误差导致模拟的失败。本文采用段文洋( 1 9 9 5 ) 提出的先积分后微 分的水动力计算公式。 在大地坐标系o - x y 中，瞬时湿表面砸) 的方程为y = y ( x ，f ) ，物面与波面 的左右两交点为b 和a ，物面上的速度势为f x ，y ( x ，f ) ，f 】，如图2 3 所示。 易知 方砂出劫 一= = - i - 即 咖y = y ；丸+ y t 对( 2 2 5 ) 式两边微分得 呶+ 妨以= 儿噍+ 咒丸+ 兵啦+ 虼 又 丢，c 5 矽，勺刃= 丢：呈矽出= 1 6 ( 2 - 2 5 ) ( 2 2 6 ) 哈尔滨工程大学硕士生毕业论文 = 芤力+ 办以+ 丢( 矽丸) 出 = j 1 力+ 慨+ 晏( 触) 防 丑l “。 j a = 正谚+ 方咒+ 彩+ 丸九以+ 丸咒 出 = 豇破+ ( v 矽) 2 出一，点未办刃 鲁息唿刃= 丢：呈妒儿斑= 一鼽嘞以+帆+丢(批)卜b = 一j 1 识以+ 办- 以+ 帆+ 晏( 丸儿) i 出 l” j ( 2 2 7 ) zz = 一争t 匆一时，y | y l + 毋y 味+ 旌y x 七母y 记+ 牵多。y x + 咖正+ 孛咖x y 0 出 aaaa - - k , + - 忙d y - 彬方一p ( 九+ 以屯+ z 站+ 虼) 出 b8b8 = 豇绣+ ( v ) 2 咖一，f ：矽豪丸刃 丢，c c 翮2 瓤织+ ( v ) 2 翮一点杀( v ) 刃 所以受力公式为 僻一尸l ( 肌t ) t - 抑2 翮 = 一p 丢，c 5 菘刃+ o c 5 圭( v ) 2 菇一杀( v ) 刃 易知 ( 2 2 8 ) ( 2 2 9 ) ( 2 3 0 ) 8 以 石y +出 、l i ， m 办 + 石 爿八b = 哈尔滨工程大学硕士生毕业论文 ( v 枷和杀( v ) 卜吐矽( w 枷) 故 除州杀( v ) 卜匿v 一扣归卜d * 班例 代入( 2 3 0 ) 得到全非线性受力公式： 户( f ：i p 丢，c c 痧元刃+ p ，( | 鼍v 一三( v 矽) 2 元 刃。2 3 。， + p 矽( 一办7 + 丸歹) i 2 4 本章小结 本章介绍了数值波浪水池的建立过程，对数值波浪水池的数学模型以及 数学实现过程都做了详尽的阐述。数值水池的基本原理为高阶边界元方法， 自由表面采用混合欧拉拉格朗日( m e l ) 方法，时间步进模式采用四阶 r u n g e k u t t a 法，并且介绍了自由面光顺处理的数学过程，以及计算模式的精 度校核和全非线性受力问题。 1 8 哈尔滨工程大学硕士生毕业论文 第3 章无反射边界条件 无反射边界条件在计算流体力学的开边界计算中有着至关重要的作用。 如果无反射边界条件不能b l 徊p c 好的消波，则反射波浪将会与原来的波浪产生叠 加效果，致使计算结果失真，因此建立一个有效的无反射边界条件将是决定 数值模拟结果准确与否的主要因素之一。本章将基于多次透射公式的人工边 界条件和人工阻尼区方法分别应用到数值波浪中，并针对两种消波方法各自 局限性，提出将两种方法耦合的消除不规则波的无反射边界条件。 3 1 人工边界条件 本章将详细介绍一种基于多次透射公式的人工边界条件，该方法由廖振 鹏院士提出并在许多领域得到广泛应用。孙善春、徐蒙曾利用此方法研究了 二维及三维物体的线性水动力学问题并得到了理想的数值计算结果。 3 1 1 人工边界条件概述 用透射理论直接在人工边界( 此处即远方辐射条件) 上模拟其应该满足 的边界条件，称这种边界条件为人工边界条件( a r t i f i c i a lb o u n d a r yc o n d i t i o n ， 简记为a b c ) 。 从物理观点看来，a b c 应保证且仅保证由发自内场振源产生的外行散射 波穿过人工边界进入外场无限域，即外行散射波满足a b c ，而内行波动则不 满足a b c ，基于这一观点，国内外近数十年的有关研究工作可划分为两个方 向：其一是模拟外行波穿过整个人工边界进入无限域；其二是分别模拟外行 波穿过人工边界的每一点分别进入无限域。前者形成积分型a b c ，称为全局 ( g l o b a l ) a b c ；后者形成微分型a b c ，称为局部( 1 0 c a l ) a b c 。 全局a b c 保证穿出整个人工边界的外行波满足无限域内的所有场方程 和物理条件，包括无穷远处辐射条件。建立全局a b c 的思路可以概括为：从 1 9 哈尔滨工程大学硕士生毕业论文 这些数学物理条件或满足这些条件的解析解出发，一般通过傅立叶变换在频 域内建立a b c ，然后利用反变换得到时域a b c 。由于满足了无限域内的全 部数学物理条件，以连续形式给出的全局a b c 是精确的。在近场波动的数值 模拟中需要将其离散化，故在实际上它只具有和整个问题的离散模型相同的 精度。 局部a b c 用于保证射向人工边界任一点的外行波从该点穿出边界。局部 a b c 的连续形式是微分公式，其离散形式则为有限差分公式。局部a b c 的 主要特征是时空解耦，即一个边界节点在某一时刻的运动仅与其邻近节点在 邻近时刻的运动有关，而与所有其他节点和其他时刻无关。这一解耦特征极 大的简化了数值计算。此类a b c 因其实用性为研究者所重视。自2 0 世纪6 0 年代末以来，局部a b c 已成为a b c 研究的一个焦点。常用的局部a b c 均 建立在单向波概念之上，可视作某种单向波动方程。所谓单向波动在物理上 可以理解为从人工边界一侧单向传送到另一侧的波动。例如从内场振源向外 部无限域辐射的波动。目前局部a b c 的研究工作沿着两条技术路线进行：其 一是针对特定的无限域模型导出特定的局部a b c ；其二是对一般的无限域模 型导出具有普适性的局部a b c 。前者在国际上己流行多年，并导致一系列有 实用价值的局部a b c 的建立；后者主要在我国发展，并创建了透射理论和相 应的局部a b c 。 0 1 2 透射理论及透射公式 导出局部a b c 最简单的思路 莫过于模拟外行波在人工边界上 穿越边界的过程。为此，首先给出 外行波的一般表达式。设x 轴在所 考虑的边界点上垂直于人工边界 人工边界 计算区 li 2 l i c 口f l 出。 1 1r 1 1 r 图3 1推导多次透射公式的集合关 哈尔滨工程大学硕士生毕业论文 并指向外部无限域( 图3 1 ) ，则沿x 轴传播的单向波，其一般表达式为 舭x ) = 乃( c _ t - x ) 歹 ( 3 1 ) 该波动函数由任意多个单向波动组成，且这些单向波动的相速度c ，一般不 同，f 为未知的任意函数。考虑其中的一个分量，略去下标后可以写成 矽( f ，石) = f ( c ，f 一功( 3 2 ) 可知矽是波动自变量c ：f x 的函数，所以 妒l j f + a t ，功= q ，x c ，a t )( 3 3 ) 式中缸为任一假定的时间步距。上式表明，外行波沿x 轴的视传播可以通过 用点工一c x a t 在时刻t 的值替换点x 在f + 址时刻的值来模拟。由于c 。一般是 未知的，所以用人工波速c 。代替，则可写成 矽o + a t ，功= 4 ( t ，x c 。a t ) + o + a t ，x ) ( 3 - 4 ) 式( 3 - 4 ) 的误差项可表示为 妒9 + a t ，x ) = + f ，z ) 一妒( f ，x c 。a t ) ( 3 5 ) 可知，对于给定的巳，和址，该误差项为波动自变量c - x 的函数a 因此， 误差项也是一个相同类型的单向波。这是一个重要的观察结果。在给定的c ， 乞和址的条件下，此误差项可以写成 矿o + f ，x ) = 石( c 工t - x ) ( 3 - 6 ) 式中石是波动自变量c x t - - x 的某一函数。因此，仿效式( 3 4 ) ，误差波 矽( f + 出，x ) 可以写成 ( f + a t ，x ) = 矽o ，x - - c 。f ) + ? 矽 + a t ，x ) ( 3 7 ) 将式( 3 5 ) 中的t 和x 分别换成t 一址和x - c a t ，则得 o ，x - - c 。f ) = ( f ，x c 。f ) 一( f a t ，x - 2 c 。f ) ( 3 8 )