（水声工程专业论文）有限长双层圆柱薄壳结构振动有限元方法研究.pdf

哈尔滨工程大学硕十学位论文 a bs t r a c t t h er e s e a r c ho nt h ef e mf o rc a l c u l a t i n gt h ev i b r a t i o no ft h ed o u b l e c y l i n d r i c a lt h i ns h e l l s t r u c t u r eh a sm a i n l yb e e nd o n ei nt h ed i s s e r t a t i o n t h e e l a s t i cc y l i n d r i c a ls h e l li st h et y p i c a ls t r u c t u r eo fu n d e r w a t e rv e h i c l e s i ti s i m p o r t a n tt os t u d yt h ev i b r a t i o no fu n d e r w a t e re l a s t i cc y l i n d r i c a ls h e l li np r a c t i c e a tp r e s e n t ，t h er a d i a t e ds o u n df i e l df r o mu n d e r w a t e rs t r u c t u r ev i b r a t i o nc a l lb e c a l c u l a t e db yt h ec o m b i n a t i o no ff e ma n df e m + b e m w h e nt h es o u n df i e l di s c a l c u l a t e db yf e m ，t h ew h o l ec a l c u l a t e dr e g i o nh a sb e e nm e s h e df u l l y , w h i c h m a k e st h ec a l c u l a t i o na m o u n th u g e w h e nt h es o u n df i e l di sc a l c u l a t e db y f e m + b e m ，t h er e g i o no fe l a s t i cs t r u c t u r ec a l lb em e s h e do n l y c o m p a r e dt of e m ， f e m + b e mh a sag r e a ti m p r o v e m e n t h o w e v e r , t h e r ei sn oa c o u s t i cf i n i t ee l e m e n t t y p ei nb e m s ot h es t r u c t u r ev i b r a t i o np r o b l e m ，w h i c ht h ef l u i de x i s t si nt h e d o u b l es h e l l s c a nn o tb ec a l c u l a t e db yb e m o n l y t h e r e f o r e 3 - da n d2 - df l u i d s o l i dc o u p l i n gf e ma l g o r i t h mh a v eb e e n d e v e l o p e di nt h ed i s s e r t a t i o n m e a n w h i l e ，t h em e t h o dt os o l v et h el a r g es p a r s e m a t r i xe q u a t i o nh a sb e e ni n v e s t i g a t e d t h ed e v e l o p e dt y p eo ff i n i t ee l e m e n t si sa sf o l l o w s ：f i r s t ，3 一ds h e l le l e m e n t s ， i n c l u d i n gt h et r i a n g u l a ra n dq u a d r i l a t e r a le l e m e n t s ；s e c o n d ，s o l i d a n df l u i d e l e m e n t s ，i n c l u d i n gh e x a h e d r o n ，t e t r a h e d r o na n dp r i s me l e m e n t s ；t h i r d ，2 一d a x i s s y m m e t r i ce l e m e n t s t h e s et y p e so ft h ee l e m e n t sc a nb eu s e dt os o l v et h e v i b r a t i o np r o b l e mo fc o u p l i n gs t r u c t u r e sc o n s i s t i n go ft h ed o u b l es h e l l ，t h er i b si n t h es h e l l ，t h ef l u i db e t w e e nd o u b l es h e l l s t h en u m e r i c a lv a l i d a t i o nh a sb e e nd o n eo nt h ed e v e l o p e df i n i t ee l e m e n t a l g o r i t h mi nt h ep a p e r t h er e s u l t sa r ec o n s i s t e n tw i t ht h a ts o l v e db ya n s y s a t t h es a m et i m e t h en u m e r i c a lc a l c u l a t i o no ft h es p h e r i c a lm o d e l sv i b r a t i o nh a s b e e nd o n eb yt h es e l f - d e v e l o p e da l g o r i t h m c o m p a r e dt ot h ea n a l y t i c a ls o l u t i o n ， t h ee r r o ro ft h es e l f - d e v e l o p e da l g o r i t h mi sl e s s k e yw o r d s ：u n d e r w a t e ra c o u s t i c s ；f e m ；d o u b l ee l a s t i cs h e l l ；v i b r a t i o n 哈尔滨工程大学 学位论文原创性声明 本人郑重声明：本论文的所有工作，是在导师的指导下，由 作者本入独立完成的。有关观点、方法、数据和文献的引用已在 文中指出，并与参考文献相对应。除文中已注明引甩的内容外， 本论文不包含任何其他个人或集体已经公开发表的 乍晶成果。对 本文的研究做出重要贡献的个人和集体，均已在文中以明确方式 标明。本入完全意识到本声明的法律结果由本人承担。 作者( 签字) ：铆翟 弱期：渺譬年3 月铷e t 7 哈尔滨工程大学 学位论文授权使用声明 本人完全了解学校保护知识产权的有关规定，基i 研究生在校 攻读学位期闻论文工作盼知识产权属于哈尔滨工程大学。哈尔滨 工程大学有权保留并向国家有关部f 1 或机构送交论文的复印件尊 本人允许哈尔滨工程大学将论文的部分或全部内容编入有关数据 库进行检索，可采用影印、缩印或扫描等复制手段保存和汇编本 学位论文，可以公布论文的全部内容。同时本人保 正毕业后结合 学位论文研究课题再撰写的论文一律注明作者第署名单位为哈 尔滨工程大学。涉密学位论文待解密骺适用本声明。 本论文( 臼在授予学位后鼓瑟可圈在授予学位1 2 个月后图 解密后) 由哈尔滨工程大学送交有关部门进行保存、汇编等。 作者( 签字) ：针楚导师( 签- 7 - - ：纫丢： l 圜期：扩鼍每呈月挣e l扣尹年- 3 月妇露 l 一 哈尔滨工程大学硕十学位论文 第1 章绪论 1 1 论文背景及其理论模型的提出 水下航行器的水下声学特性一直是人们所关注的问题，因为它直接影响 其水下作战性能。随着电子科学技术的发展，声纳系统趋于低频、大功率的 方向发展，作用距离越来越远，声纳对水下目标的识别和跟踪能力大大增强。 因此，降低水下航行体的辐射噪声不仅可以提高其自身的隐蔽性，提高水声 对抗能力；同时，还可以增大航行体自身声纳系统的作用距离，从而大大提 高水中兵器的战斗力。显然，加强对我国水下航行体降噪的研究工作将有着 十分重要的现实意义。 在水下航行体整体设计阶段，如何能对航行体的整体噪声进行预报，是 获得安静型水下航行体的前提，以确定这些航行器在方案设计中就具有噪声 方面的优势，达到较好的隐身效果；水下航行体噪声源大致可分为结构噪声、 螺旋桨噪声和流噪声。其中，在低速航行时，由机械动力源激起壳体振动所 产生的结构噪声是水下航行体辐射噪声的主要部分。因此，降低结构噪声对 安静型水下航行体的设计十分重要的。 但我国在水下航行体减振降噪方面的工作还相对落后。对水下航行体的 辐射噪声精确预报是一个复杂的系统工程，因为其牵涉到水下弹性结构与水 介质的耦合振动问题及其声场特性的预报问题，涉及结构的弹性振动、声学 和结构与水介质的耦合振动等多方面的理论知识和计算技术。五六十年代， 计算机技术远没有像今天这样发达，其计算速度和计算能力很差。当时无法 对不规则复杂结构的振动和声场特性进行数值预报，只能对具有正交坐标曲 面的结构，如球、无限长柱等，在高低频近似的情况下用特殊函数展开或者 围道积分进行解析计算，这在工程上的实用性限制很大。随着计算机技术的 飞速发展，其计算速度和计算能力大大提高。差分法、有限元法、边界元法 等数值计算技术也迅速发展，利用计算机数值求解一般弹性结构的振动和声 辐射技术也日趋成熟。哈尔滨工程大学水声工程学院曾承担了“九五国防重 点预研课题的研究任务。该课题的主要目的是要搞清楚单、双壳体结构表面 振动及声场的分布特性及其相关特性，并给出单、双壳体结构受激振动时水 哈尔滨丁稗大学硕十学位论文 ；i 1i i ；葺 下辐射噪声的估算方法。从而可以为安静型水下航行体的声学设计和现役的 其它水下航行体振动与噪声治理工作提供参考。 目前，已经出现了许多较成熟、通用的大型数值计算软件，如美国a n s y s 公司开发的有限元软件a n s y s ”哪和比利时l m s 公司开发的有限元+ 边界元软 件s y s n o i s e ”引，等等。数值计算已经成为复杂结构振动和声学特性的工程计 算设计的重要技术之一。但是在使用现有的商业有限元软件和边界元软件进 行计算时，发现这两种软件各有优缺点，有限元软件只可以计算小型的双层 壳模型，边界元软件不能计算双层壳模型，因为其没有流体单元；从短期来 看，可以使用现有的软件计算研究，但从长远来看，要计算水下航行体的振 动与噪声辐射情况，尤其要预报水下航行体的实际模型的噪声辐射情况，必 须要有自己的核心的算法，而不是仅仅使用商业软件简单的计算；鉴于此， 自行开发流固耦合有限元算法，就显得至关重要。 1 2 弹性结构水下振动声辐射研究方法概述 1 2 1 解析方法 解析方法”1 只适用于平板、圆柱壳、球壳、椭球壳等具有正交几何表面 的简单结构进行水下振动和声场求解。优点是计算速度快，预报结果精度高。 但不适合于解决实际复杂结构的水下声辐射问题。 1 2 。2 数值方法 数值方法包括有限元法、边界元法、d a a 2 法吲、有限元边界元法、 统计能量法叫2 1 。 ( i ) 有限元法( f e m ) 2 0 世纪5 0 年代，有限元法一经问世，就显示出其巨大的优越性，迅速 被应用于声辐射问题的分析计算。有限元法的优点是适用于求解有限区域内 声场，对结构边界的形状适应能力比较强。 薄壳是工程中常用的一种特殊结构。薄壳是由薄膜状态和弯曲状态组成， 所以对薄壳问题的分析比单纯对薄膜问题和薄板问题的分析更为复杂，这也 使得壳体问题的有限元分析至今仍然是有限元分析领域的重要课题，国内外 学者对薄壳的有限元分析做了大量工作。薄壳有限元分析中的薄膜有限元部 哈尔滨：r 稗大学硕十学位论文 分是有限元领域较难处理的问题，在单独作薄膜单元分析时不需要考虑绕着 膜单元法线转动的自由度，但在作壳体有限元分析时需要考虑这个自由度， z i e n k i e w i c z ，o c 1 9 1 在其所著的书中提出直接在单元刚度矩阵中的对应位置 处加入一个数值作为虚假的转动刚度；y u n u s ，s h a hm ，p a w l a k 、t i m o t h yp 和 c o o k ，r d b 4 1 在他们发表的文献中用弹簧类型的刚度模拟转动自由度所对应 的刚度，此种方法在实际中得到广泛应用，但其存在另一个致命缺点，那就 是引入这种弹簧类型的转动刚度时使单元矩阵奇异，奇异的后果是方程的解 不唯一，k a n o k n u k u l c h a i 瞄副在他发表的文章中提出在单元刚度矩阵对角线位 置加上一个小的罚刚度可以消除矩阵的奇异性；w i l s o n ，e l p s - 3 7 1 等人提出能 够改善常规线性单元性能，特别是对其用于弯曲应力状态性能改善具有重要 意义的w i l s o n 非协调元，等参元具有良好的适应性和表达格式的简明性，因 而得到广泛应用，但其精度和效率仍是不够高，w i l s o n 非协调元的引入对改 善等参元的计算精度和提高计算效率是很有意义的，特别是对三维问题的有 限元分析，但w i l s o n ，e l 所提出的非协调元只有在单元形状是平行四边形 时雅克比矩阵行列式是常数，此时单元能够通过分片试验p 1 ，如果在实际的应 用中如果只能是平行四边形单元，那是很不方便的，且严重限制了这种单元 的使用，随后w i l s o n ，e l 提出对雅克比行列式的改进方法，改进后的w i l s o n 非协调元虽然可以在不限制单元形状的情况下可以通过分片试验，这种方法 是强迫不规则单元的等参变换以一变为规则单元时的情形，雅克比逆阵是将 整体坐标系下物理空间中真实畸变单元，映射到自然坐标系下，采用改进后 的w i l s o n 非协调元后，得到的局部坐标系下的单元已不再是原畸变单元的真 实映射，这有损于单元公式表征畸变单元真实外形拓扑性质的能力，另外该 法使雅克比变换为常量，已不能真实反映原畸变单元的真实微分面积，因此， 采用该方法时，虽然不规则单元可以通过分片试验，但不管是从物理上( 等参 变换) 还是从数学上( 数值积分) 讲，都以单元公式表征畸变单元真实特征的能 力受到损失为代价的；因此鹿晓阳旺副等人提出了对w i l s o n 非协调元中的雅克 比变换的改进方法，在计算雅克比矩阵时采用了求逆犯规法，即所谓的“求逆 犯规胆习”，这实质是放松了w i l s o n 提出的满足分片试验的条件，经验证该改 进单元有良好的精度；对于膜单元中转动自由度的处理a l l m a n ，d j 1 2 刨提出了 一种独特的方法，该法更能模拟实际的平面内的转动刚度，在a l l m a n ，d j 哈尔滨j r 程大学硕十学位论文 所提出的三角形单元的基础上c o o k ，r do 2 7 1 提出了捌有a l l m a n 转动自由度类 型的四边形膜单元，但是拥有a l l m a n 类型转动自由度的单元具有两个多余的 零能模式：1 ) 在平面内的旋转自由度为常数，2 ) 由减缩积分引起的著名的沙 漏模式( 三角形单元中不存在这个零能模式) ，为了消除这两个零能模式， m a c n e a l ，r h 1 2 8 在其发表的文章中提出了消除这两个零能模式的方法，效果 良好。弯曲单元部分是壳体单元的另一重要部分，对于这部分，国内外学者 也提出了不少单元类型，如基于薄板理论的协调元与非协调元、位移和转动 各自独立插值的板单元、应力杂交板单元，由于薄板这种特殊的结构，使得 这些板单元或多或少都会存在数值上的问题，b a t o zj - l 1 2 0 l 等人提出了基于离 散k i r c h h o f f 理论( d k t ) 的三角形薄板单元，该板单元采用了独特的理论，使 单元的精度和效率均比较好，在实际中常常被采用，只要将d k t 单元简单处 理就可以获得基于离散k i r c h h o f f 理论( d k q ) 的四边形薄板单元。 与壳体单元相比，实体单元算法比较简单，一般文献p 1 上介绍的方法已 能满足工程中的实际需求。 但是有限元法在声辐射分析计算中的不足之处有：1 ) 对于二维声辐射问 题，有限元方法和有限差分法一样，需要在整个分析域内进行单元剖分、变 量插值等，分析自由度庞大；2 ) 对于工程中常见的在无限域中的外部声辐射 问题，有限元方法的剖分截止边缘难以确定，并会由此带来计算误差。 ( 2 ) 边界元法( b e m ) 边界元法的研究开始于2 0 世纪6 0 年代。它在域内采用解析表达式，只 在求解域的边界上进行离散，是半解析半数值方法。优点是：1 ) 使数值计算 的维数降低了一维，从而减少了问题的自由度和原始信息量；2 ) 适用于求解 无限区域的声场。缺点是：1 ) 存在奇异积分和特征频率处解不唯一的问题； 2 ) 适用于单频计算，由于边界元求解方程中的系数矩阵的元素是与频率相关 的，对每一个不同的激励频率，系数矩阵都需要重新进行计算，这使得在较 宽频率范围内进行结构声学分析时计算量非常大。 ( 3 ) 双重渐近近似法( d a a 2 ) 2 0 世纪7 0 年代初期，g e e r s 等人相继提出了d a a l ，d a a 2 的方法。该 方法的优点是：1 ) 在高频和低频逐步趋向精确解，中频段是一个光滑的过渡， 除谐振点附近误差较大外，整个频段的精度都比较好：2 ) 适用于描述非正规 4 哈尔滨：r 程大学硕十学位论文 几何形状表面结构的声阻抗；3 ) d a a 2 法的系数矩阵与频率无关，使数值计 算的工作量和数据存储量大大减少。国内，张敬东m 刎利用d a a 2 法研究了水 下球壳和椭球壳受激振动的声辐射问题。结果证明壳体表面形状为球形时的 结果可靠性要好于椭球壳。在本征值附近，d a a 2 法的结果是不可靠的。 ( 4 ) 有限元+ 边界元方法( f e m + b e m ) 该方法用有限元法描述复杂弹性壳体的结构振动，用边界元方法对结构 的外表面进行划分，使求解的维数降低一维，并且对结构表面的几何形状无 严格的要求。其优点是对结构形状变化适应性强，计算量较小，适应于中低 频范围有界域和无界域内的声学分析。缺点是存在频率失效问题，在处理实 际工程结构宽带或高频范围的声学分析存在一定的困难。 ( 5 ) 统计能量法 该方法主要适合于高频段，对有规则结构的振动与声学分析比较有效， 计算量也小于有限元+ 边界元法，对复杂结构的分析还可借助于有限元方法。 但统计能量分析法在实际应用中的问题在于耦合损耗因子和模态密度的确定 比较困难，还需要进一步完善。 1 3 本文研究方法及研究内容 水下结构有单壳体和双壳体两种。其典型结构形式是圆柱壳。对无肋单 壳体的水下振动和声学特性的研究文献较多，研究方法也很成熟。对于加肋 单、双壳体，特别是加肋双层壳体情况，由于其结构复杂，无论是理论计算， 水下振动与声学的实验测量都很困难，但在工程上却有重要现实意义。实际 工程设计上，一般是利用环肋、纵肋将两层壳体焊接成一体的，且双壳之间 通常含有流体。计算这种复杂结构的水下振动和声学特性一般地只能采用数 值法。 1 3 1 研究方法 本文研究的理论模型有如下特点：复杂的双层加肋壳体结构：无限域中 的流体与结构的耦合振动。这就要求研究方法对复杂结构有较强的适应性， 并且能便于求解无限域中的声辐射问题。由于有限元算法和有限元+ 边界元 算法比较成熟，并且有通用的商用计算软件，所以常采用有限元加边界元结 哈尔滨t 稗大学硕十学位论文 合计算。以前的做法是：首先利用有限元软件计算加肋双层圆柱壳与水下声 场的耦合振动。然后，将该结构的外表面位移提取出来，将其输入给边界元 软件，利用其边界元技术来计算结构的辐射声功率、近场声学量分析、结构 近场声能流分布图形、声辐射效率、远场指向性等声学量。由于有限元软件 利用有限元法计算无限区域内的弹性结构与流体介质的耦合振动和声辐射 时，结构外部的应该包上一层流体区域，在流体区域的外层再包上一层吸声 单元，对于三维问题，计算的节点过多；边界元软件不能处理双层壳间含流 体的加肋柱壳的水下振动与声场的耦合问题；现在的想法是，自行开发双层 壳间含流体的加肋柱壳的流固耦合有限元算法，在该方程中加入边界元表面 方程，这样就可以获得双层壳水下振动时的壳体表面位移，再利用边界元的 场点方程就可以计算水中的任一点的声场；使用该法不用在双层壳的外表面 包上流体区域，所以流固耦合有限元方程中的矩阵维数会变小，减小了计算 量，同时将有限元法和边界元法统一在一个方程中，使计算过程简化。本文 的主要工作是建立流固耦合有限元算法。 1 3 2 研究内容 本文利用有限元和边界元方法研究了水下双层加肋圆柱壳体受激振动和 声辐射问题，主要工作内容如下： 第1 章介绍本选题的背景、意义、研究现状、研究设想及预期结果。 第2 章介绍有限元与边界元基本理论；并介绍了如何用边界元表面方程 修正流固耦合有限元方程。 第3 章介绍具体内容有：( 1 ) 介绍具体的有限元单元算法，包括：三维壳 体单元、三维实体单元、三维流体单元；( 2 ) 介绍流固耦合矩阵的计算方法； ( 3 ) 介绍整体坐标系与局部坐标系之间变换方法；( 4 ) 介绍两种由单元矩阵形成 整体矩阵的方法；( 5 ) 介绍大型稀疏矩阵方程的快速求解方法；( 6 ) 简单介绍二 维流固耦合有限元方程的形成过程。 第4 章验证第3 章中的算法的正确性，将自行开发的算法与a n s y s 计算 结果对比，同时将自行开发的算法与球壳解析解对比。 哈尔滨t 程大学硕士学何论文 第2 章有限元和边界元基本理论 本章首先介绍流固耦合有限元理论和边界元理论，接着将有限元表面方 程加入流固耦合有限元方程，求解修正后的流固耦合有限元方程，得到结构 外表面的法向振速，再通过边界元场点方程，便可以获得水中任一点的声场。 2 1 流固耦合有限元基本理论 有限元方法计算结构振动与声耦合问题时，流固耦合系统模型示意图见 图2 1 ，对结构和流体都进行有限元网格离散，建立有限元方程进行分析。 鼠 ：流体区域：k 弹性结构；s o ：流同交界面；咒：力边界条件：s ，位移边界条件 图2 1 流固耦合系统模型示意图 真空中含有阻尼的弹性结构振动的有限元方程为p 1 ： 鸠 矽 + 【c ：】 矽 + 【k 】 u = f ) ( 2 1 ) 其中：【必】为结构的质量矩阵；【疋 为结构的刚度矩阵； c 、】为结构阻 尼矩阵； 只 为结构载荷力r e - 量。 在结构参数( 包括几何参数和材料参数) 给定，及有限元网格生成后， m 。】、【k 】、【c 】已经完全唯一确定了， 只 是载荷力向量，它的不同决定 了不同分析类型： 当 只) = 刑) 时， e ) 为时间的任意函数，属于瞬态结构分析，依赖方 程( 2 1 ) 式； 当( 只) = r p e 删) 是时间的简谐函数，则对应的分析类型为谐波分析， 依赖的方程可化为： ( 一彩2 鸠】+ 弘 c 4 - k ) 甜) - 昂g )( 2 - 2 ) 7 当 b ) = 0 ，属于模念分析，即计算 m ，】 田“e 】 西 “k 。】 “) = 0 有 非零解的本征值问题。取 c 、 - 0 ，则有 阻。 妙 + 医。】砂) = 0( 2 3 ) 利用本征向量的正交性，可以给出结构的本征模态和对应的本征值、模 态频率。 当 f ) _ e ) ，即 竽 _ o ，则属于静力分析，根据准静态原理求解力 作用终点时的状态 ) = 西) _ 0 ，依赖的方程简化为 k 甜) - 圪)( 2 - 4 ) 当弹性结构置于声学介质中时，在流体与结构的交接面s 上流体与结构 之间存在着相互作用口钆 1 ，声压对结构产生一个面力的作用。根据虚功原型3 1 可将该面力等效移置到单元节点上，于是，弹性结构与声场的耦合振动矩阵 方程为： m 。 矽) + c 】 d + 【k 。 u = 只) + ) ( 2 5 ) 式中： 弓 = 旧7 p ) 。 将( 2 - 1 ) 、( 2 - 5 ) 两式联立，可写成统一的矩阵方程如下： 搿黜刚葛) 1 黜1 铷刚黜陋6 ， 式( 2 6 ) 为流固耦合有限元方程，求解流固耦合有限元方程就可以获得弹 性结构上的位移自r a 度和流体节点上的声压自r a 度。 2 2 边界元基本理论 假设一振动体位于声速为c ，密度为p 的无限流体介质中，记表面区域 为s ，包围的区域为口，外部区域为见，如图2 2 所示。则在外部区域声场 满足的波动方程为： v 2 p c i = 吉鼍笋 ( 2 - 7 ， 其中：p ( r p ，) 一瞬时声压；f 一时间变量；c 一声速； v 2 一拉普拉斯算 子： 哈尔滨二r 程大学硕十学位论文 图2 2 辐射声场示意图 对于谐和声场，取e 删；对于非谐和声场，将( 2 7 ) 式作傅立叶变换，可 得到h e l m h o l t z 方程： v 2 p ( r ，国) + k 2 p ( r p ，缈) = 0( 2 - 8 ) 式中：k = 竺为流体介质中波数；c o 为角频率； c 方程( 2 8 ) 的解可表示成积分形式，即声场中任一点p 的声压可表示成 h e l m h o l t z 积分方程： 州驴【掣g ( h 掣姗 其中：表示s 面上外法线方向偏导数；为s 面上的点； d 行 区域的点；g 为自由空间的格林函数，表达式为： 一肚i i i i g ( 。啪2 褊 ( 2 9 ) ，。为外部 ( 2 1 0 ) 式( 2 9 ) 中口的取值与场点的位置有关，式( 2 - 11 ) 列出了当场点在s 面外、 s r 面内与s 面上等情况下的值。 4 7 r 0 2 n q f 在s 面外 ；套s s 誓空，且s 为平滑面 c 2 一， 芦在面上，且s 为平滑面 、 f 在s 面上，且s 为非平滑面 利用表面s 上声压与法向振速的关系式，则有 害= 一j p c k u 。 ( 2 1 2 ) 锄 ” 、 将( 2 1 2 ) 式代入方程( 2 - 9 ) 中有 9 哈尔溟t 程大学硕十学位论文 唰【p ( i ) 艺掣+ ( _ j p c k u 加( 两筇( i ) ( 2 - 1 3 ) 式( 2 1 3 ) 是研究结构表面声场与外部声场的基本关系式。 利用边界元技术，对式( 2 1 3 ) 离散化，可得如下矩阵表达式 m _ b r o ) ( 2 - 1 4 ) p ) = 【c p o + 【d 】舷 ( 2 一i s ) 其中： ， 圪) 为源面上声压和法向振速向量； p ) 为待求场点的 声压向量； 式( 2 1 4 ) 给出了源表面上声压和法向振速之间的关系，由式( 2 1 4 ) 得到： e ) = 彳b 】 k ( 2 1 6 ) 将式( 2 1 6 ) 代入式( 2 1 5 ) 中，得 尸 = ( 【4 】叫【b 】【c 】+ d 】) k ( 2 1 7 ) 因此如果声源表面上的法向振速向量己知，就可计算出声场中声压 户 。 2 3 有限元+ 边界元方程 双层壳模型见图2 3 ，该模型的流固耦合有限元方程可用式( 2 6 ) 表示，式 ( 2 6 ) 所表示的流固耦合有限元方程是双层壳模型( 壳间含水) 在真空中的情 况，如果将双层壳放在水中，双层壳振动时，壳外表面会对水有作用力，同 时水也会对外壳有反作用力，此时要对式( 2 6 ) 右端载荷项进行修正，即在右 端载荷项加上边界元的表面方程式( 2 - 1 4 ) ，则忽略式( 2 6 ) 中阻尼项并在式( 2 6 ) 的右端载荷上加上式( 2 1 4 ) 得式( 2 18 ) 。 4 也 螗、“、峨虬 1卜 d i l 卅1 图2 3 缩比模型横剖面图和纵剖面图 i m p 尚懈 + i 篙1 - k r 1 j 1 lf l 尸) v j ， = i + t a - l 【础圪)尺 m 舱) + m k ， j 1 尸) - p 儿 ( 2 - 1 8 ) 式( 2 1 8 ) 中的 k ) 为法向振速，而 u 为结构上的位移自由度，所以要将 k 转化为位移形式，易知法向位移与法向振速的关系为 乩) - i 1 圪) ，再 将 u 转化为 u ，假设壳体表面第i 个节点上的外法线矢量为 = ， ，则第i 个节点上的x ，y ，z 方向上的位移与该点上的法向位移 系为 = h “纱 “，：j i 船l k j ( 2 1 9 ) 则壳体外表面的所有节点x ，y ，z 方向上的位移与对应点上的法向位移关 则 = 去 = 7 i 乏 _ - ： ： l 乏j2 ，缈1 ： ：乏j ( 2 - 2 0 ) = 归 胛 ( 2 2 1 ) 其中， 甜) 为外壳与流体耦合的节点上x ，y ，z 方向上的位移向量。 咖附缸! 咖如如 h 砂 船； 血 h h； 三怠， 昱： 鬈； + 高1i 譬j ：； = l 心i f o ) i + 国c 彳】- l c b ，t 门，t “， 【- p ， r 】 m ， j 1 户 f + l 【o 】 鬈， j 1 尸) ，2 + j 国。月1 。d 1 门“ ( 2 2 2 ) 将式( 2 2 2 ) 的右端项中含有位移的部分移到左端，并将相同的项合并，就 可得到有限元加边界元方程，此时解式( 2 2 2 ) 就可得到壳体外表面的位移，再 用场点方程式( 2 1 7 ) 就可以解出水中的任一点的声场。 2 4 本章小结 本章首先介绍了流固耦合有限元的基本原理和基本关系式，接着介绍边 界元的基本理论，利用结构表面法向节点位移或法向节点振速计算场点中的 声压，最后介绍了有限元加边界元方程的形成过程。本章的主要目的是介绍 有限元加边界元法的实现过程，在本文的第三章、第四章中并没有对边界元 算法过多涉及，主要工作还是集中在流固耦合有限元算法的丌发上。 哈尔滨t 稃大学硕十学位论文 第3 章流固耦合有限元算法及其求解方法研究 在第二章中由流固耦合有限元方程和边界元表面方程导出了有限元+ 边 界元算法的最终表达式( 2 2 2 ) ，流固耦合有限元方程式( 2 6 ) 是式( 2 2 2 ) 的重要 组成部分，对于流固耦合有限元方程式( 2 6 ) ，只要能将式( 2 6 ) 中的整体刚度 矩阵和整体质量矩阵求出，其它物理量将很容易获得，本章将具体介绍如何 求出整体刚度矩阵和整体质量矩阵，并介绍大型稀疏矩阵方程求解方法。 本章的主要研究内容如下： ( 1 ) 介绍三维壳体单元、实体单元、流体单元的具体算法。 r 2 ) 介绍三维流固耦合矩阵计算方法。 ( 3 ) 介绍三维结构整体坐标系与局部坐标系之间变换方法。 ( 4 ) 介绍两种由单元矩阵形成整体矩阵的方法。 ( 5 ) 介绍大型稀疏矩阵方程的快速求解方法。 ( 6 ) 简单介绍二维流固耦合有限元算法。 3 1 三维模型单元 3 1 1 壳体单元 三维壳体单元由膜单元和弯曲单元组合而成，如果只考虑膜单元，则单 元每个节点有两个自由度：分别沿着x ，y 轴的位移自由度，如果只考虑弯曲 单元，则每个单元有三个自由度：沿着z 轴的位移，以及绕着x ，y 轴的转动 自由度；将膜单元和弯曲单元组合在一起可形成壳单元，壳单元有六个自由 度，除了上面介绍的自由度外，还应该在膜单元部分增加一个绕着z 轴转动 的自由度，也就是这个绕着三轴转动的自由度的存在；使得壳体单元的膜部 分，有多种单元类型与之对应，且给实际应用中带来诸多困难，在本文中只 详细介绍其中一种，其它两种膜单元类型只作简单介绍；3 1 1 1 和3 1 1 2 分 别介绍三角形壳单元和四边形壳单元；四边形壳单元模型和三角形壳单元模 型如图3 1 。 哈尔滨j i - 仟i i 火学硕士学位论文 四边彤单兀蜕化的三角形单元 曷为单元局部坐标系 图3 1 四边形壳单元和三角形壳单元模型 3 1 1 1 三角形壳体单元 ( 1 ) 首先介绍三角形壳体单元的膜部分，然后介绍三角形壳单元的弯曲 部分，膜部分如图3 2 ，节点自由度为“，v ，伊。 x u v 图3 2 三角形壳单元膜部分 三角形单元的膜部分位移插值函数矩阵表示形式为 r “ u = ll = n a ( 3 - 1 ) l v j 其中： a = v l “2v 2 甜3v 3 ： ( 3 2 ) 一i nm 0 2m 0 3 m oi ( 3 - 3 0 ， l 0 0 7 v ，i 7 i = 厶，n 2 = l 2 ，3 = 厶，厶，f = 1 ，2 ，3 ( 3 4 ) a 为节点位移向量，为插值形函数矩阵表示，厶，f _ 1 ，2 ，3 为三角形面 积坐标。 将上面的位移插值函数的矩阵形式代入几何方程中，可以表示为 占= its y f = 上u = l n a ( 3 - 5 ) 1 4 哈尔滨t 稗大学硕十学位论文 其中： 三= 瑚 ( 3 - 6 ) ( 3 7 ) b 为应变矩阵。 应变矩阵b 中是形函数对全局坐标的导数，为了计算刚度矩阵和质量矩 阵时可以方便的使用高斯积分，可使用等参变换法将应变矩阵b 变换为形函 变换方法如下 厶= 1 一孝一r ，掣：掣一一a n ，i = 1 , 2 , 3 却砚，弘 形函数对局部坐标的导数写成矩阵形式为 雅克比行列式为 ( 3 8 ) ( 3 9 ) 一1 1j ( 3 - 1 0 ) i 州鹏 u ，一t a 善 a 孝 l 州鹕 l - a 7 7a 7 7 ( 3 1 1 ) ( 3 1 2 ) 为了编程上的方便，把应变矩阵b 和形函数的矩阵形式分别改写成如式 1 5 m 0 0 m m 0 纠到 o m o a 一砂 7 o l a 一缸 o a一砂a一缸 o a 一砂 a 一苏 o 。，。l ，7、一ru 凯邓叭一弘 一三 |垂“叭一弘 坐=吐7，黜铲叭一西 自令 则 叉 毁 o 1 n一 盟鸳丛却鹕一曾盟却州一鹫盟却 r_，j m 儿儿毛恐恐 。l 旷p 盟西鹏一却 巩一管盟却 哈尔滨7 - 1 ：1 7 1 大学硕十学何论文 ( 3 - 1 3 ) 和式( 3 - 1 4 ) 形式 b = o n , oo 0 n 2 墩出 o o n 1oo o y o n i o n 1 0 o n e _ 砂 苏 砂 ：l l o o 2 oo o n 3oo o x o n 2 o0 o n 30 砂砂 o n 2 o o n 3 o n 30 西x瓠 巩 00 n 3 0 0 n 2 00 n 3 0 ( 3 1 3 ) ( 3 - 1 4 ) 由有限单元法理论口1 知，刚度矩阵为 l k m l = b r d b jt d孝dri(3-15) 一i i 质量矩阵为 ll 帆，= 尸ff n r n l j ld f d v ( 3 1 6 ) 一l l 其中，d 为弹性矩阵p 1 ，p 为材料的密度，f 为材料厚度。 式( 3 1 5 ) 表示的刚度矩阵的计算方法只是每个节点有两个自由度的情形， 实际上膜单元在和板单元组成壳单元时，膜单元还需要一个垂直与膜平面的 转动自由度，与该转动自由度对应的刚度称之为转动刚度，用弹簧类型口4 1 的 刚度代替；此时的单元刚度矩阵存在奇异，所以为了消除奇异所带来的数值 不稳定还要在单元刚度矩阵中加入一个小的罚刚度【25 1 。 以上介绍的是双线性插值函数没有引入w i l s o n 非协调元的情况，引入 w i l s o n 非协调元的情况见文献口卜2 孙；如过膜单元的转动部分采用a l l m a n 类型 的刚度，可见参考文献仁融引。 ( 2 ) 三角形壳单元的弯曲部分如图3 3 ，节点自由度为w ，或，0 。，该弯 曲单元部分采用d t k 胆训单元，如图3 4 。 根据基尔霍夫假设，平板弯曲问题假设为二维问题，全部应力和应变可 以用板中面的挠度w 表示p ，即 u ( x ，y ，z ) = 一z 竽 ( 3 1 7 ) ( y ( x ，y ，z ) ：一z o _ w( 3 1 8 ) 缈 1 6 哈尔滨j 厂程大学硕十学佗论文 w ( x ，y ，z ) = w ( x ，y ，0 ) 兰w ( x ，y )( 3 - 1 9 ) 图3 3 三角形壳单元的弯曲部分图3 4d k t 单元 尾用差值函数的形式表示为 66 | 3 x ：z n , f l ，l ，p y = n i $ y | i = l t = l i 1 = 厶( 2 厶一1 ) ，2 = 厶( 2 l 2 一1 ) ，3 = g ( 2 9 一1 ) 【4 = 4 ：厶，5 = 4 厶厶，6 = 4 厶厶 上l + 厶+ 厶= 1 厶，厶，厶为面积坐标 由基尔霍夫假设， ，2 ( 3 2 2 ) ( 3 2 3 ) ( 3 - 2 4 ) 每个角节点和边中节点需满足剪切应变 = 0 ( 在角点1 ，2 ，3 处)( 3 2 5 ) 屈+ _ o w = o ( 在边中节点4 ，5 ，6 处) ( 3 2 6 ) o s w 沿着单元边成二次方变化 l x 塑鲫玎嘞名矧渐 盟矿 协， 鱼谚 盟矿枷 即 挑一酽 _ ， 删 一拈零 蛐 拈 置缸 j 1lllll 自，- ，- 一雕挑一缸挑一砂 单 + + 在成 膨 哈尔滨一l ：科大学硕十学位论文 即 c 芸，。一亳w ，一言c 罢h 未，一去c 罢) ， 其中，七表示边l ，中点，t ，表示边扩的长度。 沿着边切线方向的转角成线性变化 i b n k = 三( 、p m + p n 、) 节点位移定义为 u2 l w i 或g - w 2 最z 吼z 比以，g ，j 在每个边上的转动关系为 ； 挑 a s 跏 a 门 = 匕 础 北 ( 3 2 7 ) ( 3 2 8 ) ( 3 2 9 ) ( 3 3 0 ) ( 3 - 3 1 ) 弘c o s ( x ，) ，j = s i n ( x ，2 ) ( 3 - 3 2 ) 因为忽略了剪切变形，所以节点的转动自由度满足基尔霍夫边界条件， 乱：坐，臼 = 一竺 砂 7 3 x 由式( 3 1 7 ) 一式( 3 - 3 3 ) ，展，。可表示为 刖糍钟 其中，钟( 工：，l 3 ) ，了( 上：，厶) 详细表达式见参考文献啦0 1 。 由位移应变关系为 b ( 厶，9 2 ) ：石1 z 以 忌= b u y 3 、屹? 2 + y 、2 峨b 一玛，彰上一x 。：左 一屯。磁上一五：h 互如+ y 3 。h 互岛+ m ：月互岛 其中，x 0 = x ，x ，= y j y s 。 则d k t 单元的刚度矩阵为 ( 3 3 3 ) ( 3 3 4 ) ( 3 3 5 ) ( 3 - 3 6 ) 哈尔滨t 程大学硕十学位论文 i 一厶 k d 盯= 2 a fb 7 d b d l ：蜴 00 其中：a 为三角形单元面积； 将挠度写成位移插值函数的形式为 w = n jn ：3 】u = n u 上+ l ；l 2 一，鹭一。葺 11 b 2 ( 厶茸+ i ll 。l 2 l 3 ) 一6 3 ( 茸厶+ i ll 。：厶) c ：( 厶葺+ 弘1l 2 l 3 ) 一c 3 ( l ；l 2 + 2 l 。l 2 l 3 ) b i = y 2 一y 3 ，c 1 = x 3 一x 2 ( 3 3 7 ) ( 3 3 8 ) ( 3 - 3 9 ) ( 3 - 4 0 ) ，可通过轮换下标i 一主一j 得到，d k t 单元的质量矩阵为 1 1 - 台 m 艄= 2 ai ip 。n d l 2 d 厶 ( 3 - 4 1 ) ( 3 ) 一般来说每个节点自由度的排列顺序为：“，v ，w ，臼。，曰，伊，则膜单元 每个节点上的刚度矩阵子块可写为 圹眨圳k m t l 2 k 2 i 。1 p 4 2 ， 则三角形壳体单元的每个节点上的刚度矩阵子块为 兀也 0 】2 。，【k m