行政院國家科學委員會專題研究計畫 成果報告

生醫多相流固偶合模式之發展(第 2 年) 研究成果報告(完整版)

計 畫 類 別 ： 個別型

計 畫 編 號 ： NSC 97-2221-E-216-009-MY2

執 行 期 間 ： 98 年 08 月 01 日至 99 年 10 月 31 日 執 行 單 位 ： 中華大學機械工程學系

計 畫 主 持 人 ： 牛仰堯

計畫參與人員： 碩士班研究生-兼任助理人員：湯漢威 碩士班研究生-兼任助理人員：張哲誠

報 告 附 件 ： 出席國際會議研究心得報告及發表論文

處 理 方 式 ： 本計畫涉及專利或其他智慧財產權，2 年後可公開查詢

中 華 民 國 99 年 11 月 09 日

行政院國家科學委員會補助專題研究計畫成果報告

生醫多相流固偶合模式之發展

計畫類別： „ 個別型計畫 □整合型計畫 計畫編號：NSC 97-2221-E-216-009-MY2

執行期間：2008 年 8 月 1 日至 2010 年 10 月 31 日

執行機構及系所：中華大學機械工程學系

計畫主持人： 牛仰堯 共同主持人：

計畫參與人員：碩士生(兼研究助理) 湯漢威、張哲誠

成果報告類型(依經費核定清單規定繳交)：□精簡報告 „ ^完整報告

本計畫除繳交成果報告外，另須繳交以下出國心得報告：

□赴國外出差或研習心得報告

□赴大陸地區出差或研習心得報告

□出席國際學術會議心得報告

□國際合作研究計畫國外研究報告

處理方式：除列管計畫及下列情形者外，得立即公開查詢

□涉及專利或其他智慧財產權，□一年 „ 二年後可公開查

詢

摘要

本研究為延續前期之工作由國立台灣大學附屬醫院最新進之 MRA 醫療影像設備進行 真人之主動脈掃瞄、並配合影像處理掃描所獲得資料進行主動脈計算流體力學分析，本期 將發展考慮血液動力與血管壁交互作用之流固耦合模式。以整合 Eulerian 和 Lagrangian 之統一座標下應用於人工壓縮不可壓縮流方程式，並發展可變形血管動態模擬並考慮血液 與血管壁之交互作用。在統一座標下，流體方程式和幾何方程式以守恆形式呈現並且於每 個時間歩階同步更新。根據統一座標方法判斷幾何守恆之準確性並且控制網格的速度，避 免移動血管之壁面或是靠近表面之邊界層上網格嚴重變形。本團隊研究證明在 MRA 下動脈 血液流體運動之模擬地可行性，本研究已完成三維統一座標中發展血液動力流體之流固耦 合反應模式，考慮血液與血管壁之交互作用，分析結果將與影像處理掃描所獲得結果互相 驗證、以期更精確模擬出血液流動之行為，作為提供安全血壓值及醫生臨床主動脈手術及 用藥之參考範圍。

為期兩年之計畫：

第一年，我們將推導出統一座標下，二維及三維人工壓縮法之Navier-Stokes方程式和相關 性之特徵系統。將適用的壁面移動模式和複雜血管彈性模式納入於統御方程式作為數值模 式。以MRA 影像處理掃描所獲得的實例中，驗證在統一座標下之流固耦合模式並且應用 於二維之血管流體模擬。第二年，我們將使用MRA影像技術完成人體血管之壁面移動，

之後修改現行的血管壁面模式，更準確的模擬流固耦合及相關的血液流體行為及壓力，剪 力形式。我們也將發展三維血液動力含流固耦合程式並且選擇真實的實例，驗證統一座標 下流固耦合之精確性。

Abstract

To continue our previous efforts on the simulation of a human aortic flow based on the techniques of CFD and MR phase-contrast velocimetry in National Taiwan University Hospital,

a novel unified artificial compressibility solver is initialized based on the unified Eulerian and Lagrangian coordinate transformations will be developed to simulate Hemodynamics of deformable aorta. In additions, the fluid-structure interactions will be considered. Based on the unified coordinates, the flow equations and geometry equations can be expressed in conservation form and updated simultaneously during each time step. Thus, the accurate estimation of geometry conservation and controlling the grid velocity based on the unified approach can avoid severe grid deformation caused by moving vessel walls or boundary layers (considered as slip lines) near the surfaces. Our team work has demonstrated its feasibility in the simulation of artery blood flow motions with MRA; and with the consideration of the future novel fluid-structure interaction model in the current Hemodynamics flow solvers in unified coordinates, our noninvasive simulation of blood flows would become possible.

During this 2-year project:

In the first year, we will derive mathematical formulation of the artificial 2D and 3D incompressible Navier-Stokes Equations in a unified coordinate system and the related eigensystems. Then implement a numerical modeling of governing equation with suitable wall moving model & compliant arterial model. In the validated cases, we will verify the accuracy of the current fluid-structure model in unified coordinates and application on 2D Aortic Flow Simulations. In the second year, we will use MRA imaging techniques to achieve the wall motion of a human aorta, then to modify the current vessel wall model to accurately simulate fluid-structure interaction and the related blood flow behaviors and pressure, shear stress formations. Also we will develop a three-dimensional Hemodynamics code with the fluid- structure interaction and choose validate cases to verify the accuracy of the current fluid- structure model in unified coordinatesTo implement unsteady flow calculations, a dual time stepping strategy including the LU decomposition method is used in the pseudo-time iteration and the second-order accurate backward difference is adopted to discretize the unsteady flow term. The original FORTRAN code is converted to the MPI code and tested on a 64-CPU IBM SP2 parallel computer. The test results show that a significant reduction of computing time in running the model and a near-linear speed up rate is achieved up to 32 CPUs at IBM SP2 processors. The speed up rate is as high as 31 for using IBM SP2 64 processors The test shows efficient of parallel processing to provide prompt simulation of 3D cavity, unsteady dropping airfoil and blood flows in an aortic tube with a linear elastic modelling of wall motion [2,3] is included here .

Keywords: fluid-structured interaction, incompressible flow, parallel computation, Chorin’s artificial compressibility

Introduction

Moving body simulation is one of attractive topics in the Computational Fluid Dynamics (CFD) area. Recently, several techniques of moving body simulation are applied on the study of bio-fluid dynamics problems such as blood flows through vessels and organs, also flapping wings. The key issue to achieve accurate simulation on these topics is required to reply moving grid or dynamic mesh algorithm. One common way to deal with dynamic mesh is rely on the so- called Lagrangian method. Computational cells in the Lagrangian coordinates, on the other hand, are literally fluid particles. Consequently, it is capable of producing sharp slip line resolution due to no convective flux across cell interfaces with minimized numerical diffusion. However, the disadvantage is that the computational cells exactly follow fluid particles always brings severe grid deformation, causing inaccuracy and even breakdown of the computation once the fluid velocity is used as the mesh moving velocity, To prevent this from happening, the most famous Lagrangian method in use at the present time—the arbitrary Lagrangian–Eulerian (ALE) technique—uses continuous rezoning. However, ALE requires continuous interpolations of flow variables and computational geometry that may result in unnecessary numerical inaccuracy.

Recently, to understand the connection between the Lagrangian method and the Eulerian appraoch, an unified Eulerian and Lagrangian coordinate transformation was proposed by Hui et al [2, 3] to solve the Euler equations and achieve sharp of resolution of the contact line correctly.

As we know, in the frame work of unified coordinates approach, the fluid equations and geometric evolution equations are written in a combined conservative form, which is different from the fluid equations alone in the pure Eulerian approach. The hybrid type coordinate system considers the flow variables to be functions of time and of some permanent identification of pseudo-particles which move with velocity hq, q being the velocity of fluid particles. It includes the Eulerian coordinates as special case when h = 0 and the Lagrangian when h = 1. The unified coordinate system decides the grid velocity set to be hq , where q is the fluid velocity and h is a parameter which is determined by constraint conditions, such as the mesh alignment with the slip surface, or keeping grid angle during the mesh movement. Therefore, the grid velocity can be

changed locally according to the value of h. With a prescribed grid velocity, the inviscid flow equations are written in a conservative form in the computational domain (λ, ξ, η,ζ ), as well as the geometric conservation laws which control the mesh deformation. Therefore, numerical diffusion across the slip line can be reduced to a minimum with the crisp capturing of the contact discontinuity. T

Based on the Hui’s idea, we would like to extend the previous work [4, 5] to derive three- dimensional incompressible flow equations under the Euler- Lagrangian coordinate. In the framework of Euler-Lagrangian coordinates, the unsteady artificial compressibility based incompressible flow equations are derived and the related moving geometry equations can be achieved in conservation form and updated simultaneously during each time step. Thus, the accurate estimation of geometry conservation and controlling the grid velocity based on the unified approach are expected to avoid severe grid deformation and computation breakdown caused by moving body or boundary layers (considered as slip lines in Lagrangian coordinates).

Also, a unified artificial compressibility approach is developed to simulate the moving body flows with viscous effects and fluid-structure interaction. Test cases including the three- dimensional lid-driven cavity flow, a dropping airfoil, and a pulsating aortic tube are used to verify the computations. Under this circumstance, the current FORTRAN code is converted to the MPI code tested on a 64-CPU IBM SP2 parallel computer.

Governing Equation

The governing equations of the flow considered are the time-dependent incompressible Navier-Stokes equations (1) in the Cartesian coordinate. After introducing the pseudo- compressibility to connect pressure with continuity equation based on Chorin [1], the considered equations of motion of the fluid can be compactly written in in the following nondimensional conservation form. like

[ ]

1 Re

Q E F T

N Q

t x y z

∂ +∂ +∂ +∂ = ∇

∂ ∂ ∂ ∂ ⁽¹⁾

Where

2

2

2

0 0 0 0 0 1 0 0

, , , ,

0 0 1 0 0 0 0 1

p u v w

u u p uv uw

Q E F T N

v uv v p vw

w uw vw w p

β β β

⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞

⎜ ⎟ ⎜ + ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟

= = = = =

⎜ ⎟ ⎜ ⎟ ⎜ + ⎟ ⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟

⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ + ⎠ ⎝ ⎠

Based on the transformation of unified coordinate dt d

dx hud Ad Ld Pd

dy hvd Bd Md Q d

dz hwd Cd Nd Rd

λ

λ ξ η ζ

λ ξ η ζ

λ ξ η ζ

⎧ =

⎪ = + + +

⎪⎨ = + + + ′

⎪⎪ = + + +

⎩

(2)

where h is an arbitrary function of coordinates (λ, ξ, η, ζ ) and u, v, w are x-, y- and z- component of the fluid velocity q , Re is the so-called Reynolds number, β is pseudo compressibility factor; respectively. Let

D^h

hu hv hw

Dt t x y z

∂ ∂ ∂ ∂

≡ + + +

∂ ∂ ∂ ∂ (3)

denotes the material following the pseudo-particle, whose velocity is

hq

^Dh _0, ^Dh _0, ^Dh ₀

Dt Dt Dt

ξ = η = ζ = (4)

That is, the curvilinear coordinate are material functions of the pseudo-particles, and hence are their permanent identifications. Accordingly, computational cells move and deform with pseudo- particles, rather than with fluid particles as in Lagrangian coordinates. Furthermore, the geometrical state variables satisfy the compatibility conditions as

, , , , ,

A hu L hu P hu B hv M hv Q hv C hw N hw R hw

, , ,

λ ξ λ η λ ζ λ ξ λ η λ ζ λ ξ λ η λ ζ

∂ =∂ ∂ =∂ ∂ =∂ ∂ =∂ ∂ =∂ ∂ ′=∂ ∂ =∂ ∂ =∂ ∂ =∂

∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ (5)

From the above transformation matrix, the governing equation becomes

v v v

E F G

Q E F G

λ ξ η ζ ξ η ζ

∂ ∂ ∂

∂ +∂ +∂ +∂ = + +

∂ ∂ ∂ ∂ ∂ ∂ ∂ (6)

Where _Qis a preserved variable vector,E, F and G are flux vectors and _Ev, _FvG_v are viscous terms. Like

( )

( )

( )

( )

( )

( )

( )

( )

1 1

1 1

1 1

0 0

, ,

0 0 0 0 0 0

x x

y y

z z

p q q

u uq h p uq h p

v vq h p vq h p

w wq h p wq h p

A hu

B hv

Q C E hw F

L M N P Q R

β ξ β η

ξ ξ η η

ξ ξ η η

ξ ξ η η

⎛ ∇ ⎞ ∇

⎛ ⎞ ⎜ ⎟

⎜ ⎟ ⎜ ∇ − + ⎟ ∇ − +

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ∇ − + ∇ − +

⎜ ⎟

⎜ ⎟ ⎜ ∇ − + ⎟ ∇ − +

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

=⎜ ⎟⎜ ⎟ =⎜ ⎟ =

⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟′

⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎝ ⎠ ⎝ ⎠

i i

i i

i i

i i

( )

( )

( )

( )

1 1 1

0 0

0 , 0

0 0 0 0

0 0

x y z

q

uq h p

vq h p

wq h p

G hu

hv hw

hu hv hw β ζ

ζ ζ

ζ ζ

ζ ζ

⎛ ⎞ ⎛ ∇ ⎞

⎜ ⎟ ⎜ ∇ − + ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ∇ − + ⎟

⎜ ⎟ ⎜ ⎟

∇ − +

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ =⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎝ ⎠ ⎝ ⎠

i i i i

Wall model

In this work, we consider a deformable aorta in the numerical tests. The fluid-structure interaction is required to consider in the simulation. Here, the wall compliance is modelled using an independent ring model to compute the vessel deformations. This model assumes that the structural nodes move only in the radial direction. In spite of its intrinsic limitations, the extreme simplicity of this model makes it very popular. A linear elastic model equation to describe the wall motion can be written as a damped oscillator:

m ²₂r d r kr P_w

t t

∂ + ∂ + =

∂ ∂ (7)

Where , _{2 2}, 2

(1 )

w

m h k Eh d mk

ρ a

= = ν =

− and h is the wall thickness, ^ρw the wall density, E the Young’s modulus, v the Poisson ratio, a the vessel radius, p the pressure, r the wall displacement and _pw the pressure force at the wall. The radial displacement of each structural node can be obtained by solving equation (7) by using a fourth order Runge-Kutta scheme. The fluid-structure equations (6) & (7) are solved in an uncoupled way. Both the solutions of fluid and structure equations are updated in an unsteady time marching manner. The pressure loads at the vessel wall predicted by the fluid solver are transferred to be the source terms in the structure equation at the same time step. After that, the wall displacement is updated at each grid point

along the whole tube and vessel. Also, the wall mesh velocity can be estimated based on the last two wall displacements during the previous two time steps to achieve the second-order accurate estimation. Then, the new wall position and mesh velocity are substituted into the fluid solver and the related boundary conditions. Therefore, the updated pressure load can be re-predicted and a cycle of simulation fluid-structure interaction during the same time step is completed until the fluid-structure interaction is repeated until mass conservation criteria is satisfied in the fluid solver. However, the strategy of wall motion estimation may result in the moving grid distortion to produce excessive numerical errors. To avoid numerical instability, the geometry conservation included in equation (6) and a grid re-generation may be the necessary procedure.

Numerical Results

In the numerical flux approximation, one three-order accurate upwinding flux extrapolation for the derivatives ∂_Eξ which have the form. A third-order upwind flux at the cell interface is defined by

1 1 1 1 1

2 1 2 2 2 2

1 1

[ ( ) ( )] [ ] 2 ^{j j} 6

i i i i i

E₊ = EQ₊ +EQ + Δ −Δ +Δ −ΔE⁺₋ E⁺₊ E⁻₊ E⁻₊ (8) Where The flux difference is taken as

1 1

2 ( ) 2

i i

E^±₊ A Q Q^± ₊

Δ = Δ (9) where A represents the Jacobian matrix^. The A⁺and A⁻matrices are computed first by forming a diagonal matrix of the positive eigenvalues and multiplying through by the similarity transform, and since the A⁺matrix plus the A⁻ matrix equals the original Jacobian matrix, we have

1 +

−

− +

−

= Λ

= A A A

X X

A (10)

where X is the matrix of right eigenvectors of A ,X is its inverese. The flux difference is ⁻¹ evaluated at the midpoint by using the average of Q:

_{1 ( 1 1)}

2 ^{j j}

Q= Q₊ +Q₊ (11)

The eigenvectors in (11) of equation (6) as shown in [4] linearly independent, forming a

completed basis in the state space. The system equations (8) are therefore regarded as hyperbolic for all values of h. This includes the Eulerian coordinates as a special case when h = 0 and the Lagrangian one when h = 1. As we have known from Hui’s works [2,3], it is shown that for the smooth solutions of the system of two-dimensional Euler equations of gas dynamics written in the classical Lagrangian coordinates is equivalent to the same system written in the unified coordinates with h = 1. The steps of the proof can easily be repeated for the current artificial compressibility flow equations to show its weak hyperbolicity existed on the classical Lagrangian coordinates. To avoid possible numerical difficulties arising from the lack of a completed set of eigenvectors, we shall use 0 ≦ h < 1 to keep the hyperbolic character of the artificial compressibility equations. In the following test cases, we choose h =0.99 for all the steady-state cases to keep the weak hyperbolicity for the model equations we use. In the unsteady cases, h is assumed as grid velocity for the moving body simulation to keep the geometry conservation law. Also, β is kept to be 10 in all the cases.

Parallel implementation and Discussions

The current parallelized flow solver presents the results of two-dimensional and three- dimensional numerical simulation of the flow field in the followings:

First, the lid-driven cavity flow problem is a widely used benchmark test for the incompressible Navier-Stokes code validation. With the simplicity of geometry, the driven cavity flow contains complicated flow physics driven by multiple counter rotating vortices on the corners of the cavity depending on the Reynolds number. From all our computations for Re=3200 The computations are performed on a 64x64x64 on non-uniform grids as shown in Figure 1(a) which are clustered near the wall and stretched from the wall to the cavity center.

The velocity plot on a plane and 3D streamtrace patterns in the driven cavity flow are plotted in Figure 1. It is shown that one primary vortex near the center and three corner eddies are captured.

Also one small secondary zone in the lower right corner is visible. Secondary, a two- dimensional experiment of a falling airfoil in a water tank conducted by Andersen, Persavento, and Wang [5], for the falling airfoil, the freely falling trajectory was measured based on the flow conditions with the Reynolds numbers 1100. The meshes of 400×100 are used for the calculations. Fig. 2(a) presents trajectories of the falling airfoil, where the black one is the experimental measurement in [5] and the purple one is the trajectory from the computation.

Overall, the two trajectories are close to each other and have an identical slope. The vorticity field at four instants during a full rotation of the airfoil is also presented in Fig. 2(b). Thirdly, to perform numerical simulation on an aortic tube under fluid-structured interaction, one cycle of heartbeats is 0.855 seconds according to MRI data. The equation (6) was solved for a peak Reynolds number of 5000 at the inlet of ascending of aorta and numerical boundaries were chosen based on flow conditions: (i) MRA scanned flow rate [6] at the inlet of ascending aorta (ii) Surface traction free and zero velocity gradients at the outlet of descending aorta. (iii) MRA scanned flow rates at the outlet of three branches (iv) Grid velocity as the vessel wall velocity condition. Then, the final results were achieved at the fifth cycle of the computation which was starting from the initial conditions as zero velocity. Velocity vectors and shear stress distributions on vessel walls on the cross sections of the aorta are shown in our computations as in Figure 3 and 4. The computed averaged shear stresses along aortic wall with and without elastic assumptions are observed. It is noted that a computed peak value of the wall shear stress along the aortic wall at the aortic arch and the wall shear stress drop at downstream of the aortic arch during t =2/4 T. These phenomena may be resulted from the variation of the vessel diameter and the presence of the bifurcation. The inlet flow rate approaches zero in the late diastole, so the wall shear stress distributions are approaching flat. It demonstrates that wall shear stresses were highly dynamic, but were generally high along the vicinity of the branches and low along the lesser curvature, particularly in the descending thoracic aorta. The maximum wall shear stress distribution is presented on the aortic arch in the late systole.

Elastic stenotic tube

Figure 9 demonstrates that a completed cycle of the velocity contours in an elastic stenotic tube. The computation is performed based on a 200x51 grid cells. The non-dimensional geometrical parameters are as follows: inlet diameter 1; pre-stenosis length 5, and stenosis length 1; a long post-stenosis domain 31 are chosen in order to minimize the influence of downstream boundary conditions. In order to avoid considering the turbulence effects, a mild stenosis with only 25% area reduction is considered. The flow was assumed to be laminar, incompressible and Newtonian, and the walls only deform in the radial direction with the grid velocity obtained from equation (7). Therefore the h value is floating with grid velocity at each grid

point during every physical time step. At the inlet boundary, a sinusodial incoming velocity profile as described in equation (12) is used. The fluid-structure interaction is included in the calculations. The pressures at the inlet can be obtained from the solution of the governing equation by assuming an prescribed incoming velocity distribution. At the outflow boundary, traction-free conditions is assumed.

( )

( ) 1 1 cos(2 )

u λ =2 − πλ (12)

In this case, a Reynolds number defined at the inlet of 800 is selected. As total velocity contour plots shown in Fig. 9, it is shown that the reverse flow distal to the stenosis and the re- circulation region moves to downstream from the early systolic cycle to late diastolic cycle.

Also we can find that the formation of the recirculation region is strongly effected by the compliant vessel. The separation mainly appears around the flow over the neck of the stenosis, and then it disappears around the region along the convergent wall and restarts at the beginning of the divergent wall. It is shown that the separation zone shows up periodically along the compliant vessel wall. In Figure 10, pressure distributions are shown to distribute along the stenotic wall. It is shown that the pressure drop occurs around the neck of the stenotic region in the late diastolic cycle and at the beginning of systolic cycle. Next, the pressure distributions along the stenotic region for four different Reynolds number with 100, 200, 400, and 600 are also shown in Figure 10. It is shown that the pressure drop occurs around the neck of the stenotic region no matter what the value of Reynolds number is. A negative pressure difference is increased with the increasing Reynolds number. Furthermore, the estimations of wall shear stress distributions for the whole wall and the stenotic area are depicted in Figure 10. It is very encouraging to find that the prediction of the location of maximum wall shear stress is consistent with the analytic studies which demonstrate the location of maximum wall shear stress is always upstream of the neck of the stenosis and independent of Reynolds number in the whole pusatile cycle. In addition, a strong oscillating wall shear stress distribution is shown to appear around the neck of the constriction-enlargement region for the lower wall which is depicted in Figure 11. It is noted that the prediction of the location of maximum wall shear stress is always found around the beginning of enlargement and independent of the Reynolds number. Also a positive peak of the shear stress distribution is found to appear on the lower

aortic arch during the systole and the early diastole, and then changed to be a negative peak during the late diastole.

The computation is very time consuming for time accurate and pseudo-time evolutions in the above calculations. The parallel computation technology is very necessary in the three- dimensional cases. In our parallel computation tests, the test results also display very promising potential of parallel processing as shown in Table. The original standard FORTRAN based incompressible Navier-Stokes code coupled with a linear elasticity model was converted to be a MPI based solver, also, it was tested on IBM SP2 690A parallel System. The parallel system is consisted of 415.2 GFLOPS with 96 CPUs Multiprocessor (SMP) nodes connected by High- performance switches, each node contains four POWER3 processors, four GB of main memory and 192 GB of hard disk. The clock rate of the processor is 1.9 GHz. The Floating Peak Performance of each processor is SPECcfp2000 of 1898. Each processor comprises eight execution units, a 32 KB instruction cache, 64 KB data cache and an on-board bus interface unit.

There are three fixed point units, two floating point units, two load/store units and one branch/dispatch unit in each processor. The MPI code is paralleled along both the longitude and latitude. The tested result of the model is shown in table. The model results were then carefully validated to ensure that the two versions of the model produce virtually the same result. We have made several test runs and the results are summarized in Table A. Apparently, the concept of parallel processing suited the current dual-time Navier Stokes Solver very well. The model can take advantage of the MPI code fully, since minimal amount of data transfer among CPUs is required for solving the governing equation explicitly. A significant reduction of computing time in running the model and a almost linear speed up rate is achieved up to 32 CPUs in all the different data partition. The speed up rate is as high as 32 for using 64 processors at the same time. It provides very promising potential for prompt diagnosis using modern CFD technology.

Concluded Remarks

z In this paper, an unsteady artificial compressibility solver for moving body simulation based on unified coordinate approach is proposed and developed. In the framework of unified coordinates, a unified moving body approach, include the unsteady

incompressible flow equations and moving geometry equations, with grid velocity as hq achieves conservation form and updated simultaneously during each time step.

Accurate estimation of geometry conservation and controlling the grid velocity based on the unified approach can avoid severe grid deformation and computation breakdown caused by moving body or boundary layers (considered as slip lines in Lagrangian coordinates). Also, a linear elastic modeling of wall motion is included here for the consideration of the fluid-structure interaction.

z In this study, a parallel incompressible flow-structured solver instrumented with MPI is developed on the simulation of bio-fluid flow problems. The test results show that a significant reduction of computing time in running the model and a near-linear speed up rate is achieved up to 32 CPUs at IBM SP2 platforms. The speed up rate is as high as 32 for using IBM SP2 64 processors Also, the overall accuracy on the pressure on blade surfaces are in good agreement with validated data. The test shows very promising potential of the current parallel flow code to provide prompt simulation of the current flow test cases.

Acknowledgments

We acknowledge financial support from the National Science Council under the project NSC 97-2221-E-216-009-MY2 and the computer facility supported by the National Center for High- Performance Computing, Hsin-Chu, Taiwan, ROC

Figure 1 A plot of 64 x 64 x 64 grid distribution and flow patterns for a 3D driven cavity flow (Re=3200)

(a) Numerical validation of the trajectory of a falling airfoil (the red solid line--- the measured data [6], the dashed line--- simulation)

(b) Vorticity plots at different time for the dropping airfoil during a full rotation

Figure 2 A simulation of a falling airfoil flow (Re=3200)

(a) 1/4 T 2/4 T

Figure 3 Blood flow velocity vector through the aortic vessel during a full cardic cycle (Re=5000)

(a) 1/4 T (b) 2/4 T

(c) 3/4 T (d) 4/4T

Figure 4 the shear stress distributions along the Aortic vessel during a full cardic cycle (Re=5000)

t=0/19 T

t=5/19T

t=10/19 T

t=15/19 T

Figure9 Velocity contour plots in a pulsatile cycle for a deformable tube

Figure 10 Variations of pressure at wall for Reynolds numbers 100 200 ： ； ；400 600 ； ：from X=4 to X=7 at time=3/4 T

Figure 11 Variations of shear stress at stenotic wall for Reynolds numbers=800 from X=2 to X=9 in pulsatile cycle

Table 2-D with data partition on IBM SP2

CPUs Elapsed time(sec) Speed up Efficiency 1

4 8 16 32 64

107438 35001 15803 10998 6043 3152

1.00 3.19 6.2 10.1 18.2 32.2

100.%

80.%

75.%

70.%

64%

51.%

REFERENCES

[1] Chorin A.J. “A Numerical Method for Solving Incompressible Navier-Stokes Equations”,

Journal of Computational Physics, 2 (1967) 12-26

[2] W. H. Hui, P. Y. Li, and Z. W. Li, ”A Unified Coordinate System for Solving the Two- Dimensional Euler Equations ,” J. Comput. Physics, 153,596-637, 1999

[3] W. H. Hui and S. Koudriakov . “A Unified Coordinate System for Solving the Three- Dimensional Euler Equations”, Journal of Computational Physics, 2001, 172: 235-260

[4] Yang-Yao Niu, Yi-Hsun Lin, W. H. Hui, Chien C. Chang. International Journal for Numerical Methods in Fluids 61(9):1029-1052 (2009)

[5] Yang-Yao Niu*, Chih-Hung Chang, Wen-Yih I. Tseng and Hsu-Hsia Pen , “A Transient HLLC type Incompressible Riemann Solver and Application on Aortic Flows Simulation “, Communications in Computational Physics, Vol. 5, No. 1, pp. 142-16, 2009

[6] A. Andersen, U. Persavento and Z. J. Wang, Unsteady aerodynamics of fluttering and tumbling plates, J. Fluid Mech., 541 (2005), 65-90

(ISCM II) 學術研討會心得報告

中華大學機械工程學系

報告人: 牛仰堯

本 次 參 加 的 會 議 為 The International Chinese Association for Computational Mechanics (ICACM) 每2年一度的盛大研討會，有來自各地的計 算力學專家與學者，包括大學教授、研究人員、各國家實驗室的專職人員、工

從事計算力學的研究人員，有從理論上去探討，有從實驗上去觀察，也有從 電腦上去模擬，以發掘其中的奧祕。近幾年來由於電腦的快速發展，使得研究 人員使用電腦的機會大為增加。探討理論者借助電腦來推衍方程式，實驗觀察 者仰賴電腦收集與分析大量的數據資料，電腦模擬者更是運用電腦來分析各種 問題。而計算力學在工程上更是一個日益受重視的領域，不僅因為電腦模擬技 術的精進，更由於電腦模擬快速有效的解決許多工程上的問題。

.

在十餘年前，電腦的使用尚未普及，一般工程的問題多仰賴實驗與理論 分析。不僅花費的時間冗長，費用龐大，更由於體能與實驗材料的限制，

許多複雜的問題不容易獲得理想的答案。但是在沒有更理想的方法之前，

對於未知的事項只能以相關的經驗去克服。直到 80 年代以後大型與快速的 電腦漸普及，電腦容量亦可應付各種研究的需求，因此計算模擬分析便成 為工程研究中的一個主要分支，在工程研討會中接近半數的研究報告是屬 於此一範疇，可見一斑。參加國際性的相關研討會，可說是最直接也是最 有效的方式。透過與國際上學有專精的專家學者交流並交換意見，以獲得 最新最完整的流體工程研究資訊，並了解發展趨勢，提供給國內研究人員，

作最佳的服務。尤其在今日隨電腦科技的突飛猛進，各種計算方法日新月 異，許多過去無法嘗試的技巧，都逐步進入實用階段，面對如此的環境，

端賴期刊、書籍已無法應付此種變遷，因此在民國 98 年 11 月 30 日至 8 月 3 日由本人出席了 2009 年 International Computational Mechanics (ICACM) 工程會議。本次會議我與國科會台法幽蘭多相流計畫成員(樹德科技大學胡 舉軍教授，台灣大學曾予恆教授，台灣大學薛克民教授，暨南國際大學戴 義欽教授，中央研究院郭志禹博士與法國尼斯大學 Boniface.NKONGA 教 授，南特大學 Christophe.Berthon 教授，合組一個論壇，就兩相流的數值 方法進行討論。

本次會議的地點在Hong Kong & Macau。會議於11月29日開始辦理報到註 冊，主辦單位同時展示一些會議相關論文集及專業書刊，並有部份工作人員舉 行會議前之準備集會。接下來的四天為正式的會議，包括專題演講、邀請演講、

專題討論及論文發表等議程。由於計算力學涵蓋範圍廣泛，參加人數眾多，大 會依專業細分為26大項，論文發表以20個場地同時進行，共計有200場，每場 約有5～6篇論文，會場外的大廳有軟體廠商、出版商廠商的展示。由於場次眾 多，個人僅能選擇性的參與部份論文發表。大會的子題包括了：

• CAD, CAM and CAE

• Adaptive Materials Systems, Structures and Smart Material

• Biomechanics

• Electromagnetism

• Engineering Sciences and Physics

• Environmental Science and Engineering

• Fluid Mechanics and Heat Transfer

• Geomechanics, Geographic Information Systems

• Microtribology and Micromechanics

• Methods in the Life Sciences

• Nanotechnology

• Nonlinear Dynamics

• Computer Simulation of Processes and Manufacturing

• Data and Signal Processing

• Meshless and Wavelet Methods

• Multiple-Scale Physics and Computation

• Parallel Computing

• Scientific Visualization

由這些副題可以看出，民生工業相關的力學研究，例如電腦繪圖工程、流體 機械、環境或工業應用等，這也是目前國際間研究重點並積極爭取民生工業之 計畫。也就可以解釋研討會中論文發表的方向。 本次研討會主辦單位學會是 香港城市大學與澳門大學。研討會範圍深入工程的各個層面，由於許多場次是 同時進行，因此只能選擇性的參與。對於未能參加之場次則蒐集所發表之論文 作為日後參考之用。在會場上，有機會與各國的專家學者接觸、交流，瞭解到 國際上在計算工程的發展現況、研究方向、及未來展望，也讓他們了解到在台 灣的電腦計算在計算工程研究方面的現況。 此行可謂成果豐碩，確實達到開 會之任務需求。

.

量化

成果項目 ^{實際已達成}

數（被接受 或已發表）

預期總達成 數(含實際已

達成數)

本計畫實 際貢獻百

分比

單位

備 註 （ 質 化 說 明：如 數 個 計 畫 共 同 成 果、成 果 列 為 該 期 刊 之 封 面 故 事 ...

等）

期刊論文 0 0 90%

研究報告/技術報告 0 0 100%

研討會論文 4 5 100%

論文著作 篇

專書 0 0 100%

申請中件數 0 0 100%

專利 已獲得件數 0 0 100% 件

件數 0 0 100% 件

技術移轉

權利金 0 0 100% 千元

碩士生 0 0 100%

博士生 0 0 100%

博士後研究員 0 0 100%

國內

參與計畫人力

（本國籍）

專任助理 0 0 100%

人次

期刊論文 2 4 100%

研究報告/技術報告 0 0 100%

研討會論文 3 4 100%

論文著作 篇

專書 0 0 100% 章/本 申請中件數 0 0 100%

專利 已獲得件數 0 0 100% 件

件數 0 0 100% 件

技術移轉

權利金 0 0 100% 千元

碩士生 0 0 100%

博士生 0 0 100%

博士後研究員 0 0 100%

國外

參與計畫人力

（外國籍）

專任助理 0 0 100%

人次

得獎項、重要國際合 作、研究成果國際影響 力及其他協助產業技 術發展之具體效益事 項等，請以文字敘述填 列。)

成果項目 量化 名稱或內容性質簡述

測驗工具(含質性與量性) 0

課程/模組 0

電腦及網路系統或工具 0

教材 0

舉辦之活動/競賽 0

研討會/工作坊 0

電子報、網站 0

科 教 處 計 畫 加 填 項

目 計畫成果推廣之參與（閱聽）人數 0

值（簡要敘述成果所代表之意義、價值、影響或進一步發展之可能性） 、是否適 合在學術期刊發表或申請專利、主要發現或其他有關價值等，作一綜合評估。

1. 請就研究內容與原計畫相符程度、達成預期目標情況作一綜合評估

■達成目標

□未達成目標（請說明，以 100 字為限）

□實驗失敗

□因故實驗中斷

□其他原因 說明：

2. 研究成果在學術期刊發表或申請專利等情形：

論文：■已發表 □未發表之文稿 □撰寫中 □無 專利：□已獲得 □申請中 ■無

技轉：□已技轉 □洽談中 ■無 其他：（以 100 字為限）

1.Yang-Yao Niu, Hsieh,C.T. , Chang,C.C., and Chu, C.C, ＇How does a Gurney flap enhance the aerodynamics force＇, AIAA Journal, in press, 2010

3. 請依學術成就、技術創新、社會影響等方面，評估研究成果之學術或應用價 值（簡要敘述成果所代表之意義、價值、影響或進一步發展之可能性）（以 500 字為限）

在統一座標中發展血液動力流體之流固耦合模式，使用 MRA 影像技術完成人體血管之壁面 移動，以考慮血液與血管壁之交互作用，血管模擬流固耦合及相關的血液流體行為及壓 力，剪力形式，分析結果與影像處理掃描所獲得結果互相驗證、更精確模擬出血液流動之 行為，可以作為提供安全血壓值及醫生臨床主動脈手術及用藥之參考範圍。