流固耦合模拟分析之探讨

1、MSC.Dytran 流固耦合模拟分析之探讨摘要：流场会驱动固体；而固体的运动也可能驱动流体，甚至引起流场振动。这称为流固互制或流固耦合。MSC.Dytran软件的流固耦合功能包含一般耦合、任意偶合，且采用拉格朗日法与欧拉法分别描述固体与流体的运动。Dytran中流固耦合计算种类拉格朗日的元素节点依附在材料上，节点随着材料质点作运动，故各物理量也作用在节点上随材料流动而变化。相反，除任意耦合外，欧拉元素网格与节点不随时间而变，其物理量虽也作用在欧拉元素节点上，但对于通过欧拉元素面的各时间的质量、动量与能量的进与出，加之模拟，即模拟元素面的材料流，而不模拟各材料质点的时间历程。因为对一般固体材料

2、，要模拟各材料质点的时间历程，因此大多用拉格朗日法。而对于流体不需要模拟材料质点的时间历程，故采用欧拉法，MSC.Dytran的欧拉法需用三维的计算域、三维的体元素与DMAT通用材料。此外，欧拉法容许一个元素内含有两种以上的材料，这就是模拟计算材料流的扩散与混合行为欧拉法模拟计算材料流的扩散与混合行为。欲推广应用型的计算软件，需有充分的应用范例。关于流固耦合的模拟，除MSC.Dytran中文范例手册（1999b）与Example Problem Manual(2001a)论述到一般耦合与任意耦合的应用范例之外，本文进一步探讨其他的应用。而搭配的MSC.Patran软件除有Results的后处理

3、工具外；该软件的Insight工具，能透明地看到体元素所构成的欧拉域内部，因此，更需用Insight，去展示欧拉域内部的流固耦合计算结果。一、 前言MSC.Dytran之流固耦合计算功能大致上包括一般耦合(general coupling)和任意耦合(ALE, Arbitrary Lagrangian-Eulerian)两大类。欲推动一项泛用型计算软件被广泛地应用，须有可供参考的文件及丰富的应用案例。而MSC.Dytran之中范例手册包括一个一般耦合和两个任意耦合的应用范例。MSC.Dytran 之Example Problem Manual (2001)也包括两个一般耦合、两个任意耦合的范例

4、。本文之要旨是经由案例，探讨MSC.Dytran软件对流固耦合（欧拉域）与固体（拉格朗日域）的一般耦合与任意耦合的模拟计算，以及MSC.Dytran软件对流固耦合计算之前后处理与模型结果展示。本文所叙述的流固耦合计算应用场合如下： 水下之固体物的告诉移动 造波板与之波浪水槽 海面上之高速物体撞击混凝土墙 上游侧与下游侧水深相等情况之潜堰 上游侧与下游侧水深不等情况之固定堰 上游侧与下游侧水深不等情况之阶梯式渠道 固定开度之水利闸门 隧洞内之气爆压作用在钢板上 隧洞内之水流推动固体物二、 对运动的数学描述法拉格朗日法与欧拉法是对运动现象的两类不同的数学描述，可说是分别对材料质点流与空间流之描述。

5、拉格朗日法与欧拉法之元素网格可在同一计算模型内，但拉格朗日法的元素与欧拉法的元素分别拥有节点，只采用介面(interface)，称为耦合面，才能将两者连结在一起；否则，纵使两者在空间内相互重叠，也彼此不相干，即忽视对方之存在。2.1 拉格朗日法对固定的坐标系而言，拉格朗日元素的节点可相对地运动。因节点系附在材料上，故材料连续体之节点系一起随着材料质点流而运动。各拉格朗日元素的质量是不变量(invariant)，但其元素体积可随时间而改变。此外，速度、压力强度或质量密度等物理量系作用在拉格朗日元素的节点上，因此，各物理量系随着材料流(material flow)而改变。因对固体材料之行为， 较须

6、追踪各材料质点之时间历程，故适宜采用拉格朗日法。拉格朗日法也适宜用以分析材料破坏(failure)或应变硬化(strain hardening)问题。2.2欧拉法除任意耦合(ALE)外，所有欧拉元素的网格与节点均保持固定，不随着时间或其他运动或变形。换言之，欧拉元素不随时间而变。各时间之速度、 压力强度或密度等物理量也是作用在欧拉节点上，系计算经欧拉元素面的各时间的质量、动量与能量等的进与出之量，而不追踪各材料质点之时间历程。 因对流体较不须追踪各材料质点之时间历程，故一般系采用欧拉法。MSC.Dytran之欧拉法必定使用三维的计算域及三维的体元素(solid elements)，且限于通用材

7、料 DMAT。欧拉法之特点是：须采用较大之计算域，计算结果才不会受到计算域的边界之影响.欧拉法的特点和要求由于欧拉法系仿真经过欧拉元素面的材料流；且在一个元素之内，容许两种以上之材料，因此，应用欧拉法可模拟计算空气或水等材料流之扩散与混合现象。欧拉法计算混合现象三、 流固耦合计算法之种类当固体影响流体之后，被改变后之流体反过来影响固体；另一方向，亦然，就是流固耦合。数值模拟计算是探讨分析流固耦合问题的主要方法之一。由于流场动压变化所产生之流场特性、固体之几何形状、振幅与振频等互有关系，故流固耦合本质相当复杂。纵使是单方向之流固耦合分析，对影响固体振动之水动力的附加质量也大多是估计。同时，也不易

8、准确地预测及量化固体系统内部之阻尼与流体吸收动能之效应。因此，相关之数值模拟计算的难度相当高。关于流固耦合之计算分析法可分为两大类如下：3.1 单方向之流固耦合分析此类分析为简化之流固耦合计算。即考虑固体单方向影响流体，但不考虑流体反过来影响固体；或反之。3.2 双方向之流固互制分析此类分析包括流固耦合的顺序分析法与完全的流固耦合分析法两种如下：3.2.1 流固耦合之顺序分析法此法系先分别计算流体与固体领域，每完成其中的一个领域之计算后，将计算结果作为另一个领域之荷载(loading)或边界条件，进行另一个领域的计算；反过来，也是相同的作法。当计算软件系由固力计算程序与流力计算软程序结合而成时

9、，就采用此种分析法，MSC.Dytran 就是一例。 有些文献作较严格的定义，称此法并非耦合(uncoupled)分析，仅能称为流固互制分析。此法并非耦合分析，仅能称为流固互制分析3.2.2 完全的流固耦合分析法将所有之流固耦合相关的参数、边界与荷载等均融入流场与固体所共享之控制方程式组内，再采用数值计算法，求解耦合的(coupled)联立方程式，故作较严格定义的文献认为此法才是真正的流固耦合分析。此种方式虽最完整，但难度也最高，故使用者最少。真正的流固耦合分析四、流固耦合计算功能MSC.Dytran 之流固耦合计算就是拉格朗日域(固体)与欧拉域(流场)的耦合计算，输入档内容直接相关者系 CO

10、UPLE卡或 COUPLE1卡。在计算过程内，拉格朗日域与欧拉域分别进行计算，每完成一个领域之计算后，再计算另一个领域；即前一领域之计算结果作为另一个领域计算所需之荷载或边界条件。至于 MSC.Dytran 的流固耦合计算，则分为两种如下：4.1 一般耦合通常是采用拉格朗日法模拟固体，以及采用欧拉法模拟流体。至于一般耦合(general coupling)，大多是拉格朗日的固体在欧拉的流场范围内运动，即拉格朗日域驱动欧拉域一般耦合应用在Lagrange的固体在Euler的流场范围内运动。；流场虽有速度，但代表流场的欧拉格网系固定及不受拉格朗日的固体之影响。换言之，在流固耦合过程内，欧拉格网不移

11、动，也不变形。一般耦合计算之前处理大多用封闭之假壳(dummy shell) 耦合面隔开拉格朗日域与欧拉域；于起始计算之时，拉格朗日的固体至少有一微小量(譬如，0.001 m)重叠在欧拉域的范围内，且固体、流体或两者有运动，才能启动流固耦合之计算器制。若拉格朗日域与欧拉域毫无重叠，则无法起动流固耦合计算。当然，拉格朗日的固体可完全位于欧拉域内，不因拉格朗日的固体运动而使欧拉格网移动或变形。此种流固耦合最适用在固体驱动流体现象之模拟，譬如，水上飞机降落在水面上之行为、水面下之物体的运动所引致之流体动力行为、隧道内之高速车辆引致的气动力行为、 管道内的固体活塞推动气体等。此外，一般耦合也可应用在具

12、有不规则的固定固体边界之流场模拟计算，即将不动的固体边界与流场之关联视为流固耦合现象。一般耦合须在拉格朗日的几何体之外露面(几何面段或元素表面，经由SURFACE 卡建立)，定义封闭的假壳元素群，作为流固耦合面，形成闭合之体积，完全隔开拉格朗日域与欧拉域；且对假壳元素，须进行 Equivalence 与 Verify 法线方向之前处理。耦合面虽为虚拟，但也须输入物理性质(Properties)为 2D/Summy Shell。然后，对耦合面的全部假壳元素， 使用 MSC.Patran 软件的 Coupling 等，就可完成相关的前处理。在计算过程内，耦合面是欧拉流场格网的边界，即假壳元素之耦合

13、面与仿真流体之体元素相接触(Contact)， MSC.Dytran 先计算欧拉流体施加在耦合面上之荷载，然后，耦合面使拉格朗日格网产生应力与变形。至于输入档相关的 Couple 卡上之参数 COVER，则定义被耦合面所包裹者（即不含欧拉领域者）是属内部(INSIDE)，还是外部(OUTSIDE)。譬如，石头掉入水中之模拟，其 COVER 是 INSIDE。汽车的安全气囊的模拟之COVER 是 OUTSIDE。4.2 任意耦合任意耦合主要是用来仿真欧拉域驱动拉格朗日域，在该耦合过程内，于计算模型承受来自欧拉域的荷载之后，拉格朗日的元素与欧拉的体元素之网格均可能变形或移位。任意耦合最适宜用来仿真

14、流体驱动固体，譬如，容器内之气体爆炸、空中之鸟体撞击飞机壳体、气爆压作用在壳形固体物上、 管路内的流体推动固体、因压力波之钢管的膨胀或收缩等之动力行为。任意耦合应用范围任意耦合之流固界面或耦合面不须是封闭的面，可为两片以上的不相连接之面段，也可为一个拉格朗日域对应于两个不相连接之欧拉域。因目前的 MSC.Patran 软件尚不能用来直接进行 ALE 相关之前处理， 故于实际建置任意耦合的计算模型时，须先经由 Master-slave surface 的接触过程，建立流固耦合所作用的面段，再于执行 MSC.Patran 之后，使用文本编辑器如记事本，修改 *.DAT 文件内之接触相关指令，将 C

15、ontact 卡 改写成 ALE 卡及建立 ALEGRID 卡。此外，也须包含所有相关的欧拉节点，建立 ALEGRID 组 (Group/Set)，且使用文本编辑器(如记事本)去定义ALEGRID 卡。在建置主(Master)从(Slave)接触模型的过程内，一般是被 动者为主，主动者为从。主从关系的建立ALE 之耦合面系经由输入档内的 ALE 卡产生。一般而言，ALE 耦合面之欧拉与拉格朗日两类元素节点须是一对一，位置重合。由于 ALE 耦合面之互制移动与变形，代表固体之拉格朗日格网随着时间而移动或变形，且欧拉格网也在空间移动或变形。即固体变形时，耦合面之位置与形状随之改变，并 带动邻接耦合

16、面的欧拉格网与其他部分欧拉格网作相应之移动与变形。因 此，一方面，流体材料在欧拉格网内移动；另一方面，欧拉格网本身也在移动，以致欧拉网格的位置与形状系随时间作不断地改变。惟因欧拉网格一定是体元素，故 ALE 之仿真只具有 3 个自由度(DOF)。至于欧拉元素之节点可定义为 ALE 节点之必要条件，则如下： 须为单种欧拉材料； 适用于固体变形光滑者； 欧拉材料不能具有剪切刚度； 对流场，不使用 Roe 求解器 (根据 Philip Roe 教授的理论之流力计算法)。五、 一般耦合计算之应用例5.1 固体驱动流体之仿真计算5.1.1 隧道内之高速车辆对流场的效应高速车辆从较空旷的环境进入空间受限的

17、隧道内之后，行进的车体会较显着地压缩前方与其周围的空气，以致形成明显的活塞之气动力（aerodynamic force）效应，如空气压与行车阻力均增大，以及在隧道出口产生微气压波(micropressure waves)等，而微气压波可使乘客感觉耳朵不适 。此种车体进入隧道时之周围空气压的波动称为Kelvin-Helmholtz 不稳现象。根据日本新干线之实测结果，微气压波的最大波压约与隧道内的车速的三次方成正比，且与离开隧道出口之距离成反比。针对此现象，通常的高速铁路设计基准之一系满足隧道内之最大空气压 的变动率不超过 3.0 kPa/sec。此外，隧道内之有轨列车若有显着的振动， 大多系因

18、车辆系统之阻尼不足与隧道内之气动力等引致列车之钢轮、钢轨与空气流场的耦合振动。于高速列车驶入隧道内时，因隧道内的空气压与其他气动力特性会有明显的改变，故列车之侧向振动可能加剧，使乘客感到不适，甚至造成乘客之恐惧感及影响行车之安全。若隧道内系铺设双向轨道，双向行车，列车大致占隧道之一半空间，则车体左右两侧所受之空气压会显 着不同，空气压差更使得车体有侧向振动加剧之倾向。基于上述之考虑，本例系针对高速列车进入隧道后之情况，应 用 MSC.Dytran，进行三维的流固耦合计算，探讨空气动力行为，亦即根据 模拟计算结果，探求隧道内的空气压与空气流速等之时间与空间的变化情 形，俾利于将来进一步研析隧道内

19、可能引致的高速车体振动。兹考虑由 x-、y-、z-坐标轴所构成之三维空间，z- 与 x- 轴形成水平面，y-轴之反方向是重力方向，车行方向为 z- 轴。几何形状包括隧道内的空气流场与车体两部分。计算例的隧道长 40 m，其横断面系由直径9.2 m 之半圆与高度 2 m、宽度 9.2 m 之矩形所构成。隧道内的空气流场与 进入隧道之车体均是三维的几何实体(solids)，且均由 8 个节点且 6 个面 的 CHEXA 体元素(solid elements)构成。空气流场共采用 1,120 个欧拉的体元素，材料为理想气体之通用材料 Ideal Gas (DMAT)，元素之物理性质(Properti

20、es)系选用Hydro (Peuler1)。车体之模拟系使用 8 个三维之拉格朗日的体元素，其 物理性质选用弹塑性之通用材料 ElasPlas (DMAT)。而根据高速铁路列车 之数据，车之质量密度采用 329 kg/m3。空气流场与车体之有限元素模型如图1。所考虑之车体宽度 3 m、高度 3.5 m 及长度 10m，行经隧道之左半部 内，如图2 所示。图 1 交通隧道之几何与有限元素模型 (图标之车体在右上角，车体外部包 裹着流固耦合面，车行方向系由图示之右上至左下；对隧道之上游与下游断 面均给予欧拉之 flow 边界条件)流固耦合面系密闭地包裹车体，并由 24 片二维之假壳元素(dummy

21、 shell)构成，紧附在拉格朗日的体元素上。至于该元素材料的物理 性质，也选用 “Dummy Shell”。 因流固耦合面是假壳元素，故不给予材料数据。空气之材料数据系采用伽马尔(Gamma law)之状态方程式(EOS, equation of state)。因空气压是质量密度、比热比与单位质量之内能的函数，故空气之状态方程式如下：p = (1) e在上式中，p = 空气压；图 2 隧道与车体之横断面 (自隧道进口往出口看之立视图)图 3 第 0.65 sec 时之隧道内面的流速分布云纹图 = 比热比 = 气体之等压比热对等容积比热之比值 (ratio of specific heats)

22、 =C p = 等压之比热；Cv = 等容积之比热； = 空气之质量密度 = 1.2887 gm/cm3 ；e = 空气单位质量之内能。弹塑性之车体系采用多项式状态方程式(EOSPOL)，主要性质如下： 质量密度 = 329 kg/m3 剪切模数 = 8.18×1010 N/m2 降伏应力 = 3.50×108 N/m2关于隧道内之空气流场，元素性质采用 Hydro (Peuler1)， 各欧拉的体元素内之初流速为零，各体元素内之空气质量密度为 1.2887gm/cm3，并经由 MSC.Patran 之 “Load & BC/Init. Cond. Euler.”，

23、 分别依序建立”Shape”、”Initial Values”与 ”Region Definition” 三个 option，构成空气流场的起始条件。在隧道入口断面与出口断面的空气流场均采用 flow 边界条件。该边界条件是：在欧拉的体元素与外界之接触面上，容许空气流入或流出。至于不给 flow 条件之欧拉域的边界，则预设(default)为欧拉材料不 穿透之刚性壁(rigid wall)，也就是固定边界。车体强制速度的边界条件是：于所有时候，z- 轴向之速度 60m/sec (即 216 km/hr)作用在所有之拉格朗日的体元素节点上，该速度代 表车速，并经由 MSC.Patran 之 ”L

24、oads/BCs/Velocity”菜单建立。流固耦合面是由 24 片假壳元素所构成，密闭地包裹车体，也是欧拉的空气流场之边界。空气压作用在耦合面上，再传至拉格朗日的体元 素节点上。耦合边界条件系经由 MSC.Patran 之 Load & BC/Coupling 菜单建立。本例系采用一般耦合(general coupling)计算，即欧拉的 体元素之网格不随时间而改变。至于 COVER 为 “INSIDE” ，则表示假壳元 素所包裹之空间内，无欧拉元素存在。于 开 始 执 行 MSC.Patran 之 前 ，经由 MSC.Patran 之 “Analysis/Analyze/Inpu

25、t Deck/Translate” 的 “ExecutionControl”菜单下，设定仿真计算时间步长的总个数、起始计算之时间步长(INISTEP)与容许最小之计算时间步长(MINSTEP)等参数。本例系设定时间 步长的总个数为 4 万个、起始之时间步长为 1.0×10-5 秒、最小之时间步长 为 1.0×10-6 秒。MSC.Dytran 之欧拉域的计算结果输出系以每个元素为单位。为存取隧道内的空气流场之计算结果，于计算之前，在 MSC.Patran 之 “Analysis/Analyze/Input Deck/Translate”菜单的 “OutputRequest”

26、 下，建立 Archive 及 “Element Output”型式之输出文件。本例系每 1,000个时间步长，输出计算结果一次，输出项目包括 Pressure、x-vel、y-vel、z-vel (即空气压与三个坐标轴向之流速分量)等。计算完 成后，经由 MSC.Patran 之“Analysis/Read Archive Files/Object= Model and Results/Translate”菜单，读取相关之 Archive 档(即*.ARC 文件)的数据。图3 显示第 Cycle 40,000 (第0.65 sec)之隧道内面的 空气流速分布的云纹图，该图系由隧道出口左侧往进

27、口之方向看。该云纹图所显示者包括：隧道左侧壁、隧道顶与出口断面之各元素的空气平均流速。 大致上，车头附近之空气流速较高，车尾附近之空气流速较低。在隧道出口断面上，车体左侧的空气元素(#574)与右侧的空气元素(#540)之平均空气压的历时曲线比较如图4，其显示左右两侧之空气压显着不同。因车体系靠近隧道左侧壁，故于各时间，车体左侧的平均空气压明显地较车体右侧为大。图示之空气压最大值系发生于第 0.35 sec，即车头行至隧道中点之时。因本例所考虑之车行时间与隧道长度均尚短，故计 算结果尚未得到隧道内的负空气压与压力波。图 4 隧道出口断面的车体左侧与右侧元素之平均空气压的历时曲线未来可能应用 M

28、SC.Dytran 软件，进一步研究之隧道与高速 列车的课题如下： 隧道内的微气压波之模拟计算 车体因不平衡的气动力之振动模拟计算5.1.2 水下之固体物的高速移动本计算例所用之欧拉域系长 20 m (x 向)，宽 21 m (z 向)及高10 m (y 向)的矩形立方体，在 x-y-z 向的有限元素之网格数分别为 20、20、21。负 y 轴向是重力的方向。欧拉域的下部 5 m 深为水域，上部 5 m 深为模拟空气之空 (void) 域 。因水域与空域之材料均Hydro(peuler1)，故须分别在 MSC.Patran 的 Loads/BCs 菜单下，分Shape、Initial Valu

29、e 与 Region Definition 的步骤，给予欧拉的起始条件。运动速度为 10 m/sec( velocity 边界条件)的水下高速物体系由拉格朗日元素所构成，它长 4 m，宽 1 m 及 高 1 m ，如图5 所示，对水域底面与上游边界之净距均为 1 m。欧拉域之上游面(x= 0 m)与下游面采用 flow 边界条件，针对水之材料，x、y 与 z 向的流速分量均给予零。欧拉域内的所有节点的起始流速也给予零。图 5 采用 FMAT 等于 0.002488 的等值面模拟水面本例系应用 MSC.Patran 软件，采用流固的一般耦合计算，在水内，用假壳(dummy shell)元素包裹高

30、速固体物的外部，形成封闭的流固耦合面，经 MSC.Patran 的有限元素之 Equivalence 及 Verify 假壳的法线方向，使所有之法线指向外，以致假壳所构成之耦合面无破洞；然后，藉由该耦合面及 Loads/BCs/Create/Coupling 菜单之设定，执行流固的一般耦合计算。图5 显示波动的自由水面与水下之流速分量；但未显示水下固体物之移动。换言之，Insight 工具虽显示欧拉域内部之动态状况，但不显示欧拉域范围内之拉格朗日固体物的运动状况。此为目前 MSC.Patran的限制。5.1.3 造波板与之波浪水槽本计算模型是 20m(x 向)×10m(y 向)

31、15;3m(z 向)的波浪水槽，内含 7m 高(y 向)、3m 宽(z 向)与 0.5m 厚(x 向)之造波板。波浪水槽之内部是欧拉域，分成 20×20×3 个体元素；造波板是拉格朗日域，分成 1×14×3个体元素，如图6。如图所示的圆点代表造波板之外部包裹假壳的耦合面，使拉格朗日域与欧拉域产生一般耦合计算。造波板之所有元素节点均有正弦时变性之 x 向的速(velocity BC)，最大速度为 1 m/sec，周期为 1 sec。图 6 模拟造波板与波浪水槽之有限元素格网欧拉域之下部 5 m 是水，材料是 LinFluid；上部是空气，当作空域处理。整个

32、欧拉域之上游面(x= 0 m)与下游面(x= 20 m)均给与予flow 边界条件。水与空域之性质均是 Hydro(Peler1)，搭配 MSC.Patran的Loads/BCs 菜单下之欧拉起始条件(Init. Cond. Euler)，但仅给予水域 x 向之初速度 0.1 m/sec。在代表水槽之欧拉域的下游面之 flow 边界条件，系给予流速0.1 m/sec (x 向)；若欧拉域之上游面之 flow 边界条件，也是 x 向流速0.1 m/sec，则第 1.0 sec 时的上游边界附近会有水面特别壅高之现象，如图7 所示。若上游面之边界流速改为0.1 m/sec，则该壅高水面就可消除。图

33、 7 水槽内的造波板(在中央者)引致约第 0.8 sec 时之水面壅高及流速的空间分布云纹图5.1.4 海面上之高速物体撞击混凝土墙本例之目的在显示一个计算模型可包括两组以上的一般耦合，且被包裹耦合面的固体物可有一部分在水域内与模拟空气之空域内。本例也探讨一般耦合与固体对固体之撞击同在一个计算模型内的状况。计算模型之几何体是三个较小的长方体(solid)在一个较大长方体内，如图8。较小的长方体由左至右，分别称为 maker、body 与 struc，均是用拉格朗日法描述之固体。尺寸 0.5m×7m×3m 的 maker 是造波板，用来模拟海洋环境；1m×1m

34、15;1m的body，用来模拟高速运动之水上物体；尺寸 1m×9m×3m的struc用来模拟混凝土墙。尺寸为 20m×10m×3m 之较大的长方体称为 euler，系用来模拟包括水与空气之欧拉域，其范围的 x 向系由 0 m 至 20 m。计算起始时的 maker 与 body 之 x 向净距是 5.5 m，body 与 struc 之 x 向净距是2.0 m。水上物体 body 的初速度是 60 m/sec，朝 x 坐标轴的方向。图 8 含自由水面的水域、造波板、高速固体物与混凝土墙构成之模型注：最大之长方体是包含水与空气之欧拉域，较小长方体是拉革朗日

35、域。图8 由左至右，依序为：由 1×14×3 个三维体元素所构成造波板 maker、由 2×2×2 个三维体元素所构成之高速物体 body、由 4×18×6 个三维体元素构成之混凝土墙 struc(详如图9)。虽 maker、body 与 struc均是弹塑性的通用材料 ElasPlas (DMAT)，但 struc 的材料参数与前两者不同且含材料破坏模型。对包含水与空域之欧拉域，其在 x 值等于 0 m 与 20 m 之边界面分别全面给予水流进与流出的 flow 边界条件，且均有 x 向之流速 0.1m/sec，材料对象仅是水。对其

36、它未给予边界条件之四个面，则预设为固定边界。同时，对所有的欧拉元素，给予欧拉之起始条件，含 “shape”、 “initial values”与 “region definition”三部分。水域内之初流速为 0.1 m/sec。至于空域，则不给予初流速。图 9 高速固体物与混凝土墙之有限元素格网本例对仿真造波板之 maker 各元素节点，给予正弦运动之velocity 边界条件，最高速度为 1 m/sec，周期为1秒。正弦运动速度数据输入过程系经由 MSC.Patran 的 “Loads/BCs/Create/Velocity”菜单，配合 Field 菜单所建置之各时间的正弦函数值。对高速固

37、体物之各元素节点，系给予初速度为 x 向 60 m/sec之荷载条件(Loads/BCs)。对混凝土墙的最底面之所有元素节点，给予 x、y 与 z 三个方向之位移量为零的 displacement 边界条件(相当于节点固定之条件)；其前处理过程，可针对最底面之整个几何面(即 solid 3.3)，程序会自动地选定相应的元素节点。对高速固体物与混凝土墙之接触及撞击的模拟，系给予主面与从属面接触(master-slave surface contact)之条件。主面系位于混凝土墙上，从属面则位于高速运动的固体物上。此外，也在 “Execution Controls/Inertial Loads”之

38、菜单下，给予重力加速度之大小(9.81m/sec2)与其作用方向。在 MSC.Patran 的 “Execution Controls/Execution Control Parameters”下之菜单内，系给予计算终止时间(或终止步数End Step) 0.07 sec、起始时间步长t 为 3.9×10-5 sec、容许之最小时间步长(Min. t) 为 1.0×10-14 sec。于计算之时间步长将短于 Min. t时，计算即告终止。图10 与图11 系应用 MSC.Patran 的 Results 工具绘出计算模型的外围各欧拉元素之FMAT值的云纹图，它呈现高速物体引

39、致之自由水面位置的变动与水面波。两图之时间分别为开始计算后之第 0.07 sec与 0.20 sec 时，图上约中央处是造波板，靠右者是混凝土墙；因高速物体系位于欧拉域的内部，故不能经由 Results 工具看到它。于第 0.06 sec时，初速度 60 m/sec之固体物已贯穿混凝土墙。图12 显现混凝土墙受到固体物接触及贯入后，混凝土墙之材料破坏与塑性应变的分布情形。因水与混凝土之阻尼效应，在接触及撞击过程内，固体物的速度会随时间而降低。图 10 于第 0.07 sec 时之 FMAT 云纹图显示自由水面波动情形图 11 于第 0.20 sec 时之 FMAT 云纹图显示自由水面的波动图

40、12 高速固体物贯穿混凝土墙与混凝土墙之塑性应变分布图 13 Insight 工具显示计算模型内部与 FMAT 等值面所代表之自由水面图13 是采用 MSC.Patran 的 Insight 工具所产生之影像，它显示第 0.04 sec 时的自由水面之波动与水下的流速向量。产生该图像文件的过程是：于使用 MSC.Patran 的 “Insight Application”菜单之 “Action/Tool” 制 作 Insight 工 具 之 后 ， 选 取 menu bar 上之 “insight control/animation control” ( 动 画 控 制 ) 下的 “anima

41、tion setup” (动画建置)，然后，点选 “save frames to files”，再点选 “animate”，如此，若有 8 幅动画之画面(frames)，就会建立且储存 8 个图像文件。如图13 所示，Insight 工具不显示欧拉域内之拉格朗日的固体物之移动，此为目前之 Insight 工具的限制。于第 0.03 sec 时，高速固体物尚未接触混凝土墙。于第 0.04sec时，高速固体物约贯入混凝土墙厚度之一半。于第 0.06 sec 时，高速固体物已完全贯穿混凝土墙，且离开混凝土墙。固体物之速度引致流速变化。在 0.07 sec 的过程内，流场内元素得最高流速变化如表1。表

42、 1 各时间之流场内的元素最高流速时间(sec)0.000.010.020.030.040.050.060.07流速(m/sec)0.136.660.538.128.523.722.018.15.2 应用流固耦合功能形成流场较复杂之固体边界由于含不规则之固体边界之流场空间不能划分成矩形格网(iso-mesh)如图14，不利于给予含自由液面之欧拉的起始条件，为避免此项限制，可将流体空间先采用 simple solid (即四面体、五面体或六面体)定义成欧拉域；然后，在该欧拉域内，迭上不规则之拉格朗日 solid，形成不规则之流场的固体边界。如此，对流场空间仍可使用矩形格网。图 14 有限元素的矩

43、形格网之平面5.2.1 上游侧与下游侧水深相等情况之潜堰包含水与空气的欧拉域之矩形体系长 20 m(x 向)、高 3 m(y向)及宽 3 m(z 向)，分成 20×12×3 个元素。欧拉域的下部 1.5 m 是水域，上部之 1.5 m是代表空气之空域。堰是主轴与流向垂直之结构物，用以抬高水位，俾藉重力取水及引水。拉格朗日的混凝土堰体尺寸是 1m×1m×3m，其上游面位于 x= 9 m 处，堰体被划分成 2×2×3 个元素，如图15，堰体上的饼图样代表它被一般耦合的接口所包裹着经由 MSC.Patran 的 Loads/BCs 菜单，在

44、欧拉域的上游与下游边界面，分别给予 flow 边界条件；针对水之流进与流出，给予流速10 m/sec。欧拉域下部 1.5 m 深之水域也均给与初流速 10 m/sec。另由Analysis 之 “Inertial Loads”表，给予重力加速度 9.81 m/sec2，作用在负 y 轴向。相关的是须由 Loads/BCs，对堰体底面节点给予位移量为零之边界条件；否则，堰体会因重力作用而往下掉落至欧拉域的范围外。图 15 拉格朗日域(含耦合面之堰体)及欧拉域(水与空气)的元素格网图 16 第 0.026762 sec 时的越过潜堰的水域之流速空间分布图 17 采用 FMAT 为 0.002488

45、 之等值面显示越过潜堰之水面形状本例之计算起始时间步长为 1×10-4 sec，容许最小时间步长为 1×10-14 sec。根据计算结果，图16 是 Cycle 240 (第0.026762 sec)时之流速空间分布，在堰体附近之流速提高，最高达 13.1 m/sec。由于水下潜堰的效应，越过潜堰处之水面会起伏，如图17 所示；该图显示每个元素内之 FMAT 值等于 0.002488 的等值面，其接近自由水面。应用 MSC.Dytran 之优点是不但可模拟分析流况，且可同时模拟计算混凝土堰体之应力与应变量。惟 MSC.Dytran 之缺点是针对由上游往下游之冲击水流，适宜用

46、在瞬时(transient)行为之模拟计算；想由瞬时计算至稳态(steady)，往往会遇到流速变得很高(即计算发散)与计算时间步长t 变得很小，以致计算中止的情况。本例于 Cycle 570(第 0.063588 sec)以后，就有计算显著发散的现象，故该时间以后之计算结果不可信。图 18 水域与空域之欧拉元素及堰体之拉格朗日元素格网5.2.2 上游侧与下游侧水深不等情况之固定堰关于本例的整个计算模型之材料，除水与混凝土两种外，还有空气是性质 Hydro (Peuler1)及欧拉体(solid)的空域。本例的欧拉域系长 20 m(x 向)、宽 3 m(y 向)及高 10 m(z 向)之长方体。

47、而采用拉格朗日法描述之堰体系高 1.5 m(y 向)、厚为 1 m(x 向)及长 3 m(z 向)。堰体完全在欧拉域之内，且堰底位于欧拉域之底面。欧拉域之元素数目为 20×20×3，共有1,200个矩形元素，拉格朗日之元素数目为 2×3×3，如图18 所示。经由 MSC.Patran 之 Loads/BCs 菜单下的 flow 边界条件，给予整个欧拉域上游边界面(x = 0 m) 之流速 4 m/sec 与下游边界面(x = 20 m)之流速为零，flow 边界条件之材料均仅针对水。至于图18 下部中央之堰体外部符号，则显示堰体被一般耦合面包裹着，使堰体

48、与欧拉域有流固的一般耦合效应。图 19 起始 Euler1 群组(图左)与 Euler2 群组(图右)的水域与空域定义图 20 欧拉域于第 0.58 sec 时之 FMAT 数值图显示瞬时的水面形状在 MSC.Patran 的 Analysis 菜单下，对整个计算模型给予重力荷载 9.81 m/sec2，重力方向是负 y 轴向。以堰厚中央位置(x= 10 m)为界线，将将欧拉域分成 Euler1与 Euler2 两个群组(group)，Euler1 由 x= 0 m 至 x= 10 m，Euler2由 x= 10 m 至 x= 20 m，俾给予两部分不同之欧拉起始条件(Init. Cond.E

49、uler.)，如图19。Euler1 的下部 6 m 是水域，起始流速为 4 m/sec；上部 4 m 是空域。Euler2 的下部 1 m 是水域，起始流速为零；上部 9 m 是空域。图20 是 Cycle 6300(第 0.58 sec)时之水流剖面，那时候，混凝土堰体之最大有效应力为 1.58×104 N/m2。本例所用之计算起始时间步长t 为1×10-5 sec，容许之Min.t 为 3×10-16 sec。计算结果，于时间超过 Cycle 6300(第 0.576261sec)，计算时间步长小于 Min. t，因而计算中止。本例之重点在于： 使用 Eul

50、er1 与 Euler2 两个群组，去定义欧拉域的上游部分与下游部分的不同水深之边界条件与起始条件。对此例，值得思考是： 下游部分的 Euler2 之水域的起始流速若不是零，而是 4 m/sec或其他数值，则所代表之意义及计算结果之差异程度如何。5.2.3 上游側與下游側水深不等情況之階梯式渠道如同前节之固定潜堰，所用欧拉域系长 20 m(x 向)、宽 3 m(y向)及高 10 m(z 向)之长方体。惟有四个长方体连在一起及构成拉格朗日域。每个长方体均厚 2 m(x 向)及宽 3 m(z 向)，而高度(y 向)分别为 4 m、3m、2 m与1 m，形成阶梯状，如图21 之左下方。欧拉域类似前节

51、之固定潜堰，将 x等于0 m至10 m之范围定义为 Euler1 群组，且定义下部分 5 m 深为水，6 m/sec 之上游面 flow边界条件与该水域初流速；上部分之5 m深为代表空气之空域。图 21 四个几何长方体构成一个阶梯状之有限元素体与欧拉域之阶梯状边界类似地，欧拉域 x 等于 10 m 至 20 m 之范围定义为 Euler2群组，且定义下部分 1 m 深为水，3 m/sec 之下游面 flow 边界条件与该水域初流速；上部分之 9 m 深为代表空气之空域。前述之构成拉格朗日域的四个长方体，材料使用混凝土，它们一起建置有限元素格网，直接用所形成之阶梯状体元素群之外露的元素面作为耦合

52、面，进行 Loads/BCs 菜单下之一般耦合建置。此方式省去在拉格朗日之体元素外部建构 dummy shell 的耦合面、equivalence 与 verifynormal 等程序，仍可完成流固之一般耦合，如图21；然而，此方式是否到处行得通(universal)仍待进一步探讨。如图21 所示者是 Cycle 780 时的水域范围图，大于 Cycle 780 者之时间步长小于容许之 Min.t，因而计算中止。本例之重点目的在于用阶梯状的拉格朗日之体元素，形成欧拉域之固定边界，而欧拉域仍可采用矩形之元素格网，利于设定欧拉域的 flow边界条件与内部之起始条件。此外，如同前节，对欧拉域内的不同

53、范围可分为两个以上之群组，俾设定不同之起始条件与 flow 边界条件。5.2.3 固定开度之水利闸门固定开度之水利闸门(hydraulic gate)也是欧拉域的固定边界之一种。本例也是流固之一般耦合分析的应用。流体是水与代表空气的空域，固体是闸门。垂提式闸门开度固定为 2 m，闸门上游侧之水深为 5 m，水之流速是在 x 向为 20 m/sec。本例之欧拉域的几何是 x 向为 20 m 长，y 向为 10 m 高，z向为 3 m 宽，采用四节点六面体之有限元素共 20×20×3 个，性质为欧拉之体元素 Hydro (Peuler1)。在几何体之上游与下游边界，分别给予 f

54、low 边界条件，具有 x 向之流速 20 m/sec。该欧拉域在 x 向，分为各长 10 m 之两个群组(group)，称为 euler1 与 euler2。模拟闸门之几何体是 x 向为 0.3 m 长，y 向为 8 m 高，z 向为 3 m 宽，采用四节点六面体之有限元素共 1×16×3 个，性质为拉格朗日之体元素。闸门之几何体外部包裹假壳(dummy shell)元素所构成之闭合面，经有限元素之 verify，让各法线(normal)之方向均朝外，确保无破洞，作为流固耦合面。对 euler1 群组之上部 5 m 深，定义为 void1 域之空隙；下部5 m深，定义为w

55、ater1域之水，且初始有x向之流速20 m/sec。euler2 群组，则全部定义为 void2 域之空隙，不给予初始流速。图 22 第 0.027737 sec 时之闸门附近的水流纵剖面本例目的之一是验证：不但 water1 域之水面(即 water1 域与 void2 域之接口)会随着时间变动，且 water1 域之水可横向地往 void2域扩展。如图22 是 cycle 1300 (时间为第 0.027737 sec)之水体立剖面。本例提出群组(group)之应用；若非应用群组，则无法对上、下游两部分之欧拉域设定不同之起始条件，即上游之 euler1 域的水深为5 m，下游之 eule

56、r1 域的水深为零。重点之一是： 闸门体须与 water1 域重叠一部分，而与 void2域重叠数。本例系将 0.3 m 厚之闸门体完全与 water1域重叠之下游面重合。本例的难处是：因用显式的时间积分法，计算之时间步长愈来愈小，以致算不下去，无法算至流场达稳态之时。譬如，本例计算迄 cycle 300 (时间仅第 0.027737 sec)，接着之计算时间步长小于 1.0×10-14sec，而中止计算。一般而言，MSC.Dytran 对固体不动，但流体运动之模拟，其流速计算较易发散(diverge)，譬如本例，流场上游与下游的边界流速与初始的 x 向之流速均为 20 m/sec，

57、计算达 cycle 130 (时间仅第0.002639 sec)，最大流速即发散至 7.73×103 m/sec。六、 任意耦合計算之應用例大致上，任意耦合分析可应用在气体或液体的流场驱动固体物之互制状况。执行任意耦合计算时，拉格朗日域与欧拉域不须迭接，只须相连接。本文所述之两个例子包括隧洞内之气爆压作用在钢板上、隧洞内之水流推动固体物等两项如下。6.1 隧洞内之气爆压作用在钢板上本例是进行包裹在 3m×3m×3m 几何立方体(solid 1)外部的壳元素之拉格朗日域如图23 之左半部，以及模拟空气的 3m×3m×3m 几何立方体(solid

58、2)之欧拉域如图23 之右半部，两者之任意耦合分析。两者的相连接处是任意耦合面，即 solid 1 的左侧面与 solid 2 的右侧面之 x 坐标分别为 x 等于零与 9 m。初始的爆炸源之半径为 0.05 m，中心点系位于 solid2 几何体的中央，即 x= 4.5 m、y= 1.5 m 与 z= 1.5 m 处；该坐标数值系采用文本编辑器在 *.dat 档内的 TICEUL 卡上，予以键入。壳元素是弹塑性的通用材料 ElasPlas(DMAT)之钢，欧拉域的材料是符合理想气体律之空气。图 23 第 0.008099 sec 时的壳形固体物与爆炸之空气外围的位移量分布图 24 第 0.0

59、08099 sec 时之壳体变形与位移量之分布任意耦合面之建置系经由 MSC.Patran 的主面与从属面之接触(Contact/Element Uniform/Option= Master-Slave Surface)过程。主面是二维壳形体之右侧面段，可选相应之几何面solid 1.2，予以设定。从属面则是三维欧拉体元素群之左侧外露面所构成之面段，可选相应之几何面solid 2.1，予以设定。根据历时的流固互制模拟计算结果，Cycle 200 (第 0.008099 sec)时，壳形体接触爆炸压的面之最大位移量是 0.372m，如图24 之云纹图所示。6.2 隧洞内之水流推动固体物本例系应用任意耦合功能，仿真隧洞内之水流推动固体物的行为。如图25，左边是计算模型的上游侧，右边是下游侧。计算模型系由两个几何长方体与一个立方体所构成，即近乎在中央之立方体