Plenary Lectures Merge B5

Plenary Lectures - Abstracts

Residual-driven Variational Multiscale Turbulence Modeling for Large Eddy Simulation of Incompressible Flow Y. Bazilevs, V.M. Calo, J.A. Cottrell, T.J.R. Hughes Institute for Computational Engineering and Sciences, The University of Texas at Austin, USA e-mail: {bazily,victor,jac3,hughes}@ices.utexas.edu

A. Reali Structural Mechanics Department, University of Pavia, Italy e-mail: [email protected]

G. Scovazzi 1431 Computational Shock and Multi-physics Department, Sandia National Laboratories, USA e-mail: [email protected]

Abstract The objectives of recent variational multiscale work in turbulence have been to capture all scales consistently and to avoid use of eddy viscosities altogether. This holds the promise of more accurate and efficient LES procedures. In this work, we describe a new variational multiscale formulation, which makes considerable progress toward these goals [1].

Keywords: variational multiscale methods, large eddy simulation, turbulence modeling, Isogeometric Analysis, NURBS, incompressible flows, homogeneous isotropic turbulence, turbulent channel flows

1

Summary

We begin by taking the view that the decomposition into coarse and fine scales is exact. For example, in the spectral case, the coarse-scale space consists of all Fourier modes beneath some cut-off wave number and the fine-scale space consists of all remaining Fourier modes. Consequently, the coarse-scale space has finite dimension whereas the fine-scale space is infinite dimensional. The derivation of the coarse- and fine-scale equations proceeds, first, by substituting the split of the exact solution into coarse and fine scales into the Navier-Stokes equations,

then, second, by projecting this equation into the coarse- and fine-scale subspaces. The projection into coarse scales results in a finite dimensional system for the coarse-scale component of the solution, which depends parametrically on the fine-scale component. In the spectral case, in addition to the usual terms involving the coarse-scale component, only the cross-stress and Reynolds-stress terms involve the fine-scale component. In the case of non-orthogonal bases, even the linear terms give rise to coupling between coarse and fine scales. The coarse-scale component plays an analogous role to the filtered field in the classical approach, but has the advantage of avoiding all problems associated with homogeneity, commutativity, walls, compressibility, etc. The projection into fine scales results in an infinite-dimensional system for the fine-scale component of the solution, which depends parametrically on the coarse-scale component. We also assume the cut-off wave number is sufficiently large that the philosophy of LES is appropriate. For example, if there is a well-defined inertial sub-range, then we assume the cut-off wave number resides somewhere within it. This assumption enables us to further assume that the energy content in the fine scales is small compared with the coarse scales. This turns out to be important in our efforts to analytically represent the solution of the fine-scale equations. The strategy is to obtain approximate analytical expressions for the fine scales then substitute them into the coarse-scale equations which are, in turn, solved numerically. If the scale decomposition is performed in space and time, the only approximation in the procedure is the representation of the fine-scale solution. To provide a framework for the fine-scale approximation, we assume an infinite perturbation series expansion to treat the fine-scale nonlinear term in the fine-scale equation. By virtue of the smallness of the fine scales, this expansion is expected to converge rapidly under the circumstances described in many cases of practical interest. The remaining part of the fine-scale Navier-Stokes system is the linearized operator which is formally inverted through the use of a matrix Green’s function. The combination of a perturbation series and Green’s function provides an exact formal solution of the fine-scale Navier-Stokes equations. The driving force in these equations is the Navier-Stokes system residual computed from the coarse scales. This expresses the intuitively obvious fact that if the coarse scales constitute a good approximation to the solution of the problem, the coarse-scale residual will be small and the resulting fine-scale solution will be small as well. This is the case we have in mind and it provides a rational basis for assuming the perturbation series converges rapidly. Note that one cannot use such an argument on the original problem because in this case the perturbation series would almost definitely fail to converge. (If we could have used this argument, we would have solved the Navier-Stokes equations analytically! Unfortunately, this is not the case.) The formal solution of the fine-scale equations suggests various approximations may be employed in practical problem solving. We are tempted to use the word “modeling” because approximate analytical representations of the fine scales constitute the only approximation and hence may be thought of as the “modeling” component of the present approach, but we want to emphasize that this is very different from classical modeling ideas which are dominated by the addition of ad hoc eddy viscosities. We will present numerical results that demonstrate

that eddy-viscosity terms are unnecessary in the present circumstances. There are two aspects to the approximation of the fine scales: 1) Approximation of the matrix Green’s function for the linearized Navier-Stokes system; and 2) approximation of the nonlinearities represented by the perturbation series. The first and obvious thought for the latter aspect, nonlinearity, is to simply truncate the perturbation series. This idea is pursued in conjunction with some simple approximations of the Green’s function. It turns out there is considerable experience in local scaling approximations of the Green’s function based on the theory of stabilized methods; Hughes [4], Hughes et al. [5], Hughes and Sangalli [6], Hughes, Scovazzi and Franca [7]. The Green’s function is typically approximated by locally defined algebraic operators (i.e., the “τ ’s” of stabilized methods) multiplied by local values of the coarse-scale residual. An outline of the presentation is summarized as follows: we begin by presenting the mathematical details of the variational multiscale theory. This represents our general approach to LES-style turbulence modeling and is independent of the specifics of the discrete spaces utilized to represent the coarse scales. The relationship between this version of the variational multiscale method and classical stabilized methods is delineated. It is noted that that the variational multiscale method includes additional terms. Both conceptually and from the point of view of actual implementation, stabilized methods may be viewed as historical stepping stones leading to the more coherent variational multiscale formulation. We then present our numerical studies of forced isotropic turbulence at Reλ = 165 and Reλ = ∞. (Reλ is the Taylor microscale Reynolds number.) We begin with a description of the approximation spaces consisting of NURBS elements (non-uniform rational B-splines, see, e.g., Rogers [14], Piegl and Tiller [13], Farin [3], and Cohen, Riesenfeld and Elber [2]). In the case of the rectilinear geometry considered, NURBS reduce to B-splines, which have been advocated for turbulence calculations previously (see Kravchenko, Moin and Moser [8], Shariff and Moser [15], Kravchenko, Moin and Shariff [9], and Kwok, Moser and Jim´enez [10]). We employ trivariate linear, quadratic, and cubic NURBS with periodic boundary conditions. Linear trivariate NURBS turn out to be identical to trilinear hexahedral finite elements, but the higher-order NURBS are different than classical higher-order finite elements. We perform a dispersion error analysis for NURBS versus classical finite elements on simple, linear, one-dimensional advective and diffusive model problems, and conclude that NURBS have better approximation properties than classical finite elements. We employ meshes of 323 , 643 , 1283 , and 2563 to explore convergence with mesh refinement (h-convergence). We also examine the behavior of increasing order from linear to cubic on fixed meshes (k-convergence). In the case of Reλ = 165, we compare with the DNS spectral results of Langford and Moser [11]. Energy spectra and third-order structure functions are presented. Sample energy spectra results are presented in Figure 1. In the case of Reλ = ∞ we also clearly see the development of an inertial subrange. We present results for turbulent channel flows at Reτ = 395. (Reτ is the wall-friction Reynolds number.) We employ meshes of 323 and 643 . This time the mesh is graded in the wall-normal direction to better capture the boundary layer. Again, we consider convergence from the h- and k-refinement perspectives. A

striking result is how much better quadratic elements are than linear elements. For a mesh of 643 , the quadratic and cubic results are essentially identical to the DNS results of Moser, Kim and Mansour [12] for first- and second-order statistics (see Figure 2), and for a mesh of 323 they are in close agreement. We close with conclusions and suggested future directions for research.

References [1] Y. Bazilevs, V.M. Calo, J.A. Cottrell, T.J.R. Hughes, A. Reali, and G. Scovazzi. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Computer Methods in Applied Mechanics and Engineering, 197:173–201, 2007. [2] E. Cohen, R. Riesenfeld, and G. Elber. Geometric Modeling with Splines. An Introduction. A. K. Peters Ltd., Wellesley, Massachusetts, 2001. [3] G.E. Farin. NURBS Curves and Surfaces: From Projective Geometry to Practical Use. A. K. Peters, Ltd., Natick, MA, 1995. [4] T. J. R. Hughes. Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Computer Methods in Applied Mechanics and Engineering, 127:387–401, 1995. [5] T. J. R. Hughes, G. Feij´oo., L. Mazzei, and J. B. Quincy. The variational multiscale method–A paradigm for computational mechanics. Computer Methods in Applied Mechanics and Engineering, 166:3–24, 1998. [6] T. J. R. Hughes and G. Sangalli. Variational multiscale analysis: the fine-scale Green’s function, projection, optimization, localization, and stabilized methods. SIAM Journal of Numerical Analysis, 45:539–557, 2007. [7] T. J. R. Hughes, G. Scovazzi, and L. P. Franca. Multiscale and stabilized methods. In E. Stein, R. De Borst, and T. J. R. Hughes, editors, Encyclopedia of Computational Mechanics, Vol. 3, Computational Fluid Dynamics, chapter 2. Wiley, 2004. [8] A. G. Kravchenko, P. Moin, and R. Moser. Zonal embedded grids for numerical simulation of wall-bounded turbulent flows. Journal of Computational Physics, 127:412–423, 1996. [9] A. G. Kravchenko, P. Moin, and K. Shariff. B-spline method and zonal grids for simulation of complex turbulent flows. Journal of Computational Physics, 151:757–789, 1999.

[10] W. Y. Kwok, R. D. Moser, and J. Jim´enez. A critical evaluation of the resolution properties of B-spline and compact finite difference methods. Journal of Computational Physics, 174:510–551, 2001. [11] J. A. Langford and R. D. Moser. Optimal LES formulations for isotropic turbulence. Journal of Fluid Mechanics, 398:321–346, 1999. [12] R. Moser, J. Kim, and R. Mansour. DNS of turbulent channel flow up to Re=590. Physics of Fluids, 11:943–945, 1999. [13] L. Piegl and W. Tiller. The NURBS Book (Monographs in Visual Communication), 2nd ed. Springer-Verlag, New York, 1997. [14] D. F. Rogers. An Introduction to NURBS With Historical Perspective. Academic Press, San Diego, CA, 2001. [15] K. Shariff and R. D. Moser. Two-dimensional mesh embedding for B-spline methods. Journal of Computational Physics, 145:471–488, 1998.

10

k −5/3 1

E

DNS 32 0.1

64 128

0.01

256 1

10

k (a) C 0 -continuous linear NURBS 10

k −5/3 1

E

DNS 32 0.1

64 128

0.01

1

10

k (b) C 1 -continuous quadratic NURBS 10

k −5/3 1

E

DNS 32 0.1

64 0.01

1

10

k (c) C 2 -continuous cubic NURBS Figure 1: Energy spectra for h−refinement. Reλ = 165.

25

64

U+

20

P=2

3

15

P=1 DNS

10

P=3

5

0 0.1

1

10

100

y+ (a) Mean stream-wise velocity

w+

1.5 1 0.5

DNS

P=1

P=2

v+

1.5

P=3

1 0.5 0 4

64 3

P=1

3.5

u+

3 2.5 2

P=3

1.5 1

DNS

0.5 0

P=2 100

200 +

y (b) Velocity fluctuations

300

400

Figure 2: Turbulent channel flow at Reτ = 395 computed on a mesh of 643 elements: krefinement interpretation of results.

Geometric Flow for Quality Surface/Volumetric Modeling (Extended Abstract) Chandrajit L. Bajaj ∗ Department of Computer Science, Center for Computational Visualization, Institute for Computational Engineering and Sciences University of Texas, Austin, TX 78712 Email: [email protected]

Guoliang Xu†, Qin Zhang LSEC, Institute of Computational Mathematics, Academy of Mathematics and System Sciences, Chinese Academy of Sciences, Beijing, 100080 Email: {xuguo, zqyork}@lsec.cc.ac.cn

Abstract We present a general variational framework for a higher-order spline level-set (HLS) method and apply this to smooth surface constructions. Starting from a first order energy functional, we derive the general level set formulation, and provide an efficient solution of a second order , time-dependent, geometric partial differential equation (termed geometric flow), using a C 2 B-spline basis. We also present a fast cubic C 2 B-spline interpolation algorithm based on convolution and the Z-transform, which exploits the local relationship of interpolatory cubic spline coefficients with respect to given function data values. We provide two demonstrative smooth surface construction examples of our HLS method. The first is the construction of a smooth surface model (an implicit solvation interface) of bio-molecules in solvent, given their individual atomic coordinates and solvated radii. The second is the smooth surface reconstruction from a cloud of points generated from a 3D surface scanner.

1

Introduction

Level set methods are increasingly being used in the solution of time-dependent (evolutionary) partial differential equations [6]. Here we consider its application to smooth surface construction ∗ Supported in part by NSF grant CNS-0540033 and NIH contracts P20-RR020647, R01-GM074258, R01GM07308,R01-EB004873. † Supported in part by NSFC grant 10371130 and National Key Basic Research Project of China (2004CB318000).

1

(a)

(b)

(c)

C2

cubic B-spline Molecular Surface Reconstruction[2, 10, 4]: (a) shows Fig 1.1: Variational the van der Waals surface of a molecule. (b) shows the corresponding solvent excluded molecular surface constructed using our C 2 tri-cubic spline level-set method. (c) illustrates that the smooth solvent excluded surface constructed tightly encloses the van der Waals surface (a). by presenting a higher order spline level set generalization and solution of an appropriate geometric evolutionary partial differential equation. The level set geometric formulation is derived by minimizing an energy functional defined with respect to the surface and its first derivatives. Given a non-negative function g(x) over a domain Ω ∈ R3 . Find a surface Γ in Ω, such that the energy functional   E(Γ) = g(x)dx +  h(x, n)dx (1.1) Γ

Γ

is minimal, where x and n is a surface point and its surface normal, respectively. Further, h(x, n) is another given non-negative function defined over R3 × R3 which is used for regularizing the constructed smooth surface. Finally,  ≥ 0 is a constant. Many solid and physical modeling problems, such as surface (solid boundary) reconstruction, and physically based simulation of deformable interfaces could be formulated as minimizing an energy in the form of (1.1). By minimizing the energy functional (1.1), a partial differential equation (PDE) in the level-set formulation can be generated. The PDE is solved using the higher-order spline level-set (HLS) method that we present here. Fig. 1.1 and 1.2 show two examples of smooth surface construction, an interface (surface) that separates a molecule’s atoms from the solvent atoms (typically water), and a smooth surface reconstruction fit to a point set generated by a 3D surface scanner. Why Use a Level-set Method? In shape deformation simulations, topology changes may occur. This topology change makes parametric form surface tracking difficult. However, implicit form surface deformation could overcome this difficulty. Implicit surface splines, such as tetrahedral A-patches and prism A-patches, have been successfully used in computer graphics and surface modeling in the past decades (see [1, 5, 12] for references), however mostly used in static surface modeling. The level-set method described here allows one to dynamically deform and track an implicit surface spline using a governing PDE, which describes various laws of motion depending on geometry, external forces, or a desired energy minimization (see [6, 7, 9, 11]). Fur2

(a)

(b)

(c)

Fig 1.2: Variational C 2 cubic B-spline Point - Sampled Surface Reconstruction[]: (a) shows the point set data produced by a 3D surface scanner. (b) shows the the reconstruction result using our C 2 tri-cubic spline level-set method. (c) shows the reconstruction result for a C 0 tri-linear level-set method (also generated from our implementation). thermore, the underlining data structure is simple and topological changes are handled easily, with the computation being restricted to a thin shell (traditionally called a narrow band for evolution) surrounding the level-set. Why Use a Higher-order Spline Method? The level-set surfaces obtained from classical level-set methods are generally bumpy due to the use of piecewise tri-linear interpolation from the discrete function data computed on a rectilinear grid. To produce a better quality surface, a denser grid needs to be used. However, the increased grid resolution substantially increases the computation costs. Another drawback of using discrete data over grids is the non-trivial requirement of estimating derivatives for smooth interpolation. In many surface construction problems, such as the construction of molecular surface, the underlining surface is at least C 1 smooth. Therefore, a smooth level-set function is highly desirable. In this paper, we are solving second-order geometric partial differential equations. In the solutions of these PDE’s we utilize accurate estimates of mean surface curvature Therefore, we use C 2 tri-cubic spline as the level-set function basis. Note that tri-cubic is the lowest order B-spline that could achieve C 2 continuity in 3D. The advantages using C 2 spline function bases include: 1. Since the level-set function is C 2 smooth, the level-set surface is G2 smooth. There do exist a finite number of critical level-set values where the level-set may have a finite number of isolated singular points (i.e the gradient of the level-set function vanishes). However working in a finite precision numerical domain one automatically avoids these finite set of critical level-set values. 2. Derivatives up to the second order and curvatures, which appear in the governing geometric partial differential equations, are easily and exactly computed from the C 2 level-set function. 3. Using smooth level-set functions implies that larger and higher-order spline iso-surface 3

patches could be directly generated.

2

Mathematical Notation and Definitions

We first introduce some useful notation and definitions for geometric quantities on level-sets M in terms of the corresponding level-set function φ. Let φ : Ω → R be some smooth function on a domain Ω ⊂ R3 . Suppose Mc := {x ∈ Ω : φ(x) = c} is a level-set of φ for the level value c. For the sake of simplicity,we simply write M = Mc and assume that ∇φ = 0 on M. Hence by the implicit function theorem, Mc is a smooth surface and the normal n=

∇φ ∇φ

on the tangent space Tx M is defined for every x on M. Using a co-area formula, the energy functional can be defined as   E[φ] := e[Mc ]dc = ∇φh(x, n)dx. R

Ω

Next we compute  d

E [φ], ϑ = E[φ + ϑ]=0 d     ∇(φ + ϑ) d  ∇(φ + ϑ)h x, = dx  d ∇(φ + ϑ) =0  Ω  d = (∇(φ + ϑ))=0 h(x, r)dx d Ω  ⎞  ⎛ ∇(φ+ϑ)   d ∇(φ+ϑ) dx ⎠  dx + + hn ∇φ ⎝hx  d d

(2.1)

(2.2)

=0

Ω

The remaining task is to calculate  d ∇φ · ∇ϑ (∇(φ + ϑ))=0 = , d ∇φ  d

∇(φ+ϑ) ∇(φ+ϑ)

d

    

= =0

where operator P = I − identity mapping.

(2.3)

∇ϑ∇φ − (∇φ ⊗ ∇φ)∇ϑ/∇φ ∇φ2

= ∇φ−1 P∇ϑ,

∇φ ∇φ

⊗

∇φ ∇φ

(2.4)

is the projection onto the tangent space and I indicates the

4

Next we compute x (0) =

dx() d |=0 .

Since

φ(x()) + ϑ(x()) = c, Taking differential with respect to , and then taking  to be zero, we have (∇φ)T x (0) + ϑ(x) = 0.

(2.5)

Let x (0) = α(x)

∇φ + β(x)T (x), ∇φ

(2.6)

whereT (x) is a tangent vector at x. substituting (2.6) into (2.5), we obtain α(x) = −

ϑ(x) . ∇φ

Since the motion in the tangential direction does not alter the shape of the surface, we ignore the tangential movement and take x (0) = −ϑ(x)

∇φ ∇φ2

(2.7)

Substituting (2.3), (2.4) and (2.7) into (2.2), we obtain  ∇φ · ∇ϑ (2.2) = h(x, n) ∇φ Ω   (∇h)T ∇φ ϑ (∇n h)T P∇ϑ) + ∇φ − + dx ∇φ2 ∇φ

(2.8)

Equation (2.8) is the weak formulation of Euler-Lagrange equation. If we utilize the formula of integrate by parts taking into account ϑ ∈ C∞ 0 (Ω), Euler-Lagrange equation   T ∇φ (∇h) ∇φ −div h(x, n) + P∇n h − =0 (2.9) ∇φ ∇φ is deduced. Obviously, this equation is a second-order nonlinear partial differential equation. If we write the left hand side of (2.9) as an operator acted on function φ, we can construct geometric (gradient) flow as follows: ∂t φ = −∇φL(φ). For further details, please see [3].

5

3

Algorithm Outline

From minimizing the energy (1.1), we obtain the following evolution equation (see [3] for details). ∂φ ∂t

 = (g + h)div

∇φ ∇φ

 ∇φ +  div(P∇n h)∇φ

+ 2[∇(g + h)]T ∇φ = L(φ) + H(∇φ), where

 L(φ) = (g + h)div

∇φ ∇φ

(3.1)

 ∇φ,

H(∇φ) =  div(P∇n h)∇φ + 2[∇(g + h)]T ∇φ, ∇φ ∇φ P = I − ∇φ ⊗ ∇φ is a projection operator onto the tangent space and I indicates the identity mapping. ∇ and ∇n denote the usual gradient operator with respect to x and n. Note that L(φ) is a parabolic term and H(∇φ) is a hyperbolic term. Hence, in solving equation (3.1) in the following, the first order term H(∇φ) is computed using an upwind scheme (see [8] for the reason of using an upwind scheme) over a finer grid, the higher order term L(φ) is computed using a spline presentation defined on a coarser grid. Consider the solution of equation (3.1) over the domain Ω = [a, b] × [c, d] × [e, f ] ∈ R3 . For simplicity, we assume b − a = d − c = f − e > 0. We suppose that the domain Ω is uniformly partitioned with vertices G0 = {xijk }nijk=0 := {xi }ni=0 × {yj }nj=0 × {zk }nk=0 , where

xi = a + iΔx, yj = c + jΔy, zk = e + kΔz, and Δx = Δy = Δz = (b − a)/n. Let Gl be the set of vertices of the grid which is generated by binary subdividing G0 uniformly l times. Let φ be a piecewise tri-linear level-set function over the grid Gl , Φ be a tri-cubic spline approximation of φ over the grid G0 . In general, l is chosen as 0 or 1 or 2. In our implementation, we take l = 1. If l = 0, Φ and φ are defined on the same grid G0 . This is the simplest case. The aim of the following algorithm is to compute the spline level set function Φ. Algorithm Steps. 1. Initialization. Given a initial Γ, construct a piecewise tri-linear level-set function φ over the grid Gl . If necessary, apply a re-initialization step to set φ to be a signed distance function to Γ (see [3] for details). Convert φ to Φ (see [3]). 2. Evolution. Resampling Φ to obtain a new φ over the grid Gl . Compute L(Φ) and H(∇φ) in the thin shell N = {(xi , yj , zk ) ∈ Gl : |φ(xi , yj , zk )| < H}. Update φ in N for one time step to get φ˜ by an ODE time stepping method (see [3] for details).

6

3. Re-initialization. Applying re-initialization step to φ˜ in the shell N = {(xi , yj , zk ) ∈ Gl :

min

−1≤μ,ν,λ≤1

|φi+μ,j+ν,k+λ | ≤ H}.

to get a new φ (see [3] ). Convert φ to Φ (see [3] for details). Go back to step 2 if the termination condition is not satisfied. 4. Iso-contouring. Extract 3-sided or 4-sided iso-surface patches (vertices with normals) of Φ = c, where c is a given iso-value. Remark 2.1. For the problem of molecular surface construction, the grid size G0 should be less than the radii of atoms so that atoms are distinguishable from the level set surface. In our implementation, the grid size is chosen to be one-half of the minimal value of the atom radii. Remark 2.2. The aim of using l > 0 is to make φ a more accurate approximation of the signed distance function. The larger the value of l we use, the better approximation of the signed distance function we have. Since the scanned data to be approximated in general suffers from noise, we use the approximation Φ over a coarse grid G0 for denoising. Furthermore, for generating larger level-set surface patches, again a coarse grid needs to be used. Acknowledgment. A substantial part of this work in this paper was done when Guoliang Xu was visiting Chandrajit Bajaj at UT-CVC. His visit was also supported in part by the J. T. Oden ICES visitor fellowship.

References [1] C. Bajaj, J. Chen, and G. Xu. Modeling with cubic A-patches. ACM Transactions on Graphics, 14(2):103–133, 1995. [2] C. Bajaj, V. Pascucci, A. Shamir, R. Holt, and A. Netravali. Dynamic maintenance and visualization of molecular surfaces. Discrete Applied Mathematics, 127:23–51, 2003. [3] C. Bajaj, G. Xu, and Q. Zhang. A Higher Order Level Set Method with Applications to Smooth Surface Constructions. ICES Report 06-18, Institute for Computational Engineering and Sciences, The University of Texas at Austin, 2006. [4] C. Bajaj, G. Xu, and Q. Zhang. A Fast Variational Method for Construction of Smooth Molecular Surfaces. ICES Report 08-19, Institute for Computational Engineering and Sciences, The University of Texas at Austin, 2008. [5] C. Bajaj, G. L. Xu, S. Evans R. J. Holt, and A. N. Netravali. NURBS Approximation of A-Splines and A-Patches. Intl. J. on Computational Geometry and Applications, 2001. [6] S. Osher and R. Fedkiw. Level Set Method and Dynamic Implicit Surfaces. Springer, New York, 2003. 7

[7] S. Osher and J. Sethian. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79:12–49, 1988. [8] S. Osher and C.W. Shu. High-order essentially nonoscillatory schemes for hamilton-jacobi equations. SIAM Journal of Numerical Analysis, 28(4):907–922, 1991. [9] D. Peng, B. Merriman, S. Osher, H. Zhao, and M. Kang. A PDE-based fast local level set method. Journal of Computational Physics, 155:410–438, 1999. [10] Y. Zhang, G. Xu, and C. Bajaj. Quality meshing of implicit solvation models of biomolecular structures. Computer Aided Geometric Design, 23:510–530, 2006. [11] H. Zhao, S. Osher, B. Merriman, and M. Kang. Implicit nonparametric shape reconstruction from unorganized points using a variational level set method. Computer Vision and Image Understanding, 80(3):295–319, 2000. [12] W. Zhao, G. Xu, and C. Bajaj. An Algebraic Spline Model for Molecular Surfaces. In Proceedings of the 2007 ACM symposium on Solid and physical modeling, pages 297–302, Beijing, China, 2007. IEEE pub.

8

Advances in LES-based Turbulence Modeling Venkat Raman Dept. of Aerospace Engineering & Engineering Mechanics, The University of Texas at Austin, USA. [email protected]

Robert D. Moser Dept. of Mechanical Engineering, The University of Texas at Austin, USA. [email protected]

Abstract Recent research progress in turbulence modeling are discussed with focus on two topics: optimal LES and optimal theory based performance estimation. Optimal LES is an approach in which the subgrid model is formulated as minimum mean square error estimate. It has the advantage of being perfectly general, but requires information about the statistics of the small-scale turbulence. New advances in representing the multi-point correlations in both isotropic and wall-bounded flows are discussed, as is the performance of LES simulations based on these models. In the second part, some recent progress in analyzing the performance of models used to described turbulent combustion are discussed. Based on the concept of optimal error estimation, conventional models for the sub-filter variance of mixture-fraction are analyzed. A new dynamic procedure that provides improved performance is also discussed. Finally, the interaction of numerical errors with sub-filter models is studied in an effort to identify the more suitable formulations for LES-based combustion simulations.

Keywords: large eddy simulation, optimal estimation, dynamic modeling procedure. Large eddy simulation (LES) is now considered an attractive tool for studying turbulent flows. While many applications of LES have shown very good prediction of the flow field, many lingering questions regarding sub-filter modeling and the interaction of numerical and modeling errors still remain. One approach to describing these errors in the optimal LES procedure. Optimal LES is is based on the observation that the large scale fields being simulated do not provide sufficient information to reconstruct the small, scales, or even the evolution of the large scales [1]. The unknown small scales and therefore the LES evolution thus need to be treated statistically. Optimal LES models are formulated by postulating a model dependency and then minimizing the mean square error in representing the exact model term. Such models

Figure 1: Three-dimensional energy spectra from finite-volume optimal LES of infinite Reynolds number isotropic turbulence using a range of resolutions from 163 to 1283 . Also shown is the k −5/3 slope and the result of filtering a k −5/3 spectrum. can be formulated in terms of small-separation multi-point velocity correlations. The problem of LES modeling is thus explicitly reduced to the problem of modeling these correlations. Given this information, optimal models can be constructed that account for the errors of the numerical scheme used to solve the equations[2], and that are valid even in the presence of strong anisotropy and inhomogeneity[3]. Optimal LES is thus one approach which can address the shortcomings of current LES models. The required multi-point correlations include the 2-point second order, 3-point third order and 4-point fourth order correlations. For high Reynolds number isotropic turbulence, where a Kolmogorov inertial range exists, models for these correlations are available or have recently been developed[4]. They have been used to construct optimal LES models, which yield remarkably good results. For example, an isotropic LES based on a finite-volume discretization (filter) and the correlation models produces spectra that are consistent with the finite volume filtering of a Kolmogorov k −5/3 spectrum (figure 1). In a wall bounded flow, however, modeling the correlations is more difficult. A new formulation for the anisotropy and inhomogeneity of the two-point second correlation based on the structure tensors of Kassinos et al [5] has been developed. A comparison of thse model correlations with those determined from a DNS is shown in figure 2. Details of the optimal LES and multi-point correlation modeling approaches will be discussed, as will their appliaction in LES simulation.

r2

fig1

fig2

fig4

0.05

0.05

0.05

0

0

0

0

−0.05 −0.05

0

0.05

−0.05 −0.05

fig5

r2

fig3

0.05

0

0.05

−0.05 −0.05

fig6

0

0.05

−0.05 −0.05

fig7

0.05

0.05

0.05

0.05

0

0

0

0

−0.05 −0.05

0 r1

0.05

−0.05 −0.05

0 r1

0.05

−0.05 −0.05

0 r1

0

0.05

fig8

0.05

−0.05 −0.05

0 r1

0.05

Figure 2: Comparison of the two-point correlation tensor components from the model and DNS in the x-y plane with no separation in z. The correlations are centered at y + = 100 in a turbulent channel flow, with Reτ = 940 In the second part of the lecture, we discuss the performance of models used to describe turbulent combustion. Most combustion models use a passive scalar, termed mixture fraction, to describe the thermochemical state of the gas-phase. In LES, the filtered gas-phase properties can be obtained if the sub-filter variance of mixture fraction is known. This measure of sub-filter scalar energy has to be modeled and several models are available in literature. Recently, the optimal error estimation procedure was used to evaluate sub-filter models [6, 7]. It was found that the dynamic models, not surprisingly, provided the least error for a range of filter-widths. However, simple Taylor’s series based analysis of the dynamic procedure found that certain key terms are being neglected in the model formulation [7]. When included, the new procedure was found to provide lower errors compared to the conventional procedure. To understand the impact of numerics on model performance, apriori tests were conducted using different discretization schemes. It was found that the numerical error is of the same order as modeling error. Further, numerical errors have a “benign” effect on certain models leading to reduced overall error. These interesting findings also indicate that mathematical structure of the model is very important for reducing the inaccuracies due to numerical discretization [8, 9].

References [1] J. Langford and R. Moser, “Optimal LES formulations for isotropic turbulence,” Journal of Fluid Mechanics 398, 321 (1999). [2] P. S. Zandonade, J. A. Langford, and R. D. Moser, “Finite volume optimal large-eddy simulation of isotropic turbulence,” Physics of Fluids 16, 2255 (2004). [3] S. Volker, P. Venugopal, and R. D. Moser, “Optimal large eddy simulation of turbulent channel flow based on direct numerical simulation statistical data,” Physics of Fluids 14, 3675 (2002). [4] H. Chang and R. D. Moser, “An intertial range model for the three-point third-order velocity correlation,” Physics of Fluids 19, 105111 (2007). [5] S. Kassinos, W. Reynolds, and M. Rogers, “One-point turbulence structure tensors,” Journal of Fluid Mechanics 428, 213 (2001). [6] A. Moreau, O. Teytaud, and J. P. Bertoglio. Optimal estimation for large-eddy simulation of turbulence and application to the analysis of subgrid models. Phys. Fluids, 18, 1-10, 2006 [7] G. Balarac, H. Pitsch and V. Raman, “Modeling of sub-filter scalar variance using the concept of optimal estimator”, Physics of Fluids, 20 035114, 2008 [8] G. Balarac, H. Pitsch and V. Raman, “Modeling of the sub-filter scalar dissipation rate using the concept of optimal estimators”, Submitted to Physics of Fluids, 2008 [9] V. Raman, G. Balarac, and H. Pitsch, “Minimizing numerical errors in the computational sub-filter scalar variance”, To be Submitted to Physics of Fluids, 2008

Modeling Two Phase Flow Dynamics for Deformable Interfaces (Extended Abstract) Chandrajit Bajaj, Albert Chen, Richard Hankins, Bong-Soo Sohn Department of Computer Science, Center for Computational Visualization, Institute for Computational Engineering and Sciences University of Texas, Austin, TX 78712 Email: [email protected]

Abstract We describe a Stokesian (slow viscous) flow model for producing interacting deformable surfaces. This captures several phenomena in their nativity, such as cell membranes interacting, or soft docking of flexible molecules, etc. Starting from any initial configuration of closed, compact ”interfaces”, the interfaces are continually evolved, as well as deformed based on the relative viscosities and the interfacial tension. The velocity computation and the interfacial dynamics are achieved via a Boundary Element formulation of the governing Stokesian flow equation, while the interface evolution and topology maintenance utilizes level set representation and underlying function updates. Effects such as coalescence, break-up, and additional near interface interactions, can also be accurately captured. These last three effects in particular require adaptive refinement of meshed geometry and controlled coupling of the numerical errors in computation to yield topologically realistic looking phenomological modeling.

1 Introduction We present a technique to model the interaction of deformable surfaces using two-phase Stokesian (slow viscous) flows. These interfaces can represent air bubbles in a viscous liquid, oil droplets in a suspension, or cellular membranes subjected to hydrodynamics forces. Two-phase fluid simulations have been a topic of interest in computer graphics. In [5] Foster and Metaxas describe a technique for simulating free-surface water flows by solving the full Navier-Stokes equations in three dimensions via a finite difference scheme and representing the free surface with massless marker particles. In [4] Foster and Fedkiw implement a solution to the full Navier-Stokes equations using a modified version of the semi-Lagrangian technique introduced by Stam in [10] to capture the complex behavior of free water surfaces. To represent the air-water interface, they used a level set method as introduced in [9]. In [3] a particle level set method was introduced for improved interface capturing. Massless marker particles are advected with the level set data and used to repair the level set in regions of degradation due to the use of a coarse grid for animation. The combination of a semi-Langrangian 1

Navier-Stokes solver coupled to a particle based level set method allows for complex modeling of the dynamics at the air-water interface. In all of the above treatments, the density of the air is assumed to be zero. This is a reasonable approximation as the density of air is 1000 times less than that of water. A constant pressure boundary condition is applied everywhere at the interface and the Navier-Stokes equations are solved only within the liquid region. Using these techniques, while visually realistic simulations can be achieved which capture the convective (and turbulent) interaction of gases and liquids, many interesting features of air-water interfacial dynamics arising from non-convective, non-turbulent, and slow viscous flows cannot be observed. These interfacial deformations and dynamics arise from interfacial tension and the curvature of the interface. In order to capture such interfacial deformations in a simulation, for even visual realism, the Navier-Stokes equations must be accurately solved in both fluid regions, and pressure boundary conditions must be applied at the two phase fluid interface, while carefully orchestrating the accuracy of interfacial geometry and the interfacial velocities. Such an accurate two phase fluid simulation to produce realistic visual animations of deforming interfaces, is the main contribution of this paper. In [7], Hong and Kim with perhaps similar goals to ours, modified the semi-Lagrangian scheme [10] coupled with the volume-of-fluid method (VOF) introduced in [6] to simulate bubbles in liquids. The authors calculate the interfacial tension, and thereby are able to capture certain interfacial deformations, the VOF method has accuracy limitations for effects such as bubble flattening, coalescence, bubble necking, and break-up, and additional subtle near bubble interaction. We achieve this accuracy through our precise representation of the interface geometry, coupled to a topology tracking technique and a stable and accurate boundary element fluid solver, allowing us to observe these phenomena at significantly higher resolution. Our solution can be broken down into three main sub-areas: accurate boundary element mesh representation of the interface, careful error bounded calculation of physical quantities (velocities, surface tension) on the interface by regularization and adaptive mesh refinement, and topological tracking for precise near-bubble interactions. For the calculation of physical quantities we use the boundary element method (BEM) on adaptive geometries. Given an initial configuration of bubbles, our adaptive BEM solver, estimates the velocities on the interfacial boundary with greater accuracy, due to several advantages it has over competing methods. It reduces the dimensionality of the problem by one, and focuses computational effort on the boundary which, for two phase simulations, is the region of interest thereby yielding superior accuracy of the BEM over both finite element methods (FEM) and finite difference methods (FDM), where the interface is represented indirectly using a discretized volume domain. We also present a dynamic remeshing algorithm for smooth evolving interfaces of an objects such as bubbles. The interface is discretized to triangular or quadrilateral mesh where the velocity of each vertex is computed by the BEM. Since an interface changes its geometric shape including surface area and curvature distribution, meshes with fixed vertex count and connectivity cannot represent the evolving interfaces accurately. Therefore, the number of vertices and vertex connectivity need to be adjusted according to given geometric properties of an interface for each timestep. The quality of triangular or quadrilateral elements, often measured with element shape, also needs to be good enough for accurate BEM calculation. The main difficulty occurs when the topology of interfaces change. For instance, two bubbles may coalesce or a single bubble may break up into two bubbles. We utilize up-sampling / down-sampling methods and maintain a dynamic octree for tracking and controlling the mesh topology. 2

2 Deformable Interfaces A sketch of our method of accurate two-phase simulations is below. Details of each individual computational step is given in [1]. Algorithm Sketch 1. Initial interface : specify an initial configuration and level set representation of bubbles, as well as initialize the fluid physical parameters 2. Interface parametrization : quality mesh extraction for the level set 3. Interfacial velocity computation : use an accurate BEM to calculate interfacial velocities and analyze its associated errors 4. Oracle : check topology change condition through the Oracle and modify underlying functions that define the bubble interface. This process involves channel identification and modification of interfaces. 5. Interface evolution : evolve the interface using a level set method 6. Interference check : check geometric interference and adjust timestepping, go back to mesh extraction and BEM calculation for continuing the simulation.

2.1 Initialization of bubbles Initial geometry data is specified and given as inputs to the boundary element solver. The sign distance field to the given initial geometry is computed using an implicit C2 cubic B-Spline level set function [2]. Further normals and curvatures of the interfaces are then estimated from this piecewise analytic representation.

2.2 Stokes flow and boundary integral equations Let us denote ρ as the fluid density, μ as the fluid viscosity, u as the fluid velocity field, p as the fluid pressure. In addition, let g represent some external field, such as gravity, that acts on the fluid. The full Navier-Stokes equations for incompressible flow are given by, ρ ( ∂∂ut + u · ∇u) = −∇p + μ ∇2 u + ρ g, and ∇ · u = 0. By specializing to flows with low Reynolds number and low acceleration parameter we ignore the inertial and convective terms in the momentum equation. In this limit, viscous forces dominate and we have a linearized version of the Navier-Stokes equations called Stokes equation, −∇p + μ ∇2 u + ρ g = 0 Our choice of numerical solution technique requires that we reformulate the governing equations as integral equations over the boundary interfaces. Since we are interested in updating the interface between the two fluids we only need an integral relation for the interfacial velocity. The most general expression for the velocity at a point x0 on a bubble interface Γ is given by [14], uk (x0 ) =

1 1−λ 4π 1 + λ

 Γ

ui (x)Ti jk (x, x0 )nˆ k (x)dΓ(x) + Fj (x0 ) 3

(2.1)

Fj (x0 ) =

2 1 (u∞ (x0 ) − 1+λ j 8πμs

 Γ

f (x)nˆ i (x)Gi j (x, x0 )dΓ(x))

(2.2)

2.3 Numerical Technique In the above relations we assume N bubbles with varying viscosities immersed in a Stokes fluid. The viscosity ratio of the bubble to the fluid it is embedded in is, λ = μμs , and u∞j is the asymptotic Stokes velocity that the bubbles are immersed in. The function f is the boundary condition specifying the pressure difference at the interface and is given by, f = 2γκ + (ρ f luid − ρ bubble )z. In the above, γ is the surface tension of the drop, nˆ is the normal pointing into the suspending fluid and κ = 12 ∇ · n is the extrinsic mean curvature of the boundary. For our calculations we use the free space Green’s functions for Stokes flow given by, Gi j (x − x0 ) =

δi j xˆi xˆj + 3 , xˆ = x − x0 r r

and Ti jk (x − x0 ) = −6

xˆi xˆj xˆk r5

We generate a quadrilateral mesh from our B-spline level set approximation[11], to represent the interface separating the two fluids and discretize the integral equation. We break up the integrals over each surface into a sum of integrals over each quadrilateral face of the mesh. The integration over most faces may be done using standard quadrature techniques. However, the Green’s functions appearing in these expressions show divergent behavior as the evaluation point approaches the surface over which we are integrating. These singular and near-singular integrals must be handled carefully in order to obtain accurate results from a BEM when two droplets are close. Details are given in [1].

2.4 Error Analysis Proper error analysis is important for our algorithm. We use this information to determine the validity of the boundary element calculation and how to improve it by surface mesh refinement. We describe two different error tools that we use for feedback during the simulation. The first one is a global error measure and derives from the incompressibility condition of the governing equations we are using to model the fluids. This condition in integral form becomes,  Γ

u · ndΓ ˆ =0

(2.3)

We calculate this value and use it to determine which octree level to mesh the geometry for each timestep. The other method we implement is a local measure of the error. We use bilinear interpolation when calculating values on the surface from the values at the vertices which were obtained through solving the boundary element system. We implement a technique to estimate the error at each face of the mesh and add more geometry if the error is outside acceptable tolerance. The error analysis implemented is that

4

of [8]. The technique is chosen based on computational speed and ease of implementation. The error on the m-th face is estimated to be 

Eu (m)2 =

Γm

| u0 − uˆ |2 dΓm

(2.4)

where u0 is the predicted exact solution which we approximate by higher order interpolation and uˆ is the numerical solution for the velocity. Finally we measure the relative error by calculating, Erel,u (m)2 = 

Eu (m)2 0 2 Γ | u | dΓ

(2.5)

While rigorous mathematical bounds do not exist for local errors in collocation methods this technique serves its purpose in quantifying the error regions of sparse geometry where the calculation could be improved and have been implemented by authors in boundary element methods as well as finite element methods [12] [13]. This error analysis is done at each time step following the boundary element calculation of the interfacial velocities. Given a user defined tolerance ε1 > 0 we check that the velocity error computed over each face of the mesh satisfy Erel,u < ε1 . If this condition is satisfied then we proceed with the interface evolution. If it is not satisfied then we refine the faces of the mesh where the tolerance is exceeded and recalculate the interfacial velocities using this refined mesh.

2.5 Topology Control Assume f t is a function at time t where its level set represents deformable interfaces.  f t is a piecewise t t trilinear function that approximates f . M is a mesh that approximates the level set. The level set topology defined in  f t is preserved during mesh extraction process. Boundary Element Method (BEM) is applied to computing velocities which can be used for updating t+1 at time the function  f t at time t to evolve the interface in viscous flows. This generates the function f t, where level set mesh Mt+1 can be extracted. Topology changes of bubbles may occur under various conditions such as minimum distance or contact surface area between two bubbles. The moment of topology changes can be also chosen manually. We introduce an oracle that is an independent procedure to decide whether the topology change occurs or not based on the user-specified conditions. We also need a remeshing procedure to actually change the topology of meshes for the bubbles if the oracle decides that. Our algorithms for the oracle and remeshing support coalescence and breakup of bubble interfaces. Details are given in [1].

2.6 Interface Update The interface is updated using our higher order level set method [2]. The evolution of the level set is governed by the level set equation,

∂ ϕi  + v · ∇ϕi = 0 ∂t

5

(a) t = 1

(b) t = 40

(c) t = 51

(d) t = 66

(e) t = 100

(f) t = 200

Fig 3.1: Bubble coalescence. Several timesteps are shown here depicting bubble coalescence. Bubbles don’t tend to coalesce easily in pure Stokes flow. For the simulation a term was added to the fluid solver to simulate the intermolecular forces at work that cause coalescence. We see that the bubbles tend to flatten out as they approach each other. A sudden joining occurs after they have been in close proximity for enough time. The joined bubble returns to a spherical shape due to surface tension.

3 Results 3.1 Implementation We have an implementation of the above technique that runs on a Linux platform. The implementation consists of two main libraries: a meshing library, and a boundary element calculation library. These two sets of code are packaged in a graphical user interface where initial data can be input and physical parameters are defined. The package is also capable of outputting mesh data for analysis and rendering. Figure 3.1 shows results of a simulation of two bubbles coalescing. The input data are two uniform spheres separated by a small distance. The viscosity ratio is λ = 0.5 and a small gravity field is enabled so that the bubbles rise slightly due to buoyant forces. There is no asymptotic flow in the suspending fluid. The simulation shows that bubbles tend to exist in a flattened state before intermolecular forces take over and the bubbles ultimately join. Figure 3.2 shows two droplets deforming as they pass each other in a shearing flow given by u∞ = (γ z, 0, 0) where γ is a user defined parameter that controls the strength of the flow that can be used to 6

(a) t = 1

(b) t = 20

(c) t = 30

Fig 3.2: Bubbles deforming in shear flow. This simulation shows bubbles deforming as they move past each other in a shearing flow. The viscosity ratio for the simulation is 0.5 and a small gravity field is applied to make the bubbles rise slightly. The shear flow is keyframed to slowly dissipate. We observe the shearing and deformations as the bubbles interact and then return to their spherical shapes as the speed of the flow decreases. control the strength of the field. For this animation we set γ = 1 and add a small gravity field. We observe the bubbles slowly rising and moving past each other. As the lower bubble rises it begins to move to the right as it moves to a height where the asymptotic flow is positive. The flow is keyframed to gradually dissipate over the coarse of the simulation and we see the bubbles slow down and return to their spherical shapes. Within our simulations there are several parameters which may be adjusted to create a desired simulation based animation. The user can specify an initial configuration of droplets and a Stokes flow to embed the droplets in. This flow has parameters which control the speed and variation in the pressure gradient along the different spatial axes. In addition, the strength of the gravity field or other long range body forces may also be specified by the user and short range forces may be added to simulate intermolecular interactions. Adjustment of these parameters leads to very different results. The viscosity ratio parameter is particularly important. As a special case we set λ = 1 in the boundary integral equation eliminating the integral over the double layer potential. The result is that there is no linear system to solve so interfacial velocities are computed simply by evaluating, uk (x0 ) =

2 1 (u∞ (x0 ) − 1+λ j 8πμs

 Γ

f (x)nˆ i (x)Gi j (x, x0 )dΓ(x))

It should be noted that all phenomena of bubble interactions mentioned here can be observed in this special case, but the computational demands are much less than the λ = 1 solution. Acknowledgment. This research was supported in part by NSF grants IIS-0325550, CNS-0540033 and NIH contracts P20-RR020647, R01-EB00487, R01-GM074258, R01-GM07308

7

References [1] C. Bajaj, A. Chen, R. Hankins, and B. Sohn. A C2 cubic B-Spline Level Set Approach to Simulating Deformable Interfaces. TICAM Report 08-xx, Texas Institute for Computational and Applied Mathematics, The University of Texas at Austin, 2008. [2] C. Bajaj, G. Xu, and Q. Zhang. A higher order level set method with applications to smooth surface constructions. TICAM Report 06-18, Texas Institute for Computational and Applied Mathematics, The University of Texas at Austin, 2006. [3] J. Ferziger D. Enright, R. Fedkiw and I. Mitchell. A hybrid particle level set method for improved interface capturing. J. Comp. Phys., 183:83–116, 2002. [4] N. Foster and R. Fedkiw. Practical animation of liquids. In ACM Press/ACM SIGGRAPH, E. Fiume, Ed., Computer Graphics Proceedings, Annual Conference Series, pages 23–30, 2001. [5] N. Foster and D. Metaxas. Realistic animation of liquids. Graphical Models and Image Processing, 58(5):471–483, 1996. [6] C.W. Hirt and B.D. Nichols. Volume of fluid (vof) method for the dynamics of free boundaries. J. Comp. Phys., 39:201–255, 1981. [7] J.-M. Hong and Kim C.-H. Animation of bubbles in liquid. In In Proceedings of Eurographics 2003, pages 253–262, 2003. [8] E. Kita and N. Kamiya. A new adaptive boundary element refinement based on simple algorithm. Mech Res Commun., 18(4):177–86, 1991. [9] S. Osher and J.A. Sethian. Fronts propagating with curvature dependent speed: algorithms based on hamilton-jacobi formulations. J. Comp. Phys, 79:12–49, 1988. [10] J. Stam. Stable fluids. In In Proceedings of SIGGRAPH 99, ACM SIGGRAPH/Addison Wesley Longman, Computer Graphics Proceedings, Annual Conference Series, ACM, pages 121–128, 1999. [11] Y. Zhang and C. Bajaj. Adaptive and quality quadrilateral/hexahedral meshing from volumetric data. Computer Methods in Applied Mechanics and Engineering (CMAME), 195:942–960, 2006. [12] O.C. Zienkiewicz and J.Z. Zhu. A simple error estimator and adaptive procedure for practical engineering analysis. Int. J. Numer. Methods. Engng., 24:337–57, 1987. [13] O.C. Zienkiewicz, J.Z. Zhu, and N.G. Gong. Effective and practical h-p-version adaptive analysis procedures for the finite element method. Int. J. Numer. Methods. Engng., 28:879–91, 1989. [14] Z.A. Zinchenko, M.A. Rother, and R.H. Davis. A novel boundary integral algorithm for viscous interaction of deformable drops. Phys. Fluids, 9:1493–1511, 1996.

8

Isogeometric Analysis of Fluid-Structure Interaction Y. Bazilevs, V.M. Calo, T.J.R. Hughes Institute for Computational Engineering and Sciences, The University of Texas at Austin, USA e-mail: {bazily,victor,hughes}@ices.utexas.edu

Y. Zhang Department of Mechanical Engineering, Carnegie Mellon University, USA [email protected]

Abstract Isogeometric analysis is a recently developed methodology based on technologies that were originated in the field of computational geometry and widely used in design, graphics and animation. It includes standard finite element analysis as a special case, but offers other possibilities that are unique and powerful. It allows more precise and efficient geometric modeling, it simplifies mesh refinement, and it possesses superior approximation properties. Isogeometric analysis has been applied to numerous problems in solid and fluid mechanics. This paper describes an isogeometric formulation for fluid-structure interaction of incompressible fluid flow and nonlinear solids. The fluid discretization derives from a residual-based variational multiscale formulation, applicable to laminar and turbulent phenomena. Both fluid and solid domains may undergo large motions, and the geometry and kinematics are fully compatible across fluid-structure interfaces. A strongly coupled solution algorithm is adopted to preclude instabilities that often afflict weaklycoupled procedures. Applications to the human cardiovascular system are emphasized.

Keywords: Isogeometric Analysis, NURBS, Fluid-Structure Interaction, Vascular Modeling, Navier-Stokes Equations, Elastic Arterial Wall, Mesh Movement, Blood Flow

1

Introduction

Isogeometric Analysis based on NURBS (non-uniform rational B-splines) was first introduced in [1] as an attempt to improve on and generalize the standard finite element method. Further study of isogeometric analysis showed that results superior to standard finite elements are obtained in the context of structural vibrations [2]. Mathematical analysis of the isogeometric

approach was performed in [3]. Optimal approximation estimates in p, the polynomial order used to define NURBS functions, were obtained for h-refined meshes. Stability and optimal convergence was proved mathematically and verified numerically for problems of compressible and incompressible elasticity, Stokes flow, and scalar advection-diffusion. In this paper, NURBSbased isogeometric analysis is applied to fluid-structure interaction (FSI) problems with particular emphasis on arterial modeling and blood flow (see [4] for a more detailed exposition). It is believed that the ability of NURBS to accurately represent smooth exact geometries, that are natural for arterial systems, but unattainable in the faceted finite-element representation, and the high order of approximation of NURBS, should render fluid and structural computations more physiologically realistic. This work adopts the arbitrary Lagrangian-Eulerian (ALE) framework. The arterial wall is treated as a nonlinear elastic solid in the Lagrangian description governed by the equations of elastodynamics. Blood is assumed to be a Newtonian viscous fluid governed by the incompressible Navier-Stokes equations written in the ALE form. The fluid velocity is set equal to the velocity of the solid at the fluid-solid interface. The coupled FSI problem is written in a variational form such that the stress compatibility condition at the fluid-solid interface is enforced weakly. The ALE equations require the specification of the fluid region motion. This motion is found by solving an auxiliary static linear-elastic boundary-value problem for which the fluid-solid boundary displacement acts as a Dirichlet boundary condition. Galerkin’s method is employed for the structural and the fluid subdomain motion parts of the formulation, while a residual-based multiscale method for the fluid equations. The resultant semi-discrete equations are advanced in time using the generalized-α algorithm. The kinematic constraint is enforced strongly by requiring basis functions to be C 0 -continuous across the fluid-solid interface. The coupled nonlinear system resulting from the NURBS discretization of the FSI equations is solved monolithically, that is, the fluid, the structural, and the mesh solution increments in the Newton iteration are obtained simultaneously. The effect of the structural and the mesh motion on the fluid equations is included in the left-hand-side matrix for robustness. The coupled system is solved iteratively by the GMRES procedure with simple diagonal scaling.

2 2.1

Numerical examples Blood flow in an idealized aneurysm

In this test case, taken from [5], we examine pulsatile flow in an idealized aneurysm. The problem setup is shown in Figure 1. A time-periodic velocity waveform, specified at the inflow plane, is parabolically distributed over the circular surface. The domains proximal and distal to the aneurysm region are assumed to have rigid walls, while the aneurysm wall is elastic. A resistance boundary condition is applied at the outflow.

Figure 1: Idealized aneurysm problem setup. Figure 2 shows the inflow and outflow waveforms. Note the outflow lags the inflow due to the distensibility of the aneurysm wall. This well-known phenomenon was also observed in practice as well as computations of other researchers. Figure 2 also shows excellent agreement with reference results of [5].

2.2

Blood flow in a patient-specific abdominal aorta

This computational example makes use of patient-specific geometry. Data for this model was obtained from 64-slice CT angiography of a healthy male over 55 years of age. Some preprocessing, including contrast enhancement, denoising, and segmentation, was necessary in order to find the lumenal surface of the blood vessel. Figures 3(a) - 3(c) show the geometrical model, the control mesh, and the NURBS mesh of the abdominal aorta. As in the previous example, we specify a time-periodic velocity waveform at the inflow, while all outflows are assigned a resistance boundary condition. Figure 3(d) shows contours of the arterial wall velocity magnitude plotted on a deformed configuration during early systole.

References [1] T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194 (2005) 4135-4195.

Figure 2: Idealized aneurysm. Inflow and outflow waveforms. Notice the time lag attributable to the distensibility of the wall. [2] J.A. Cottrell, A. Reali, Y. Bazilevs, and T.J.R. Hughes. Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering, 195 (2006) 52575296. [3] Y. Bazilevs, L. Beirao da Veiga, J.A. Cottrell, T.J.R. Hughes, and G. Sangalli. Isogeometric analysis: Approximation, stability and error estimates for h-refined meshes. Mathematical Methods and Models in Applied Sciences, 16 (2006) 1031-1090. [4] Y. Bazilevs, V.M. Calo, Y. Zhang, and T.J.R. Hughes. Isogeometric fluid-structure interaction analysis with applications to arterial blood flow. Computational Mechanics, 38 (2006) 310-322. [5] A.-V. Salsac, M.A. Fernandez, J.-M. Chomaz, and P. Le Tallec. Effects of the flexibility of the arterial wall on the wall shear stress and wall tension in abdominal aortic aneurysms. Proceedings of the 58th Annual Meeting of the Division of Fluid Dynamics, Chicago, IL, November 2005.

Figure 3: Patient-specific abdominal aorta. (a) Geometrical model; (b) Control mesh; (c) NURBS mesh; (d) Contours of arterial wall velocity magnitude plotted on a current configuration during early systole.

Modeling Complexities in Turbulent Spray Combustion Venkat Raman Dept. of Aerospace Engineering & Engineering Mechanics, The University of Texas at Austin, USA. [email protected]

Abstract Turbulent spray combustion is a common process found in aircraft and automobile engines as well as chemical reactors. Spray combustion is a multilscale multiphysics problem involving the nonlinear interaction of spray evolution and evaporation, gas-phase turbulent mixing, and combustion. A detailed description of this complex process requires models for each of the physical processes involved. In this study, we focus on the development of the multiphase combustion models. Currently, spray combustion models are directly extended from corresponding single-phase combustion models. However, combustion in the presence of evaporating droplets has many unique features that render invalid many of the assumptions underlying the single-phase models. The objective of this work is to introduce a novel probability density function (PDF) based approach, which addresses many of the multiphase combustion modeling challenges. Numerical algorithms and preliminary results are presented here.

Keywords: Spray combustion , large eddy simulation, probability density function, combustion regime. Liquid fuel based combustion is widely encountered in aircraft and automobile engines as well as chemical reactors. In typical combustors, the liquid fuel is sprayed using an atomizer that produces a fine mist of droplets. The droplets evolve in the background turbulent gasphase flow and evaporate. This fuel vapor then mixes with the gas-phase oxidizer and reacts in a high-temperature environment. The proper dispersion of fuel inside the combustor is critical in maintaining stability and in reducing emissions. The complete description of the spray combustion system requires models for tracking the spray droplets, gas-phase turbulent flow, and combustion. Since the focus of this work is the modeling of the combustion process, standard and state-of-the-art approaches for the other two components will be used. The spray droplets will be evolved using a Lagrangian approach. The gas-phase turbulent flow will be described using the large eddy simulation (LES) method.

Single-phase combustion models are derived based on the nature of the combustion regime. For instance, the fuel and oxidizer could be molecularly mixed before entering the combustion chamber leading to premixed combustion. If the fuel and oxidizer mix insider the combustor, a non-premixed combustion process is supposed to exist. The controlling parameters and the flame structure are vastly different in the two cases giving rise to very different combustion models. Since the combustion regime is determined by the boundary conditions, the combustion model can be chosen apriori. In spray combustion, the fuel is released in vapor form through droplet evaporation, which in itself depends on the physical dispersion of droplets and the rate of liquid evaporation. Both these processes are dependent on the gas-phase turbulence and combustion. Consequently, the combustion regime can change within the reactor due to a number of reasons including the relative rates of evaporation and mixing, effect of droplet inertia on the gas-phase turbulence, and the spatial distribution of the reaction zone [1]. Fig. 1 shows the flame structure for two different droplet structure in a coflowing oxidizer stream. When the droplets are close enough, the flame is pushed outside and a flame front typical of premixed combustion is formed. If the droplets are dispersed far enough, the flame front advances into the inter-droplet spacing leading to a partially-premixed combustion process. While these cases demonstrate the extreme effects, droplet structure and its evolution play a critical role in determining the combustion process. For this reason, a combustion model that presumes the combustion regime will not be fully valid in spray combustion. In this work, a novel method termed the probability density function (PDF) approach is used [2]. The main advantage in this approach is that the chemical source terms appear closed and do not require modeling. This, in turn, implies that the combustion regime is not fixed in a simulation. While this is certainly advantageous, mixing of scalars has to be modeled and poses a tremendous challenge. In the lecture, models recently proposed by the author to address this issue will be discussed. In the PDF approach, the transport equation for the joint-PDF of all chemical species and other scalars that define the thermochemical state of the system is solved directly. This transport equation is high dimensional and cannot be solved using conventional finite-volume or finite-difference discretization scheme. Typically, a Lagrangian Monte-Carlo method is used. In practical simulations, this PDF solver is coupled with the spray and turbulence models and evolved in a temporally accurate manner. This coupled solver was used to simulate canonical flows in an effort to verify and validate this computational tool (Fig. 2). The nature of information exchange between the solvers is very important in the context of spray combustion. It will be demonstrated that the current techniques for simulating such flows change the underlying physics of the combustion process.

References [1] V. Raman and H. Koo, “Effect of Droplet Inertia on Evaporation and Spray Combustion”, Submitted to the The 32nd International Combustion Symposium, 2008 [2] V. Raman and H. Koo, “LES/Filtered-Density Function Approach for Turbulent Spray Combustion”, Submitted to Combustion and Flame, 2008

Figure 1: Multiple droplet combustion with detailed resolution of the droplet interface. Droplets are dispersed closely (right) and (right) far apart.

Figure 2: Spray droplet population superimposed on the instantaneous contours of gas-phase temperature. The inset shows more details about the reaction zone.