（船舶与海洋结构物设计制造专业论文）用边界元方法对构件附加水质量的研究.pdf

中文摘要 结构振动特性除受其刚度影响外，与其质量也存在着很大关系。准确的计算 附加水质量对于分析构件在水中的运动规律具有重要意义。一直以来人们普遍认 为水中物体引起的附加水质量仅与物体的形状有关。本文运用三维边界元方法， 通过数值计算，证明了在一定条件下流体边界、构件长度及构件周围管节点对附 加水质量均存在较大影响。本文以方形和圆形构件为例，分析了构件附加水质量 受这些因素影响的变化规律。 本文的数值模型采用高阶三维边界元方法，边界的几何形状和流场变量采用 八节点的四边形单元，用三次形函数插值，这种单元不仅在单元内具有高阶的精 度，同时又在单元之间具有很好的连续性。本文对单元内的数值积分进行了有效 的处理以保证其具有较高精度。使用速度势连续性条件妥善处理了对数值解的稳 定性有重要影响的自由表面和边壁交界处的边界条件。 通过数值计算得出附加水质量除受构件形状决定外，流体边界、构件长度和 其周围的管节点对构件附加水质量均有很大的影响。本文给出了这些因素对附加 水质量的影响范围。孤立柱体的计算主要考虑流体边界和构件长度对其附加水质 量的影响规律。复杂构件的附加水质量的计算主要是计算其附加水质量的数值， 并与孤立柱体相比较，得出一些基本的规律。计算的结果是通过编制的边界元程 序进行的，经过与以往结论的比较，计算结果显示了良好的精确性和稳定性，可 以为工程和理论研究提供借鉴。 关键词：附加水质量，边界元，管节点，流体边界，长度 a b s t r a c t t h ed y n a m i cp r o p e r t i e so fas t r u c t u r ea r er e l a t e dn o to n l yt oi t ss t i f f n e s s ，b u ta l s o t 0i t sm a s s t h e r e f o r e ，i t sv e r yi m p o r t a n tt oc a l c u l a t ec o r r e c t l yt h ea d d e dm a s sf o r u n d e r w a t e rs t r u c t u r e s i th a sb e e nl o n ga c c e p t e dt h a tt h ea d d e dm a s si sr e l a t e do n l yt o t h es h a p eo fb o d i e s t h r o u g hn u m e r i c a le x a m p l e su s i n g3 - db o u n d a r ye l e m e n t m e t h o d ，i t sp r o v e di nt h i sp a p e rt h a tt h ef l u i db o u n d a r y , t h em e m b e rl e n g t ha n dt h e t - j o i n tm a ya l s oh a v eb i ga f f e c t so nt h ev a l u eo fa d d e dn l a s si nm a n yc a s e s t h e e f f e c t so ft h e s ef a c t o r sa r ea n a l y z e di nt h i sp a p e rt h r o u g he x a m p l e sw i t hs q u a r e m e m b e r sa n dr o u n dm e m b e r s 3 一db o u n d a r ye l e m e n tm e t h o di su s e di nt h i sp a p e rt oc a l c u l a t et h ea d d e dm a s s f o ru n d e r w a t e rs t r u c t u r e s q u a d r i l a t e r a le l e m e n t sw i t h8n o d e sa r eu s e df o ra l l e l e m e n t s ，w h i c hn o to n l yl e a dt oh i g h e ra c c u r a c y , b u ta l s om a k eas m o o t hc o n n e c t i o n a m o n gd i f f e r e n te l e m e n t s e f f e c t i v em e t h o d sa r eu s e dt oe n s u r ea na c c u r a t en u m e r i c a l i n t e g r a t i o nw i t h i ne a c he l e m e n t c o n t i n u i t yc o n d i t i o n sf o rv e l o c i t yp o t e n t i a la r c c o n s i d e r e do nt h ef r e es u r f a c e ，o nt h ei n t e r f a c el i n eo ft w ob o u n d a r ys u r f a c e sa n da t c o m e rp o i n t s ，w h i c hh a v eas i g n i f i c a n ta f f e c to nt h es t a b i l i t yo ft h en u m e r i c a l p r o c e d u r e f r o mt h en u m e r i c a lr e s u l t si nt h i sp a p e r , i t sc o n c l u d e dt h a tb e s i d e sb o d ys h a p e s f o rs t r u c t u r em e m b e r s ，t h ef l u i db o u n d a r y , t h em e m b e rl e n g t ha n dt h er - j o i n tm a ya l s o h a v eb i ga f f e c t so nt h ev a l u eo fa d d e dm a s si nm a n yc a s e s n ee f f e c t so ft h e s ef a c t o r s o nt h ea d d e dm a s sa r es t u d i e da n dg i v e ni nt h i sp a p e rt h ee f f e c t so fm e m b e rl e n g t h f l u i db o u n d a r yo na d d e dm a s sa r es t u d i e df o ri s o l a t e dp o l e s w h i l ef o rt j o i n t s ， n u m e r i c a lr e s u l t sf o ra d d e dm a s sa r eg i v e na n dc o m p a r e dw i t ht h o s ef r o ms i m p l e p o l e s t og e tt h en u m e r i c a lr e s u l t s a3 - db o u n d a r ye l e m e n tm e t h o dp r o g r a m h a sb e e n u s e d t h r o u g hc o m p a r i s o nw i t hp r e v i o u sr e s u l t s ，b e t t e ra c c u r a c ya n ds t a b i l i t ya r e s h o w nf o ro u rr e s u l t sw h i c hc o u l db eu s e df o rb o t he n g i n e e r i n ga n df u n d a m e n t a l r e s e a r c h e s k e yw o r d s ：a d d e dm a s s ，b o u n d a r ye l e m e n tm e t h o d ，t j o i n t ，f l u i db o u n d a r y , m e m b e rl e n g t h - l v 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得的 研究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他人已经发表 或撰写过的研究成果，也不包含为获得墨鲞盘茎或其他教育机构的学位或证 书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中 作了明确的说明并表示了谢意。 学位论文作者签名：咿仓 签字日期： 谝7 年，月，9 日 学位论文版权使用授权书 本学位论文作者完全了解苤鲞盘堂有关保留、使用学位论文的规定。 特授权鑫壅盘堂可以将学位论文的全部或部分内容编入有关数据库进行检 索，并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。同意学校 向国家有关部门或机构送交论文的复印件和磁盘。 ( 保密的学位论文在解密后适用本授权说明) 学位论文作者签名： 呃7 仓 j 签字日期：瑚7 年月修日 导师签名：- 7 冈也 签字日期：t 孓年月夕日 天津大学硕士学位论文用边界元方法对构件附加水质量的研究 1 1 研究工作的目的及意义 第一章绪论 随着世界科学技术的发展和工业制造规模的日益扩大，陆地上的资源已经不 能满足需求。浩瀚的海洋中有着无尽的财富，当人民意识到海洋是另一块未开发 的领地时，便把目光投向了富含各种资源的海水和海底，例如海洋波浪的利用， 海底石油和天然气的开发，希望从海洋和海洋覆盖下的陆地中获取人类需要的资 源。 石油和天然气是不可再生的一次性资源，是现代工业的命脉，而海底蕴藏着 丰富的石油天然气资源。据1 9 9 5 年的估计世界近海已探明的石油资源储量为3 7 9 亿顿，天然气为3 9 亿立方米。据不完全统计，海底蕴藏的油气资源储量占全球 油气储量的1 3 。世界海洋石油的绝大部分存在于大陆架上。关于海洋石油的储 藏量，由于勘测资料和计算方法的限制，得出的结论也是不同的。但是，海洋石 油有巨大的开采价值，这是勿容质疑的。中东地区的波斯湾，美国墨西哥之间的 墨西哥湾，英国挪威之间的北海，中国近海，包括南沙群岛海底，都是世界公认 的蕴藏海洋石油最丰富的地区。 海洋油气的开采需要船舶与海洋工程工业的发展，出现了如巨型油船，大型 集装箱船，化学品船，液化气船和张力腿平台，f p s o ，v l f s 等新船型和新海洋 平台形式。而这些水中的航行体或结构物随着海洋波浪和冲击波等流体动力激励 或其他非流体的激励而发生振动时，其周围的流场也发生变化，这种流场的变化 反过来使结构所受的动力发生变化，形成反馈的流体结构相互作用的流固耦合 问题这是一个非常复杂的跨学科问题，涉及到固体力学，流体力学，动力学和计 算数学等，并与许多工程技术领域有关。 2 l 世纪是海洋的世纪，海水中的建筑物受波浪等的冲击要发生振动，而这 些振动会使建筑物或船体受到损害甚至摧毁。所以研究海洋建筑物在水中的振动 是势在必行的，振动中结构的附加水质量在有的时候是不可忽略的，所以研究结 构的附加水质量对研究结构的振动有非常重要的意义。对整个海洋工业的发展也 有很重要的意义。 第一章绪论 1 2 目前的研究现状及理论综述 刑景棠掣l j 曾对流固耦合力学进行了阐述。人民对流固耦合问题的早期认识 源于飞机工程中的气动弹性研究。人们对飞机气动弹性问题以及土木工程中的流 固耦合问题进行了大量的研究在船水耦合动力学问题的早期研究，考虑船是一刚 性浮体探讨其由于船的刚体运动引起的扰动。h a s k i n d e 2 。3 】对这类问题作了深入理 解的分析。他应用g r e e n 定理构造了由于船的存在及运动引起的速度势，并推导 了点源g r e e n 函数的形式。至于考虑船体弹性变形的水弹性理论的创造性工作应 归功于b i s h o p 和他的学生p r i c e ，他们在其专著1 中对其研究及理论作了详述。 有些学者巧7 1 问题的屈曲问题进行了探讨。 由于流固耦合的复杂性，尽管可假设流体是无旋无粘的理想流体，在大多数 情况下仍不能求得解析解。因此，一般情况下要求助于数值分析，如有限元方法， 边界元方法和有限差分方法，并且对流体作充分简化。具体实现时，由于流体和 固体采用的数值分析不同，主要有两种方法：一种是对固体和流体都采用有限元 方法；另一种时对固体采用有限元方法，对流体采用边界元方法。 对于线性流固耦合系统的有限元分析，通常采用的有限元格式有两类，一种 是结构和流体以位移矢量为场变量的位移位移格式醒母1 另一类是以结构位移矢 量和流体场变量的混合型格式，如位移压力格式，位移位移势格式和位移 速度势n 1 。1 2 1 格式。对于位移位移格式，其有限元与结构动力学方程一样，因而 可照搬结构动力学的模态分析理论。但是，这个格式的流体节点自由度太多，并 使流体刚度矩阵奇异，分析时出现大量伪模态，从而给实际应用带来困难。对于 混合型格式，由于流体场变量己满足运动无旋的条件，因而不会出现流体矩阵奇 异的问题，此外，流体节点自由度也只有位移位移格式的三分之一。但是，位 移压力型和位移位移势型有限元方程是非对称的频域分析将涉及非对称问题， 必须采用一些非标准算法或对称化技术。戴大农位移势型有限元方程出发， 进一步建立流固耦合系统的模态分析理论，证明了系统具有实频率，进而到处系 统的一对具有三个实对称矩阵的广义特征值问题。吴一红劲位移和流体速度势， 从加权余量的伽辽金方法出发，导出了系统的有限元方程，专著”的流固耦合 问题的理论和有限元分析方法作了系统的论述。 以上可以看出对于流固耦合问题的研究己经取得了很大的成绩，但是在实际 的工程领域中，能够实际应用的工程计算软件却比较少。我们针对本课题的特点， 将流体设为均匀的，无粘，无旋的理想流体，并讨论线性的小挠度情况。在有限 天津大学硕士学位论文用边界元方法对构件附加水质量的研究 元方法中，对于不可压缩流体，最常用的是用附加水质量解耦钔。对于有界的 流体域我们可以根据实际流体范围划分单元，但是对于无限流体域，其范围比较 难确定，通常取足够大的范围来计算，但这样做工作量大，费时费力。能求解析 解的情况是极少数的， 莫里森( m o r i s o n ) ”公式对一些规则的结构给出了一些 近似公式，但对大多数工程结构就缺乏一种定量的计算方法。我们利用位移压 力格式，编制相应的程序，并对规则的结构进行了计算。在对结果分析的基础上， 进而提出了一般结构物的流固耦合系统的附加水质量计算的流体范围的确定方 法，使我们在实际计算时能够划取适当的流体范围，对实际工程计算有一定的借 鉴，也使我们对实际模型进行计算时有了明确的计算范围。 航行体或结构物在水中运行或受到水的冲击时，总是处在一定的水深中，它 除了受到流体的冲击或者动力激励外，一定水深产生的压力对其运动或者振动都 有影响。水的深度对结构的作用是不随结构运动与否而改变的，它对结构振动的 影响相当于结构具有一定的初始压力，初始压力对结构的影响主要是由于结构的 儿何非线性引起的。在有限元分析中找们可以看出主要体现在结构的“初应力刚 度矩阵”上，也称为“儿何刚度矩阵”。 由以上可以看出，在大的方向的研究有比较成熟的方法与理论，但是在一些 细微的小部分中还没有具体的研究与分析。附加水质量的研究现在还处在初级阶 段，结构振动特性除受其刚度影响外，与其质量也存在着很大关系【1 6 】。准确地计 算附加水质量对于分析构件在水中的运动规律具有重要意义。理论推导只能给出 无限流场内无限长孤立柱体的附加水质量【1 7 j ，而实际生活中的流体及位于其中的 构件均存在边界和其它构件的影响。而现阶段构件的附加水质量在工程应用中往 往是由实验得出的一个系数，但是随着边界的变化和其他构件的影响，这个系数 是不同的。在这种条件下，研究附加水质量的问题势在必行。 1 3 本文的主要工作 本文以势流理论为基础，对流体与任意三维物体相互作用问题进行理论和数 值研究。对势流理论的计算已经有很多研究，但也存在很多没有完全解决的问题， 流体中物体运动和载荷的三维频域计算方法不是很完善。本文的主要工作集中在 计算流体中结构物的附加水质量问题。 目前边界元方法仍然是研究流体与物体相互作用问题的最常用的方法，本文 使用边界元方法计算流场。采用线性方法计算运动，并建立了三维边界元的数值 模型，以解决工程中的实际问题。本文的主要工作有： 第一章绪论 1 基于势流理论，建立了计算海工建筑物上的三维数值模型。 2 编写程序进行了求解。为了更好地处理动边界，模型求解采用了边界元法。 3 计算了简单形状结构物的附加水质量，并计算结构物附加水质量随流体边 界和自身形状变化而其附加水质量的变化趋势。 4 计算一些复杂结构物的附加水质量，并总结一些基本的问题。并对结果进 行验证。 天津大学硕士学位论文用边界元方法对构件附加水质量的研究 第二章数学描述与数值方法 势流理论假设流体为不可压缩、无粘性和无旋的，因此可应用速度势来表达 流场内的流体速度分布。对于物体和波浪相互作用的水动力问题，流场中的速度 势满足拉普拉斯方程，还须满足自由表面条件、底部条件、物面条件和辐射条件 等边界条件以及适当的初始条件。虽然拉普拉斯方程是线性的，但是边界条件是 非线性的，这就给理论研究造成了极大的困难。由于这一模型求解的困难，人们 依据不同的假设，建立了多种不同的计算理论。主要可以根据对时间项的不同处 理方法分为频域方法和时域方法两类。在这两类方法中，根据物面条件和自由面 条件阶数的不同，又可分为线性和非线性方法。线性方法是指自由面条件和物面 条件的摄动阶次均为一阶的形式，而非线性方法则采用二阶甚至更高的阶次。根 据求解流场所采用的数学模型的维数不同，可分为二维理论和三维理论。根据求 解流场应用的数值方法不同，可分为有限元法，边界元法，有限差分法等。本文 采用边界元方法。 2 1 控制方程及边界条件 在波浪与物体相互作用问题研究中，通常假设水为不可压缩、无粘性和无旋 的，同时也忽略水的表面张力影响。这样的流场可以用势流理论进行分析。流场 中的速度势除应满足拉普拉斯方程外，还应满足非线性自由表面条件、底部条件、 瞬时物面条件以及适当的初始条件和辐射条件。边界元方法是求解势流问题的常 用方法。本章对势流理论和边界元方法的基本原理进行阐述。 带自由表面的势流的控制方程如下所示。在笛卡尔坐标系( x ，y ，z ) t 用速度势 矽( x ，0 来表达不可压无粘无旋三维流动( 见图2 1 ) ，坐标原点( 0 ，o ，o ) 为水中 心位置。 第二章数学描述与数值方法 f i g u r e 2 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 ff n p fe q u a t i o n s t a n g e n t i a lv e c t o r sa tp o i n tr ( t ) o ft h ef r e es u r f a c es la r ed e f i n e da s ( s ，m ) a n d o u t w a r dn o r m a lv e c t o ra sn 图2 1 非线性势流方程组3 - d - b e m 求解的计算域示意图。在自由表面s l 上任意一 点r ( t ) 处的切向矢量定义为( s ，m ) ，与它们垂直指向外的矢量定义为n 速度定义为 材= v 痧= ( “，v ，w ) ( 2 1 ) 本文假设流体为势流，采用边界元法进行求解 1 8 】。对于三维位势问题，控制 方程为l a p l a c e 方程 v 2 矽= 0f 2 2 1 考虑混合边界条件： 在s i 上，已知 a 彩 在是上， 茜已知 其中，n 为边界s 的外法矢量，s = 焉+ 是。 对于任一点q ，其相应的边界积分方程可写成： 洲+ f g ( 川) 掣钌= f 如) 掣劣 对于三维问题，式( 2 3 ) 的基本解为 g ( m ) k 一击 式中，为作用点p 至观测点q 之间的距离函数。为边 的参考点的距离。 ( 2 3 ) ( 2 4 ) 面上的点到边界面上 天津大学硕士学位论文用边界元方法对构件附加水质量的研究 窃的取值与边界光滑度有关，一般可用刚体位移等方法间接求解【1 9 1 。 边界s 按照不同的边界条件划分成不同的部分。在自由表面s 上，矽满足非线性 运动学和动力学边界条件 i d r ：“：v 矽 ( 2 5 ) d f 、 尝= 一g z + j 1v 印矽一告 ( 2 6 ) 这里r 是自由表面上流体质点的位置矢量，g 是重力加速度，是大气压， p 是流体密度，由于我们只关心动力的附加水质量，所以只考虑( 2 6 ) 式即可。 为了简化起见，本文未考虑由于波动而引起的水面变化，即水面将始终保持水平， 为了摆脱速度的非线性项的影响，我们仅考虑初始时刻，即水表面v = 0 ，且认为 附加水质量与水的速度大小无关。则由方程( 2 6 ) ，水表面的边界条件可写成： 丝：一笠 但7 ) 一= 一一 ，- o t p p 。为大气压。如忽略大气压的影响，可选择相对大气压儿= 一p 譬= o 。 由十本又仪计冗由十水的加运度向严生的附加水质量，且认为构件尺寸与拇 水深度相比很小，所以选取离构件足够远处水域的边界条件为 ”祟：o ( 2 8 ) 即离构件足够远处水体运动可以忽略。 海底的边界条件为 铲塞= o ( 2 9 ) 如构件沿y 轴负方向以单位加速度运动，则其边界条件可写成： 迎水面：q ：昙兰：一1 ( 方形) ( 2 1 0 ) 铲意= - - c o s o n ( 圆形) ( 2 1 1 ) 背水面：口。：昙挲：1 ( 方形)( 2 1 2 ) d t o n ：晏兰：e 。s 幺( 圆形) ( 2 。1 3 ) o t ( n 第二章数学描述与数值方法 2 2 边界元方法 对于波浪载荷的计算，人们提出了各种各样的计算方法，其中最著名的是有 m o r i s o n 等在1 9 5 0 年提出的计算公式。m o r i s o n 方程形式简单，便于工程应用，是 目前海岸工程设计中计算波浪载荷的最常用的方法。但是，m o r i s o n 方程的应用 有一个限制条件：物体的尺度相对于波长必须是较小的，即小尺度问题，以保证 入射波在物体周围基本上是未受扰动的。而当物体较大时，即大尺度问题，这一 点是不能保证的。入射波将在物体附近产生明显的畸变，即产生绕射现象。随着 大型结构物在海岸工程的广泛应用，绕射问题的意义越来越重要。目前，国内外 在波浪对大尺度物体作用力方面有较多的研究。当波浪高度相对于波长较小时， 可用线性理论使问题进一步简化。迄今为此，线性理论的发展已经比较成熟，对 于海岸工程中大量大直径圆柱可以很方便地求解。当波幅或波陡很大时，必须考 虑非线性问题。目前，对非线性问题的处理，就数学方法而论，主要有摄动法、 边界元法和直接进行数值求解( 主要用有限元法或有限差分法) 。 有限差分法的原理最简单、最直观，它是将计算域离散后，在原来的微分方 程中用差分代替微分，构成差分方程。有限差分法的理论严密，差分方程和它的 解必须收敛于微分方程和它的解，对于发展方程还要求满足稳定性条件，除此之 外还要限制山于数值计算引起的伪物理效应。早期有限差分网格都是矩形的，因 而限制了这一方法在复杂边界条件下的使用。近年来发展起来的网格生成技术可 以生成贴体坐标，可以在需要的地方进行网格加密，或采用自适应网格，在计算 过程中变量梯度较大时自动将网格加密。尤其是最近发展起来的非结构网格，更 贴近实际计算域的形状，这样有限差分法有了较大的优势。 有限单元法和边界元法都属于试函数方法。它们分别是将古典迎辽金法和有 限基本解法与有限单元技术相结合，在单元内设试解演化而来的。有限单元法最 大的优点是单元大小、形状可根据计算的需要确定。在所求变量梯度较大处单元 小一些，密一些，相反梯度小的地方单元可大一些，疏一些。还有一个优点是最 后形成的线性代数方程组的系数矩阵为稀疏矩阵，经过一定处理后可节约计算时 间。 边界积分方程方法的基础的出现远早于有限元法。早在十九世纪，人们就得 到了势流问题的积分方程：本世纪初，人们就利用分布电源、分布面源和分布涡 等方法解决了一些流体力学问题。尽管如此边界元法也只是在六十年代，有限元 充分发展后，才得到了迅速地发展。边界元法首先是在声波问题中，弹性力学问 题中，稳态弹性波地传播问题中得到了发展。不幸的是，由于当时有限元法的出 现及寻求问题的基本解的困难，大大的降低了这些工作的重要性。直到七十年代， 天津大学硕士学位论文用边界元方法对构件附加水质量的研究 受有限元思想的影响和由于有限元法在解决问题中遇到了困难，边界元法才作为 一种独自的方法得到了充分地发展。 由于边界元法可以使问题简化为边界积分问题，把问题降低一维，从而使未 知量大大减少，只解相当小的线性方程组。边界元法的另一个特点是，它能够处 理无穷远的边界条件，从而使它对外问题特别适用。正是由于边界元法具有这些 优点，才使它在近年来得到了迅猛地发展。目前，边界元法己经在有关弹性力学、 塑性力学、流体力学、断裂力学、结构动力学、地质力学、热传导、电场和磁场 等工程实际问题中得到了广泛地应用。边界元法解决水波问题是七十年代后才发 展起来的。七十年代，由于有限元法的应用正处于颠峰阶段，再加上有限元法比 边界元法具有更强的适应性，从而使有限元法受到了更多科技工作者的青睐。一 般地水波问题所涉及的域为无限区域，因此要处理无穷远处的辐射条件。由于边 界元法有可能解析地满足无穷远辐射条件，从而克服了有限元在处理辐射条件中 遇到地困难。尤其是在等水深情况下，用边界元法处理水波问题比较简单、适用。 因此，边界元法也受到了许多学者地重视。h w a n g 和t u c k ( 1 9 7 0 ) 解决了入射波在 直海岸地任意形状港湾内所诱发的振动问题。g a r r i s o n ( 1 9 7 2 ) 用边界元法计算了 任意形状三维物体所受到的波浪载荷。m e i ( 1 9 7 8 ) 综述了边界元法在水波绕射问 题中的应用。近些年来，边界积分方程方法( b i e m ) 逐渐用来计算重力波问题。 2 2 1 边界元概述 速度势矽的边界积分方程( 2 3 ) 和其对应的娑的边界积分方程通过边界元来 d f 求解。边界被离散成从个节点，用m 。个高阶边界单元在这些节点中的m 个中插 值。这样，在第k 个边界单元，里，边界的几何位置和物理量( 为了简单起见， 令u = ，q = o # i o n ) 都通过形函数来离散。这里，形函数由一个解析的表达式给 出，为在参考单元足。上的高阶多项式。这个参考单元是这m ，个任意形状的边 界单元通过一个变量转换得到的。参考元素上的局部坐标表示成( 孝，刁) 一l ，1 。 每个单元k 上的几何量和物理量的变化由单元的节点值 工，甜f ，q j k 和形函数，( 善，7 7 ) 来描述 x ( 孝，r ) = ，( 告，r ) x j ( 2 1 4 ) “( 古，r ) = ，( 毒，r ) u ， ( 2 1 5 ) g ( 孝，r ) = n j ( 告，r ) q ，。 ( 2 1 6 ) 这蜀= 1 ，2 ，m ，m 是每个单元s 。的节点个数，重复的下标表示求和。 第二章数学描述与数值方法 形函数为( 孝，r ) 的多项式，其系数由要求方程中的甜( 毒，7 7 ) 在节点# 上取节点 的值u ，即在方程( 2 1 5 ，2 1 6 ) 中 “( 毒( 。) ，刁( 。) ) = n j ( 专，r ，) u ，= 材， ( 2 1 7 ) 这样，对于一个有m 个节点的参考单元的第i 个节点，形函数必须满足 以( 缶，臻) = 磊，i ，j = 1 ，2 ，嬲在参考单元。上 ( 2 1 8 ) 这里点，是k r o n e c k e r 符号。对于给定的单元元素形状和阶数，方程组( 2 1 8 ) $ 3 解给 出形函数的近似表达式( 见下一部分) 。 2 2 2 三维高阶边界元 等参单元能在单元内部提供高阶的近似，但是在单元之间的公共节点上只具 有几何量和物理量的c o 阶连续性。所以等参单元单元间的连续性不好，不是十分 理想的三维边界单元。在这里，本文避免了求解位移和面力表达式的困难，而只 求解边界上某些特定点的函数值，边界上其他点的位移和面力值则通过适当的插 值公式来确定，这样就将一个求解未知函数的边界积分方程离散化为求解有限个 边界点上的未知边界量。所使用的方法，定义一种新的单元。这种单元既有高阶 的精度，又有很好的连续性。 边界单元是八节点的四边形，用三次形函数插值，插值函数用到本单元的节 点。这样，在计算方程( 2 3 ) 的边界单元积分时用单元的节点离散化就可以得到离 散化的方程。一个边界元上的节点按图示编号。单元上的节点数记做 l ，甜；，t ；( 1 = 1 2 l ) 表示单元上任意点的坐标，位移，和面力分量。令单元上任 意点石；的局部坐标为( 善，7 7 ) 则，“；，t ；可用如下的插值公式来表示： 三 = ，( 善，叩) 工； ( 2 1 9 ) i = 1 土 材，= n ，( 善，叩) “； ( 2 2 0 ) 1 = 1 t j = ，( 孝，州 ( 2 2 1 ) i = 1 式中，( 孝，刁) 为插值函数。 四边形八节点二次元的插值函数为 破= 一丢( 1 一孝) ( 1 7 7 ) ( 毒+ 7 7 + 1 ) 欢= 三( 1 誓) ( 1 刊 ( 2 2 2 ) ( 2 2 3 ) 天津大学硕士学位论文 用边界元方法对构件附加水质量的研究 珐= 丢( 1 + 孝) ( 1 一刁) ( 孝一露一1 ) 以= 吉( 1 + 孝) ( 1 一刁2 ) 九= 三( 1 + 孝) ( 1 + 7 7 ) ( 孝+ 7 7 1 ) 九= 圭( 1 一卵( 1 + 7 7 ) 庐7 = - 1 ( 1 一。( 1 + 功( 孝一7 7 + 1 ) 识= 圭( 1 一孝) ( 1 一，7 2 ) ( 2 2 4 ) ( 2 2 5 ) ( 2 2 6 ) ( 2 2 7 ) ( 2 2 8 ) ( 2 2 9 ) f i g u r e 2 2 s k e t c ho f 8 - n o d ec u b i cb o u n d a r ye l e m e n ts a n dc o r r e s p o n d i n gr e f e r e n c e e l e m e n t 名矿 图2 2 节点三次插值的边界单元的示意图和其对应的参考单元长刀。 2 2 3 边界积分的离散 边界积分方程( 2 3 ) 的积分转化成沿边界单元的积分之和，在每个单元上， 积分都转化到参考单元s 。上计算。这个转化关系为 嘉：鼍掣_ 些掣) 箬( 砌o ) ) - t ( 2 3 0 ) a 芒8 芒| 、 a l j a 芒 u u l 、“i u 。jj 、 骞= 掣 学) 考“( 僦) ) ) ( 2 3 ，) a 7 7a 7 7 、 a 7 a 刀 叭川” 、 7 这里j - l ，2 m ，在单元s 。上( 可= 1 ，2 m ) 。由以上分析可知 出一苏 8 考两 第二章数学描述与数值方法 蝣2 击毒，所c 善= 古高 c 2 3 2 ， ”斟驴矧 亿3 3 ， 第三个矢量定义为同一点上的垂直方向上的矢量 善：罢i o x ( 2 3 4 ) 一= 一一 8 a 鼍a q 、 对应的单位矢量定义为 孵刀) - 去善郜煳 ( 2 3 5 ) 这个矢量指向计算域的外部，只要定义所考虑的单元的( s ，m ) 的叉积指向外法 向。 单兀转化的j a c o b i a n 矩阵定义为 j = 袅，兰伽丁 ( 2 3 6 ) 在边界积分方程o p 多g k 个单元积分中要用到为j a c o b i a l l 矩阵的行列式的大小 l o r ( 考，7 ) i = 啊办2 ，k = 1 , 2 ，m s ( 2 3 7 ) 由以上方程可知这个值可以在单元的任何一点都能解析的计算出来。 通过转化，方程( 2 3 ) 中的积分可以表达成以下的离散形式 f 。，髻g ，钌= 兰j = l 兰k = l 【g e ，( 善刃) g ( x ( 善，矾工，) l j ( 善，刀) l d 弘7 7 ) 鼍( z ，) 2 氍n sm s 啪- 触- i f - = 和等 亿3 8 、 o 鲁办= 琵p 彤荆掣眦荆恸她，) = - o - - 。毛。 以= ”力 2 1 矧j 5 1 ( 2 3 9 ) 这里1 = 1 ，2 n s ，d 。和k 。分别代表各个单元和整体的d i r i c h l e t 矩阵，e o 和k d 为 n e u m a n 矩阵。这里的下标j 为边界上的整体节点的编号，表示单元k 上的节点值。 方程( 2 3 8 ) 和( 2 3 9 ) 中的g r e e n 函数，形函数和j a c o b i a n 行列式分别由以上的 分析方程中给出。 利用方程( 2 3 8 ) 和( 2 3 9 ) ，边界积分方程( 2 3 ) 最终表达成 天津大学硕士学位论文 用边界元方法对构件附加水质量的研究 c c | 弘l = k i q j k 寸链0 d ( 2 4 0 ) 这里i = 1 , 2 ，n s ， 方程( 2 4 0 ) 对应的边界条件为：( i ) 甜= 矽或a 铆的d i r i c h l e t 条件( 2 6 ) ； ( 2 ) q 2 a 锄或a 2 矽l a t a 的n e u m a n 条件( 2 8 ) ，( 2 9 ) 。通过将未知的节 点变量移到左边，将已知的移到右边，最终得到的线性代数方程组为 c + k d ，”) “户一k p t d q g = k p l d g p 一 c 一+ k g t ”心 1 ) 这里l = 1 ，2 n s ，表示节点的总体编号；g = l ，2 n g ，代表自由表面s 止的具有d i r i c h l e t 条件的点；p = l 2 一n 。，代表边界5 tu 5 2 上具有n e l l m a n n 条件的点。c 是一个由系 数口，组成的对角矩阵。 2 2 4 刚体模型方法( r i g i dm o d em e t h o d ) 方程( 2 4 1 ) 的系数c 可以由纯几何的方法得到，即计算节点处的实体角。 这些系数可以间接地由一种叫做刚体模型法的方法更为精确和有效地得到。 考虑一个均匀的d i r i c h l e t 问题，这里在整个边界上都己知u ，而且 甜= c o n s t 0 。这时法向的量q 必须在每个节点为零。方程( 2 4 1 ) 简化为 c f l + k j j ”) = 0 ( 2 4 2 ) 要求对所有的l ，花括号内的值均为零。这样，通过分离方程左边的对角元素， 我们得到 _ e c f ，+ ”) = 一砀“ 如。) “ ，1 = 1 ，2 ，n s( 2 4 3 ) 由这个方程可以解出方程( 2 4 2 ) 的每一行方程的对角元素项，其值即它的非对角 系数的和的相反数。在离散形式的方程组( 2 4 1 ) q b ，这些对角项被直接替换掉。 这个方法可以提高数值解的精度；而且从物理意义上讲，对于势流，这也对 应于使离散的问题在特定条件下精确满足零流量条件。 2 2 5 离散后边界条件在角点的处理 边界条件和法向方向在边界之间的交界处( 例如自由表面或底部和边壁间的 交界处) 通常是不同的。这些相交的地方被称做边缘，对应的离散的节点称做角 点。要在数值模型中反映这些不同，本文采用双节点方法来表达这些点，它们具 有相同的坐标，但是不同的法向。这样，在每个节点上，有两个不同的离散的边 第二章数学描述与数值方法 界积分方程来表达。 对于d i r i c h l e t - n e u m a m a 条件，我们有方程，例如( i ) l = p ，在边界s ，上；( i i ) l = f ， 在自由表面s ，上。因为速度势必须是单值的，即在给定的点上唯一，在最后的 方程组中，这两个方程中的一个必须改为满足条件丸= r ( 即速度势连续) ，方程 的形式变为 k 刍= o ，w p ，k 品= m ，l p = m 办， m = i k n l 嘶 ( 2 。4 4 ) 这里三。表示方程的右端的值，m 为权数，等于方程左端系数绝对值的最大值， 使用这个参数可以保证方程组求解的性能。 对于n e u m a r m n e u m a n n 条件，我们有，例如( i ) 1 = p ，在边界s ，上；( i i ) l - b ，在 底部s 上。这是速度势连续条件写成丸一九= 0 ，这两个量都是未知的。在最终的 方程组里，将两个方程中的一个换成表达这个条件，另一个方程形式不变。改变 后的方程为 x 三= 0 ，j p ，b ，k p b = 一m ，灭二= m ，l p = 0 m = 附i 懈 ( 2 4 5 ) 类似的连续性条件同样使用于相应的边界积分方程在角点处的a a t 在三个边界面的交界处，用三角点来表示这些点。在最终的方程组中，对应 于这些点的三个边界积分方程中的两个替换成表示速度势及其对时间取偏导的 连续性条件。 2 2 6 网格的生成 网格的生成采用十分简化的办法。初始的网格等距划分。在自由面上，网格 节点固结在流体质点上，网格节点随着流体质点的移动而移动。在其它的面上， 网格节点在z 方向上均匀分布。由于所采用的高阶边界元十分有效，在计算中自 由表面无需进行数据光滑或网格的重划。注意到切向矢量的求法( 见后) ，在每个 方向上单元的数量不得少于8 个。 天津大学硕士学位论文 用边界元方法对构件附加水质量的研究 2 3 数值积分 由于解析形式十分复杂，方程( 2 3 8 ) ，( 2 3 9 ) q b 的积分相对于每个节点由数值 积分得出。对于每个单元，这些积分表达成当地矩阵d ，。，e ，。 当参考的节点不属于要积分的单元时，使用一般的g a u s s 积分来积分即可。 当参考节点属于要积分的单元时( 图2 3 ，j - 1 ) ，g r e e n 函数和其导数里的r 在单元的 其中一点变成零( 方程( 2 4 ) ) 。当就j 1 时，积分d ，是弱奇异性的，而尾是非奇 异性的。对于前者的积分，本文采用子单元法，详见后面。对于后者，实际上， 当，一0 时出现的强奇异性已经消除掉了，成为系数a ，的一部分。当使用刚体模 型方法后，项玩已经通过方程( 2 4 3 ) 0 p 它的总体形式e ，”间接计算得到。 0 一1 l 。一 i ：+ 1 一 f 9 ； 1 1 2 f i g u r e 2 3s k e t c ho fc o o r d i n a t et r a n s f o r m a t i o n sf o rw e a k l ys i n g u l a ri n t e g r a l si na q u a d r i l a t e r a le l e m e n ts n 图2 3 四边形边界单元量，。内弱奇异积分的坐标转换示意图。源点与节点1 重合。 2 3 1 常规积分 方程( 2 3 8 ) 和( 2 3 9 ) 里的常规积分用二维的高斯积分来计算，每个方向上取 个高斯点。这些积分的形式为 + l + j。 厶= ，n ( 孝，, 1 ) a c d r z = 艺艺w g ( 以， ) ( 2 4 6 ) 一一l5 r 一- - i1一11 1h = l 这里乃代表积分d o 。或乞的核函数，( w g ，) 和( 以，九) 分别为高斯权数和高 斯点。 第二章数学描述与数值方法 2 3 2 弱奇异积分 弱奇异积分对应于方程( 2 3 8 ) 中当j = l 时的d o 。对于这种情况，本文采用子 单元法 2 0 】l 来求积分( 参见图2 3 ) 。将四边形单元拆成两个三角形单元i 和i i ，对于 每个子单元，先将坐标转化为极坐标，消除掉1 r 的弱奇异性，然后再将极坐标变 换成可以应用高斯积分的坐标，最后用高斯积分得到结果。下面以子单元i 为例 说明。 子单元内极坐标与母单元的坐标关系如下 孝= p c o s o l r = p s i n 矽一1 d 喜d r l = p a p a o( 2 4 7 ) p ，0 的变化范围为 0s 9 霓| 冬 0 p 2 c o s t 9 定义高斯积分坐标( 孝，d 刁) ，极坐标与孝7 ，r 的关系为 秒= 詈( 1 + 孝7 ) p = ( 1 + - ) _ l ( 2 4 8 ) c o s 专( 1 + 善) 经过变化，孝，r 的变化范围为一】孝1 ，一1sr l ，这样就可以利用高斯积分 进行数值积分。将( 2 4 8 ) 代入( 2 4 7 ) 可得最终的坐标转化关系 孝= r r =