流固耦合模拟分析之探讨

./MSC.Dytran流固耦合模拟分析之探讨摘要：流场会驱动固体；而固体的运动也可能驱动流体,甚至引起流场振动。这称为流固互制或流固耦合。MSC.Dytran软件的流固耦合功能包含一般耦合、任意偶合,且采用拉格朗日法与欧拉法分别描述固体与流体的运动。Dytran中流固耦合计算种类拉格朗日的元素节点依附在材料上,节点随着材料质点作运动,故各物理量也作用在节点上随材料流动而变化。相反,除任意耦合外,欧拉元素网格与节点不随时间而变,其物理量虽也作用在欧拉元素节点上,但对于通过欧拉元素面的各时间的质量、动量与能量的进与出,加之模拟,即模拟元素面的材料流,而不模拟各材料质点的时间历程。因为对一般固体材料,要模拟各材料质点的时间历程,因此大多用拉格朗日法。而对于流体不需要模拟材料质点的时间历程,故采用欧拉法,MSC.Dytran的欧拉法需用三维的计算域、三维的体元素与DMAT通用材料。此外,欧拉法容许一个元素内含有两种以上的材料,这就是模拟计算材料流的扩散与混合行为欧拉法模拟计算材料流的扩散与混合行为。欲推广应用型的计算软件,需有充分的应用X例。关于流固耦合的模拟,除MSC.Dytran中文X例手册〔1999b〕与ExampleProblemManual<2001a>论述到一般耦合与任意耦合的应用X例之外,本文进一步探讨其他的应用。而搭配的MSC.Patran软件除有Results的后处理工具外；该软件的Insight工具,能透明地看到体元素所构成的欧拉域内部,因此,更需用Insight,去展示欧拉域内部的流固耦合计算结果。Dytran中流固耦合计算种类欧拉法模拟计算材料流的扩散与混合行为前言MSC.Dytran之流固耦合计算功能大致上包括一般耦合<generalcoupling>和任意耦合<ALE,ArbitraryLagrangian-Eulerian>两大类。欲推动一项泛用型计算软件被广泛地应用,须有可供参考的文件与丰富的应用案例。而MSC.Dytran之中X例手册包括一个一般耦合和两个任意耦合的应用X例。MSC.Dytran之ExampleProblemManual<2001>也包括两个一般耦合、两个任意耦合的X例。本文之要旨是经由案例,探讨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>,大多是拉格朗日的固体在欧拉的流场X围内运动,即拉格朗日域驱动欧拉域一般耦合应用在Lagrange的固体在Euler的流场X围内运动。；流场虽有速度,但代表流场的欧拉格网系固定与不受拉格朗日的固体之影响。换言之,在流固耦合过程内,欧拉格网不移动,也不变形。一般耦合计算之前处理大多用封闭之假壳<dummyshell>耦合面隔开拉格朗日域与欧拉域；于起始计算之时,拉格朗日的固体至少有一微小量<譬如,0.001m>重叠在欧拉域的X围内,且固体、流体或两者有运动,才能启动流固耦合之计算器制。若拉格朗日域与欧拉域毫无重叠,则无法起动流固耦合计算。当然,拉格朗日的固体可完全位于欧拉域内,不因拉格朗日的固体运动而使欧拉格网移动或变形。此种流固耦合最适用在「固体驱动流体」现象之模拟,譬如,水上飞机降落在水面上之行为、水面下之物体的运动所引致之流体动力行为、隧道内之高速车辆引致的气动力行为、管道内的固体活塞推动气体等。此外,一般耦合也可应用在具有不规则的固定固体边界之流场模拟计算,即将不动的固体边界与流场之关联视为流固耦合现象。一般耦合应用在Lagrange的固体在Euler的流场X围内运动。一般耦合须在拉格朗日的几何体之外露面<几何面段或元素表面,经由SURFACE卡建立>,定义封闭的假壳元素群,作为流固耦合面,形成闭合之体积,完全隔开拉格朗日域与欧拉域；且对假壳元素,须进行Equivalence与Verify法线方向之前处理。耦合面虽为虚拟,但也须输入物理性质<Properties>为2D/SummyShell。然后,对耦合面的全部假壳元素,使用MSC.Patran软件的Coupling等,就可完成相关的前处理。在计算过程内,耦合面是欧拉流场格网的边界,即假壳元素之耦合面与仿真流体之体元素相接触<Contact>,MSC.Dytran先计算欧拉流体施加在耦合面上之荷载,然后,耦合面使拉格朗日格网产生应力与变形。至于输入档相关的Couple卡上之参数COVER,则定义被耦合面所包裹者〔即不含欧拉领域者〕是属内部<INSIDE>,还是外部<OUTSIDE>。譬如,石头掉入水中之模拟,其COVER是INSIDE。汽车的安全气囊的模拟之COVER是OUTSIDE。4.2任意耦合任意耦合主要是用来仿真欧拉域驱动拉格朗日域,在该耦合过程内,于计算模型承受来自欧拉域的荷载之后,拉格朗日的元素与欧拉的体元素之网格均可能变形或移位。任意耦合最适宜用来仿真「流体驱动固体」,譬如,容器内之气体爆炸、空中之鸟体撞击飞机壳体、气爆压作用在壳形固体物上、管路内的流体推动固体、因压力波之钢管的膨胀或收缩等之动力行为。任意耦合应用X围任意耦合应用X围任意耦合之流固界面或耦合面不须是封闭的面,可为两片以上的不相连接之面段,也可为一个拉格朗日域对应于两个不相连接之欧拉域。因目前的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"菜单的"OutputRequest"下,建立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工具虽显示欧拉域内部之动态状况,但不显示欧拉域X围内之拉格朗日固体物的运动状况。此为目前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用来模拟混凝土墙。尺寸为20m×10m×3m之较大的长方体称为euler,系用来模拟包括水与空气之欧拉域,其X围的x向系由0m至20m。计算起始时的maker与body之x向净距是5.5m,body与struc之x向净距是2.0m。水上物体body的初速度是60m/sec,朝x坐标轴的方向。图8含自由水面的水域、造波板、高速固体物与混凝土墙构成之模型注：最大之长方体是包含水与空气之欧拉域,较小长方体是拉革朗日域。图8由左至右,依序为：由1×14×3个三维体元素所构成造波板maker、由2×2×2个三维体元素所构成之高速物体body、由4×18×6个三维体元素构成之混凝土墙struc<详如图9>。虽maker、body与struc均是弹塑性的通用材料ElasPlas<DMAT>,但struc的材料参数与前两者不同且含材料破坏模型。对包含水与空域之欧拉域,其在x值等于0m与20m之边界面分别全面给予水流进与流出的flow边界条件,且均有x向之流速0.1m/sec,材料对象仅是水。对其它未给予边界条件之四个面,则预设为固定边界。同时,对所有的欧拉元素,给予欧拉之起始条件,含"shape"、"initialvalues"与"regiondefinition"三部分。水域内之初流速为0.1m/sec。至于空域,则不给予初流速。图9高速固体物与混凝土墙之有限元素格网本例对仿真造波板之maker各元素节点,给予正弦运动之velocity边界条件,最高速度为1m/sec,周期为1秒。正弦运动速度数据输入过程系经由MSC.Patran的"Loads/BCs/Create/Velocity"菜单,配合Field菜单所建置之各时间的正弦函数值。对高速固体物之各元素节点,系给予初速度为x向60m/sec之荷载条件<Loads/BCs>。对混凝土墙的最底面之所有元素节点,给予x、y与z三个方向之位移量为零的displacement边界条件<相当于节点固定之条件>；其前处理过程,可针对最底面之整个几何面<即solid3.3>,程序会自动地选定相应的元素节点。对高速固体物与混凝土墙之接触与撞击的模拟,系给予主面与从属面接触<master-slavesurfacecontact>之条件。主面系位于混凝土墙上,从属面则位于高速运动的固体物上。此外,也在"ExecutionControls/InertialLoads"之菜单下,给予重力加速度之大小<9.81m/sec2>与其作用方向。在MSC.Patran的"ExecutionControls/ExecutionControlParameters"下之菜单内,系给予计算终止时间<或终止步数EndStep>0.07sec、起始时间步长∆t为3.9×10-5sec、容许之最小时间步长<Min.∆t>为1.0×10-14sec。于计算之时间步长将短于Min.∆t时,计算即告终止。图10与图11系应用MSC.Patran的Results工具绘出计算模型的外围各欧拉元素之FMAT值的云纹图,它呈现高速物体引致之自由水面位置的变动与水面波。两图之时间分别为开始计算后之第0.07sec与0.20sec时,图上约中央处是造波板,靠右者是混凝土墙；因高速物体系位于欧拉域的内部,故不能经由Results工具看到它。于第0.06sec时,初速度60m/sec之固体物已贯穿混凝土墙。图12显现混凝土墙受到固体物接触与贯入后,混凝土墙之材料破坏与塑性应变的分布情形。因水与混凝土之阻尼效应,在接触与撞击过程内,固体物的速度会随时间而降低。图10于第0.07sec时之FMAT云纹图显示自由水面波动情形图11于第0.20sec时之FMAT云纹图显示自由水面的波动图12高速固体物贯穿混凝土墙与混凝土墙之塑性应变分布图13Insight工具显示计算模型内部与FMAT等值面所代表之自由水面图13是采用MSC.Patran的Insight工具所产生之影像,它显示第0.04sec时的自由水面之波动与水下的流速向量。产生该图像文件的过程是：于使用MSC.Patran的"InsightApplication"菜单之"Action/Tool"制作Insight工具之后,选取menubar上之"insightcontrol/animationcontrol"<动画控制>下的"animationsetup"<动画建置>,然后,点选"saveframestofiles",再点选"animate",如此,若有8幅动画之画面<frames>,就会建立且储存8个图像文件。如图13所示,Insight工具不显示欧拉域内之拉格朗日的固体物之移动,此为目前之Insight工具的限制。于第0.03sec时,高速固体物尚未接触混凝土墙。于第0.04sec时,高速固体物约贯入混凝土墙厚度之一半。于第0.06sec时,高速固体物已完全贯穿混凝土墙,且离开混凝土墙。固体物之速度引致流速变化。在0.07sec的过程内,流场内元素得最高流速变化如表1。表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,不利于给予含自由液面之欧拉的起始条件,为避免此项限制,可将流体空间先采用simplesolid<即四面体、五面体或六面体>定义成欧拉域；然后,在该欧拉域内,迭上不规则之拉格朗日solid,形成不规则之流场的固体边界。如此,对流场空间仍可使用矩形格网。图14有限元素的矩形格网之平面5.2.1上游侧与下游侧水深相等情况之潜堰包含水与空气的欧拉域之矩形体系长20m<x向>、高3m<y向>与宽3m<z向>,分成20×12×3个元素。欧拉域的下部1.5m是水域,上部之1.5m是代表空气之空域。堰是主轴与流向垂直之结构物,用以抬高水位,俾藉重力取水与引水。拉格朗日的混凝土堰体尺寸是1m×1m×3m,其上游面位于x=9m处,堰体被划分成2×2×3个元素,如图15,堰体上的饼图样代表它被一般耦合的接口所包裹着经由MSC.Patran的Loads/BCs菜单,在欧拉域的上游与下游边界面,分别给予flow边界条件；针对水之流进与流出,给予流速10m/sec。欧拉域下部1.5m深之水域也均给与初流速10m/sec。另由Analysis之"InertialLoads"表,给予重力加速度9.81m/sec2,作用在负y轴向。相关的是须由Loads/BCs,对堰体底面节点给予位移量为零之边界条件；否则,堰体会因重力作用而往下掉落至欧拉域的X围外。图15拉格朗日域<含耦合面之堰体>与欧拉域<水与空气>的元素格网图16第0.026762sec时的越过潜堰的水域之流速空间分布图17采用FMAT为0.002488之等值面显示越过潜堰之水面形状本例之计算起始时间步长为1×10-4sec,容许最小时间步长为1×10-14sec。根据计算结果,图16是Cycle240<第0.026762sec>时之流速空间分布,在堰体附近之流速提高,最高达13.1m/sec。由于水下潜堰的效应,越过潜堰处之水面会起伏,如图17所示；该图显示每个元素内之FMAT值等于0.002488的等值面,其接近自由水面。应用MSC.Dytran之优点是不但可模拟分析流况,且可同时模拟计算混凝土堰体之应力与应变量。惟MSC.Dytran之缺点是针对由上游往下游之冲击水流,适宜用在瞬时<transient>行为之模拟计算；想由瞬时计算至稳态<steady>,往往会遇到流速变得很高<即计算发散>与计算时间步长∆t变得很小,以致计算中止的情况。本例于Cycle570<第0.063588sec>以后,就有计算显著发散的现象,故该时间以后之计算结果不可信。图18水域与空域之欧拉元素与堰体之拉格朗日元素格网5.2.2上游侧与下游侧水深不等情况之固定堰关于本例的整个计算模型之材料,除水与混凝土两种外,还有空气是性质Hydro<Peuler1>与欧拉体<solid>的空域。本例的欧拉域系长20m<x向>、宽3m<y向>与高10m<z向>之长方体。而采用拉格朗日法描述之堰体系高1.5m<y向>、厚为1m<x向>与长3m<z向>。堰体完全在欧拉域之内,且堰底位于欧拉域之底面。欧拉域之元素数目为20×20×3,共有1,200个矩形元素,拉格朗日之元素数目为2×3×3,如图18所示。经由MSC.Patran之Loads/BCs菜单下的flow边界条件,给予整个欧拉域上游边界面<x=0m>之流速4m/sec与下游边界面<x=20m>之流速为零,flow边界条件之材料均仅针对水。至于图18下部中央之堰体外部符号,则显示堰体被一般耦合面包裹着,使堰体与欧拉域有流固的一般耦合效应。图19起始Euler1群组<图左>与Euler2群组<图右>的水域与空域定义图20欧拉域于第0.58sec时之FMAT数值图显示瞬时的水面形状在MSC.Patran的Analysis菜单下,对整个计算模型给予重力荷载9.81m/sec2,重力方向是负y轴向。以堰厚中央位置<x=10m>为界线,将将欧拉域分成Euler1与Euler2两个群组<group>,Euler1由x=0m至x=10m,Euler2由x=10m至x=20m,俾给予两部分不同之欧拉起始条件<Init.Cond.Euler.>,如图19。Euler1的下部6m是水域,起始流速为4m/sec；上部4m是空域。Euler2的下部1m是水域,起始流速为零；上部9m是空域。图20是Cycle6300<第0.58sec>时之水流剖面,那时候,混凝土堰体之最大有效应力为1.58×104N/m2。本例所用之计算起始时间步长∆t为1×10-5sec,容许之Min.∆t为3×10-16sec。计算结果,于时间超过Cycle6300<第0.576261sec>,计算时间步长小于Min.∆t,因而计算中止。本例之重点在于：使用Euler1与Euler2两个群组,去定义欧拉域的上游部分与下游部分的不同水深之边界条件与起始条件。对此例,值得思考是：下游部分的Euler2之水域的起始流速若不是零,而是4m/sec或其他数值,则所代表之意义与计算结果之差异程度如何。5.2.3上游側與下游側水深不等情況之階梯式渠道如同前节之固定潜堰,所用欧拉域系长20m<x向>、宽3m<y向>与高10m<z向>之长方体。惟有四个长方体连在一起与构成拉格朗日域。每个长方体均厚2m<x向>与宽3m<z向>,而高度<y向>分别为4m、3m、2m与1m,形成阶梯状,如图21之左下方。欧拉域类似前节之固定潜堰,将x等于0m至10m之X围定义为Euler1群组,且定义下部分5m深为水,6m/sec之上游面flow边界条件与该水域初流速；上部分之5m深为代表空气之空域。图21四个几何长方体构成一个阶梯状之有限元素体与欧拉域之阶梯状边界类似地,欧拉域x等于10m至20m之X围定义为Euler2群组,且定义下部分1m深为水,3m/sec之下游面flow边界条件与该水域初流速；上部分之9m深为代表空气之空域。前述之构成拉格朗日域的四个长方体,材料使用混凝土,它们一起建置有限元素格网,直接用所形成之阶梯状体元素群之外露的元素面作为耦合面,进行Loads/BCs菜单下之一般耦合建置。此方式省去在拉格朗日之体元素外部建构dummyshell的耦合面、equivalence与verifynormal等程序,仍可完成流固之一般耦合,如图21；然而,此方式是否到处行得通<universal>仍待进一步探讨。如图21所示者是Cycle780时的水域X围图,大于Cycle780者之时间步长小于容许之Min.∆t,因而计算中止。本例之重点目的在于用阶梯状的拉格朗日之体元素,形成欧拉域之固定边界,而欧拉域仍可采用矩形之元素格网,利于设定欧拉域的flow边界条件与内部之起始条件。此外,如同前节,对欧拉域内的不同X围可分为两个以上之群组,俾设定不同之起始条件与flow边界条件。5.2.3固定开度之水利闸门固定开度之水利闸门<hydraulicgate>也是欧拉域的固定边界之一种。本例也是流固之一般耦合分析的应用。流体是水与代表空气的空域,固体是闸门。垂提式闸门开度固定为2m,闸门上游侧之水深为5m,水之流速是在x向为20m/sec。本例之欧拉域的几何是x向为20m长,y向为10m高,z向为3m宽,采用四节点六面体之有限元素共20×20×3个,性质为欧拉之体元素Hydro<Peuler1>。在几何体之上游与下游边界,分别给予flow边界条件,具有x向之流速20m/sec。该欧拉域在x向,分为各长10m之两个群组<group>,称为euler1与euler2。模拟闸门之几何体是x向为0.3m长,y向为8m高,z向为3m宽,采用四节点六面体之有限元素共1×16×3个,性质为拉格朗日之体元素。闸门之几何体外部包裹假壳<dummyshell>元素所构成之闭合面,经有限元素之verify,让各法线<normal>之方向均朝外,确保无破洞,作为流固耦合面。对euler1群组之上部5m深,定义为void1域之空隙；下部5m深,定义为water1域之水,且初始有x向之流速20m/sec。euler2群组,则全部定义为void2域之空隙,不给予初始流速。图22第0.027737sec时之闸门附近的水流纵剖面本例目的之一是验证：不但water1域之水面<即water1域与void2域之接口>会随着时间变动,且water1域之水可横向地往void2域扩展。如图22是cycle1300<时间为第0.027737sec>之水体立剖面。本例提出群组<group>之应用；若非应用群组,则无法对上、下游两部分之欧拉域设定不同之起始条件,即上游之euler1域的水深为5m,下游之euler1域的水深为零。重点之一是：闸门体须与water1域重叠一部分,而与void2域重叠数。本例系将0.3m厚之闸门体完全与water1域重叠之下游面重合。本例的难处是：因用显式的时间积分法,计算之时间步长愈来愈小,以致算不下去,无法算至流场达稳态之时。譬如,本例计算迄cycle300<时间仅第0.027737sec>,接着之计算时间步长小于1.0×10-14sec,而中止计算。一般而言,MSC.Dytran对固体不动,但流体运动之模拟,其流速计算较易发散<diverge>,譬如本例,流场上游与下游的边界流速与初始的x向之流速均为20m/sec,计算达cycle130<时间仅第0.002639sec>,最大流速即发散至7.73×103m/sec。六、任意耦合計算之應用例大致上,任意耦合分析可应用在气体或液体的流场驱动固体物之互制状况。执行任意耦合计算时,拉格朗日域与欧拉域不须迭接,只须相连接。本文所述之两个例子包括隧洞内之气爆压作用在钢板上、隧洞内之水流推动固体物等两项如下。6.1隧洞内之气爆压作用在钢板上本例是进行包裹在3m×3m×3m几何立方体<solid1>外部的壳元素之拉格朗日域如图23之左半部,以与模拟空气的3m×3m×3m几何立方体<solid2>之欧拉域如图23之右半部,两者之任意耦合分析。两者的相连接处是任意耦合面,即solid1的左侧面与solid2的右侧面之x坐标分别为x等于零与9m。初始的爆炸源之半径为0.05m,中心点系位于solid2几何体的中央,即x=4.5m、y=1.5m与z=1.5m处；该坐标数值系采用文本编辑器在*.dat档内的TICEUL卡上,予以键入。壳元素是弹塑性的通用材料ElasPlas<DMAT>之钢,欧拉域的材料是符合理想气体律之空气。图23第0.008099sec时的壳形固体物与爆炸之空气外围的位移量分布图24第0.008099sec时之壳体变形与位移量之分布任意耦合面之建置系经由MSC.Patran的主面与从属面之接触<Contact/ElementUniform/Option=Master-SlaveSurface>过程。主面是二维壳形体之右侧面段,可选相应之几何面solid1.2,予以设定。从属面则是三维欧拉体元素群之左侧外露面所构成之面段,可选相应之几何面solid2.1,予以设定。根据历时的流固互制模拟计算结果,Cycle200<第0.008099sec>时,壳形体接触爆炸压的面之最大位移量是0.372m,如图24之云纹图所示。6.2隧洞内之水流推动固体物本例系应用任意耦合功能,仿真隧洞内之水流推动固体物的行为。如图25,左边是计算模型的上游侧,右边是下游侧。计算模型系由两个几何长方体与一个立方体所构成,即近乎在中央之立方体solid1、上游侧之长方体solid2与下游侧之长方体solid3。solid2与solid1之上游面相连接,且solid1之下游面与solid3相连接。在几何体solid1之外部6个面各建置2×2片壳元素,即24片壳元素,材料是钢,厚度为5mm,性质是弹塑性的通用材料ElasPlas<DMAT>,用vonMises降伏模型搭配最大塑性应变MPS之破坏模型。建置几何体solid2成为2×2×25个体元素,材料是水。建置solid3成为2×2×20个体元素,是空域。均给予几何体solid2与solid3欧拉体与Hydro<Peuler1>性质,且搭配欧拉之初始条件。在几何体solid2上游面给予flow边界条件,流速为10m/sec。solid3的下游面给予flow边界条件,流速为零。经由MSC.Patran的Loads/BCs菜单下的主面与从属面的接触,设定流固之任意耦合面。该主面是solid1的上游与下游面,从属面包括solid2的下游面与solid3的上游面。而在执行MSC.Patran所形成之*.dat输入文件内,则增加两组ALEGRID指令,分别针对solid2与3从属面上之元素节点,使它们承受荷载之后能移动。图25水流隧洞与位于隧洞中央之固体物的有限元素格网图26推动隧洞内的固体物之水流速度分布云纹图图26所显示者是,Cycle120<第0.003212sec>时之水流元素流速<elementvelocity>分布,最大流速为10.2m/sec,而该时之固体物已有的位移量是0.105m。图27第0.008218sec时之隧洞内的水压空间分布从此例发现,任意耦合面与flow边界条件之作用面不能在同一面上。此外,任意耦合面可为拉格朗日之壳元素的面搭配欧拉之体元素的面,一组是主面,另一组是从属面,经由MSC.Patran之Loads/BCs菜单下之Contact,于形成*.dat文件之后,再采用「记事本」程序修改该文件。在本例内,主面包括分开之solid1的上游面与下游面,从属面包括分开之solid2的下游面与solid3的上游面。本例值得进一步探讨者是：本文系在几何体solid3的下游面给予流速为零的flow边界条件；另一作法是不给flow边界条件,使它变成固定边界。两者之物理意义的差异如何。此外,几何体solid3当作空域或真正考虑成空气域,会有多大差异,也值得加以探讨。七、结论与建议1.任何软件均有其限制,纵使是泛用型的计算软件亦然。MSC.Dytran软件之特点是对时间采用显式之积分法。因显式法于模拟分析延时太长时,往往有流场发散之困扰,其主要征兆是流速变得很高,计算时间步长<或称为计算时间距>变得很小,以致小于容许之min.∆t而使得计算中止。换言之,显式积分法使MSC.Dytran局限在短时间之瞬时<transient>或冲击行为为主的模拟计算上,而不适用在较长时间之稳态<steadystate>行为的模拟。2.MSC.Dytran软件对「流体驱动固体」之流固耦合的模拟,往往使得流场速度的计算易快速发散,计算时间步长变得短于所定之Min.∆t,以致计算中止。因此,MSC.Dytran的流固一般耦合较适宜用在「固体趋动流体」之互制,而较不适宜用在「流体趋动固体」之互制。3.由于Min.∆t的限制,在一般状况下,不是想算到第几cycle,就能算到第几cycle。因不知算到第几cycle之计算时间步长会短于Min.∆t,故须经由多次试算,相当费时。希望将来新版之MSC.Dytran软件可预先提供可算到第几个cycle的相关信息。4.对MSC.Dytran的计算模型而言,只要是固体与固体有接触,就须使用接触<contact>指令。相似的,只要是流体与固体有接触,就须使用流固耦合<generalcouplingorALEcoupling>的指令。5.一般仿真固体的拉格朗日元素与仿真流体的欧拉元素,经由耦合面,才能连结在一起,产生互制作用；否则,纵使流固两者在空间内相连接或相重叠,彼此也相互视若无睹,互不相干,无互制效应。6.MSC.Dytran之文献系混合称呼「流固互制」与「流固耦合」分析。惟严格而论,其计算只宜称为流固互制分析。爰此,所称的一般耦合与任意耦合,宜分别改称为一般互制与任意互制。文献内的其他耦合名词也几乎宜改为互制。惟为避免混乱