流固耦合模拟分析之探讨.docx

流固耦合模拟分析之探讨.docx

《流固耦合模拟分析之探讨.docx》由会员分享，可在线阅读，更多相关《流固耦合模拟分析之探讨.docx（40页珍藏版）》请在冰点文库上搜索。

流固耦合模拟分析之探讨

MSC.Dytran流固耦合模拟分析之探讨

摘要：

流场会驱动固体；而固体的运动也可能驱动流体，甚至引起流场振动。

这称为流固互制或流固耦合。

MSC.Dytran软件的流固耦合功能包含一般耦合、任意偶合，且采用拉格朗日法与欧拉法分别描述固体与流体的运动。

拉格朗日的元素节点依附在材料上，节点随着材料质点作运动，故各物理量也作用在节点上随材料流动而变化。

相反，除任意耦合外，欧拉元素网格与节点不随时间而变，其物理量虽也作用在欧拉元素节点上，但对于通过欧拉元素面的各时间的质量、动量与能量的进与出，加之模拟，即模拟元素面的材料流，而不模拟各材料质点的时间历程。

因为对一般固体材料，要模拟各材料质点的时间历程，因此大多用拉格朗日法。

而对于流体不需要模拟材料质点的时间历程，故采用欧拉法，MSC.Dytran的欧拉法需用三维的计算域、三维的体元素与DMAT通用材料。

此外，欧拉法容许一个元素内含有两种以上的材料，这就是模拟计算材料流的扩散与混合行为。

欲推广应用型的计算软件，需有充分的应用范例。

关于流固耦合的模拟，除MSC.Dytran中文范例手册（1999b）与ExampleProblemManual（2001a）论述到一般耦合与任意耦合的应用范例之外，本文进一步探讨其他的应用。

而搭配的MSC.Patran软件除有Results的后处理工具外；该软件的Insight工具，能透明地看到体元素所构成的欧拉域内部，因此，更需用Insight，去展示欧拉域内部的流固耦合计算结果。

一、前言

MSC.Dytran之流固耦合计算功能大致上包括一般耦合（generalcoupling）和任意耦合（ALE,ArbitraryLagrangian-Eulerian）两大类。

欲推动一项泛用型计算软件被广泛地应用，须有可供参考的文件及丰富的应用案例。

而MSC.Dytran之中范例手册包括一个一般耦合和两个任意耦合的应用范例。

MSC.Dytran之ExampleProblemManual（2001）也包括两个一般耦合、两个任意耦合的范例。

本文之要旨是经由案例，探讨MSC.Dytran软件对流固耦合（欧拉域）与固体（拉格朗日域）的一般耦合与任意耦合的模拟计算，以及MSC.Dytran软件对流固耦合计算之前后处理与模型结果展示。

本文所叙述的流固耦合计算应用场合如下：

♦水下之固体物的告诉移动

♦造波板与之波浪水槽

♦海面上之高速物体撞击混凝土墙

♦上游侧与下游侧水深相等情况之潜堰

♦上游侧与下游侧水深不等情况之固定堰

♦上游侧与下游侧水深不等情况之阶梯式渠道

♦固定开度之水利闸门

♦隧洞内之气爆压作用在钢板上

♦隧洞内之水流推动固体物

二、对运动的数学描述法

拉格朗日法与欧拉法是对运动现象的两类不同的数学描述，可说是分别对材料质点流与空间流之描述。

拉格朗日法与欧拉法之元素网格可在同一计算模型内，但拉格朗日法的元素与欧拉法的元素分别拥有节点，只采用介面（interface），称为耦合面，才能将两者连结在一起；否则，纵使两者在空间内相互重叠，也彼此不相干，即忽视对方之存在。

2.1拉格朗日法

对固定的坐标系而言，拉格朗日元素的节点可相对地运动。

因节点系附在材料上，故材料连续体之节点系一起随着材料质点流而运动。

各拉格朗日元素的质量是不变量（invariant），但其元素体积可随时间而改变。

此外，速度、压力强度或质量密度等物理量系作用在拉格朗日元素的节点上，因此，各物理量系随着材料流（materialflow）而改变。

因对固体材料之行为，较须追踪各材料质点之时间历程，故适宜采用拉格朗日法。

拉格朗日法也适宜用以分析材料破坏（failure）或应变硬化（strainhardening）问题。

2.2欧拉法

除任意耦合（ALE）外，所有欧拉元素的网格与节点均保持固定，不随着时间或其他运动或变形。

换言之，欧拉元素不随时间而变。

各时间之速度、压力强度或密度等物理量也是作用在欧拉节点上，系计算经欧拉元素面的各时间的质量、动量与能量等的进与出之量，而不追踪各材料质点之时间历程。

因对流体较不须追踪各材料质点之时间历程，故一般系采用欧拉法。

MSC.Dytran之欧拉法必定使用三维的计算域及三维的体元素（solidelements），且限于通用材料DMAT。

欧拉法之特点是：

须采用较大之计算域，计算结果才不会受到计算域的边界之影响.

由于欧拉法系仿真经过欧拉元素面的材料流；且在一个元素之内，容许两种以上之材料，因此，应用欧拉法可模拟计算空气或水等材料流之扩散与混合现象。

三、流固耦合计算法之种类

当固体影响流体之后，被改变后之流体反过来影响固体；另一方向，亦然，就是流固耦合。

数值模拟计算是探讨分析流固耦合问题的主要方法之一。

由于流场动压变化所产生之流场特性、固体之几何形状、振幅与振频等互有关系，故流固耦合本质相当复杂。

纵使是单方向之流固耦合分析，对影响固体振动之水动力的附加质量也大多是估计。

同时，也不易准确地预测及量化固体系统内部之阻尼与流体吸收动能之效应。

因此，相关之数值模拟计算的难度相当高。

关于流固耦合之计算分析法可分为两大类如下：

3.1单方向之流固耦合分析

此类分析为简化之流固耦合计算。

即考虑固体单方向影响流体，但不考虑流体反过来影响固体；或反之。

3.2双方向之流固互制分析

此类分析包括流固耦合的顺序分析法与完全的流固耦合分析法两种如下：

3.2.1流固耦合之顺序分析法

此法系先分别计算流体与固体领域，每完成其中的一个领域之计算后，将计算结果作为另一个领域之荷载（loading）或边界条件，进行另一个领域的计算；反过来，也是相同的作法。

当计算软件系由固力计算程序与流力计算软程序结合而成时，就采用此种分析法，MSC.Dytran就是一例。

有些文献作较严格的定义，称此法并非耦合（uncoupled）分析，仅能称为流固互制分析。

3.2.2完全的流固耦合分析法

将所有之流固耦合相关的参数、边界与荷载等均融入流场与固体所共享之控制方程式组内，再采用数值计算法，求解耦合的（coupled）联立方程式，故作较严格定义的文献认为此法才是真正的流固耦合分析。

此种方式虽最完整，但难度也最高，故使用者最少。

四、流固耦合计算功能

MSC.Dytran之流固耦合计算就是拉格朗日域（固体）与欧拉域（流场）的耦合计算，输入档内容直接相关者系COUPLE卡或COUPLE1卡。

在计算过程内，拉格朗日域与欧拉域分别进行计算，每完成一个领域之计算后，再计算另一个领域；即前一领域之计算结果作为另一个领域计算所需之荷载或边界条件。

至于MSC.Dytran的流固耦合计算，则分为两种如下：

4.1一般耦合

通常是采用拉格朗日法模拟固体，以及采用欧拉法模拟流体。

至于一般耦合（generalcoupling），大多是拉格朗日的固体在欧拉的流场范围内运动，即拉格朗日域驱动欧拉域；流场虽有速度，但代表流场的欧拉格网系固定及不受拉格朗日的固体之影响。

换言之，在流固耦合过程内，欧拉格网不移动，也不变形。

一般耦合计算之前处理大多用封闭之假壳（dummyshell）耦合面隔开拉格朗日域与欧拉域；于起始计算之时，拉格朗日的固体至少有一微小量（譬如，0.001m）重叠在欧拉域的范围内，且固体、流体或两者有运动，才能启动流固耦合之计算器制。

若拉格朗日域与欧拉域毫无重叠，则无法起动流固耦合计算。

当然，拉格朗日的固体可完全位于欧拉域内，不因拉格朗日的固体运动而使欧拉格网移动或变形。

此种流固耦合最适用在「固体驱动流体」现象之模拟，譬如，水上飞机降落在水面上之行为、水面下之物体的运动所引致之流体动力行为、隧道内之高速车辆引致的气动力行为、管道内的固体活塞推动气体等。

此外，一般耦合也可应用在具有不规则的固定固体边界之流场模拟计算，即将不动的固体边界与流场之关联视为流固耦合现象。

一般耦合须在拉格朗日的几何体之外露面（几何面段或元素表面，经由SURFACE卡建立），定义封闭的假壳元素群，作为流固耦合面，形成闭合之体积，完全隔开拉格朗日域与欧拉域；且对假壳元素，须进行Equivalence与Verify法线方向之前处理。

耦合面虽为虚拟，但也须输入物理性质

（Properties）为2D/SummyShell。

然后，对耦合面的全部假壳元素，使用MSC.Patran软件的Coupling等，就可完成相关的前处理。

在计算过程内，耦合面是欧拉流场格网的边界，即假壳元素之耦合面与仿真流体之体元素相接触（Contact），MSC.Dytran先计算欧拉流体施加在耦合面上之荷载，然后，耦合面使拉格朗日格网产生应力与变形。

至于输入档相关的Couple卡上之参数COVER，则定义被耦合面所包裹者（即不含欧拉领域者）是属内部（INSIDE），还是外部（OUTSIDE）。

譬如，石头掉入水中之模拟，其COVER是INSIDE。

汽车的安全气囊的模拟之COVER是OUTSIDE。

4.2任意耦合

任意耦合主要是用来仿真欧拉域驱动拉格朗日域，在该耦合过程内，于计算模型承受来自欧拉域的荷载之后，拉格朗日的元素与欧拉的体元素之网格均可能变形或移位。

任意耦合最适宜用来仿真「流体驱动固体」，譬如，容器内之气体爆炸、空中之鸟体撞击飞机壳体、气爆压作用在壳形固体物上、管路内的流体推动固体、因压力波之钢管的膨胀或收缩等之动力行为。

任意耦合之流固界面或耦合面不须是封闭的面，可为两片以上的不相连接之面段，也可为一个拉格朗日域对应于两个不相连接之欧拉域。

因目前的MSC.Patran软件尚不能用来直接进行ALE相关之前处理，故于实际建置任意耦合的计算模型时，须先经由Master-slavesurface的接触过程，建立流固耦合所作用的面段，再于执行MSC.Patran之后，使用文本编辑器如记事本，修改*.DAT文件内之接触相关指令，将Contact卡改写成ALE卡及建立ALEGRID卡。

此外，也须包含所有相关的欧拉节点，建立ALEGRID组（Group/Set），且使用文本编辑器（如记事本）去定义ALEGRID卡。

在建置主（Master）从（Slave）接触模型的过程内，一般是被动者为主，主动者为从。

ALE之耦合面系经由输入档内的ALE卡产生。

一般而言，ALE耦合面之欧拉与拉格朗日两类元素节点须是一对一，位置重合。

由于ALE耦合面之互制移动与变形，代表固体之拉格朗日格网随着时间而移动或变形，且欧拉格网也在空间移动或变形。

即固体变形时，耦合面之位置与形状随之改变，并带动邻接耦合面的欧拉格网与其他部分欧拉格网作相应之移动与变形。

因此，一方面，流体材料在欧拉格网内移动；另一方面，欧拉格网本身也在移动，以致欧拉网格的位置与形状系随时间作不断地改变。

惟因欧拉网格一定是体元素，故ALE之仿真只具有3个自由度（DOF）。

至于欧拉元素之节点可定义为ALE节点之必要条件，则如下：

•须为单种欧拉材料；

•适用于固体变形光滑者；

•欧拉材料不能具有剪切刚度；

•对流场，不使用Roe求解器（根据PhilipRoe教授的理论之流力计算法）。

五、一般耦合计算之应用例

5.1固体驱动流体之仿真计算

5.1.1隧道内之高速车辆对流场的效应

高速车辆从较空旷的环境进入空间受限的隧道内之后，行进的车体会较显着地压缩前方与其周围的空气，以致形成明显的活塞之气动力（aerodynamicforce）效应，如空气压与行车阻力均增大，以及在隧道出口产生微气压波（micropressurewaves）等，而微气压波可使乘客感觉耳朵不适。

此种车体进入隧道时之周围空气压的波动称为Kelvin-Helmholtz不稳现象。

根据日本新干线之实测结果，微气压波的最大波压约与隧道内的车速的三次方成正比，且与离开隧道出口之距离成反比。

针对此现象，通常的高速铁路设计基准之一系满足隧道内之最大空气压的变动率不超过3.0kPa/sec。

此外，隧道内之有轨列车若有显着的振动，大多系因车辆系统之阻尼不足与隧道内之气动力等引致列车之钢轮、钢轨与空气流场的耦合振动。

于高速列车驶入隧道内时，因隧道内的空气压与其他气动力特性会有明显的改变，故列车之侧向振动可能加剧，使乘客感到不适，甚至造成乘客之恐惧感及影响行车之安全。

若隧道内系铺设双向轨道，双向行车，列车大致占隧道之一半空间，则车体左右两侧所受之空气压会显着不同，空气压差更使得车体有侧向振动加剧之倾向。

基于上述之考虑，本例系针对高速列车进入隧道后之情况，应用MSC.Dytran，进行三维的流固耦合计算，探讨空气动力行为，亦即根据模拟计算结果，探求隧道内的空气压与空气流速等之时间与空间的变化情形，俾利于将来进一步研析隧道内可能引致的高速车体振动。

兹考虑由x-、y-、z-坐标轴所构成之三维空间，z-与x-轴形成水平面，y-轴之反方向是重力方向，车行方向为z-轴。

几何形状包括隧道内的空气流场与车体两部分。

计算例的隧道长40m，其横断面系由直径9.2m之半圆与高度2m、宽度9.2m之矩形所构成。

隧道内的空气流场与进入隧道之车体均是三维的几何实体（solids），且均由8个节点且6个面的CHEXA体元素（solidelements）构成。

空气流场共采用1,120个欧拉的体元素，材料为理想气体之通用材料IdealGas（DMAT），元素之物理性质（Properties）系选用Hydro（Peuler1）。

车体之模拟系使用8个三维之拉格朗日的体元素，其物理性质选用弹塑性之通用材料ElasPlas（DMAT）。

而根据高速铁路列车之数据，车之质量密度采用329kg/m3。

空气流场与车体之有限元素模型如图1。

所考虑之车体宽度3m、高度3.5m及长度10m，行经隧道之左半部内，如图2所示。

图1交通隧道之几何与有限元素模型（图标之车体在右上角，车体外部包裹着流固耦合面，车行方向系由图示之右上至左下；对隧道之上游与下游断面均给予欧拉之flow边界条件）

流固耦合面系密闭地包裹车体，并由24片二维之假壳元素（dummyshell）构成，紧附在拉格朗日的体元素上。

至于该元素材料的物理性质，也选用“DummyShell”。

因流固耦合面是假壳元素，故不给予材料数据。

空气之材料数据系采用伽马尔（Gammalaw）之状态方程式（EOS,equationofstate）。

因空气压是质量密度、比热比与单位质量之内能的函数，故空气之状态方程式如下：

p=（γ－1）ρe

在上式中，p=空气压；

图2隧道与车体之横断面（自隧道进口往出口看之立视图）

图3第0.65sec时之隧道内面的流速分布云纹图

γ=比热比=气体之等压比热对等容积比热之比值（ratioofspecificheats）=

Cp=等压之比热；

Cv=等容积之比热；

ρ=空气之质量密度=1.2887gm/cm3；

e=空气单位质量之内能。

弹塑性之车体系采用多项式状态方程式（EOSPOL），主要性质如下：

•质量密度=329kg/m3

•剪切模数=8.18×1010N/m2

•降伏应力=3.50×108N/m2

关于隧道内之空气流场，元素性质采用Hydro（Peuler1），各欧拉的体元素内之初流速为零，各体元素内之空气质量密度为1.2887gm/cm3，并经由MSC.Patran之“Load&BC/Init.Cond.Euler.”，分别依序建立”Shape”、”InitialValues”与”RegionDefinition”三个option，构成空气流场的起始条件。

在隧道入口断面与出口断面的空气流场均采用flow边界条件。

该边界条件是：

在欧拉的体元素与外界之接触面上，容许空气流入或流出。

至于不给flow条件之欧拉域的边界，则预设（default）为欧拉材料不穿透之刚性壁（rigidwall），也就是固定边界。

车体强制速度的边界条件是：

于所有时候，z-轴向之速度60m/sec（即216km/hr）作用在所有之拉格朗日的体元素节点上，该速度代表车速，并经由MSC.Patran之”Loads/BCs/Velocity”菜单建立。

流固耦合面是由24片假壳元素所构成，密闭地包裹车体，也是欧拉的空气流场之边界。

空气压作用在耦合面上，再传至拉格朗日的体元素节点上。

耦合边界条件系经由MSC.Patran之Load&BC/Coupling菜单建立。

本例系采用一般耦合（generalcoupling）计算，即欧拉的体元素之网格不随时间而改变。

至于COVER为“INSIDE”，则表示假壳元素所包裹之空间内，无欧拉元素存在。

于开始执行MSC.Patran之前，经由MSC.Patran之“Analysis/Analyze/InputDeck/Translate”的“ExecutionControl”菜单下，设定仿真计算时间步长的总个数、起始计算之时间步长（INISTEP）与容许最小之计算时间步长（MINSTEP）等参数。

本例系设定时间步长的总个数为4万个、起始之时间步长为1.0×10-5秒、最小之时间步长为1.0×10-6秒。

MSC.Dytran之欧拉域的计算结果输出系以每个元素为单位。

为存取隧道内的空气流场之计算结果，于计算之前，在MSC.Patran之“Analysis/Analyze/InputDeck/Translate”菜单的“Output

Request”下，建立Archive及“ElementOutput”型式之输出文件。

本例系每1,000个时间步长，输出计算结果一次，输出项目包括Pressure、x-vel、y-vel、z-vel（即空气压与三个坐标轴向之流速分量）等。

计算完成后，经由MSC.Patran之“Analysis/ReadArchiveFiles/Object=ModelandResults/Translate”菜单，读取相关之Archive档（即*.ARC文件）的数据。

图3显示第Cycle40,000（第0.65sec）之隧道内面的空气流速分布的云纹图，该图系由隧道出口左侧往进口之方向看。

该云纹图所显示者包括：

隧道左侧壁、隧道顶与出口断面之各元素的空气平均流速。

大致上，车头附近之空气流速较高，车尾附近之空气流速较低。

在隧道出口断面上，车体左侧的空气元素（#574）与右侧的空气元素（#540）之平均空气压的历时曲线比较如图4，其显示左右两侧之空气压显着不同。

因车体系靠近隧道左侧壁，故于各时间，车体左侧的平均空气压明显地较车体右侧为大。

图示之空气压最大值系发生于第0.35sec，即车头行至隧道中点之时。

因本例所考虑之车行时间与隧道长度均尚短，故计算结果尚未得到隧道内的负空气压与压力波。

图4隧道出口断面的车体左侧与右侧元素之平均空气压的历时曲线

未来可能应用MSC.Dytran软件，进一步研究之隧道与高速列车的课题如下：

♦隧道内的微气压波之模拟计算

♦车体因不平衡的气动力之振动模拟计算

5.1.2水下之固体物的高速移动

本计算例所用之欧拉域系长20m（x向），宽21m（z向）及高10m（y向）的矩形立方体，在x-y-z向的有限元素之网格数分别为20、20、21。

负y轴向是重力的方向。

欧拉域的下部5m深为水域，上部5m深为模拟空气之空（void）域。

因水域与空域之材料均Hydro（peuler1），故须分别在MSC.Patran的Loads/BCs菜单下，分Shape、InitialValue与RegionDefinition的步骤，给予欧拉的起始条件。

运动速度为10m/sec（velocity边界条件）的水下高速物体系由拉格朗日元素所构成，它长4m，宽1m及高1m，如图5所示，对水域底面与上游边界之净距均为1m。

欧拉域之上游面（x=0m）与下游面采用flow边界条件，针对水之材料，x、y与z向的流速分量均给予零。

欧拉域内的所有节点的起始流速也给予零。

图5采用FMAT等于0.002488的等值面模拟水面

本例系应用MSC.Patran软件，采用流固的一般耦合计算，在水内，用假壳（dummyshell）元素包裹高速固体物的外部，形成封闭的流固耦合面，经MSC.Patran的有限元素之Equivalence及Verify假壳的法线方向，使所有之法线指向外，以致假壳所构成之耦合面无破洞；然后，藉由该耦合面及Loads/BCs/Create/Coupling菜单之设定，执行流固的一般耦合计算。

图5显示波动的自由水面与水下之流速分量；但未显示水下固体物之移动。

换言之，Insight工具虽显示欧拉域内部之动态状况，但不显示欧拉域范围内之拉格朗日固体物的运动状况。

此为目前MSC.Patran的限制。

5.1.3造波板与之波浪水槽

本计算模型是20m（x向）×10m（y向）×3m（z向）的波浪水槽，内含7m高（y向）、3m宽（z向）与0.5m厚（x向）之造波板。

波浪水槽之内部是欧拉域，分成20×20×3个体元素；造波板是拉格朗日域，分成1×14×3个体元素，如图6。

如图所示的圆点代表造波板之外部包裹假壳的耦合面，使拉格朗日域与欧拉域产生一般耦合计算。

造波板之所有元素节点均有正弦时变性之x向的速（velocityBC），最大速度为1m/sec，周期为1sec。

图6模拟造波板与波浪水槽之有限元素格网

欧拉域之下部5m是水，材料是LinFluid；上部是空气，当作空域处理。

整个欧拉域之上游面（x=0m）与下游面（x=20m）均给与予flow边界条件。

水与空域之性质均是Hydro（Peler1），搭配MSC.Patran的Loads/BCs菜单下之欧拉起始条件（Init.Cond.Euler），但仅给予水域x向之初速度0.1m/sec。

在代表水槽之欧拉域的下游面之flow边界条件，系给予流速0.1m/sec（x向）；若欧拉域之上游面之flow边界条件，也是x向流速0.1m/sec，则第1.0sec时的上游边界附近会有水面特别壅高之现象，如图7所示。

若上游面之边界流速改为－0.1m/sec，则该壅高水面就可消除。

图7水槽内的造波板（在中央者）引致约第0.8sec时之水面壅高及流速的空间分布云纹图

5.1.4海面上之高速物体撞击混凝土墙

本例之目的在显示一个计算模型可包括两组以上的一般耦合，且被包裹耦合面的固体物可有一部分在水域内与模拟空气之空域内。

本例也探讨一般耦合与固体对固体之撞击同在一个计算模型内的状况。

计算模型之几何体是三个较小的长方体（solid）在一个较大长方体内，如图8。

较小的长方体由左至右，分别称为maker、body与struc，均是用拉格朗日法描述之固体。

尺寸0.5m×7m×3m的maker是造波板，用来模拟海洋环境；1m×1m×1m的body，用来模拟高速运动之水上物体；尺寸1m×9m×3m的struc用来模拟混凝土