 

 

 

LIBRARY ﬂ

Michigan State 4
University J

This is to certify that the
thesis entitled

APPLICATION AND DEVELOPMENT OF NUMERICAL
METHOD FOR MOVING BOUNDARY AND COMPLEX
GEOMETRY PROBLEMS

presented by
Yanbing Li

has been accepted towards fulﬁllment
of the requirements for the

 

 

 

 

. Department of Mechanical
M. 8. degree In Engineering
’2
f Major Professor's Signature
I 1/3 /0 2,
Date

MSU is an Afﬁrmative Action/Equal Opportunity Institution

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 cJCIFIC/DateDuepss-p. 1 5

THE APPLICATION AND DEVELOPMENT OF NUIVIERICAL METHOD FOR
MOVING BOUNDARY AND COMPLEX GEOMETRY PROBLEMS

By

Yanbing Li

A THESIS
Submitted to
Michigan State University

in partial fulﬁllment of the requirements
for the degree of

MASTER OF SCIENCE
Department of Mechanical Engineering

2002

ABSTRACT

THE APPLICATION AND DEVELOPMENT OF NUMERICAL METHOD FOR
MOVING BOUNDARY AND COMPLEX GEOMETRY PROBLEMS

By
Yanbing Li
Various methods have recently been introduced to alleviate the difficulties associated
with simulating moving boundary problems. Among them, the ﬁctitious domain method
is explored and methods to mitigate the errors associated with this technique are
presented in this work. A summary of moving mesh strategies currently in use in
commercial computational ﬂuid dynamics soﬁware is ﬁrst provided. Signiﬁcant issues
associated with moving gn'd techniques such as the interpolation errors, the human cost
for mesh generation, and the quality of the resulting mesh are discussed by solving, with
the help of a commercial code, a variety of ﬂuid dynamics problems assocrated with in-
cylinder ﬂows of internal combustion engines. The Lagrange multiplier/ﬁctitious domain
method is then discussed by solving simple heat transfer problems. The ﬁctitious domain
method is based on the use of Lagrange multipliers that do not match with an underlying
mesh. Such an approach introduces errors on the adjacent nodes and a parametric study
of various factors affecting the quality of the results is performed. Two approaches for
reducing the errors are proposed. A ﬁrst method uses modiﬁed boundary conditions to
reduce the errors on the adjacent nodes. The method is based on a predictor/corrector
scheme and is called the “ﬁctitious constraint” method. A second method consists of
simply modifying the shape of the boundary by matching the complex boundary with the
underlying mesh. Improvements in the solutions by using these techniques are illustrated

with the help of one-dimensional and two-dimensional problems.

To my wife, my parents, and my grandmother in-law

iii

TABLE OF CONTENTS

LIST OF TABLES vii
LIST OF FIGURES viii
NOMENCLATURE xvi
1. Introduction . l
2. An Introduction to Moving Grid Technique 5
2.1 The Governing Equations 6

2.2 Transformed Governing Equations in a Body-Fitted Coordinate System 6
2.3 The Geometric Conservation Law (GCL) 8
2.4 Solution Interpolations 10

3. Application of Moving Grid Technique to In-Cylinder Flow Simulation 12

3.1 In-Cylinder Flow Simulations of Internal Combustion Engine 12
3.2 Moving grid strategy in current CFD codes 14
3.2.1 Spring Deforming 14
3.2.2 Local Regridding 16
3.2.3 Layering 17
3.3 Turbulence Modeling 18
3.4 Description of the Problem 21
3.5 Numerical Setups of the Simulation 23
3.6 Experiment Veriﬁcations 24

3.7 Simulation Results Discussion 27

iv

3.8 Issues with Moving Grid Technique for In-cylinder Flow Simulation 38
4. An Introduction to Lagrange Multiplier/ Fictitious Domain Methods,
Methods for Reducing Errors in LM/FDM Applied to Heat Transfer
Problems 47
4.1. Introduction 47

4.2. Lagrange Multiplier/Fictitious Domain Methods and the Dirichlet

Problem 48
4.2.1. A Model Problem 48
4.2.2. Lagrange Multiplier/Fictitious Domain Formulation 49
4.2.3. Finite Element Approximations 50

4.3. Effects of the Presence of a Multiplier Within an Element 52
4.4. Parametric Study on the Sources of the Errors 53
4.4.1. Auxiliary Boundary Condition Inﬂuence 57

4.4.2. Mesh Size and Relative Position of Lagrange Multiplier Inside

Element Inﬂuence 59

4.5. Reducing the Error 64
4.6. Fictitious Constraint Methods 68
4.7. Shape Reconstruction Method 76
4.8. LM/FDM in Two Dimensions 80
4.8.1. Boundary Integrations 80

4.8.2. Implementation of Fictitious Constraint LM/FDM in 2-D 82

4.8.3. Implementation of Shape Reconstruction LM/FDM in 2-D 86

5. Conclusions 89
APPENDIX 91

Turbulence Boundary Layer

1. The Concept of y+

2. Wall Function

3. Near Wall Mesh Generation for Wall Functions

BIBLIOGRAPHY

vi

91

91

92

92

93

LIST OF TABLES

Table 1 Engine Speciﬁcations and Calculation Conditions
Table 2 Values of standard k—e model constants
Table 3 Numerical model data

Table 4 Computing the ﬁctitious constraint

vii

22

23

24

85

LIST OF FIGURES

Figure 1 Discretization of the original problem using a body ﬁtted mesh
Figure 2 Spring deforming on a simple cylinder at the end of compression
Figure 3 Local regridding applied to the valve region during intake stroke
Figure 4 Layering applied to the valve region during intake stroke

Figure 5 Geometry of the generic engine and the boundary deﬁnition
Figure 6 Computational Mesh of the generic engine

Figure 7 Average in-cylinder pressure comparision

Figure 8 Flow pattern comparison at CA 450 ATDC

Figure 9 Flow pattern comparison at CA 540 ATDC

Figure 10 Flow pattern comparison at CA 630 ATDC

Figure 11 Flow pattern comparison at CA 660 ATDC

Figure 12 Selected plans for the generic engine

Figure 13 Velocity distributions at 370 ATDC (x-plane and y-plane)
Figure 14 Velocity distributions at 380 ATDC (x-plane and y-plane)
Figure 15 Velocity distributions at 390 ATDC (x-plane and y-plane)
Figure 16 Velocity distributions at 400 ATDC (x-plane and y-plane)
Figure 17 Velocity distributions at 410 ATDC (x-plane and y-plane)
Figure 18 Velocity distributions at 420 ATDC (x-plane and y-plane)
Figure 19 Velocity distributions at 430 ATDC (x-plane and y-plane)
Figure 20 Velocity distributions at 440 ATDC (x-plane and y-plane)

Figure 21 Velocity distributions at 450 ATDC (x-plane and y-plane)

viii

l7

18

22

24

25

26

26

27

27

28

31

31

31

32

32

32

33

33

33

Figure 22 Velocity distributions at 460 ATDC (x-plane and y-plane) 34

Figure 23 Velocity distributions at 480 ATDC (x-plane and y-plane) 34
Figure 24 Velocity distributions at 510 ATDC (x-plane and y-plane) 34
Figure 25 Velocity distributions at 540 ATDC (x-plane and y-plane) 35
Figure 26 Velocity distributions at 550 ATDC (x-plane and y-plane) 35
Figure 27 Velocity distributions at 580 ATDC (x-plane and y-plane) 35
Figure 28 Velocity distributions at 610 ATDC (x-plane and y-plane) 36
Figure 29 Velocity distributions at 640 ATDC (x-plane and y—plane) 36
Figure 30 Velocity distributions at 660 ATDC (x-plane and y-plane) 36
Figure 31 Velocity distributions at 680 ATDC (x-plane and y-plane) 37
Figure 32 Velocity distributions at 700 ATDC (x-plane and y-plane) 37

Figure 33 Velocity distributions at 550(left) and 660 (right) ATDC (z-plane) 37
Figure 34 Velocity distributions at 680(leﬁ) and 700 (right) ATDC (z-plane) 38
Figure 35 Solution interpolation from “old” grid to “new” grid introduces error 39
Figure 36 Computational mesh at selected plans for a 50k generic engine numerical
model 40
Figure 37 Velocity absolute error distributions at selected plans for a 50k generic engine
numerical model 40
Figure 38 Velocity relative error distributions at selected plans for a 50k generic engine
numerical model 40
Figure 39 Computational mesh at selected plans for a 50k generic engine numerical

model (Method I) 41

ix

Figure 40 Computational mesh at selected plans for a 50k generic engine numerical
model (Method 11) 41
Figure 41 Computational mesh at selected plans for a 50k generic engine numerical
model (Method IH) 42
Figure 42 Streamline contours comparison (left: method 1, middle: method H, right:
method III) at selected plan (x-plan) for a 50k generic engine numerical model at CA 460
ATDC 42
Figure 43 Streamline contours comparison (left: method 1, middle: method H, right:
method IH) at selected plan (y-plan) for a 50k generic engine numerical model at CA 460
ATDC 42
Figure 44 Streamline contours comparison (left: method I, middle: method 11, right:
method IH) at selected plan (z-plan) for a 50k generic engine numerical model at CA 460
ATDC 43
Figure 45 Streamline contours comparison (left: method 1, middle: method H, right:
method IH) at selected plan (x-plan) for a 50k generic engine numerical model at CA 630
ATDC 43
Figure 46 Streamline contours comparison (left: method 1, middle: method H, right:
method H1) at selected plan (y-plan) for a 50k generic engine numerical model at CA 630
ATDC 43
Figure 47 Streamline contours comparison (left: method 1, middle: method H, right:
method H1) at Selected plan (z-plan) for a 50k generic engine numerical model at CA 630

ATDC 44

Figure 48 Average predicted in-cylinder turbulence kinetic energy comparison for three
moving mesh strategies 44
Figure 49 Turbulence kinetic energy distribution comparison (left: method 1, middle:
method H, right: method III) at selected plan (y-plan) for a 50k generic engine numerical
model at CA 460 ATDC 45
Figure 50 Turbulence kinetic energy distribution comparison (left: method 1, middle:
method H, right: method III) at selected plan (z-plan) for a 50k generic engine numerical
model at CA 460 ATDC 45
Figure 51 Turbulence kinetic energy distribution comparison (left: method 1, middle:
method 11, right: method 111) at selected plan (y-plan) for a 50k generic engine numerical
model at CA 630 ATDC 45
Figure 52 Turbulence kinetic energy distribution comparison (left: method I, middle:
method H, right: method 111) at selected plan (z-plan) for a 50k generic engine numerical
model at CA 630 ATDC 46
Figure 53 Illustration of the geometry associated to the Dirichlet problem presented in
(4.1) 49
Figure 54 Discretization mesh set, discretization of ﬁctitious domain (I and embedded
boundary 7 51
Figure 55 A linear approximation of the solution enforces the embedded boundary
constraint but cannot capture a discontinuity in the ﬁrst derivative at other location but
the nodes. 53
Figure 56 Absolute error /Relative Error introduced by the LM/FDM for 2-D steady state

heat conduction problem 53

xi

Figure 57 The extended problem possess discontinuity across the embedded boundary
54

Figure 58 LM/FDM solutions coincide with the exact solution at nodes when no

discontinuity locates inside the corresponding element 56
Figure 59 The inﬂuence of the auxiliary boundary condition 57
Figure 60 Deﬁnition of 6 58
Figure 61 a - I9 relation 58
Figure 62 The inﬂuence of the local gradient jump 59
Figure 63 L2 norm oscillates with the mesh reﬁned 60
Figure 64 Energy norm oscillates with the mesh reﬁned 60
Figure 65 Xkrel oscillates with the mesh reﬁned (a = 1) 61
Figure 66 Plots of the L2 norm of error versus element size 62
Figure 67 Plots of the energy norm of error versus element size 63

Figure 68 Plots of the L2 /energy norm of error versus the relative location of Lagrange

multiplier xkrel . 63

Figure 69 Solution comparison of the LM/FDM methods and BFM methods 64
Figure 70 Mesh reﬁnement may get a better solution 65
Figure 71 Mesh reﬁnement may produce a worse solution 65

Figure 72 By adjusting auxiliary boundary condition adjustment a better solution can be

obtained 65
Figure 73 High order element can get a better solution 66
Figure 74 L2 norm oscillates with the mesh reﬁned for quadratic elements 67

Figure 75 Energy norm oscillates with the mesh reﬁned for quadratic elements 67

xii

Figure 76 A LM/FDM Solution which satisfying a speciﬁc embedded boundary condition
(“Fictitious Constraint”) coincides with the exact solution at nodes 68
Figure 77 Linearity assumption of the boundary-neighboring zones: Variable
distributions in each separated boundary-neighboring zones (such as zone 1 and zone 2
which belong to different domain-original domain or auxiliary domain) have similar
slope 70
Figure 78 Derivative similarity assumption: Variable derivative of the LM/FDM solution
in each region (region 1 or region 2) is similar to the exact solution 71
Figure 79 Implementation of ﬁctitious constraint method in l-D (The solid line is the
exact solution; the dotted line is the LM/FDM solution and dashed line the ﬁctitious
constraint solution) 71
Figure 80 “Predict-Correct Fictitious Constraint” methods obtaines better solution 73
Figure 81 L2 norm comparison for three methods : Standard LM/FDM, Fictitious
Constraint LM/FDM, Body Fitted methods 74

Figure 82 Energy norm comparison for three methods: Standard LM/FDM, Fictitious

Constraint LM/FDM, Body Fitted methods 74

Figure 83 L2 norm is not very sensitive to relative location of the Lagrange multiplier

inside elements for Fictitious Constraint LM/FDM methods 75
Figure 84 Energy norm is not very sensitive to relative location of the Lagrange
multiplier inside elements for Fictitious Constraint LM/FDM methods 75
Figure 85 Plots of the L2 norm of error versus element size 76

Figure 86 L2 norm and energy norm is not very sensitive to the local gradient jump for

Fictitious Constraint LM/FDM methods 76

xiii

Figure 87 Shifting the location of Lagrange multiplier from xk to the nearest node x1 in

the corresponding interfacial element can obtain a better solution 77

Figure 88 Energy norm and L2 norm comparisons for three methods: Standard

LM/FDM, Fictitious Constraint LM/FDM, Shape Reconstruction LM/FDM 78

Figure 89 Energy norm and L2 norm are not very sensitive to local gradient jump

(compare to the standard LM/FDM methods) but much more sensitive (compare to the
ﬁctitious constraint LM/FDM) for Shape Reconstruction methods 79

Figure 90 Energy norm and L2 norm are not very sensitive to relative location of the

Lagrange multiplier inside elements (compare to the standard LM/FDM methods) but
much more sensitive (compare to the ﬁctitious constraint LM/FDM) for Shape
Reconstruction LM/FDM methods 80
Figure 91 Absolute errors and relative error for the stand LM/FDM using Gauss
Legendre quadrature 82
Figure 92 Absolute errors and relative error for the stand LM/FDM using Collocation
method 82
Figure 93 Illustration of the implementation ﬁctitious constraint LM/FDM 83
Figure 94 Absolute errors and relative error for the Fictitious constraint LM/FDM using
Gauss Legendre quadrature 86
Figure 95 Absolute errors and relative error for the Fictitious constraint LM/FDM using
collocation method 86
Figure 96 Illustration of the Shape Reconstruction in 2-D 87
Figure 97 Absolute errors (left) and relative error (right) for the shape reconstruction

LM/FDM using Gauss Legendre quadrature 88

xiv

Figure 98 Absolute errors (left) and relative error (right) for the shape reconstruction

LM/FDM using collocation method 88

XV

'0
“1-

”GI

Pstatic 3

T static 1

NOMENCLATURE

Roman Symbols
pressure, Psi
x, y, 2 component of velocity, m / s

x, y component of body force, N
turbulent kinetic energy, m2 /SZ

constants used in the a balance equation

constant used in mixing length turbulence model

velocity component, m/ .9

mean velocity component, m / s

ﬂuctuation from the mean velocity component 17,- , m / 3
mean pressure, Psi

ﬂuctuation from the mean pressure, Psi

static pressure, Psi

temperature, K

static temperature, K

turbulence intensity,

xvi

lit:

0k:

~(DI

mean velocity magnitude, m / 3

weak solution of the problem over the entire domain
weak solution of the problem over ﬁnite element
ﬁnite dimensional spaces for Q

ﬁctitious domain element

embedded boundary element

Greek Symbols

ﬂuid density, kg / m3

ﬂuid viscosity, Pa-s
turbulent ﬂuid viscosity, Pa -s
turbulent dissipation rate, m2 /s3

turbulent Prandtl number for k

mean ﬂuid density, kg / m3

ﬂuctuation from the mean ﬂuid density, kg / m3

spatial domain

boundaries of spatial domain

weak solution of the problem on the boundary

xvii

ATDC:

BDC:

BFM:

CA:

CFD:

GCL:

LM/FDM:

weak solution of the problem over the boundary elements

ﬁnite dimensional spaces for y

discretization of the ﬁctitious domain Q

discretization of the embedded boundary y

Abbreviations
after top dead center
bottom dead center
body-ﬁtted mesh method
crank angle
computational ﬂuid dynamics
geometric conservation law
Lagrange multiplier/ﬁctitious domain method

top dead center

xviii

Chapter 1
Introduction

Computational ﬂuid dynamics (CFD) still encounters difﬁculties when dealing with
geometrically complex and dynamically complex problems [4]. Among them, moving
boundary problems have been arising as a continuing challenging issue in a variety of
important engineering application areas [2]. Typical examples are heat transfer problems
and chemical reaction problems where there may be signiﬁcant changes of material
properties, physrcal-chemical properties or ﬂow features. In these situations, the
interfaces move under the interface of the ﬂow ﬁeld by internal or external forces and in
turn affect the behavior of the ﬂow, and may subject to instabilities under certain
conditions [2,3]; the prediction of the mechanism in such moving boundary systems is ‘
very important but very difficult to analyze “numerically” due to existing of the
complicating mechanisms around the interface region [2].

Numerous techniques exist for tracking moving boundaries [2,3]. Among them,
Lagrangian methods, which are based on a moving grid technique, are widely used in
commercial soft wares. The use of a body-ﬁtted grid system that conforms to the actual
geometry of the problem and is updated (redistribution and reﬁnement) allows to track
the interface (moving boundary) explicitly. The boundary condition can be applied at the
exact location of the interface at each time step. Application of moving grid techniques to
problems with complex geometry and extensive interfacial activities (e.g. large
deformation and topological change) requires a strong remeshing; this corresponds to an

increase of the human cost and thus limits its applications (such as design evaluations

with various geometry conﬁgurations). The formulation of a moving grid technique is
reviewed and related issues such as geometric conservation law of the coordinate

mapping and solution interpolation are discussed in chapter 2.

In the context of transient IC engine simulations, moving grid techniques are widely
adopted in today’s most commercial (FIRE, FLUENT, Star-CD) and no-commercial
(Kiva) codes. Numerical simulations with such CFD platforms are performed to study the
in-cylinder ﬂow activities and help the design evaluation [11, 13, 18-26]. In chapter 3, a
summary of the moving mesh strategies is provided; Simulations result using FIRE for an
Generic IC engine intake and compression stroke process are then preSented.
Experimental data are also employed to verify the results. Finally the major issues that
limit the application of moving gird technique are discussed based on the work of the in-

cylinder ﬂow simulations for Generic 1C engine.

Fictitious domain methods, which can be term as one of the “combined Eulerian-
Lagrangian” methods, has been proposed as a way to reduce the complexity of meshing
for moving boundary problems. In this method, the original complex geometry associated
to the problem at hand is embedded into an extended simpler ﬁxed domain (the ﬁctitious
domain), the interface is tracked as an independent entities (either same order of
dimension or low order of dimension of the extended domain), the computations are
performed on the ﬁxed ﬁctitious domain mesh whose topology is independent of that of
the interface and the original boundary conditions on the interface are enforced in the
extended domain with the help of Lagrange multiplier [27-31] or by using a penalty

method [32]. Therefore the solution is transformed to a one-domain approach and the

meshing task is simpliﬁed into two independent geometries, the discretizations then can
be done separately and need not to be conform to each other. This obviates the need for

repetitive remeshing.

The use of Lagrange multiplier to enforce the boundary conditions however may
introduce errors [39]: the multipliers do not necessarily enforce the nodal boundary
conditions at nodal locations and may fall within elements. This can introduce signiﬁcant
errors on the adjacent nodes when the extended solution involves a discontinuous
derivative across the embedded boundary, e.g. enforcing a Dirichlet boundary condition
within the extended domain. The errors are due to an inappropriate approximation of the
local gradient at the location of a multiplier. A one-dimensional model problem is ﬁrst
used to study possible improvements. To reduce the errors, it is ﬁrst possible to reﬁne the
mesh. It is observed that a uniform reduction in the element size changes the relative
location of the Lagrange multiplier within an element and the errors are found to Oscillate
as the mesh is reﬁned uniformly for a given extended domain. The errors are also
functions of the relative locations of the Lagrange multipliers within the elements for a
ﬁxed element size. For ﬁxed relative Lagrange multiplier position, the Lagrange
multiplier/ﬁctitious domain methods (LM/FDM) result in a lower-order convergence rate
in L2 norm compared to a solution based on a body conforming mesh. The magnitudes
of the errors are related to the jump of the local solution gradient and the relative location

of the multiplier within the corresponding element. The errors are minimized as the jump

of the local gradients disappears. This is illustrated in chapter 4.

Finally in chapter 4, ﬁctitious constraint method and shape reconstruction method are
proposed to mitigate the errors: With ﬁctitious constraint methods, instead of satisfying
the original embedded Dirichlet boundary condition, a two-step predictor-corrector
procedure is adopted to adjust the local boundary conditions to values that reduce
signiﬁcantly the errors. These local values are called the ﬁctitious constraint and are
obtained based on the computed solution obtained from the standard Lagrange
multiplier/ﬁctitious domain method. With the shape reconstruction methods, the
boundary conditions are adjusted to be applied on the nearest adjacent nodes (edges) to
the Lagrange multiplier, the LM/FDM computations are then performed with the
reconstructed boundaries of the problem. Computation of the errors shows that the
solutions obtained with these methods are less sensitive to the jump of the local gradients
and the relative location of the Lagrange multipliers within the elements than the
solutions of standard Lagrange multiplier/Fictitious domain method. A higher-order

convergence ratio for L2 norm can be achieved when the relative Lagrange multiplier

position is ﬁxed. This is tested by solving simple one- and two-dimensional problems.

Chapter 2

An Introduction to Moving Grid Technique

In moving grid methods, the motion of the moving boundary is explicitly known at each
instant. a boundary ﬁtted mesh (as shown in ﬁgure 2.1) is generated at each time step to
explicitly track the interface’s evolution. The corresponding grid motion can then be
incorporated into the numerical scheme via the geometric conservation law [10]. When
the motion of the interfaces begin to distort the mesh, regridding (grid
redistribution/reﬁnement) is needed and consequently the solution obtained on the old

mesh needs to be interpolated onto the new mesh.

 

 

 

 

 

 

 

Figure l Discretization of the original problem (left) using a body ﬁtted mesh (right)
There are three important issues related to moving grid technique. First the governing
equations need to be transformed into a body-ﬁtted coordinate system. Second, the
geometric conservation law needs be satisﬁed to derive an equivalent differential relation
that can take the mesh movement into account. Finally, solution interpolations are needed

to transfer the data from the old grid to the new grid during a mesh regridding.

In this chapter, the fundamental ﬂuid dynamics conservation laws are ﬁrst presented for a
Newtonian ﬂuid. We then discuss these three issues: coordinate transformation using
generalized curvilinear coordinates (Thompson et al. 1985), the geometric conservation
law and interpolations of the solution.

To simplify the presentation, we will use 2D cased to highlight the issues involved, an
extension to 3D geometrical geometry can follow the same concept without qualitative

modiﬁcations.
2.1 The Governing Equations

The governing equations in Cartesian coordinates for two-dimensional, compressible

ﬂow can be written in dimensional form as:

 

Continuity:
_a_P_+ 600“) + 6C0") = 0 (21)
at 6): 6y .

x-momentum:

6,014 Bpuu apuv 6P 6 ( an) 6 Bu
=-—-+ —— — — — 2.2
at + 6x + 6y 6x [6x ﬂax +6y ﬂay +pgx ( )

 

y-momentum:

0/” 6W apw___6_’: 2. Q: 3 2
at + 8x + ay " ay+Iax[”axI+ay[”ay)I+pgy (2‘3)

 

For simplicity, we did not consider the source terms in the above equation.
2.2 Transformed Governing Equations in a Body-Fitted Coordinate System

Following W. Shyy’s work in [2,3], the coordinate transformations are discussed below:

In order to transform the governing equation to a moving, body-ﬁtted coordinate system
using generalized curvilinear coordinates, we introduce a time-dependent invertible

mapping transformation (for a two-dimensional space case):

x =X(6,7],T)
y = y(€.77,r) (2-4)
I: 2'

The time varying irregular physical domain is then mapped to a ﬁxed uniform
computational space. The ﬂow equations are then recast in a body-ﬁtted curvilinear

coordinate system (6,77) and are:

 

3(Jp)+a(pU)+a(pV)=0

 

 

6t 65 677 (2'5)

and

2 .6Wu) «ML- 6_P_ 6_P u _

atUpul 6: + a” {”75}; ygan}+a§[J(qlug qzuq)]+ (2.6)
a
5BR?” - (1214.13)] + pgx°J

2 WW) amp. (1’. 9: 1.41. _

64m” a: + an {Fan ””aéii'arid‘m‘ WWII“ (2.7)
6 p
‘3;[7(43V77 “ (12"; I] + Pg y 'J

where the subscripts 5 and 7] denote 6/66 and 6/677 , respectively; J is the Jacobian of
the coordinate transformation J = xgy” — xnyé: and represents the volume element in the

transformed coordinate.

and U and V are the contravariant velocity components

U =(u-5r)y,7 -(v—y)x,]

2.8
V=(v—y)x§-(u—5c)yé: ( )

where x§,x,7, yg, y” are the metrics of the coordinate transformation, and 5c, yare the

Cartesian components of the grid velocity vector deﬁned as:

n n—l n_ n—l
x:_x__._x—, y=Z_—L_

2.9
At At ( )

where the superscripts (n) and (n —1) denote the current time step and the previous time

step, respectively. The metrics ql , q2 , q3 are deﬁned as:

m=%+ﬁ

42 = xgxn +y§y77 (2.10)
2

‘13 = x; +y§

The contravariant velocity component and Cartesian velocity components at the boundary
are then computed to enforce mass conservation; and the kinematics condition can be
enforced at the interface (moving boundaries) via the conformed curvilinear coordinate
system.

After discretization, the physical domain 0 then is covered by a grid consisting of a
ﬁxed number of nodal points distributed over Q and along its boundary.

2.3 The Geometric Conservation Law (GCL)

For moving boundary problems, part of the boundaries of the physical domain moves in
time, the boundary-conforming character of the mapping transformation 2.4 implies a
corresponding motion of the grid points in physical domain. Ignoring this grid motion can
introduce errors in the computational ﬂow ﬁelds with ﬁnite-difference method [5]. This
issue can be addressed via an auxiliary equation derived via the concept of “geometric

conservation law” (GCL) as discussed below:

Let Aebe an arbitrary, ﬁxed element in the computational domain 06 enclosed by a
smooth boundary 6A , and let A(t) ={}|}=}(E,z)v5 e Ae}be the corresponding

element in the physical domain 0 under the time-dependent coordinate transformation
it: ﬂat), then the change in volume of A(t) equals the total ﬂux through the surface

6A(t) , which can be expressed as:
i j d; = Ia x_t-d6A(t) (2.11)
dt A(t) AG)

where E; is the element velocity).This is the integral form of GCL [5,6],satisfying the

differential GCL scheme that governs the special volume element under an arbitrary
mapping will eliminate numerical oscillations and instabilities for solutions on moving
grid [5].

In the differential scheme, grid motion and geometric conservation are handled in a
natural way through the contravariant velocities and the Jacobian evaluations [2,3], A
Jacobian transport equation is derived by considering a uniform density and velocity
ﬁeld,under a time dependent coordinate transformation. The following identity derived

from the mass continuity equation results:

6.1 6 . . 6 . .
E+E£(-xyn+yxﬂ)+$(-yxg+xy§)=0 (2.12)
Integrating the above equations using the same time integration scheme over the same

control volume used for mass conservation leads to the following equation:

(Jn _Jn-1)
T+(—jyn +jrx”)e —(—5cy,, +yxn)w+(-yx; +595)" -(—yx; +595)? :0 (2.13)

This formulation ensures that the Jacobian is updated in every computational mesh
according to the time dependent grid movements to guarantee geometric conservation in
the discrete form of the conservation laws.

An alternative approach is using ﬁnite-volume method to interpret the metric and
Jacobian, which establishes a direct physical relationship as following [9]:

If A§=An=1,then:

_ area of element in 2—D (2 14)
_ volume of element in 3—D -
and
N§=J§xT+J§xF+J§xi€
(2-15)

= area vector of the element face(3—D)

This implies that with ﬁnite volume method, geometric conservation law is automatically
satisﬁed.

2.4 Solution Interpolations

Moving grid techniques require grid rearrangement (regridding) when the boundary
movement causes the grid to be skewed and unevenly distributed. At these situations,
solution interpolations are necessary to transfer the data from the old grid to the new grid.
The basis strategy for solution transfer is search and interpolate [10]: Each new grid is
ﬁrst located within an element of the old mesh; the solution value is then interpolated
locally using element basis functions and element nodal values from the old grid.

For example, Assuming that the element ethat containing node p with coordinates

(x p: y p)has been identiﬁed, then the variable value a at p can be interpolated using the

10

nodal values and element basis functions in local reference element coordinates

(emeﬁ:
Ne ,.
“(xp’yp)= Z“;V’j(5p”7p) (2°16)

F1

where Ne is the number of nodes for element e, the local reference coordinates

(5p ,7] p ) can be determined using the global value (x p , y p) , the element basis function

61- and nodal values (x; i) by the solving the following equation:
Ne A
xp " Z xiV’j(5p”’p) =0
1:1 (2.17)
Ne A
yp ‘ Zlyiy’j (5pi’7p) = 0
J:

11

Chaper 3
The Application of Moving Grid Technique to the In-Cylinder Flow Simulation

3.1 In-Cylinder Flow Simulations of Internal Combustion Engine

Gas ﬂow in a cylinder of an internal combustion engine has a profound inﬂuence on the
performance of the engine [17,18,19]: the ﬂow into the cylinder through the inlet valve or
valves forms a impinging jet, and establishes organized motions in the cylinder (swirl -
about the cylinder axis- and tumble-orthogonal to the cylinder axis). These in-cylinder
ﬂow characteristics are key issues to ensure stratiﬁcation requirement (efficient
convective transport of fuel to the plug; high turbulence level to initiate combustion; a
stable and reproducible generated motion) and ﬁmdamental considerations for the
exhaust emissions of an IC engine. To design an environmentally ﬁiendly (low emission)
IC engine, which would also operate smoothly and produce high power with low fuel
consumption, an improved understanding of these thermo and ﬂuid dynamics of the in-

cylinder process is thus very important.

For decades, the investigation of the in-cylinder ﬂow patterns for IC engines have been
achieved mame by traditional methods based on extensive experiments: Engineers
developed new combustion systems by making variations in previously successful
conﬁgurations. Given the IC engine’s high state of reﬁnement and the physical
complexity of the in-cylinder processes, this cut-and-try process is not sufﬁcient to create
the signiﬁcant improvement now sought, plus the experiment set-up is time consuming

and expensive.

12

The advent of computers and the possibility of performing ”numerical” experiments
based on computational ﬂuid dynamics analysis provide a new way of the simulating and
designing internal combustion engines. By easily examining the various tradeoffs that
must be made to move current design toward the optimum, CFD enables the transient
analysis in the operation condition, which can yields detailed information on the ﬂow
ﬁeld in a timely and cost-effective manner, thus meets today’s performance targets

requirement.

However the three-dimensional ﬂow numerical simulation about a helical/direct intake
port-vertical/canted valve-cylinder system is remaining as a challenge to the CFD group

The main reasons are as following [18, 19, 21]:

o It is extraordinarily difﬁcult to generate computational grids of high quality due
to complicated port shape, chamber shape, valve position and the valve lift and

piston motion strategy.

0 The ﬂow ﬁeld in an internal combustion engine is turbulent and comprises many
time and length scales (from eddies that are essentially large enough to ﬁll the
available engine cylinder space down to eddies often substantially below a
millimeters in size). It is impossible with present techniques and facilities to
obtain a detailed numerical solution of the Navier-Stokes equations that can
account for all the in-cylinder turbulence time and length scales, and proper

turbulence models need to be introduced.

Fortunately, there are a number of commercial and in-house CFD codes have been

developed today, each of employing one or more numerical methods and techniques to

13

circumvent these difﬁculties. Among them, FIRE, Star-CD, FLUENT and KIVA-3V are

the most widely used codes in industries and academics.

3.2 Moving Grid Strategy in Current CFD Codes

Almost all these CFD software (FIRE, Star-CD, FLUENT,KIVA-3V) use the Lagrange
based moving grid technique to track the valve and piston motion for the engine
simulation [13,14, 15,16, 25], that is, the mesh attaches to the moving boundaries (valve
surface and piston head) and deforms with the computational domain, and is subjected to
the topology changes in the course of its evolution. There is three mesh motion strategies
to accommodate the volume deformation associated with these in-cylinder motions:

spring deforming, local regridding and layering.
3.2.1 Spring Deforming

In the spring-based deforming method, a block of mesh cells (single layer or multi-layers)
deforms like springs due to the movement of the mesh boundary: the edges between any
two mesh nodes are idealized as a network of interconnected a springs. A displacement at
a given boundary node will generate a force proportional to the displacement along all the
springs connected to the node, and so the displacement of the boundary node is
propagated through the whole block [26]. The displacement of each node can be
controlled by specifying different “spring coefﬁcient” of each “springs” to get a better
control over some speciﬁc region (especially the boundary layer regions). A simple
example of the spring-based deforming is shown in Figure 2 for a cylindrical volume

where one end of the cylinder is moving.

14

In the Spring Deforming Method, the connectivity of the mesh cells and nodes in the
corresponding block remains the same, the topology and the resolution of the grid system
do not change, only the nodes position may change due to the boundary motion. It is well
suited for the situation where the motion of the interested region (piston head or valve
surface) usually keeps the same (or reverse) direction and the displacement is not large

enough to cause some ill-shaped (such as high aspect ratio) cells.

Spring deforming allows the use of a low grid resolution during the early time steps thus
can speed the analysis through the initial transient behaviors. When the ﬂow variables
settle, layering (as discussed in 3.2.3) to a higher grid resolution can obtain the necessary

accuracy.

Since this method always involves grid motions, it requires the Jacobin to be upgraded at

each time step to satisfy the GCL.

IIIIIII
I IIIIII IIIIIIIIII IIIIII
‘ IIIIIIIIIIIIIIIIIIIIIIIIIII

 

Figure 2 Spring deforming on a simple cylinder at the end of compression (Left: piston
head at BDC, right: piston head at TDC)

3.2.2 Local Regridding

When the boundary displacement is large, the moving boundaries can distort the cells too
much (especially for the gap region between valve surface and seat), which will create
bad quality cells (e.g. high aspect ratios, twisted faces, negative volumes, negative
volumes, etc.), which lead to a high degradation of a convergence rate and a poor
accuracy as the solution is advanced to the next time step. Local regridding is thus
necessary to adjust the mesh to the level with acceptable quality. This can be done either
by changing the local mesh topology (ﬁgure) or using smoothing algorithms and others
similar techniques to redistribute the grid nodes without changing the connectivity
information. An interpolation is thus needed to maps the corresponding boundary
conditions and ﬂow variable data that is generated during the analysis from the old grid
to the new grid between two adjacent time steps. This method may lead to the change of
connectivity and grid resolution as topology adjustment is employed. Figure 3 gives such

an example.

Local regidding method is widely used in the regions (such as the upper portion of the
combustion chamber and port region around valve seat) that involve the valve movement
and extensively meshing effort is needed to match the geometries at the position where
the combustion deck is complex or a canted valve is involved in the motion. However,
the application of hybrid mesh techniques (such as hexahedron elements couples with
tetrahedral elements or pyramid element in FLUENT and Star-CD) can simplify this

process.

16

 

 

Figure 3 Local regridding applied to the valve region during intake stroke (Left: valve
opening, right: valve closing)

3.2.3 Layering
When the boundary displacement is large, by simply adding or removing a cell layer in

the boundary motion direction, the grid resolution can be adjusted to obtain a block of
cells (single layer or multi-layers) with good aspect ratio and acceptable layer-thickness
(which is very import for the boundary layer regions). The layering method is particularly
suited in the cylinder region that is directly above the piston head, and can also be

utilized in the regions above the valve, including the valve seat region.

Layering changes the resolution of the grid system and it can maintain the topology of
grid system in the direction where there is no motion involved. An interpolation is also
needed to map the data between different data set at different time steps. Consequently,
the application of layering to a single layer block or multi-layers block determines the

level of interpolation. A simple layering example is illustrated in Figure 4.

17

 

Figure 4 Layering applied to the valve region during intake stroke (Left: valve closing,
right : valve opening)

The spring deforming technique and local regridding technique can be viewed as grid
redistribution techniques while layering technique can be viewed as grid reﬁnement
method. The reason that we term them separately is that each technique is employed at

speciﬁc situations for engine simulation.

The combination of these three techniques coupled with other meshing techniques (block
structure mesh, unstructured mesh, hybrid mesh, automatic mesh generation) to deal with
the meshing task in different portion of the computational volume enables the users to
generate an acceptable computational mesh. The meshing time usually varies from 2-day
to 2 weeks depending on the grid generator’s automatic level and user’s familiarity with

the codes.
3.3 Turbulence modeling

For low speed (laminar) ﬂows without heat transfer, the equations governing the
conservation of mass and momentum can be used to describe the ﬂow exactly for

incompressible ﬂows. Turbulence however, leads to rapid velocity ﬂuctuations in both

18

space and time. So although these equations can properly describe the details of turbulent
motions, it’s too costly and often time consuming to obtain a solution with detailed
information about both the time and space variations of ﬂow variables. This leads to the

concept of“ turbulence modeling” [7, 8].

Most engineering models of turbulent ﬂow are based on the use of the Reynold’s
averaging technique [7, 8], which assumes that the quantities (that appears in the N-S
equations) at a given point in space and time are described as a superposition of some
mean part, which may vary slowly with time and a random component, which varies

rapidly. Mathematically, they can be expressed as:

at =17,“ +11,"
p=ﬁ+p'

with u—,'- =,—o-' =27 = 0
Here the quantities with bar denote the mean values and those with primes are

ﬂuctuations.

Then the ensemble-averaged mass and momentum balance equations are given as (for

incompressible ﬂuids):

0 The ensemble-averaged continuity equation:

ga+%[m+ﬁ]+-§y—[p—v+ﬂ]=o (3.2)

o The ensemble -averaged momentum equation:

19

 

 

 

—— “1—2—
Du 2__ +6u'v' +'6uw
— = —— 6—+ Vu
p Dt 'u 2vp6u[ ax 62
57 pu—_6 u"v +6v—'2 +'6vw
— - —— :y—+ ,uV (3.3)
Dt 6x ++6y 62
D»? p 2 6a 6v'w' 6w'2
p— - -— + #V W p
Dt a 6): 6y az

6_6 66

where—=r7— +v—+w—+—
Dt

ax+ay azat

The additional unknown terms in the momentum equations are know as the Reynolds

stress tensor expressed as:

f __

 

 

 

 

au'2 5?? am
0x 6y 62

—,7 .72- 'ﬁ

Reynolds stress tensor = — p 6u v 6v 6" W
5x 6y 62

8“,“) a;— aw'z
6x 62

t 6” ,

 

(3.4)

The modeling of Reynolds stress tensor with other variables introduces different

turbulence models (standard k—r: model, RNG model, k—co model..

.), among them

standard k -3 model is widely used in industry due to its numerical robustness and

economy.

In this method, the Reynolds stress tensor terms can be expressed in terms of the mean

rate of strain S j and turbulence viscosity pt :

 

, r 2
-pui uj =2ﬂtSzj-3P5g‘k

(3.5)

with :

S" -1 §£+ﬁ ' -C 53 and 6" is the Kroneche delta function
U26xj aye-”W ”ps’ ’1 '

The turbulence kinetic energy k and its dissipation s can be solved from the following

equation:
———= —£ +— +— —
th PU} ) 6xj [[ﬂ 07‘]ij
_ (3.6)
D8 614k 8 a [It 65
—= C +C k—- a —+— +— —-
th p[ 511’): .93 M £2 )1: axj[[# 01(13ij

where the production of turbulence energy is:

 

I r a ' . .
I), =—u,~ u j 319— and C51, C52, C63, 0k are empirical constants.
x o

It should also be noted that the application of k — e to engine simulations still has serious
limitations [18]: the model assumes that the turbulent transport is in the same direction as
the mean ﬂow gradients and does not consider the effects of pressure-velocity
correlations on the k equation and 3 equations. And also the .9 equation does not have a

strong physical foundation.

3.4 Description of the Problem
A traditional direct port diesel engine (generic engine) was employed in the study; the

geometric shape and engine parameters are listed in Figure 5 and Table l:

21

 

Figure 5 Geometry of the generic engine and the boundary deﬁnition

Table 1 Engine Speciﬁcations and Calculation Conditions

 

Bore x Stroke

80mm x 100mm

 

 

 

 

 

 

 

 

Squish 5.0 mm
Engine RPM 1500 Rpm
Compression Ratio 19.3
Connecting Rod 198 mm
Maximum Intake Valve Lift 8.68 mm
. 353.3 ATDC
Intake Valve Opening (after top dead center]
Intake Valve Closure 545.6 ATDC
Pstatic =1 an",

Boundary Conditions

Inﬂow Boundary

1:91am = 300K
,u, = Cllpk2 ls =100,u

I=J2k/3/U=0.l

 

 

 

U = f (Rpm)
Wall Boundary adiabatic wall, no slip
Symmetry Boundary Symmetry

 

 

Initial Conditions

 

u=v=w=0

Pstatic =1 atm,
Tstatic = 300K

p, =Cﬂpk2/£=100a

I=J2k/3/U=O.l

U = f (Rpm)

 

 

22

We assume a symmetry boundary across the engine symmetry plane thus there is no swirl
generation during the simulation and we can use half engine simulation to reduce the
computational time.

3.5 Numerical Setups of the Simulation

FIRE, a ﬁnite volume code that solves the ensemble averaged fully compressible
conservation equations for mass, momentum, and energy, was used to perform the
simulations. The Reynolds stresses are linked with the mean ensemble averaged
properties with the standard k —.9 model. The standard “ wall function” is used to bridge

the viscosity-affected region between the wall and the fully turbulent region.
The k — a model constants were left to their default settings as listed in Table :

Table 2: Values of standard. k — a model constants

 

Cy Cal Caz C53 O'k

 

 

 

 

 

 

0.09 1.44 1.92 -0.373 1

 

The discretisation equations are solved by iteration following the SIMPLE velocity-

pressure coupling algorithm.

The second order central total variation-diminishing scheme (CTVD) with minmod
limiter was used for space differencing of all the variables. The ﬁrst order fully implicit

scheme is adopted for the temporal differencing of all the variables [14].

The numerical model was made based on the geometry CAD data (surface mesh) from
another commercial code- Star-CD. The using of arbitrary multi-layer based spring
deforming and layering technique were used to handle the valve and piston motion, and

arbitrary local regriddings were employed to handle the valve seat region when the mesh

23

 

screwed during the valve motion. The resulting numerical models are shown in Figure 6

Table 3shows a summary of the model data.

 

Figure 6 Computational Mesh of the generic engine

Table 3 Numerical model data

 

 

 

 

 

 

 

Total number of cells (maximum) 404,184
Number of cells in port and valve (maximum) 99,504
Number of cells in cylinder (maximum) 304,680

 

 

3.6 Experiment Veriﬁcations

A validation based on the experiment data (provided by Dr. Shock) had been performed

to check the code quality.

The average in-cylinder pressure shows good agreement with the experimental data; the

difference is below 5%. (Figure 7)

24

45KB

 

Pavg
i

IIrlTTrTTIIII VIII

 

 

 

P A k o
q 1"r—j—J'ﬁg—rT'L—I—rﬁ'T'T'I—l-F—I‘T'T—r—L gLJ I r r r r I

420 4&1 540 60) $0 72)

0

Figure 7 Average in-cylinder pressure comparision: dash line-numerical data, solid line-

experimental data
The ﬂow vector distribution at the symmetry plan are also compared (Figure 8-11, due to
the limitation of the measurement, the experimental data only reﬂect a portion ﬂow
distribution in the whole ﬂow ﬁeld): the agreements of numerical results with
experimental result at intake stroke (CA 450, 540 ATDC (after top dead center» are
satisfactory, but clear differences remain at compression stroke (CA 630, 660 ATDC):
the ﬂow pattern is almost completely wrong. This behavior may be caused by the
inaccuracies or errors in either the numerical modeling or the experiments; currently there

is no solid explanation for this behavior of the code.

25

 

 

    

 

 

 

° ° '- -- "fry-$2.:
' '1 5’5". :1... J, I l.‘ ‘5
-1 -1 It; [5],; 5:2: :.1’-/' '- 3%“-
// I \ \ \ ' T I57,” :95: "ﬁrst? \“I‘ I'III .. I ‘1’:
‘2 \ l , ‘2 “’0’, I’oir.‘ -~:,\‘ II IIHII'KI: {$8.
/ j \ \ \ I fab-1:: "'2'” [IIIIIIIIIIII IL,“ \ \\\\:\\‘.
.3 \ \ J 1‘ ‘ ‘\ \‘C :- I’M, ‘3'“ II" [IIW‘ 1I“\ 0‘
I I,Ik~..I,I I] II I \\
4 / \ \ 4 |‘\ I‘ “ ‘\‘\\“‘- I". '9‘ ”M1 ”‘1“ [III ‘\:‘\\\\\\ :\\\\\‘|
l I ‘ ‘\ \‘ ‘\\\\‘\~ o I ["1” I” “ \\\ \\\ \\'
a >. i,\ ~_.,I‘ 'I11‘\\\\\ \\\1
/ I i I I I “‘~\‘-’IIII III M“ \‘\\I
'5 '5 I “ \ ‘s " I“ ‘1‘ \\\\\\ \\\\ [I
/ ’ I l I I I I -“.‘:\“ ‘ i’xii\ \\ ‘I‘ \\\\‘I
'.. . §\ \\ \§\ \\\\I
4 I I I l 1 I j I «a - . . g ‘ .,. ~~§Jt~ .\-

\ ‘ ‘ I l I l
.7 \ ,7
-8
'00 1 2 3 4 5 6 7 8
x

 

 

 

 

Figure 8 Flow pattern comparison at CA 450 ATDC (left-experimental data , right-

numerical data)

 

 

 

 

 

 

 

 

Figure 9 Flow pattern comparision at CA 540 ATDC (left-experimental data, right-

numerical data)

26

 

-‘N

. “ ‘ ~ ' .3", . s s ‘ s -=: _ '.
\\\\\ ‘1. \\\ 17”" SC»? 1"."
" “ -‘t ‘
IIIII‘IIII‘II “'\ ‘ "

- - . \ \ __
II I, “If "' r \ >-:- I: '. .
-= ~ ‘ ‘ ‘ ‘ ‘ -2 I II I "II 'm' -.
- \ \ \ \ \ ‘6‘ ‘II Iﬂ 1;,” :"I’ ”III”? - . ‘ . _.:’//,’ 1”], r.
\I1IIIIIIM’I ‘-, 01’:-
'3 \ \ \ \ \ ‘3 \I I II "I", "III, ’ '- , ’I, ’, r
I. ‘I “I \I “I" ,"r’, ’I’,’ III/,3?” I ' a 1":I”'," I" I '1
4 ‘ ‘ ‘ ‘ \ \ \ .4 \\ ‘I ‘I "II I If], If,” [’2’ I o ’I ' “'1' I:'
>- t I I \ \ \ \ \ >~ \“\\\\\‘\\II‘ III” :I], W,” ”iii”; ' :I',"':"' "II" '1 1'1
'5 . . I l 1 I I l '5 1‘0 whiz/I III III,
\\\ \II 1III III] I
.6 I r I I I I I I -o \\\\ \\§\‘\\\\\ \‘\\ lIII I/éhiIip‘ I“‘\ ““5311!
I / / / / / / /

 

 

 

 

 

 

 

Figure 10 Flow pattern comparision at CA 630 ATDC (left-experimental data, right-

numerical data)

 

0 .. .

a». -2. ":QCZE'

1 “an“. : ’I‘. 1 S -..::3r.,:’,".
- III "'1, 91-! -/ /, I,

31* :1th “I311", " 1:95,; 3:3-" . /,’/ 21,715" I

:‘i IIIIII’II//;%/w 'IR‘IIM III III:

 

   

 

 

 

 

Figure 11 Flow pattern comparision at CA 660 ATDC (left-experimental data, right-

numerical data)

3.7 Simulation Results Discussion

Simulation results are presented here for ﬂow analysis. Three plans are chosen here ( two
vertical plans: x-plan and y-plan and one horizontal plan: z—plan) to reveal the ﬂow

pattern. These plans are deﬁned in Figure 12: The x-plan cuts through the axis of the two

27

intake valves and normal to the horizontal plan. The y-plan cuts through the axis of the
intake and outlet valves and normal to the horizontal plan. The z-plan is a horizontal plan

and cuts through the middle of the cylinder.

 

gig i a;

y_—pl_an§: through

both valve centers

3:21am: through the

intake valve center

 

2421mm: at the
cylinder mid-height

 

Figure 12 Selected plans for the generic engine

Figures 13-34 show the ﬂow ﬁeld in the form of velocity vector distribution at different
degree ATDC along the vertical plans (Figuresl3-32) and horizontal plan (Figures 33,
34).

At the beginning of the intake stroke (Figure 13,14, CA 370,380 ATDC), with the valve
opening and piston moving down, the ﬂow is induced into the cylinder and forms an
impinging jet around the valve seat region. As the valve liﬁ and piston speed are low, the
turbulence intensity level is not strong, the ﬂow remains attached to the valve head and

seat and then separates ﬁom the valve at the inner edge of the valve seat, part of the flow

28

will hit the cylinder wall, this results in the formation of obvious recircuilation motions
under the cylinder head and valve (Figure 14, CA 380 ATDC): Two couples of opposite
tumble motions are created-one couple that formed at the position between left side of
valve seat and the top-left cylinder corner and one formed in the opposite direction. Later
on we can see that the tumbles formed under the valve head have dominant inﬂuence on

the ﬂow pattern.

As the piston accelerates downside and the valve lift increases, (Figure 15-22, CA 390-
460 ATDC), much of the directed energy in the jet is converted into turbulence, the
turbulence level increases. The tumble in the right upside of cylinder corner can not ﬁnd
enough space to grow and ﬁnally breaks up. The other three tumbles motions are
strengthened continuously as the distances from cylinder wall are far enough for them to
grow (e. g., the tumble at the left corner can shift it position downside as the piston moves
down). The structures of the two tumbles motions under the valve also changes as the
turbulence intensity evolves. The left side one (counterclockwise) continues to grow and
ﬁnally dominates the upper side in-cylinder ﬂow ﬁeld while the right side one does not
grow and ﬁnally looses most of its energy (break up in to turbulence). Another tumble
motion whose rotation direction is opposite the one under the valve is also created above
the piston head at around CA 450 ATDC (Figure 21,22, CA 450, 460 ATDC) and
dominates the lower part of the in-cylinder ﬂow. At CA 460 ATDC, the turbulence level

is strongest as the engine approaches its maximum valve lift and piston speed.

During the second half of the inlet stroke (Figure 23-25, CA 480-540 ATDC), as the
piston slows down and the valve is going to close (which means the jet is coming to an

end with fewer directed energy introduced), most of the tumble motions decays: the

29

intensity of the tumble that above the top of piston decreases markedly and ﬁnally
attenuates into two small tumbles around the center position of the cylinder (Figure 25

CA 540 ATDC). Viscosity will also have some effects on these phenomena.

During the compression stroke (Figure 26-32, CA 550-700 ATDC), the tumble motions
are ampliﬁed due to the increase of density and the changes in length scales (as engine
compressed, the charged geometry is changed). The two tumbles located at the center of
cylinder combine into one tumble motion that rotates in clockwise direction, and try to
move upside to ﬁnd more space to grow, the intensity of the upper tumble also increases
and it shifts its location to the center position right below the valve head to get more

tumble axis-cylinder distance.

As the piston approaching TDC (Figure 31,32, CA 680,700 ATDC), the vortexes created
by the tumble cannot ﬁnd sufﬁcient room to maintain their form. They ﬁrst combine into
one vortex at CA 680 ATDC and then break up into turbulence at CA 700 ATDC. A
crudely homogeneous condition can be observed at that time. The break up of the swirl
motions (strictly speaking, these motions are recirculations in the cylinder radius
direction since we assume a symmetry boundary condition across the cylinder symmetry
plane) in the compression stroke can also be observed (Figure 33, from CA 550 ATDC to
CA 660 ATDC). The breaking up of the tumble motions at CA 680 ATDC will also
result in a changing over to small swirl motions (Figure 34 CA 680 ATDC to CA 700

ATDC).

3O

 

 

 

Figure 15 Velocity distributions at 390 ATDC (x-plane and y-plane)

31

 

Figure 16 Velocity distributions at 400 ATDC (x-plane and y-plane)

 

Figure 18 Velocity distributions at 420 ATDC (x-plane and y-plane)

32

 

Figure 20 Velocity distributions at 440 ATDC (x-plane and y-plane)

 

Figure 21 Velocity distributions at 450 ATDC (x-plane and y-plane)

33

 

 
 

  
 
 

“I;IIIIIIIIINUIHWMF ‘

 

Figure 24 Velocity distributions at 510 ATDC (x-plane and y-plane)

34

 

Figure 27 Velocity distributions at 580 ATDC (x-plane and y-plane)

35

 

 

Figure 30 Velocity distributions at 660 ATDC (x-plane and y-plane)

36

 

 

 

 

 

 

 

 

 

 

 

Figure 31 Velocity distributions at 680 ATDC (x-plane and y-plane)

 

 

 

 

 

 

 

 

 

 

Figure 32 Velocity distributions at 700 ATDC (x-plane and y-plane)

 

 

 

distributions at 550(1eﬁ) and 660 (right) ATDC (z—plane)

Figure 33 Veloci

37

 

 

 

Figure 34 Velocity distributions at 680(left) and 700 (right) ATDC (z-plane)

3.8 Issues with Moving Grid Technique for In-cylinder Engine Simulation

There are three issues closely related to the moving grid technique: interpolation errors,
human cost for mesh generation and the quality of grid system. Each of them will affect
the numerical modeling process and the corresponding solution.

The interpolation error is introduced by the solution mapping process: once there is a
need for local regridding or layering, the solution on old grid needs to be transformed into
the new grid. As we have seen in section 2.4, the interpolated solution on the new grid is
obtained from the local approximations using the nodal values of the corresponding
element on the old grid, since the accuracy of the local approximation is highly depended
on the chose of local basis function, it may not be a good representative of the actual
solution and thus introduce errors. Figure 35 gives a simple one dimensional example of
such interpolation: The numerical solution at the nodes on “old” grid agrees very well

with the exact solution, while the interpolated one on the “new” grid lose accuracy.

38

 

— Exact Solution
.6— Numerical Solution at Old Mesh
0.25 - *— lnterpolated Solution at New Mesh

 

 

  

 

 

1 . ____L

"12‘ 1.3 1i4’”1f5x1i6 1.7 1.8 1.9 2

Figure 35 Solution interpolation from “old” grid to “new” grid introduces error

The second issue is related to the generation of computational grid to handle the
valve/piston motion in the context of complicated engine geometry. As we have
discussed earlier, moving grid technique requires the grid system to conform to the real
engine geometry at the corresponding time step as well as preserve suitable topology to
handle the valve/piston motion, it is usually done with a multiblock unstructured grid
system and the quality of the grid is very difficult to control: the grid system may be
skewed and unevenly distributed, especially in the region related to the valve motion. The
bad quality grid will introduce error to the numerical solution and the error will transport
and polluted others regions. Figure 36-38 give the example of a bad grid system and the
corresponding errors are obtained based on the comparison of the numerical solution for
a low-resolution grid (50k) with the one for a high-resolution grid (500k, we take it as the
exact solution here)-Figure 36 shows the gird system for a 50K engine at CA 460 ATDC
at selected x,y 2 plans and the corresponding velocity point-wise relative error and

absolute error contours were draw in Figure 37 and 38 respectively.

39

  

i

ii nmmmu luumumm» ,
ll llllllllllllll lllllllllllll Ill
llllllllllllllll lllllllllllll Ill
Illlllllllllllll llllllllllllllll

l|| Illlllllllll‘ ' llllllllllllllll‘
111 Illllllllllll mmnmumn

    

iii. "
’lll Illlllllllllillllll‘
.lll Ill ll llllllllllll
.Il Illl llllllllllllllll
5‘ nlt ll llllllllllll

    

    

      
     
 

Figure 36 Computational mesh at selected plans for a 50k generic engine numerical
model

   

Figure 37 Velocity absolute error distributions at selected plans for a 50k generic engine
numerical model

   

Figure 38 Velocity relative error distributions at selected plans for a 50k generic engine
numerical model

And also because such a limitation (geometrical conﬁrmation) exists, we ﬁnd that the

solution is highly sensitive to the grid topology and regridding strategy that were used in

40

the numerical model-to illustrate this, we employed three methods here to generate a 50K

grid for the generic engine, each grid uses a different mesh strategy. In method I and II

the grid topology in horizontal plan is exactly the same while the moving mesh part for

the near-valve region are generated using different remeshing strategy (the time and

topology for regridding are different); Method 11 and III use the same remeshing strategy

for the near-valve region and differs from each other in the horizontal plan (different

topology). Figure 39-41 shows the grids’ difference of these three methods and Figure

42-47 show the different of the streamline contours for these methods (at CA460 and CA

630 ATDC).

111qu
ll [Mill

 

11111111111111

 

   

      
  
 

lllllllll IIIIIIIIIIMI
Ill lllll Illlllllllllll
m um Illllllllllll
IllllllII Illlllllllllll
Illllllll Illlllllllllll
111qu mnmmm
IIIIIIII Illllll

    
  
  

    

Figure 39 Computational mesh at selected plans for a 50k generic engine numerical

model (Method I)

 

illlllllllllllll Ill l

 

Figure 40 Computational mesh at selected plans for a 50k generic engine numerical

model (Method 11)

 

Figure 41 Computational mesh at selected plans for a 50k generic engine numerical
model (Method 111)

     

Figure 42 Streamline contours comparison (leﬁ: method 1, middle: method 11, right :
method IH) at selected plan (x-plan) for a 50k generic engine numerical model at CA 460
ATDC

/

”ll

Figure 43 Streamline contours comparison (leﬁ: method 1, middle: method 11, right:

) a
method H1) at selected plan (y-plan) for a 50k generic engine numerical model at CA 460

ATDC

 

42

  
    
  
  

, ’ 4; '31,
/All./‘S 9)wa

Figure 44 Streamline contours comparison (leﬁ: method I, middle: method 11, right:

method H1) at selected plan (z-plan) for a 50k generic engine numerical model at CA 460

   

 

 

ATDC

 

Figure 45 Streamline contours comparison (left: method 1, middle: method II, right:
method 111) at selected plan (x-plan) for a 50k generic engine numerical model at CA 630

ATDC

 

   

Figure 46 Streamline contours comparison (left: method I, middle: method II, right:
method III) at selected plan (y-plan) for a 50k generic engine numerical model at CA 630

ATDC

43

 

Figure 47 Streamline contours comparison (left: method 1, middle: method H, right:
method 111) at selected plan (z-plan) for a 50k generic engine numerical model at CA 630
ATDC

We also compare the predicted turbulence kinetic energy for these three methods since
turbulence kinetic energy is also one of the most important issue for engine simulation:
an overproduction of turbulence kinetic energy changes spreading rate of the jet and
caused different mixing in the cylinder. The averaged turbulence kinetics energy

comparison is plotted in Figure 48, and the detailed TKE distributions at CA 460 and 630

are compared in Figure 49-52.

 

 

 

 

Figure 48 Average predicted in-cylinder turbulence kinetic energy comparison for three

moving mesh strategies

 

Figure 49 Turbulence kinetic energy distribution comparison (left: method 1, middle:
method H, right: method IH) at selected plan (y-plan) for a 50k generic engine numerical
model at CA 460 ATDC

   

Figure 50 Turbulence kinetic energy distribution comparison (leﬁ: method 1, middle:
method H, right: method HI) at selected plan (z-plan) for a 50k generic engine numerical
model at CA 460 ATDC

 

Figure 51 Turbulence kinetic energy distribution comparison (left: method 1, middle:
method H, right: method HI) at selected plan (y-plan) for a 50k generic engine numerical
model at CA 630 ATDC

45

 

Figure 52 Turbulence kinetic energy distribution comparison (left: method 1, middle:
method H, right: method H1) at selected plan (z-plan) for a 50k generic engine numerical
model at CA 630 ATDC

This tremendous difference in the prediction of turbulent kinetic energy will surely affect
main engine parameters during a design process.

The third issue is well known: moving grid technique has a strong requirement for
repetitive remeshing when intensive interface (moving boundaries) activity involved,
which is time consuming and costly. This bad feature limits its application in engine
model evaluation since every geometrical (port, valve, combustion chamber)

configuration and possible liﬁ strategies, we need to generate a separated set of grids to

do the simulation, and such a human cost is huge.

46

Chapter 4

An Introduction to Lagrange Multiplier/Fictitious Domain Methods (LM/FDM),

Methods for Reducing Errors in LM/FDM Applied to Heat Transfer Problems

4.1 Introduction

Fictitious domain methods offer the possibility of signiﬁcantly reducing the difﬁculty of
meshing complex domains. This is achieved extending a complex domain into a much
simpler one that can be discretized easily (the extended domain is called the ﬁctitious
domain). The original problem is preserved in the extended domain by enforcing the
original boundary conditions with the help of Lagrange multiplier [27-31] or by using a
penalty method [32].

The ﬁctitious domain method is especially useful when dealing with moving boundary
problems. Two independent meshes can be then used: one ﬁxed Cartesian mesh and an
independent boundary mesh that “follows” the boundary and enforces the boundary
condition. This obviates the need for repetitive remeshing. R.Glowinski et. al have
developed recently various approaches based on the use of Lagrange multiplier
distributed on the boundary of the moving object [27-29] and Lagrange multiplier
distributed over the entire domain associated to a moving object [30, 31, 33-35].

Given its easy-implementation capability that which enables a fast turn around time in
engineering modeling, Lagrange Multiplier/Fictitious domain methods can come to be a
powerful technique for the design evaluation and engineering analysis if it accuracy can
also be guaranteed. The use of multipliers to enforce an essential boundary condition

however may result in a discontinuity in the ﬁrst derivative of the solution (e.g.

47

0.)!

temperature or velocity ﬁelds). This can introduce signiﬁcant errors if the multipliers are

not coincident to nodes of the underlying mesh [39]. This situation arises since most

elements commonly used are C 0 continuous i.e. a discontinuity of the ﬁrst derivative is

allowed only at the boundaries of the elements. Errors are thus introduced on the adjacent

nodes due to the impossibility of obtaining a CO solution within an element. The use of
higher order polynomials as shape functions somewhat reduces the errors but such errors
are nonetheless still considerable compared to the use of a mesh ﬁtted to the boundary
(see below in section 4.3 and 4.4)..

In this chapter, the Lagrange multiplier/ﬁctitious domain method and its use in problems
with Dirichlet boundary conditions is discussed. A parametric study for the introduced
errors is conducted. Two possible approaches, Fictitious Constraint methods and Shape
Reconstruction method, are proposed to mitigate the errors associated to the use of
Lagrange multiplier with the corresponding element. For simplicity, we only consider
one-dimensional and two-dimensional steady state conduction heat transfer problems and
discuss the inﬂuence of the jump of local solution gradient and the presence of the

multiplier within an element.
4.2 Lagrange Multiplier/Fictitious Domain Methods and the Dirichlet Problem

4.2.1 A Model Problem

The following Dirichlet boundary value problem is ﬁrst considered:

48

I

Given feH_1(Q), T0€L2(}’),T1EL2(7),
ﬁnd a function T such that

—AT=f in (Ma), r (4.1)
T=T0,on F;

T=T1,on 7;

 

 

 

 

F

 

Figure 53 Illustration of the geometry associated to the Dirichlet problem presented in
(4.1)

In this problem, (2 is a bounded domain in Rd (d22) including an inclusion a) ,I‘ and y

are the boundaries of Q and a).
4.2.2 Lagrange Multiplier/Fictitious Domain Formulation
The Lagrange multiplier/ﬁctitious domain formulation consists of:

I Extending the domain Q\a) to the larger square domain (2 , the original
boundary condition on the 7 becomes an embedded constraint for the extended

problem
- Extending fin H1(Q\w) to fin 111(0), and Tin H1(Q\a)) to f in

H1 (a).

49

“I!

.{..i_.

I Introducing Lagrange multipliers to enforce the embedded boundary condition on
7 and making the augmented functional stationary.
The following equivalent Lagrange multiplier/ﬁctitious domain formulation is thus

obtained

Find TeH1(Q), 461.2(7) such that
a(T, v)+Iy/1vdy= Llfvdﬂ \7’ veH(l)(Q)
Lindy: Irfipdy V ,ueL2(7) (4-2)

~

T=T0 on F

 

where a(T, v)s bVT-Vv d0

4.2.3 Finite Element Approximations

The approximation spaces for the discrete variables Th ,1}, , are chosen as following:

Th ={filfi e[C°(a>]d ,T‘ilE elm" ’VEESh}
(4.3)

Ah ={AhI/lhlE7 =const.,VE7 63;}

Where d=1,2, or 3. 3h is a regular discretization of the ﬁctitious domain (2 , 3; is the

discretization of the embedded boundary 7 , E and E 7 are the corresponding ﬁctitious
domain element/ embedded boundary element. 7?}; is assumed globally CO continuous

and locally (within each ﬁctitious domain element) Cl continuous.

Problem (4.2) can be posed as:

50

Find YieTh, ,1}, EA}, such that:
61(2),, v)+Iy/1},vd7= £1]de V VET},
Lﬁﬂd7= Irder V #EAh

”T240 on r

(4.4)

 

J
R.Glowinski in 1994 [27] made the remark that the spaces T}, and A}, can be chosen to
be independent and suggested to deﬁne A}, from the intrinsic geometrical properties of
7. This approach enables the uses of non-matching mesh sets. For problem (4.4), the

discretization can be chosen as shown in Figure 54 i.e. two meshes are used, one
underlying mesh and one to capture the boundary on which to enforce the boundary

conditions.

   

Figure 54. Discretization mesh set (left), discretization of ﬁctitious domain 9 (middle)
and embedded boundary y (right)

This treatment provides substantial simpliﬁcations to the meshing task, and is particularly
well suited to the situations where w is subjected to a rigid body motion [27,28]. Both R.
Glowinski et. al [27-31, 33-35] and Betrand et. a1 [32].employed this approach to
construct the approximation spaces and simulate the flow around moving rigid bodies of

various geometries.

51

4.3 Effects of the Presence of a Multiplier Within an Element

The use of non-matching meshes to solve a problem with an embedded boundary
constraint often results in the Lagrange multipliers that are not located at the nodes of the
underlying Cartesian mesh. When the extended solution of the problem is discontinuous
(with respect to the ﬁrst derivative of the variables) across the embedded boundary, the
original governing equations (Eq. 4.1) deﬁned in (Na) are not valid for the extended
domain 0 even if the formulation of Eq. 4.2 still holds.

In a ﬁnite element approximation, the extended domain (2 is discretised into a collection
of preselected ﬁnite elements, the solution over each element then can be approximated
by a set of approximation functions derived from the interpolation theory [37]. These
approximation functions are often algebraic polynomials which are smooth and process
no discontinuities, which means when applied with LM/FDM, the solutions are still
assumed to be smooth in the interfacial regions where the corresponding Lagrange
multipliers is located, the approximation functions then can not track such a jump of the
local gradient and thus introduces local errors in the numerical solution (as the adjacent
nodes are adjusted to enforce the Dirichlet boundary conditions within the corresponding
element).

Figure 55 illustrates a LM/FDM implementation for one-dimensional steady state heat
conduction problem with non-matching mesh set. Obviously, the linear approximation of

the unknowns cannot capture a discontinuity of the ﬁrst derivative within an element.

52

   

— Exact Sch-lion

 

   

7 Eucr Solurion ‘

-- LM/FDM Solution (linen)

 

Figure 55 A linear approximation of the solution enforces the embedded boundary
constraint but cannot capture a discontinuity in the ﬁrst derivative at other location but
the nodes. The right ﬁgure is an enlargement of the solution near the location of the
multiplier.

In Figure 56, another example for a two dimensional case is provided. In this ﬁgure, the

absolute and relative errors are shown and it can be seen that the errors are very high on

the nodes surrounding the collocation points corresponding to the Lagrange multipliers.

 

Figure 56 Absolute error (Leﬁ) /Relative Error (right) introduced by the LM/FDM for 2-
D steady state heat conduction problem. Eabs and Ere] is deﬁned as the point wise

errors with Eabs =|T —T},|,E,.e} =|T—T},I/ T , the LM/FDM solution is obtained using

gauss quadrature integration.

4.4 Parametric Study on the Sources of the Errors

Second order problems with single unknown and constant coefﬁcients should have
almost no error between the exact solution and the body ﬁtted ﬁnite element solution at

the nodes [37]. Clearly this is not the case for the Lagrange multiplier/ﬁctitious domain

53

methods as observed from ﬁgure 55 and 56. The ﬁnite element solutions do not coincide

with the exact solutions at the nodes.

To study the effect of a Lagrange multiplier located within an element, a simple one-

dimensional problem is considered. The problem is shown in Figure 57.

 

 

Figure 57 The extended problem possess discontinuity across the embedded boundary

xk: Original problem: —AT = x in (xa,xk)with boundary condition

T Ixa = a and Tlxk =0; Extended problem: —AT =0 in (xa,x},)with embedded
boundary constraint: T I xk = 0 and boundary condition T lxa = Ta ’Tlxb = T},;
The exact solution to the above problem is given by:

3
x + x+d xe x ,x
T(x)= a1 3 01 1 (a k) (4.5)
azx +c2x+d2 xe(xk,xb)

where:

54

01:02 :-——;

1 3 3 1 3 3
Ta+g(xa‘xk) Tb+g(xb_xk)
cl: ; 6‘2: ; (4.6)
xa"xk xb‘xk

 

 

1 3 1 3
dlzgxk_qu; dzzgxk—czxk;

This solution will be used below for studying the quality of the LM/FDM solution.

Two measures are used to estimate the quality of the solution: the energy norm and the

L2 norm. These two measures are deﬁned as:

l

2 5
Energy norm: ||e||1=||T—T},"l ={E[%xl—%] {it} (4-7)
1
L, norm: ||e|lo=||T-T},"O ={ E(T—T},)2 dx}2 (4.8)

Clearly, other factors other than the element size that can affect the accuracy of the
LM/FDM solution. It can be observed that by adjusting the relative position of the
Lagrange multiplier inside the element or the value of the applied Dirichlet boundary
condition at the auxiliary part, better LM/FDM solutions can be obtained: in Figure 58,
when the Lagrange multiplier locates on the nodes’ position of the corresponding element
(left) or the exact solution for the extended problem does not possess discontinuity across
the embedded boundary (right), there is no requirement of having a discontinuous slope
inside the corresponding element and the LM/FDM solutions coincide with the exact

solution at the nodes.

55

 

»———~ Exact Solution l
_-~- LM/FDM Solution (lineariL‘

 

   

 

—— Exact Solution
_,_. LM/FDM Solution (linear) j

 

 

(2). 7; =0.1, 7; =0 7;, =—1
Figure 58 LM/FDM solutions coincide with the exact solution at nodes when no
discontinuity locates inside the corresponding element: Lagrange multipliers locate on the

adjacent element nodes (case-l); no discontinuity exists across the interface region (case-
2)

Three inﬂuence factors are studied here: the element size (h), the relative location of the

multipliers within the corresponding element (xkrd , where xkre} = (xk—x1)/ h

with xk the location of the Lagrange multiplier and x1 the adjacent node’s location), and

the value of the applied auxiliary Dirichlet boundary conditions (here a nondimensionless

56

parameter—the ratio of both the Dirichlet boundary conditions, a = Tb ITC, is ﬁrst used

to study the boundary condition’s effect) which determines the jump of the local gradient

of the solution.
4.4.1 Auxiliary Boundary Condition Influence

By adjusting the value of the auxiliary boundary condition (Figure 59), it can be observed
that both the errors reach the minimum values at the speciﬁc condition a z —10 , which is
compatible to the results of the body ﬁtted ﬁnite element method (BFM methods, we
employ this method in this article for the comparison, both methods use the linear
elements here), and then gradually increase when a(correspond to the auxiliary

boundary condition) changes in both directions.

 

0.4 __ Ilello (LM/FDM Solution. h=0.2) 4

0.35 ||e||1(LM/FDMSolution,h=O.2)
, __ "eno (BFM Solution, h=0.22)
_. _ ||e||1 (BFM Solution, h=0.22)

.0
o»

 

 

Iclaello.llell,°
a a.- s 1:.»

p
5?

 

 

 

 

Figure 59. The inﬂuence of the auxiliary boundary condition (xkrel = 0.5 ).

As already mentioned in section 4.3, through domain extension the exact solution may

possess a jump of the local gradient across the embedded boundary 7, it usually cannot

be avoided in the 2-D and 3D case. By adjusting the auxiliary boundary condition, the

magnitude of the jump varies. To make this clearer, another parameter 0 (Figure 60) that

57

represents the magnitude of such a jump is introduced and its relationship with the errors

is studied.

 

 

Figure 60. Deﬁnition of 49
Figure 61 shows the relationship between (9 and the auxiliary boundary condition. At

a = —10 , the exact solution smoothly across the embedded interface (6 = 180°).
Changing auxiliary boundary condition in both direction(a < —10 and a > —10), cross-

boundary discontinuity exists and varies.

 

 

 

 

 

Figure 61. a — 0 relation (where a = —10 corresponds to 6 = 180°)
In LM/FDM methods, when the boundary mesh is not conformal to the underlying

Cartesian mesh, this discontinuity will locate inside the corresponding ﬁctitious domain

58

element (the interfacial element), which means the Lagrange multiplier is located within
that element. For this local element region, the use of smooth weighting function
approximations is inadmissible. It is very difﬁcult to capture features of such a

discontinuity thus introduces error. From ﬁgure 62 it can be seen that (in the case of the

Lagrange multiplier locates inside element) when 6? = 180° (corresponds to a = —10),

which means the exact solution is smooth in the local interfacial region, minimum errors

can be achieved. The errors increase gradually when the angle becomes sharp (0 < 180°)

 

 

 

 

 

  

 

or blunt (0 > 180° ).
”5 __ ||e||0 (LMIFDM Solution. h=0.2)
0,3. - - - ||e||1(LMIFDM Solution. h=0.2)
. _ ueuo (BFM Solution. h=0.22)
0.25» __ _ no"1 (BFM Solution, h=0.22)
E 0.2 ‘.
.—_-90 151
g
0.1 ~
0.05»
0'- 1 _l . l. ._ I l -1
60 120 140 160 240

 

0 (degree)
Figure 62. The inﬂuence of the local gradient jump ( xkre} = 0.5 ).

4.4.2 Mesh size and relative Lagrange multiplier inner element position inﬂuence

A decrease in the element size should signiﬁcantly improve the precision. But from
ﬁgure 63,64, it can be observed that the absolute error does not linearly change but
exhibits some “oscillation” as the mesh is reﬁned. A continuous reduction in the element
size (here we assume uniform discretization) will change the relative location of the
Lagrange multiplier within the corresponding interfacial element. This affects the

accuracy of the local ﬁnite element approximations and the errors are greater than the

59

ones for the body-ﬁtted ﬁnite element solutions except under some speciﬁc condition for

which they are equal.

For a ﬁxed ﬁctitious domain, a change in the element size (here we assume uniform
discretization) will cause the change of the relative location of the Lagrange multipliers
within the corresponding interfacial element, it also oscillates with the element size

(Figure 65). This oscillation is synchronous to the error-element size oscillation.

 

 

 

 

 

- ° - LM/FDM .
4.5 ——.—— BFM J
.2» .
9 Y
9 ﬁx 'Hfs’fi'b'V‘I
-2.5’ ttg”t,a17'\',",'ﬂin,".1,ni,'~".“I11"," .le I» 1‘}| 1 ‘ ‘3]; :1: I,
[:5 1"; 61,1!" II 3 ‘1"; ¢ :5 1! 1’ ' ; 1‘: . l 3
£ ‘3’ 1 ‘11 :1 1; i. 11 ‘1 :1 'l1
2 ‘ l' 11 " " ,'
— -3.51‘ ll 11 11 1' l "
g 11 ll 11 ll 1:
_. ll 1: 1
41 z I l,
‘ ' I
-4.5 l .
-5r
-5.5

4.35 4.3 4575 4.7 4.65 4:6 4.35 41.5 4.35 4:4 435
logh

Figure 63 L2 norm oscillates with the mesh reﬁned (compared with BFM solution,

 

 

--_.ﬂ

4.9 _._ BFM j . g
t. I 4‘ 1"
-2 91 Yl111"1"1,

l r\,"ll|" 11.

7 11
f | H [[1 h
'2-1 9 'If‘ltll ““11“,": . ‘l .1” §

 

 

 

 

; .9 -1.8 4.7 43 4T5 4:4

Figure 64 Energy norm oscillates with the mesh reﬁned (compared with BFM solution,

a=1)

 

t
[I
ll '1‘
1|

04r ”Huh” ""1
. . 1
- 'hi lf‘lll‘f‘l’ ‘1'. . l
1"“ . 11.1 I.’ ‘ ll
| l'Alllr .‘
:1. l 1. ,1 111 1.

03* ,1'“ """1’ 1
0.2» 1'1 ““53”
0.1 1" +11 . 1' H l

 

 

 

__ 1 . 1 1 ' n 1
09015 0.02 0.025 0.03 0.035 0.04 0.045 0.05
h

Figure 65 xkre} oscillates with the mesh reﬁned (a = 1)

Thus the relative position of Lagrange multiplier also plays a key role in affecting the
solution accuracy. It is then necessary to study the inﬂuence of element size and relative

location of Lagrange multiplier within an element separately.

For ﬁxed relative location of the Lagrange multiplier (ﬁgure 66,67), it can be seen that

error linearly changed with element size. Plots of log||e||0 and log||e||l versus logh

show that:

1. When Lagrange multiplier located inside elements (xkrel e(0,l)), the errors can

be expressed as:

log N e ”0: logh + log Co + log f0 (xkrd) (4.9)

log ll elll= 10gh+10g61 +logf1(xkre1) (4.10)

which means convergence ratio of the LM/FDM solution is 1 in the L2 norm and 1

in the energy norm. In this case, the errors are greater than the ones in the body-
ﬁtted ﬁnite element methods and the convergence ratio is smaller than in the

body-ﬁtted ﬁnite element methods (which equals 2).

61

2. When Lagrange multiplier located on the boundary of the correspond elements

(Xkrel = 0 ), then:

log ||e||0=210gh+logc0 (4.11)

log ||e|l1=logh+logcl (4.12)
which is exactly the same errors as the body-ﬁtted ﬁnite element methods: the
convergence rate is 2 in the L2 norm and 1 in the energy norm. Generally, such

situation is very difﬁcult to obtain in 2-D and 3-D when LM/FDM is employed

with independent discretization for the ﬁctitious domain and embedded boundary.

 

xk =0.5 _
0 xk =0.25, 0.75 l
xk =o.10, 0.90 1
-—o- ”( =0 ‘

13.13.

 

 

 

.4 - /
2 109 “
b9 “O“Ogbg C+

15.0 -2f4 -2f2

log Melt,
1
t
1
11
11
u
1
01

 

.5

 

 

 

 

 

.2 '09 h4.0 -1:6 41.4 4.2
Figure 66. Plots of the L,norm of error versus element size. The log-log plots give the
rates of the convergence in the L2 norm. The rates of convergence are given by the slopes

of the lines (the plots shown are for linear elements).

62

      

 

 

 

“i074?” #5 2’i i' ‘
rel '
1n
.2 I /:’
V9 // /’ I ,0
_2.2 w,‘/»’( ,’/
__‘_- Mk ,' ,' 1f //
'5 0.99 /. r '.,/ |
=-2.4 #9 ,’.,r f 1/ 0
3' \\o\\'\ ’1’; ,K/a/
_ ‘9 ’ I.‘ 72/ V9“
’25- (.204 {/0 0'
/ a ,‘ o
o’ I ’
2 8 52" '3”\OQ\\°\\‘
a”!
-3' 3:3;
'%26 24 22 2 -18 16 -L4 12
bgh

Figure 67. Plots of the energy norm of error versus element size. The log-log plots give
the rates of the convergence in the energy norm. The rates of convergence are given by
the slopes of the lines (the plots shown are for linear elements).

If the element size get ﬁxed (ﬁgure 68), it can seen that the error is a parabolic-like
function of the Lagrange multiplier relative location: It will reduce to the minimum value
as it approaches the nodes of the corresponding element but will gradually increased and

ﬁnally reach the peak value as it approaches the center of the element.

 

||e||0= f0(h)+ f(xk,e}) = f0(h)+Za,-(xkre} —0.5)‘ (4.13)
i
llel|1=f1(h)+f(xkre1)=f1(h)+2bj(xkrel -0-5)1 (4-14)

1'
.23: _ 1:31:00 1 ' .000.
0.016’_;—1’l
0.014» ’_'____\ l

=-o.012» ’/ I ‘ ‘.\ J

i 0.00. /, ’ 111311,=f(h)+2bj (mm-051‘

E 0.00% , "
“W ///“\x \
°-°°4 / ||e||°=r(h)+):aI (mm-0.5) \

0.002- \ .

 

 

 

1 0:1 0L2 of: 0:4 0:15 0:0 0:7 0:8 019
Figure 68. Plots of the L2 /energy norm of error versus the relative location of Lagrange
multiplier xkre} .

63

The jump of the local gradient of the solution across the immersed boundary in the
extended domain and the presence of the Lagrange multiplier within an element inﬂuence
the accuracy of the LM/FDM solution. As a result, the numerical solution is very
sensitive to relative location of the Lagrange multiplier within an element, the LM/FDM
solution then cannot achieve the same level accuracy as the body ﬁtted method’s

(BFM’s) solution if the same type of element is used (ﬁgure 69).

 

—-—— ExactSolution
-«-—- LMIFDM Solution
-w— BFM Solution

 

 

 

 

A
r ‘r v v 1 T"
‘I

Figure 69. Solution comparison of the LM/FDM methods and BFM methods (both use
linear element) : BFM solution possess no error at the node while LM/FDM introduces
errors

4.5 Reducing the error

There are several ways to reduce the errors, obviously include mesh reﬁnement (ﬁgure
70), and adjusting the Dirichlet boundary condition for the auxiliary domain in the
extended problem (ﬁgure 72). But as shown already, the LM/FDM solution is very
sensitive to the location of the Lagrange multiplier inside element, mesh reﬁnement for
the ﬁxed ﬁctitious domain may not produce an improved solution (ﬁgure 71). And also,
as we have discussed in section 4.4, even the relative location of Lagrange multiplier

within an element is ﬁxed, LM/FDM methods converge slower than the body ﬁtted

methods. Regarding the auxiliary boundary condition adjustment, it only works for the

external domain-extending problem and the correct direction of the adjustment is not

known in most cases.

Exact Solution

 

_-_.LM/FoM Somﬂon(h=02,xh“=05)
_a- LM/FDM Somuon(h=008.xgﬂ=015)

 

 

Figure 70. Mesh reﬁnement may get a better solution

ExactSolutlon
_ .__ LMIFDM Solution (h=0.2.xkn'=0.5)

—o— LMIFDM Solution (h=0.22. xk"'=0.95)

 

 

 

Figure 71. Mesh reﬁnement may produce a worse solution: with a coarser mesh (h=0.22),
we can even get a better solution (compared with h=0.2)

_ ExactSolution
_,_ LM/FDM Somuon(a=1)
_o- LMIFDM Solution (a=-1)

 

 

Figure 72. By adjusting auxiliary boundary condition adjustment from
Tb = 0.1 (corresponds to a =1) to T}, = —0.1(corresponds to a = —l ), a better solution

can be obtained
High-order element can also provide a better solution because it use high-order

polynomials which can approximate the discontinuity much more accurate than linear

65

element (ﬁgure 73). But the error analysis shows it also is sensitive to the relative inner
element location of the Lagrange multiplier and the L2 norm does not show much
improvement compared to the BFM solution and LM/FDM linear element solution. Thus

in LM/FDM employing high-order elements is also not a featherable way for reducing

error(ﬁgure74 and 75), plus, both mesh reﬁnement and high order element is expensive.

Exact Solution

 

-«--- LM/FDM Solution (linear elements)
—o— LM/FDM Solution (quadratic elements)

 

 

 

(2) local zoom in
Figure 73 High order element can get a better solution (1) because it can approximate the

discontinuity much more accurate (2).

66

 

 

O Y
- ' - LMIFDM (linear)
4 '7— BFM (linear)
—' - LM/FOM (quadratic)
.2 ”t 1°'\.-'ﬁ'~,\,1",.x(1f
<40 '~‘°271’“-‘i’°”i ‘21:" 1..., i'liﬂ "‘ ~“ Lwténa .11» .04.;
=9 0 1 it 1:-
; 4" “ .i "j 1‘11 ,"i
6’ ”k " '1 0. .1
-6~ 1: ll l1 2 ‘1
1 1. z . l
7 l l 1
l l '
-e .

 

 

 

""7105 41.8 4.75 4.7 4.05 4:6 4.55 71:5 4.315 431 4.35
logh

Figure 74 L2 norm oscillates with the mesh reﬁned for quadratic elements

”'5 — LMIFDM (linear) -

— . - BFM (linear)
_. < LMIFDM (quadratic)

 

I
J

 

 

-1 s
'2" . 1 «.2 r ' h-\’ 7‘ I'V‘fLLL .CL’LL»
. 15' ' 1W ' \‘11 ’.‘ W, l 75": ”Ivy" A. ALP”,
,_-2.51 {JILL—4' W 0 rm.
K4” wm\/+JM ‘7 j+ J T
.3—
45.4
'1

’v'. [l
4” J

(“1’11” V" “All”
.1 ' 1' 0.

Il ‘ l

1 {I [l .l

1 l
4.5 1 J l '
i 1 1 . 1

-5;§

'09 Hell

-‘

 

 

 

-1 .6 -1.5 -1.4

log h
Figure 75 Energy norm oscillates with the mesh reﬁned for quadratic elements

.9 4:0 4.7

67

4.6 Fictitious Constraint Methods

The multipliers enforce speciﬁed constraints within the elements but induce errors on the
adjacent nodes. A possible approach to mitigate the errors on the adjacent nodes might
consist of adjusting the embedded boundary condition to a value that reduces the error
locally. This approach thus introduces artiﬁcial boundary conditions that are referred to
here as the ﬁctitious constraints. In this method, the embedded boundary condition is

adjusted to an arbitrary value that can be predicted from the standard LM/FDM solution.

This approach can be illustrated by considering a linear approximation T ', which
satisﬁes a speciﬁc embedded boundary condition T xk :17: (T; is the “Fictitious
Constraint ”), that coincides with the exact solution at the nodes (Figure 76), then

obviously T }: = T1 *(1— xkrel) + T 2 *kaI (T1,T2 are the nodal values for the interfacial

element [X 1, X 2 ]) and may not be the actual embedded boundary condition T x

k =Tk if

a discontinuity in the ﬁrst derivative of the solution exists.

Exact Solution

-4“ LMIFDM Solution T'

 

 

Figure 76. A LM/FDM Solution which satisfying a speciﬁc embedded boundary
condition (“Fictitious Constraint”) coincides with the exact solution at nodes

68

Such a bridge up with a proper value for “Fictitious Constraint ” I]: modiﬁes the
extended problem as following:

Find TeH1(Q), lieEL2 (7) such that

a('i‘, v)+Ir2vd7= [)fvdn v veH(1,(o)

LTpdy=IlTk*,ud7 V ,ueL2(7) (4-15)

~

T=T0 on r

 

where 61(2', v) 2' LIVT-Vv dQ

The solution for this modiﬁed equation can coincide with the solution for equation 4.2 in

most regions except the interfacial region.

Knowledge of the exact solutions for the two different non-interfacial regions [X 01X 1]

and [X 2,X b] of Figure 76 makes it simple to construct a smooth solution deﬁne on

~
5

[X1,X2] that satisﬁes the governing equation a( T v)+J;’/1vd7= bfvdﬂ with

boundary condition [T1,T2] and the corresponding ﬂuxes at position Xl , X 2. The
combined solution of these three regions [Xa,X1],[X1,X2] and [X2,Xb] also

satisﬁed the governing equation 4.2 and the corresponding boundary condition, only the
embedded constraint is different from the original problem. Its value can be adjusted
depending on the selection of the element type.

In this case the “Fictitious Constraint ” can be adjust to the theoretical condition

T x = T }: because the exact solution of [Ti , T2] is already known, but in practice, a two-

k
step predictor-corrector procedure based on the solution obtained for standard LM/FDM

allows to ﬁnd the appropriate value of the ﬁctitious constraint.

69

The procedure consists of ﬁrst assuming that the exact solution in the neighborhood
(must cover at least 2 elements) of the interfacial region (in both the original domain side
and auxiliary domain side) is smooth, almost linear and does not possess any other

discontinuity (Figure 77). It can then be assumed that the solution near the multiplier is

linear (such as in zone 1 deﬁned in [X 0, X k] and zone 2 deﬁned in [X k,X 3]).

Linear
Zone 2
/

 

 

 

 

 

 

Linear
Zone 1
ﬂ“‘\
\\ .....
3f?- .............. f 1.---Xk x2 it-

 

 

 

 

Figure 77. Linearity assumption of the boundary-neighboring zones: Variable
distributions in each separated boundary-neighboring zones (such as zone 1 and zone 2
which belong to different domain-original domain or auxiliary domain) have similar

slope.

In second it is assumed there are sufﬁcient number of elements in each zone so that the

variable distribution in the interfacial region such as [X1,Xk] ( [X k,X 2] can be
predicted from the neighboring element region [X0,X1] ( [X 2,X 3] in the same linear

zone.
Third, the derivative of the LM/FDM solution is assumed to approximate the derivative

of the exact solution in regions except near the interfacial element region (Figure 78).

70

The slope computed near the multiplier is extrapolated from each sides to the slope

within the element associated to a multiplier (such as regions [X 0, X1] and [X 2 , X 3]).

 

;— ExactSolution
_ + — LM /F D M s o lu tio n_ _(_linewa_r)__

Regiok

 

 

 

 

 

 

 

 

I---------------‘---—'----‘

Figure 78. Derivative similarity assumption: Variable derivative of the LM/FDM solution
in each region (region 1 or region 2) is similar to the exact solution.

The predicted solution from standard LMIFDM then can be used to obtain more accurate
nodal values for the interfacial elements! and the “Fictitious Constraint” can be obtained
by the corrections based on these nodal values. This is described in detail below with

Figure 79 for the problem described in Figure 57 as following:

 

 

Figure 79. Implementation of ﬁctitious constraint method in 1-D(The solid line is the
exact solution; the dotted line is the LM/FDM solution and dashed line the ﬁctitious
constraint solution).

71

1. First computing, using LMIFDM method (enforcing the original boundary

condition) a solution T * (which takes value If]; at nodal X 1,X 2 respectively).

For this solution the derivative in the non-interfacial region is similar to the one of
the exact solution (assumption 3). Hence the derivatives (in the LMIFDM solution)

for regions [X 0, X1] (corresponds to the derivative k1) and [X 2,X 3] (corresponds
to the derivative k2) can be estimated.

2. Once k1 and k2 are known, from assumptions 1 and 2, the derivative of the
solution in region [X 1,X k] ([Xk,X2]) is similar to that in region
[X0,X1]([X2,X3]), with assumption 3, it can be treated as the same value
obtained from the LMIFDM solution in region [X 0,X 1]( [X 2,X 3]), which implies
the slopes in the interfacial regions [X1,Xk] ([ X k,X2]) is k1(k2 ). The nodal
values T1,T2 for this interfacial element then can be approximated by a linear

extrapolation from the boundary condition to the adjacent nodes as:

T1 =k1xkre} +Tk (416)
T1 =k2(1-xkre})+Tk .

3. The ﬁctitious constraint thus can be computed based on the adjusted nodal value
and the relative location of the Lagrange multiplier inside [Xl , X2] :
a:
Tk =Tk+T1(l—xkrel)+T2xkrel (4.17)

4. The computed ﬁctitious constraint is ﬁnally used as boundary condition and the

problem is solved again.

72

A much more accurate can then be obtained (Figure 80).

 

 

 

 

——- ExactSolution g
_.- BFM Solution 1
---- Standard LMIFDM Solution

—o- Fictitious Constraint LMIFDM Solution

 

 

 

 

Figure 80. “Predict-Correct Fictitious Constraint” methods obtaines better solution. Note
1
the boundary constraint is adjust from original position A to A.

The error analysis also reﬂects such improvements: The L2 norm and energy norm is
smaller than that obtained from the standard LM/FDM methods (ﬁgure 81, 82) but still
shows some oscillations. They do approach nonetheless the results for standard body
ﬁtted ﬁnite element methods. When compare to the standard LM/FDM method, the
errors are not very sensitive to relative location of the Lagrange multiplier inside
elements (ﬁgure 83 84); the convergence ratio for L2 norm is also improved (ﬁgure 85),
and the errors are also not sensitive to the adjustment of the boundary condition (ﬁgure
86).

We also need to point out that when the discontinuity angle satisfy some speciﬁc

condition (6 e[160°,195°] for this problem), the standard LM/FDM solution is more

accurate than the ﬁctitious constraint method. The solution is then almost smooth across

the boundary so LM/FDM can approximate well the exact solution.

73

 

 

 

 

 

 

-1_5 Standard LMIFDM
- « - Fictitious Constraint LMIFDM
-2 —-—~ BFM
.N' ’l
i -3 “if .liklni.’ “"1““! 12".” I
g l l U L” 11 I; H 1
3.3.51» . l u “l , ”‘0 1
ii i ‘l rut w\' '1" funr‘lL'A'n ”I
4* " Hrrtﬁdﬂiflk‘o‘p‘ ‘t’t' ' /
t,~t mil 1"} [Aim {1 1 ,/
4.51 1‘ M
l
4.0" 4.8 43 4ft; 43 4T4 J
log h

Figure 81. L2 norm comparison for three methods : Standard LMIFDM, Fictitious
Constraint LMIFDM, Body Fitted methods

 

 

"‘7 —~ - Standard LMIFDM .
4.8 - . - Fictitious Constraint LMIFDM
BFM

 

 

 

 

 

Figure 82. Energy norm comparison for three methods: Standard LMIFDM, Fictitious
Constraint LMIFDM, Body Fitted methods

74

x10“ _ . . _
3 —— Standard LMIFDM h_0 04
:_;fi§£“2‘1§_099§ﬂi01 FM/EPM ' '

 

r

F
.

 

 

 

 

Figure 83. L2 norm is not very sensitive to relative location of the Lagrange multiplier
inside elements for Fictitious Constraint LMIFDM methods

 

— Standard LMIFDM i T “:0 04 '
0.014_- — - Fictitious Constraint LMIFDM - '

 

 

 

 

 

Figure 84. Energy norm is not very sensitive to relative location of the Lagrange
multiplier inside elements for Fictitious Constraint LMIFDM methods

75

 

 

- . - «Jo-5 I
o — “rel-0.2. 0‘8
_.- Xk =01,0.91
2 ’9' J l T (‘kmo 4 1
FOM' losillﬁllo"og C293 ' -r 3
smatde q—"+ 0".."— ’ .4
0'?” o" O _ 4r
£‘3 ”ca-H‘H" "H‘.'- Tn’ ,QW“
g .00..“‘M. -' GM‘IK A
o, e\\°5\09 v4,v:‘
'2 '4 “‘09“ v “C": I /

 

 

 

 

Figure 85. Plots of the LGorm of error versus element size. The log-log plots give the

rates of the convergence in the L2 norm. The rates of convergence are given by the slopes
of the lines (the plots shown are for linear elements).

0.35 r

 

‘ _- Hello (Standard LMIFDM)

_ _ . "an1 (Standard LMIFDM)

__.- ||e||o (Fictitious Constraints LMIFDM)
l __ ||e||1 (Fictitious Constraints LMIFDM)

0.3

025* ‘

 

‘ h=0.2

.0
N

O
.5
0|

Ilello. llell1

_O
.a.
r

0.05 *

   

.—

 

 

 

(80.845d10r0'QL—120 11 1! 180’ 200 "220 240
el‘éegre?)

Figure 86. L2 norm and energy norm is not very sensitive to the local gradient jump for
Fictitious Constraint LMIFDM methods
4.7 Shape Reconstruction Method
Since the embedded boundary condition on 7 imposes a very strong constraint on the
solution, properly relaxation of this constraint may also reduce error. One of the easiest
way to realize this idea is shifting the location of Lagrange multiplier to the nearest nodal

points in the corresponding interfacial element. Figure 87 illustrates this idea. Since we

76

assume that the mesh is ﬁne enough and no great variations near the interface region, the

modiﬁed solution should approach the original solution.

 

— Original Solution (exact) W
— — — Modified Solution (exact) ] ,' ' \\

 

  

 

Exact Solution

- + — Standard LMIFDM Solution

#0.. Shape Reconstruction LMIFDM Solution
/,._

 

(b)
Figure 87. Shifting the location of Lagrange multiplier from xk to the nearest node x1 in

the corresponding interfacial element (a) will obtain a better solution (b).

The errors analysis also shows such improvement under some speciﬁc conditions: The

L2 norm and energy norm still show the oscillation but they are smaller than that

obtained from the standard LM/FDM methods and greater than the ﬁctitious constraint
methods (ﬁgure 88); Compare to the standard LM/FDM method, the errors are not very
sensitive to relative location of the Lagrange multiplier within elements (ﬁgure 89) and
the discontinuity angle (ﬁgure 90), but much more sensitive if compared with the one

from Fictitious Constraint methods. When the discontinuity angle satisfy some speciﬁc

77

condition (06 [1300,2100] for this problem), the solution of standard LMIFDM is more

accurate than that of Shape Reconstruction method.

 

"'7 —~— Standard LMIFDM
4.8 ._._ Fictitious Constraint LMIFDM
— - — Shape Reconstruction LMIFDM

 

-1.

.2J

-2.1~

I00 llall1

-2. r

-2,3—- I . IX “X5 ,1 (‘1‘ .
-24- ° I't. ‘ ' 1'" 0
. . ' , “- 1'1 ,
11" i. 4
'2. W t
0

‘2;?_9 ' 4.8 4.7 -1.6 4.5 4.4

 

N
-3é3
O
O

 

 

 

 

. - - — Shape Reconstruction LMIFDM

W111 114).”) 111.11! 1111,11) .. .),-J];

   
 

 

 

 

-1.8 -1.7 -1.6 -1.5 -1.4

(b)
Figure 88. Energy norm (a) and [Q norm (b) comparisons for three methods: Standard

LMIFDM, Fictitious Constraint LMIFDM, Shape Reconstruction LMIFDM

78

 

0.35

0.3 ~

025*

Figure 89. Energy norm (3) and L2 norm (b) are not very sensitive to local gradient
jump (compare to the standard LM/FDM methods) but much more sensitive (compare to

 

1 ——- Standard LMIFDM
e - Fictitious Constraints LMIFDM
La — ~ Shape Reconstruction LMIFQM_1

 

 

 

 

“120’ 160 180 200 250 240
ree)

"811139
(a)

 

 

— Standard LMIFDM
—— Fictitious Constraints LMIFDM
1 — — — Shape Reconstruction LMIFDM

 

 

 

 

 

 

:\1 ”l L 7 "— 1-—F l ‘ —
80 100 120 140 160
0 (degree)

the ﬁctitious constraint LM/FDM) for Shape Reconstruction methods

 

0. 18
0.015
0.014

— Standard LMIFDM
- - - Fictitious Constraint LMIFDM
— Shape Reconstruction LMIFDM

=0.04

 

 

0.013 r
0.012 >

Ilell,

0.01 ~
0.009 .
0.008 -
0.007 -

0.011 ~

 

0.006 ——

 

 

 

x 10‘3

 

~— Standard LMIFDM 11-1304 '
8 — — - Fictitious Constraint LMIFDM ' ‘
7, _’ Shape 53°99§W°ti9n LWEPM

Ilello

 

 

 

rel
(b)
Figure 90. Energy norm (a) and L2 norm (b) are not very sensitive to relative location of
the Lagrange multiplier inside elements (compare to the standard LM/FDM methods) but
much more sensitive (compare to the ﬁctitious constraint LMIFDM) for Shape
Reconstruction LM/FDM methods

4.8 LMIFDM in Two Dimensions
The essential concepts about stand LM/FDM, Fictitious Constraint LM/FDM and Shape
Reconstruction LM/FDM that are outlined in the above sections do not change in two or

three dimensions. The only difﬁculties that emerge are related to the computation of the

boundary integrations associated with the embedded boundary 7 ( ledy , Liudr and

LTlpdy in equation 4.2 for example) and the implementation of the Fictitious

Constraint and shape reconstruction ideas in a higher dimension space.
4.8.1 Boundary integrations
Typically these embedded boundary integration can be evaluated through Gauss-

Legendre Quadrature [27] which is summarized as following (2-D case):

Suppose the function for boundary 7 take the following forms:{x=x(s) , then the

y = y(s)

integration I, f (x, y)d7 can be evaluated through the formula:

80

 

17f (X, ”‘0’ = £912 f (X(S), y(s))\/x2 (s) + y2 (s)ds , which changes to the one dimensional

problem E F (x)dx which then can be approximated using Gauss-Legendre quadrature

pending on the degree of function F (x) .
Another alternative approach developed by Betrand et a1 [32] in 1997 is the collocation

methods which employs the Dirac delta function to enforce the embedded boundary

constraint T = T01), pointwisely, in this method, the space for A}, is deﬁned as a

collection of control points {7,} which discretized the embedded boundary 7, then

i=1

A h can be expressed as:

Ah ={A},

where 6(-) is the Dirac delta function. The boundary integration then can be easily

 

N
11}, = Z&§(i-E),Il~1,lz,mljv 6 R2} , (4.18)
i=1

obtained as:
1M1X1dr=121161X—X11F1X1dr=2mxn (4.19)
i=1 i=1

This treatment simplify the computation process, but introduces errors since it only
enforce the boundary condition at the collocation points. Figure 88 ,89 show the different

between these two approach.

81

 

Figure 91 Absolute errors (leﬁ) and relative error (right) for the stand LMIFDM using
Gauss Legendre quadrature

 

Figure 92 Absolute errors (leﬁ) and relative error (right) for the stand LM/FDM using
Collocation method

4.8.2 Implementation of Fictitious Constraint LMIFDM in 2-D

The application of ﬁctitious constraint methods in 2-D is straightforward. A heat transfer
example is provided to illustrate the implementation in the context of linear element
combined with the using of collocation method ( the integration using Gauss-Legendre

quadrature method is similar):

Suppose we have an interfacial element E0 with a collocation point xk locates inside

(Figure 93). A numerical solution 731d is already obtained from standard LMIFDM.

82

 

 

 

 

 

 

 

 

 

Figure 93 Illustration of the implementation ﬁctitious constraint LM/FDM
Then the ﬁctitious constraint T * at location xk can be obtained from the following

formula:

I! 4 E0
T lxk = TO + Z Yi—pre¢i
i=1

 

 

Y=Xkre1 . (4.20)

 

4
E0 _ l Ek Ek __
where Ii—pre';§ 7; —ZITJ- ¢j|X=Xkre1
J:

Where the following are deﬁned:

83

7111‘“: Solution obtained from the standard LM/FDM for the ith (i=1,2,34)

node of element E k ;

T1131)”, e: Predicted solution for the ith (i =1 ,2,34) node of interfacial element E0;

m : the collocation points x1 ’3 relative position in element’s local
coordinate;
n: the total number of similarity neighboring elements;
k = N1, N 2,...N n is the corresponding similarity neighboring element.
For example, node 1 in interfacial element E0 has four neighbors (E0 , E1, E2 , E4),

among these neighbors only the one located in the non-interfacial region can be used for

computations, these elements are referred as the similarity neighbors (E2 ,.E4 ). A more

accurate nodal value (for node 1 in E0) then can be predicted from following:

 

4

E0 _1 Ek Ek ._

71-1.”,- 2 T1 —ZT,~ 1112141,, (421)
k=2,4 j=l

Table 4 lists all the relative information for the correction of the ﬁctitious constraint at

location xk in E0.

84

Table 4 Computing the ﬁctitious constraint

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Node Similarity
# o o
. Neighbors Neighbors Nodal value Correction
rn E0
150 -1TEk Etc
1 E0,E1,E2,E4 E2,E4 Il-pre-5k§4l 21]:]. ¢J|X=--Xk Xkrel
EEEE EEE TEO -1 TEk 4TEk-—
2 o, 2, 3, 5 2, 3. 5 2_p,e-3 2 -2 j 1’] X: re,
k=2,3,5 1:]
E E E E 750 -1 TEk_£':TEk¢._
3 0’ 5’ 7’ 8 155,57 3—pre-2 3 , j 1 X: rel
k=5,7 1:
50 -1 Ek Ek _.
k=4,6,7 1=
Fictiti 4
ous * _ E0 _.
Const T lxk -To+§7}-P’e¢i X=Xkrel ’
raint '-

 

 

 

Figure 94 and 95 show the errors obtained for the Fictitious Constraint method with the
collocation approach and Gauss-Legendre quadrature approach, respectively. We can see
the ﬁctitious constraint method for these two approaches provides signiﬁcantly more

accurate results than previous LM/FDM.

85

 

Figure 94 Absolute errors (leﬁ) and relative error (right) for the Fictitious constraint
LM/FDM using Gauss Legendre quadrature

 

Figure 95 Absolute errors (leﬁ) and relative error (right) for the Fictitious constraint
LMIFDM using collocation method

4.8.3 Implementation of Shape Reconstruction LMIFDM in 2-D

The implementation of shape reconstruction LM/FDM in 2-D is much more easier. This
can be shown in Figure 96. For collocation approach, the location of the Lagrange
multiplier at xk can be simply moved to the nearest nodal position x2. For the Gauss-

Legendre quadrature approach, the edge that connects the two neighboring Lagrange

multipliers needs to be reconstructed, one simple way to do that is just use the edges that

86

connect the nearest node in the corresponding interfacial element. Other approaches may

also be feasible.

 

 

 

 

*

x3
.1 ‘3 X3

 

7.1x

 

."

 

 

 

 

Figure 96 Illustration of the Shape Reconstruction in 2-D, where xk represents the

location of the Lagrange multipliers and x]: represents the reconstructed location (the

nearest nodal position), the thick-solid line represents the original boundary, and the
thick-dashed line represents the reconstructed boundary.

87

Figure 97 and 98 show the errors obtained for the Shape reconstruction method with the
collocation approach and Gauss-Legendre quadrature approach, respectively. The results
with gauss Legendre quadrature approach provide good results than the one from

standard LMIFDM. With the collocation approach, the results are not good.

 

Figure 97 Absolute errors (left) and relative error (right) for the shape reconstruction
LM/FDM using Gauss Legendre quadrature

 

Figure 98 Absolute errors (left) and relative error (right) for the shape reconstruction
LM/FDM using collocation method

88

Chapter 5

Conclusions

This project explored the applicability and use of moving grid technique in simulating in-
cylinder ﬂows for internal combustion engine. Current CFD codes use three strategies:
spring deforming, local regridding and layering to handle the moving mesh for engine
simulation. Simulations for a generic IC engine intake and compression stroke process
have been performed using FIRE- an engine simulation code based on the above moving
mesh strategies. The in—cylinder ﬂow process are studied and analyzed. Three major
issues that limit the application of moving gird technique: interpolation errors, human
cost for mesh generation and the quality of grid system, are discussed based on the work
of the in-cylinder ﬂow simulations for Generic IC engine.

Fictitious domain methods reduce the complexity of meshing by using a simpler auxiliary
domain and augmenting a functional to implement the original boundary conditions in the
extended auxiliary domain. A possible ﬁctitious domain approach, attractive for its
simplicity in 2-D and 3-D, is based on using Lagrange multipliers combined either with a
“collocation-like” method or the Gauss Legendre quadrature based boundary integration.
The presence of a Lagrange multiplier within an element in LM/FDM introduces
signiﬁcant errors on the adjacent nodes when the extended problems process
discontinuity across the embedded boundary, the magnitude of the errors is related to the
jump of the local gradient across the embedded boundary and the relative location of the
multiplier with respect to the adjacent nodes. Compared to the body ﬁtted methods, this

reduces the convergence ratio in L2 norm with ﬁxed relative Lagrange multiplier position

89

for the LM/FDM methods. Either adjusting the boundary condition with the factitious
constraint or changing the location of the Lagrange multiplier to the adjacent nodes can
provide improved numerical solutions for 1-D and 2-D steady state heat conduction

problem. Further investigations and improvement methods are desired.

90

APPENDIX

Turbulence Boundary Layer

1. The Concept of y+

The dimensionless symbol y+ is related to the characteristics of near-wall turbulence
ﬂows. In ﬂows along solid boundaries, there is a region of intertia-dominated ﬂow far
away form the wall and a thin layer within which viscous effects are important. Close to

the wall, viscous effects dominate the flow and the mean ﬂow velocity only depends on

the distance from the wall:
u
14+ =—=[M]=f(y+)
“r 1“

The above equation is termed law of the wall and contains deﬁnitions for the two

dimensionless group y+ and u+. u, is the friction velocity and is deﬁned as :

y+ provides a useful measure of the inﬂuence of the viscous layer (near the wall) given a

ﬂow velocity. The near-wall region can be largely subdivided into three layers. In the
innermost layer (viscous sublayer), the flow is almost laminar and viscosity plays an
important role in the momentum transfer. In the outer layer (fully turbulence layer),

turbulence plays a major role. This region is also called the log-law region as u+ is a

straight-line function of y+. There is an intermediate region between these two layers

where the effect of viscosity and turbulence are equally important.

91

2. Wall Function

Two approaches available to model the near-wall region. In one approach, the viscous-
affected inner region (viscous sublayer and buffer layer) is not resolved. Instead, semi-
empirical formulas called “wall functions” are used to bridge these regions between the
wall and the fully turbulent region. With wall functions, the need to modify the
turbulence models to account for the presence of the wall is obviated. Another approach
modiﬁes the turbulence models to enable the viscosity-affected region to be resolved with
a mesh all the way to the wall, including the viscous sublayer.

The wall function approach saves considerable computational resources because the
viscosity-affected near wall region, where the solution variables change most rapidly,
does not need to be resolved. This provides a practical option for turbulence ﬂow
simulations.

However, the wall function approach is inadequate in situations where low Reynolds
number effects prevail in the ﬂow domain and the assumptions underlying the wall
functions cease to be valid. Such situations require near-wall models that are valid in the
viscosity-affected region and hence can be integrated all the way to the wall.

3. Near Wall Mesh Generation for Wall Functions

The region near the wall is meshed ﬁner than the rest of the cross section, as it contain
the maximum amount of gradients, The distance form the wall at the wall-adjacent cells

must be determined by considering the range over which the log-law is valid. The ﬁrst

grid point away form the wall is usually placed in the log-law region with y+ 6 [30,300] .

At least ﬁve points must be placed in the boundary region to resolve the gradients

sufﬁciently for most ﬂow situations.

92

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

BIBLIOGRAPHY

Merle C. Potter and John F. Foss, Fluid Mechanics, Great Lakes Press, INC., MI,
1982

Shyy W, Udaykumar HS, Rao MM, and Smith RW, Computational Fluid
Dynamics with Moving Boundaries, Taylor and Francis, Washington DC. 1996

W Shyy and M Francois, HS Udaykumar, N N’dri and R Tran-Son—Tay, Moving
boundaries in micro-scale bioﬂuid dynamics, App Mech Rev vol 54, no 5,
September 2001

W.Shyy, S.S. Thakur, H.0uyang, J.Liu, and E. Blosch, Computational
Techniques for complex transport phenomena, New York, Cambridge, 1997

PD. Thomas and CK. Lombard, Geometric conservation law and its application
to ﬂow computations on moving grid, AIAA Vol. 17 No. 10 , pp.1030-1037, 1979

Weiming Cao, Weizhang Huang, and Robert D. Russell, A moving mesh method
based on the geometric conservation law, SIAM J. SCI. Comput. Vol. 24 No. 1,
pp. 118-142, 2002

Ching-Jen Chen, Shenq-Yuh Jaw, Fundamentals of turbulence modeling,
Washington, DC, Taylor & Francis, 1998.

Stanisic, M.M . The mathematical theory of turbulence, New York, Springer-
Verlag, 1988.

Dale A. Anderson, John C. Tannehill and Richard H. Pletcher, Computational
ﬂuid mechanics and heat transfer, Washington, DC, Taylor & Francis, 1997

Graham F. Carey, Computational Grid, Generation, adaptation, and solution
strategies, Washington, DC, London, Taylor & Francis, 1997

Hozef Arif, Application of computational ﬂuid dynamics to the modeling of ﬂow
in horizontal wells, Master thesis, Standford University, June 1999

AD. Young, Boundary layers, Washington, DC, AIAA, BSP Professional Books,
1989

FLUENT Inc. FLUENT 5 user guide, Vol 1-4, Fluent Inc., Lebanon, New
Hampshere, USA, 1998

93

[14]

[151

[16]

[171

[18]

[19]

[201

[21]

[22]

[23]

[24]

[25]

[26]

[271

AVL Inc., FIRE 7.2b user guide, AVL LIST GmbH, Graz, Austria, 2000

Amsden A.A, Kiva-3V a block structured KIVA program for Engines with
vertical or canted valves, Las Alamos National Laboratory, LA-l3313-MS, 1997

Adapco Inc. Star-CD user guide Version 3.10 B, adapco Inc. Melville, New York,
2000

John L. Lumley, Engines, an introduction, New York, Cambridge, 1999

J .1. Ramos. Internal combustion engine modeling, New York, Hemisphere Pub.
Corp., 1989.

NC. Markatos, Computer simulation for ﬂuid ﬂow, heat and mass transfer and
combustion in reciprocation engines, New York, Hemisphere Publishing
Corporation, 1989

B.Basara, Computations of automotive ﬂows using the second-moment closure,
ECCOMAS September 2000

Jiang Yong, Zhu Ning , Qiu Rong and Fan Weicheng, Transient ﬂow analysis of
the helical intake port and cylinder of a direct injection diesel engine, Japan, Fluid
and Heat Engineering study, Vol 36 No. 2, 2002

Yang-Liang Jeng, Rong-Che Chen and Chao-He Chang, Study of tumbling
motion generated during intake in a bowl-in-piston engine, Taiwan, Journal of
Marine Science and technology, Vol. 7 No.1, pp. 52-64 1999

Henk, K., Numerical simulation and experimental veriﬁcation of DI diesel intake
port designs, http://www.cyclone.nl/ports/ports.htm

Bassem Ramadan, A study of swirl generation in D1 engines using Kiva-3V, 11th
international multidimensional engine modeling user’s group meeting, MI, 2001.

P.J. Colucci, D.Lee, C.K.Lim, G.Goldin, In-Cylinder engine modeling
developments at Fluent, 12th international multidimensional engine modeling
user’s group meeting, MI, 2002.

Gisoo Hyun and Mitsuharu Oguma, 3-D CFD analysis of the mixture formation
process in an LPG D1 81 engine for heavy duty vehicles, 12th international
multidimensional engine modeling user’s group meeting, MI, 2002.

R.Glowinski, T.-W. Pan, J. Periaux, A ﬁctitious domain method for Dirichlet
problem and applications, Comput. Meth. Appl. Mech. Eng. 111, pp. 283-303,
1994

94

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[371

[381

R.Glowinski, T.-W. Pan, J. Periaux, A ﬁctitious domain method for external
incompressible viscous ﬂow modeled by Navier-Stokes equations, Comput. Meth.
Appl. Mech. Eng. 112, pp. 133-148, 1994

R.Glowinski, T.-W. Pan, J. Periaux, A Lagrange multiplier/ﬁctitious domain
method for the Dirichlet problem. Generalization to some ﬂow problems, Japan
J.Industrial Appl. Math. 12, pp. 87-108, 1995

R.Glowinski, T.-W. Pan, J. Periaux, Fictitious domain methods for
irncompressible viscous ﬂow around moving rigid bodies, in : J.R. Whiteman
(Ed), The Mathematics of Finite Elements and Applications, Highlight 1996,
Wiley, Chichester, England, pp. 155-174, 1997

R.Glowinski, T.-W. Pan, J. Periaux, Distributed Lagrange multiplier methods for
incompressible viscous ﬂow around moving rigid bodies, , Comput. Meth. Appl.
Mech. Eng. 151, pp. 181-194,1997

F. Bertrand, P.A. Tanguy, F. Thibault, A three-dimensional ﬁctitious domain
method for incompressible ﬂuid ﬂow problems, Internat. J. numer. Meth. Fluids
25, pp. 719-736, 1997

R.Glowinski, T.-W. Pan, J. Periaux, Distributed Lagrange multiplier methods for
incompressible viscous ﬂow around moving rigid bodies, , Comput. Meth. Appl.
Mech. Eng. 151, pp. 181-194,1998

R.Glowinski, T.-W. Pan, T.I. Hesla, D.D. Joseph, A distributed Lagrange
multiplier/ﬁctitious domain method for particulate ﬂows, lntemat. J. Multiphase
Flow 25, pp. 755-794, 1998

R.Glowinski, T.-W. Pan, T.I. Hesla, D.D. Joseph, , J. Periaux, A distributed
Lagrange multiplier/ﬁctitious domain method for the simulation of ﬂow around
moving rigid bodies: application to particulate ﬂow, Comput. Meth. Appl. Mech.
Eng. 184, pp. 241-267, 2000

V.Girault and R.Glowinski, Error Analysis of a ﬁctitious domain method applied
to a Dirichilet problem. Japan J. of Ind. and App.Math. 12, (1995) pp. 487-514,

J.N. Reddy, An introduction to the Finite element method, 2d ed. McGraw-Hill,
New York/Hemisphere, Washington, DC, 1993.

Antony Jameson, The present status, challenges and future developments in

computational ﬂuid dynamics, Proceedings of the 7th AGARD Fluid Dynamics
Panel Symposium, Seville, October 1995.

95

[39] Yanbing Li and Andre Benard, A method for reducing the errors in ﬁctitious
domain methods applied to heat transfer problems, accepted by The 6th ASME-
J SME Thermal Engineering Joint Conference, March 16-20, 2003

96

11111111111111