- Email: [email protected]

Antonio Posa, Marcos Vanella, Elias Balaras

PII: DOI: Reference:

S0021-9991(17)30707-6 https://doi.org/10.1016/j.jcp.2017.09.047 YJCPH 7620

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

30 March 2017 11 August 2017 24 September 2017

Please cite this article in press as: A. Posa et al., An adaptive reconstruction for Lagrangian, direct-forcing, immersed-boundary methods, J. Comput. Phys. (2017), https://doi.org/10.1016/j.jcp.2017.09.047

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An adaptive reconstruction for Lagrangian, direct-forcing, immersed-boundary methods Antonio Posaa,b,∗, Marcos Vanellaa,∗, Elias Balarasa,∗∗ a Department of Mechanical and Aerospace Engineering, The George Washington University, 800 22nd Street NW, Washington DC, 20052, USA b CNR-INSEAN, Via di Vallerano 139, 00128 Roma, Italy

Abstract Lagrangian, direct-forcing, immersed boundary (IB) methods have been receiving increased attention due to their robustness in complex ﬂuid-structure interaction problems. They are very sensitive, however, on the selection of the Lagrangian grid, which is typically used to deﬁne a solid or ﬂexible body immersed in a ﬂuid ﬂow. In the present work we propose a cost-eﬃcient solution to this problem without compromising accuracy. Central to our approach is the use of isoparametric mapping to bridge the relative resolution requirements of Lagrangian IB, and Eulerian grids. With this approach, the density of surface Lagrangian markers, which is essential to properly enforce boundary conditions, is adapted dynamically based on the characteristics of the underlying Eulerian grid. The markers are not stored and the Lagrangian data-structure is not modiﬁed. The proposed scheme is implemented in the framework of a moving least squares reconstruction formulation, but it can be adapted to any Lagrangian, direct-forcing formulation. The accuracy and robustness of the approach is demonstrated in a variety of test cases of increasing complexity. Keywords: Cartesian grids, Immersed boundary method, Direct Forcing, Lagrangian Forcing PACS:

1. Introduction Today, immersed-boundary (IB) formulations are an established alternative to classical, boundary conforming approaches, especially in ﬂuid-structure interaction problems. In most cases the equations of ﬂuid motion are solved on a ﬁxed (Eulerian) grid, and the immersed body is usually represented by a set of surface nodes and elements, deﬁning a Lagrangian grid. These two grids are almost never aligned, and therefore, a variety of strategies has been developed to enforce the boundary conditions on the surface of the body (see [13] for a review). In moving boundary problems, one needs to constantly track the Lagrangian grid, and transfer ﬂow variables and/or forcing functions from one grid to the other. As a result, the relation between the Lagrangian and Eulerian grids is important to the accuracy, robustness, and eﬃciency in the majority of immersed boundary methods. In fact, most schemes work best when the Lagrangian elements on the surface of the body are of comparable size to the local Eulerian grid. The criteria, however, for designing either one are usually very diﬀerent. In the case of the Lagrangian grid, for example, the number and size of surface elements is guided by the complexity and local curvature of the surface, while for the Eulerian grid is guided by velocity and pressure gradients. The sensitivity to the relation between the two grids also depends on the speciﬁcs of the forcing scheme. For example, direct forcing schemes, where the solution is reconstructed around an Eulerian grid node in the neighborhood of the body [5, 9, 2], are less sensitive compared to schemes, where a forcing function is computed at the Lagrangian markers and then transferred to the Eulerian nodes [21, 23, 16] (see [24] for a direct comparison). The former class ∗ Authors

have contributed equally to this work. author. Email addresses: antonio.pos[email protected] (Antonio Posa ), [email protected] (Marcos Vanella ), [email protected] (Elias Balaras) ∗∗ Corresponding

Preprint submitted to Journal of Computational Physics

September 26, 2017

f s

f

s h

l1

l1

h 2

l2

l2

s

s

(a)

(b)

Figure 1: Immersed boundary forcing stencils for single Lagrangian grid and two diﬀerent staggered Eulerian grids: Ω f and Ω s are the ﬂuid and solid regions respectively, separated by a surface Γ s . (a) Two marker points l1 , l2 are suﬃcient to marginally cover the forcing band on an Eulerian grid of cell size h. (b) Regions with no force imposition (holes) are left for the same Lagrangian markers on a ﬁner Eulerian grid of cell size h/2.

of methods deﬁne a forcing function in Eulerian mesh nodes close to the solid, employing the kinematical boundary condition on nearby body surface elements and an interpolation strategy. They produce a sharp representation of the boundary and do not require matching Lagrangian and Eulerian grids sizes. They can potentially lead to large oscillations in the solution in moving boundary problems, due to spatiotemporal discontinuities in momentum forcing as the solid transverses the grid. Several successful treatments that alleviate this issue have been proposed in the literature (see for example [26, 18, 27]). The latter class of methods is gaining popularity due to their robustness and regular force distribution, especially in cases of moving boundaries. They typically lead to small spurious force oscillations on these cases, at the expense of a smoothed boundary representation. They were introduced by Uhlmann [21] and the main idea was to incorporate Peskin’s [15] regularized Delta functions as transfer kernels between material and spatial variables, into a direct forcing formulation. This technique was implemented in a ﬁnite-diﬀerence, fractional-step context and the overall method was designed to model suspended rigid spherical particles in laminar and turbulent ﬂows. This method was extended by Vanella & Balaras [23] to unstructured surface tessellations of moving/deforming bodies of arbitrary shape, by introducing a local moving-least-square (MLS) approximation to build the transfer functions between the Eulerian and Lagrangian grids. Along the same lines Pinelli et al. [16] used the reproducing kernel particle method proposed in [11] for the interpolation and spreading. One of the main requirements in all Lagrangian forcing schemes discussed above, is that the Lagrangian and Eulerian grids need to be of comparable resolution to ensure proper transfer of forces from one grid to the other. In most cases, the adopted Lagrangian grid is at least twice as ﬁne as the Eulerian one. If this condition is not met, certain areas of the Eulerian grid may not be forced, and transpiration through the body is possible. An example where the Lagrangian grid is ﬁne enough to maintain the continuity of the forcing band across the Eulerian grid is sketched in Figure 1a, while a case where a breakdown of continuity of the forcing ﬁeld happens is shown in Figure 1b. The latter scenario is common in practical applications as the Eulerian grid is reﬁned in areas of high velocity gradients, and Lagrangian grid is reﬁned in areas with high curvature. Even in simple geometrical conﬁgurations it is diﬃcult to enforce both constraints simultaneously. Consider for example, high Reynolds number ﬂow around a sphere. The boundary layers originating from the stagnation point in the front are very thin, requiring a ﬁne Eulerian grid to resolve them. As the Reynolds number increases, the Eulerian grid needs to be reﬁned accordingly. To ensure that the forcing is properly transferred (no-slip condition correctly approximated), similar reﬁnement needs to be done to the Lagrangian grid, although the geometry is already represented by a suﬃcient number of elements. It is therefore conceivable that some form of local adaptive reﬁnement/de-reﬁnement of the Lagrangian grid, to satisfy the constraints coming form neighboring Eulerian nodes, will improve the accuracy and robustness of most immersed boundary approaches. An h-reﬁnement strategy, however, where local reﬁnement is achieved by splitting 2

existing surface elements into several smaller ones, or by locally introducing additional nodes, will signiﬁcantly increase the cost as well as the implementation complexity. In the present work we propose a strategy to locally reﬁne the Lagrangian grid that does not alter the existing data structure nor increase storage requirements. It is based on isoparametric surface element mapping and can be coupled to any IB approach. In particular we will treat the Lagrangian surface with linear triangles, but the approach can be readily extended to diﬀerent types of surface elements and polynomial approximations. The implementation of the immersed boundary technique will be done considering the Lagrangian forcing MLS scheme we proposed in an earlier paper [23], although other Lagrangian forcing schemes could be employed. In the following sections we will discuss the family of Lagrangian, direct-forcing, IB schemes in the context of a projection method for incompressible ﬂows. We will then present a new strategy to dynamically change the Lagrangian grid resolution to ensure the continuity of the forcing band across the Eulerian grid. Finally, test problems of increasing complexity with moving and stationary immersed bodies are discussed to evaluate the accuracy and robustness of the proposed method, followed by a summary and conclusions. 2. Methodologies 2.1. Overview Let us consider a solid body, ΩS , surrounded by an incompressible ﬂuid region, Ω f . Each point in a subdomain, ε, of the body is labeled by a unique set of coordinates ξ, η, ζ, which can be in general local and curvilinear. At a given time, t, the location of a particle belonging to the body is given by the vector: X = X(ξ, η, ζ, t) = X0 (ξ, η, ζ) + υ(ξ, η, ζ, t),

(1)

where, X0 , corresponds to an initial reference location and, υ, is a given displacement ﬁeld. As the material coordinates, ξ, η, ζ, are not time dependent, the particle’s velocity and acceleration are denoted as: ∂ ˙ η, ζ, t) X(ξ, η, ζ, t) = υ(ξ, ∂t

U

=

u(X, t) =

˙ U

=

∂2 ¨ η, ζ, t) X(ξ, η, ζ, t) = υ(ξ, ∂t2

(2) (3)

A location within the Eulerian three-dimensional space is described by a vector x, and the ﬁeld u(x, t) describes the velocity of the particle whose position at time t is X = x. This set of kinematical relations provides the means to combine a Lagrangian description for the body motion with an Eulerian description of the ﬂuid ﬂow, a typical feature of immersed boundary methods. For example, Peskin [15] initially used the Dirac delta function to transfer the velocity ﬁeld among material and Eulerian descriptions: U = u(X, t) = u(x, t)δ (x − X(ξ, η, ζ, t)) dx (4) This expression, together with transfer formulas for solid density and elastic internal forces, which also make use of the Dirac function, form the basis of the ﬂuid-structure interaction model for the continuum based IB method proposed in [15]. Direct-forcing, Lagrangian methods use the same strategy, once however, ﬂuid and solid regions have been discretized in space and time. The Navier-Stokes equations for incompressible ﬂow are advanced in time using an exact projection method [22]: u˜ i − uni Δt

=

L(φ) = un+1 i

=

N(˜ui , uni ) − G(φn−1 ) + 1 D(u˜i ), Δt u˜i − ΔtG(φ)

1 L(˜ui , uni ) + fi , Re

(5) (6) (7)

where i = 1, 2, 3 represents the vector component direction, u˜ i is the intermediate non-divergence free velocity, uni is the divergence free velocity at time step n, Δt is the timestep, φ is the modiﬁed pressure, and N, G, L are the discrete 3

non-linear, gradient and laplacian operators respectively. The operators include coeﬃcients that are speciﬁc to the selected time advancement scheme. In all computations below we assume explicit integration of the advective and viscous diﬀusion terms, except for the cases where cylindrical coordinates are utilized and both advective and viscous terms in the azimuthal direction are treated implicitly. Details on the fully explicit scheme can be found in [25], and on the semi-implicit scheme on cylindrical coordinate grids in [1]. The direct Lagrangian forcing scheme proposed in [21] involves the computation of a force ﬁeld on the Lagrangian markers, l, which is transferred back to surrounding ﬂuid grid points. In order to accomplish this, discrete interpolation and spreading operators are used: (8) U˜ iε,l = Ih u˜ ki ε,l k (9) fi = T h Fi , where capital letters refer to a variable located at a given Lagrangian point, l, of a subdomain, ε, (this notation will be used for the mapping discussion of the next sections). We will focus on the case of surface Lagrangian forcing methods (i.e. [15, 21]), for which Lagrangian markers are distributed on surface subdomains, requiring a local set of coordinates ξ, η. The discussion can be extended to volume forcing within the solid. The set of Eulerian ﬂuid points related to the marker, (ε, l) is deﬁned by k = 1, ..., ne. Also Ih , and T h are the discrete interpolation and spreading operators, which may be obtained by diﬀerent means. Uhlman [21] for example makes use of the regularized Delta function, δh , deﬁned in [15, 17] to construct these operations, k ε,l (10) ΔV k u˜ ki ; fik = (ε,l) ⊂ nl δh xk − Xε,l ΔVεl Fiε,l , U˜ iε,l = ne k=1 δh x − X where for each marker, (ε, l), the ﬁrst sum involves all points k = 1, ..., ne of the associated Eulerian stencil. On the other hand, the force spreading to an Eulerian point k comprises all markers, (ε, l), in the set nl related to it. The volumes ΔV k , ΔVεl are representative of the Eulerian and Lagrangian grids, and deﬁned such that the total force and moment are preserved by the transfer operation. Regardless of the deﬁnition of operators (8)-(9), the IB forcing strategy follows the following sequence: 1. Compute u˜ i from equation (5) discarding the IB force term fi . 2. Proceed marker by marker. For surface marker (ε, l), interpolate U˜ iε,l from surrounding Eulerian stencil using equation (8), and one of the methods to construct Ih (i.e. [21, 23, 16]). 3. Compute the force density Fiε,l on Lagrangian marker points as, Fiε,l =

Uiε,l − U˜ iε,l Δt

(11)

where Uiε,l is the velocity of the solid marker on direction i, deﬁned by equation (2). The previous expression is analogous to the expresion used on Eulerian direct forcing methods (see for example [5, 9, 2]). 4. Spread the force Fiε,l on Lagrangian marker back to the corresponding Eulerian stencil ( fik ) using equation (9). 5. The resulting intermediate velocity that takes into account the immersed bodies is ﬁnally computed from equation (5). When some terms on the momentum equations are treated implicitly, the ﬁrst computation of u˜ i (item 1) can be done explicitly, while its ﬁnal computation taking into account IB forces (item 5), is done implicitly [9]. This approach leads to a boundary condition approximation on the body, whose accuracy can depend on the timestep size. Schemes that alleviate this problem have been proposed in the literature [8], where an iterative procedure among implicit momentum and IB forcing may be required. 2.2. Dynamic Lagrangian Mapping In this section we propose a method to reﬁne/dereﬁne the Lagrangian grid based on polynomial isoparametric mapping. Figure 2a, shows the body surface, which is represented by a tessellation, , composed of the union of triangular subdomains {ε} ⊂ . The vertex points for this mesh are the set of nodes, ℵι , and the information for the local nodes, nι , ι = 1, 2, 3, of every element, ε, is readily known, using a connectivity list typical of ﬁnite element 4

z

hy

f

hz

1 n3

N2 = 0

n3

nˆ x2

s

n1

x1 n 2

N3 =

1 2

N1 =

1 2

N1 = 0

n1

y

0

n2 1 N2 = 2

1

N3 = 0

x

(a)

(b)

Figure 2: Triangulation representing the IB surface inside a Cartesian Eulerian grid: (a) A particular element ε has vertices n1 , n2 , n3 in local numbering, deﬁning the triangle geometric quantities (area Aε , normal nˆ ε ). (b) The natural coordinate system is used to represent the master triangle and map points within it to physical coordinates. The linear shape functions are N1 , N2 , N3 .

˙ nι , may be prescribed in time, or driven by data-structures [6, 3]. Kinematical information of these nodes, Xnι , Unι , U the general translation-rotation of the center of mass in rigid bodies, or arising from the discretization of the equations of motion of the solid in elastic structures [3]. Then, linear interpolation functions are used to deﬁne the location of marker points Xε,l within the triangle ε. In order to simplify the construction of the interpolation, a master element of unitary side length in the direction of the natural coordinates, ξ, η (see ﬁgure 2b) is considered [3]. The interpolation functions for this linear triangle are: N1 = 1 − ξ − η ;

N2 = ξ ;

N3 = η

(12)

The physical location of a particular point l within the triangle, deﬁned by natural coordinates ξ, η is given by Xε,l =

3 ι=1

Nιε,l Xnι

(13)

˙ nι ) the mapping When the same shape functions are also used to interpolate ﬁeld variables of interest (i.e. Unι , U is called isoparametric [6, 3]. The transformation from natural to planar coordinates x1 , x2 for this linear triangle is characterized by the following Jacobian matrix, ⎡ ∂x1 ∂x2 ⎤ ⎥⎥ ⎢⎢ ∂ξ ∂ξ ⎥ ⎥⎦ [J] = ⎢⎢⎢⎣ ∂x (14) ∂x ⎥ 1

∂η

2

∂η

The Jacobian J = det(J) deﬁnes the ratio of physical to natural areas dx1 dx2 = Jdξdη, and can be used for the area integration on isoparametric elements. Consider now in the three Cartesian coordinates the computation of the unit normal vector nˆ lε to the surface on the Lagrangian marker l position. In general, two vectors tangent to the element surface can be computed as follows: ε,l ε,l ε,l ∂Xε,l ∂N1 n1 ∂N2 n2 ∂N3 n3 = X + X + X = (xn2 − xn1 ) ˆi + (yn2 − yn1 ) ˆj + (zn2 − zn1 ) kˆ ∂ξ ∂ξ ∂ξ ∂ξ

(15)

ε,l ε,l ε,l ∂Xε,l ∂N1 n1 ∂N2 n2 ∂N3 n3 = X + X + X = (xn3 − xn1 ) ˆi + (yn3 − yn1 ) ˆj + (zn3 − zn1 ) kˆ ∂η ∂η ∂η ∂η

(16)

5

where ˆi, ˆj, kˆ are the Cartesian unit vectors. The normal is deﬁned as, ∂Xε,l /∂ξ × ∂Xε,l /∂η

. nˆ lε =

∂Xε,l /∂ξ × ∂Xε,l /∂η

(17)

We ﬁnd that, as expected, for the linear triangle the tangent vectors ∂Xε,l /∂ξ, ∂Xε,l /∂η and the resulting normal are constant, regardless of their location within the element ε, and therefore nˆ lε = nˆ ε . Several diﬀerent mapping methods can be devised to distribute Lagrangian markers within the triangle. In the following we describe two methods that proved simple and robust, and have shown very good behavior in numerical computations. In fact numerical experiments have shown that, for suﬃciently dense marker descriptions, both methods provide practically equal results. Consider a characteristic length hε for triangle ε, which can be the maximum side length or the diagonal of the elements bounding box. To obtain the number of marker points on both ξ, η directions, mη =mξ , we use the condition hε /mξ ≤ cΔ where Δ is the local Eulerian grid characteristic size (diagonal cell length, or more conservatively the smallest cell size on the three directions local to ε, min(h x , hy , hz )). The constant c < 1 deﬁnes the ratio of the distance between markers within the element to the local cell size. Therefore, one can set mξ =

hε + 1. cΔ

(18)

The total number of marker particles on triangle ε is then mε = mξ (mξ + 1)/2. Given that for this linear triangle the determinant J is constant, the area in physical coordinates assigned to each marker l, Alε , is assumed constant and equal to Aε /mε . This result has been used for the area computations of mapping A and B in the next sections. For surface elements deﬁned with shape functions of order higher than linear, both Alε and nˆ lε would have to be recomputed at each marker l, evaluating the interpolation functions and their derivatives at the corresponding location ξ, η. In the following, two mapping strategies are presented. 2.2.1. Mapping A: Evenly distributed points across the element In ﬁgure 3a, the marker point mapping A is illustrated on natural coordinates. In order to locate the set of points within the triangle, a variable ϕ deﬁning the distance to corners in ξ and η is used. Here, ϕ = aΔξ , where Δξ = 1/mξ . The factor a determines the marker location respect to triangle corners. A value a = 1/2 produces a fairly uniform marker paving on a regular patch of triangles in physical space (see ﬁgure 3b for a two dimensional case), and was used in this work. The method for point distribution is exactly the same for ξ or η. The ﬁrst line of markers is positioned at a distance ξ st = ϕ from the origin, and the last marker at a location ξend = 1 − 2ϕ. The distance between points is obtained as: 1 − 3ϕ ξend − ξ st Δξ p = = . (19) (mξ − 1) (mξ − 1) As expected Δηp = Δξ p . The paving algorithm consists of a double loop, ﬁrst in index j (from 1 to mη ) and then in index i (1 to mξ + 1 − j), using lower triangular limits. The index of a particular point (ε, l) can be obtained by the formula, ( j − 1)( j − 2) (20) l = i + ( j − 1)mξ − 2 and the position of the marker in physical coordinates is given by the corresponding interpolation: ξi

= ϕ + (i − 1)Δξ p

ηj

= ϕ + ( j − 1)Δηp

ε,l

X

=

N1ε,l Xn1

+

N2ε,l Xn2

(21) (22) +

N3ε,l Xn3 ,

(23)

where the shape functions associated with the mapped point are N1ε,l = 1 − ξi − η j , N2ε,l = ξi and N3ε,l = η j . Velocities ˙ ε,l are computed in the same manner. In cases of general rigid body motion, both velocity and Uε,l and accelerations U acceleration of the markers can also be deﬁned analytically, using for example, the markers location respect to the body center of mass and the rigid kinematics equations. 6

1

1

0.5

x2

0

0.5

1

0 p

1

st

1.5

end

(a)

1

0.5

0

x1

0.5

1

1.5

(b)

Figure 3: Mapping A, evenly distributed points across the element: (a) Distribution of marker points on the master element, the number of markers in the natural η direction Nη is always equal to the Nξ points on ξ. Positioning and distance variables are the same in both directions. Here Nξ = 4 for illustration. (b) A regular set of triangles in physical coordinates x1 , x2 with their corresponding mapped Lagrangian points, again Nξ = 4.

2.2.2. Mapping B: Subsection of element in regular triangular patches An alternative scheme for mapping Lagrangian markers on the surface element can be seen in Figure 4a. In this case, the master element is divided in evenly sized triangular sub-patches with side length Δξ , and Δη = Δξ in natural coordinates. The mapping is divided in two sets of markers: one corresponding to the centroids of the resulting lower triangular patches (blue in Figure 4a) of size mε1 = mξ (mξ + 1)/2. And the other composed of centroids of the upper triangular patches, which is of size mε2 = (mξ − 1)mξ /2. The total number of control points is mε = mε1 + mε2 . The paving process is now split in two double loops similar to what was discussed in the previous section. First, a double loop on j (1 to mη ) and i (1 to mξ + 1 − j) is used to map the mε1 lower triangular centroids located in, ξi = ϕ + (i − 1)Δξ

;

η j = ϕ + ( j − 1)Δη ,

(24)

where in this case, ϕ = Δξ /3. The shape functions, physical coordinates and kinematics computations are done using Eq. (23). The locations of points indexed from (mε1 + 1) to mε (markers 7 to 9 in Figure 4a) are deﬁned on a second double loop on j (1 to mη − 1) and i (1 to mξ − j), where the natural coordinates are given by Eq. (24) with ϕ = 2Δξ /3. For this particular mapping scheme regular paving is obtained in physical coordinates for a non-distorted set of triangles, as shown in Figure 4b (mξ = 3). Two separate equations of the form of Eq. (20) can be used to uniquely deﬁne each Lagrangian marker, l, on the surface element, ε. Two geometrical variables that will be used in the following sections for the construction of the IB forcing ﬁeld, are the surface area assigned to the marker, Alε , and normal outside vector, nˆ lε . As shown in Figure 4a, -in natural coordinates ξ, η- the area associated to each marker, l, is constant and equal to the total area divided by the number of markers, mε . 2.3. Adaptive MLS reconstruction The mapping schemes above enable us to practically decouple the resolution requirements of the Lagrangian and Eulerian grids. In other words, the immersed body can be described by a set of triangles that properly captures the local curvature; when the size of an element is larger than the local Eulerian grid, then the proper number of Lagrangian markers is introduced to preserve the continuity and smoothness of the forcing ﬁeld. All information related to position and kinematics of marker points does not need to be stored in memory, since it can be dynamically recomputed using the above mapping processes with minimal computational overhead. In fact, only knowledge of the kinematics of the vertex nodes in the triangulation, the elements geometrical properties (Aε ,nˆ ε ), 7

1

1

0.5

x2

0

0.5

/ 3

1

/ 3

0

3

1

3

1.5

1

0.5

0

0.5

1

1.5

x1

(a)

(b)

Figure 4: Mapping B, marker points are deﬁned in centroids of evenly distributed triangular patches: (a) Distribution of marker points on the master element and the marker numeration (we use Nξ = 3 for illustration purposes). (b) A regular set of triangles in physical coordinates x1 , x2 with their corresponding mapped Lagrangian points, Nξ = 3.

and how many markers per element side (mξ ) are present, is required to recompute the state of the Lagrangian control points l. Further, the distributed ﬂuid forces on the surface markers might be computed and stored to disk for later post processing, but are not required to be kept in memory, once total or equivalent nodal forces and moments have been computed. The adaptive MLS scheme simply adds an extra step in the process outlined in section 2.1, where step 2 is now modiﬁed as follows: 2 Proceed marker by marker. Apply mapping schemes A or B to generate data of Lagrangian marker (ε, l) (Xε,l , Uε,l , Alε at time level n + 1). For surface marker (ε, l), interpolate U˜ iε,l from surrounding Eulerian stencils using equation (8), and one of the methods to construct Ih [21, 23, 16]. Once the position of a given marker (ε, l) is computed and the associated Eulerian point stencil deﬁned, the intermediate velocities can be interpolated to the marker using the following operator: U˜ iε,l =

ne

k ε,l φε,l u˜ ki k x −X

(25)

k=1

where φε,l k is a set of interpolation functions deﬁned using moving least squares MLS (see [23] for details). Then, Eq. (11) is used to compute the forcing ﬁeld, Fiε,l , with the interpolated intermediate velocities, U˜ iε,l , and the markers velocity components Uiε,l , i = 1, 2, 3, as input. Finally this forcing ﬁeld is transferred back to the Eulerian grid: ε,l k ε,l cε,l φε,l (26) Fi fik = k x −X (ε,l) ⊂ nl

where cε,l are scaling factors deﬁned to ensure conservation of total force and torque during transfer, and nl refers to the set of mapped markers related to the Eulerian point k. In this work, cε,l is computed by requiring that: nte k=1

fik ΔV k =

mε ε⊂ l=1

8

Fiε,l ΔVεl

(27)

where nte is the total number of Eulerian nodes, mε the total number of mapped markers in element ε, ΔV k = (h x hy hz )k the volume of the cell from Eulerian point k, and ΔVεl = Alε hlε . Following [23], the local thickness, hlε , can be deﬁned as follows: ne (h x + hy + hz )k hlε = φε,l . (28) k 3 k=1 We should note that the above deﬁnition for hlε works well for isotropic grids. In the case of highly anisotropic, stretched grids, however, it can lead to forcing that lacks smoothness. To alleviate this issue one can deﬁne hlε as follows: ΔV El , (29) hlε = l (ε,l) ⊂ nl Aε where the denominator is the sum of all areas of the Lagrangian markers associated to the Eulerian node in consideration. Finally the transfer coeﬃcient, cε,l , that satisﬁes Eq. (27) can be deﬁned as: cε,l = where ΔV El =

ne k=1

ΔVεl Al hl = ε Elε , El ΔV ΔV

(30)

k φε,l k ΔV .

2.3.1. Computation of ﬂuid forces acting on the immersed body The integration of ﬂuid forces on the immersed body requires ﬁrst the computation of the distributed pressure and viscous forces on the mapped Lagrangian markers. In this work we follow the normal probe approach described in references [23, 26] to compute these distributed forces on the body. For example, in Cartesian coordinates e ∂ue ∂u fipε,l = −pε,l δi j nε j ; fiμε,l = μ ∂xij + ∂xij nε j (31) Here e refers to a point located at an outside distance hn from marker (ε, l) along the normal direction nˆ lε = [nε1 nε2 nε3 ]T . The indexes i, j = 1, 2, 3 refer to the Cartesian directions. The pressure on the marker (ε, l) is obtained by: pε,l = pe −

∂p ε,l ∂n

∂p ε,l

hn

∂n

˙ ε,l · nˆ lε = −U

(32)

˙ ε,l is computed either from the element nodes using the mapping In the present method the acceleration of the marker U interpolations previously described, or from the analytical point acceleration expression in case of rigid bodies. The viscous force computation involves spatial derivatives of velocities on the external point e, which are used to approximate the velocity gradients on the marker. Then the procedure involves the interpolation of the pressure ﬁeld and velocity derivatives from the Eulerian grid to point e. These interpolations are also done using MLS approximation, in this case around the normal probe point e [23, 26]: ne ∂φek k ∂uei e k e (33) pe = ne k=1 φk x − X pk k=1 ∂x j ui ∂x j = If the velocity derivatives are already available on the Eulerian grid their direct interpolation to the external point e can be performed in similar manner to the ﬁrst of equations (33). Another option to compute the viscous forces is to make use of the equation f με,l = −μnˆ lε × ωe , when the vorticity ﬁeld ω is already available on the Eulerian grid. Finally, the distributed forces on markers can be integrated to obtain equivalent ﬂuid forces and moments on Eulerian axes. In case of rigid bodies, for example: ε pε,l ε ε,l cm + f με,l Alε M f cm = ε⊂ m (34) × f pε,l + f με,l Alε , F f = ε⊂ m l=1 f l=1 X − X where cm stands for the solids center of mass.

9

U(t) = Uo cos( t)

nz D

D

z

y

FD/(6 π ρf ν Uo Ro)

10

5

0

í5

x ny D

nx D

í10 0

(a)

0.5

1

time/T

1.5

2

2.5

(b)

Figure 5: Transversely oscillating sphere: (a) Schematic of the simulation domain and kinematics. (b) Normalized drag force vs. time for Re = 40 and = 4. pressure drag case (i); viscous drag case (i); total drag case (i); • total drag case (ii); ◦ total drag reference computation in [12].

3. Results 3.1. Flow around a transversely oscillating sphere We will utilize this test problem to illustrate the use of both mapping schemes described above, and to stress their importance in maintaining proper marker particle density for accurate IB reconstructions. In particular, we consider an oscillating sphere of diameter, D, in an otherwise quiescent ﬂuid (see Figure 5a). The vertical motion of the sphere is given by the harmonic function U(t) = Uo cos(ωt), where Uo is a reference velocity, and ω is the oscillation angular frequency. The two non-dimensional parameters dominating this ﬂow are the Reynolds and Stokes numbers: 2 , (35) Re = Uνo D ; = ωR 2ν where R = D/2, and ν is the ﬂuids kinematic viscosity. We ﬁrst consider a case with Reynolds number, Re = 40, a non-dimensional amplitude of oscillation, Ao = 1/3.2, and period of oscillation, T = π/1.6, resulting in a Stokes number of, = 4. This conﬁguration has been studied by Mei [12], who employed a truncated Fourier decomposition of the axisymmetric streamfunction-vorticity ﬂuid equations in the frequency domain, and a time dependent ﬁnite diﬀerence method in the time domain. A domain of size 8D × 8D × 12D, in the x, y, z directions was employed with a cell size of Δ x,y = 0.0104D, Δz = 0.016D. A sphere surface triangulation with 2K triangles was used. Slip conditions where applied in the domain boundaries. Two cases are considered: case (i), where Mapping-A with a total of 49K particles is used; case (ii) where Mapping-B with 40K particles is used. We should note that both schemes were found to be insensitive to considering additional forcing for all the interior cells of the sphere. The sensitivity of direct forcing schemes to the addition of forcing inside the object is well known and has been thoroughly discussed in the literature (i.e. [16, 8]). In the following we will present results from case (i) with interior forcing, and from case (ii) without. In Figure 5b the pressure, viscous and total drag force are shown for the case (i) as a function of time. The total drag for case (ii) as well as the reference computation by Mei [12] are also included. The results clearly demonstrate that the accuracy of the predicted drag force does not depend on the mapping scheme employed, or if the interior of the sphere is forced. The comparison with the reference solution by Mei [12] is also very good. Next we consider a case with Re = 100 and = 5 (period of oscillation, T = π). The domain size is reduced to 4D × 4D × 6D in order to perform several runs for uniform and adaptive Eulerian and Lagrangian marker grids at a reasonable cost. Several uniform grid computations were conducted using the mapping scheme A, on Eulerian grids with increasing resolution of 32, 48, 64, 96, 128 cubed cells in the diameter, D, of the sphere. The same surface tessellation as in the previous case was used, and a number of marker particles, ranging from 9K to 80K for the diﬀerent Eulerian grid resolutions, was introduced using the adaptive Lagrangian mapping described above. 10

In Figure 6a, a force comparison between the present computation and the results reported in [26] is shown. The agreement is excellent for both the pressure and viscous components for the same grid resolution (64 cells along D). In Figure 6b the predicted drag force is compared for diﬀerent resolutions. As expected the viscous component converges slower when compared to the pressure contribution. 2

8 7

1.5

6 5

1

4

0.5 CD

CD

3 2

0

1 0

í0.5

í1 í2

128 cells in D 96 cells in D 64 cells in D 48 cells in D 32 cells in D

Viscous

Pressure

í1

í3 0

0.5

1

1.5 time/T

2

2.5

3

í1.5 1

(a)

1.1

1.2

time/T

1.3

1.4

1.5

(b)

Figure 6: Transversely oscillating sphere at Re = 100, = 5. (a) pressure drag; viscous drag; total drag; Symbols are corresponding reference results for the same conﬁguration and grid resolution reported in [26]. Note that force coeﬃcient curves are oﬀset by a factor of 3 for clarity. (b) Grid convergence for Lagrangian forcing and mapping A.

In many cases involving moving boundaries problems, an immersed boundary may cross regions of varying Eulerian grid resolution, or at some particular instance in time, the Eulerian grid may change due to adaptive grid reﬁnement (AMR) [25]. In such case, the adaptive Lagrangian mappings described previously provide the means of increasing (or decreasing) the number of marker particles to maintain the continuity and smoothness of the IB forcing band in a cost eﬃcient manner. To demonstrate this attractive property of the proposed method let us consider the above case of the oscillating sphere on a structured AMR grid with one or two levels of reﬁnement, as shown in Figure 7a. When the simulation starts the sphere moves within a set of grid blocks at the same reﬁnement level (Level 1), and after a certain computational time, tre f , the region of size 2D × 2D × 4D of the domain is reﬁned to Level 2, by doubling the number of blocks in each direction. As the cells per block are constant, the grid reﬁnement is increased by a factor of 2. Computations on grids of reﬁnement Levels 1-2 with 32 - 64 cells per diameter (case amr1), and 48 - 96 cells per diameter (case amr2) were performed. We set tre f = 0.814T , which coincides with a zero crossing of the pressure force. A marginal number of markers was used on Level 1 (3.7K and 9K for cases amr1 and amr2 respectively). After tre f the Eulerian grid was reﬁned and the Lagrangian grid was either kept the same, or dynamically reﬁned using the adaptive Lagrangian mapping described above. This increased the number of markers to approximately 16K and 36K for cases amr1 and amr2 respectively. It was found that if the reﬁnement of the Eulerian grid was not accompanied by local adaptive Lagrangian mapping, the computations became unstable (see Figure 7b left side). The solution obtained with the adaptive Lagrangian mapping on both AMR grids is in excellent agreement with the one on a uniform grid and the MLS forcing scheme as described in [23] (see Figure 7b right side). 3.2. Flow around a sphere Next we will consider the ﬂow around a sphere at Reynolds number, Re = 300, (based on the freestream velocity, U∞ , and the diameter, D). The sphere in this case is immersed in a cylindrical coordinates grid, which is highly anisotropic and complicates the design of the Lagrangian grid. It is composed of 222×66×482 nodes along the radial, azimuthal and axial directions respectively. Near the immersed body the radial and axial grid steps were set at 0.005D (approximately 200 nodes across the diameter), and the grid was stretched away from the body. The inﬂow, outﬂow and lateral boundaries of the domain were placed at 10D, 20D and 10D respectively, from the center of the sphere. 11

2

2 Pressure

CD

1

1

0

0 Viscous

í1 í2 0.7

CD

í1 0.8

0.9

1

í2 0

2

2

1

1

0

0

í1

í1

í2 0.7

0.8

time/T

0.9

(a)

1

0.5

í2 0

0.5

1

time/T

1

(b)

Figure 7: Oscillating sphere on a static AMR grid. (a) The Eulerian grid is reﬁned at tre f to increase the resolution around the sphere from Level 1 (top) to Level 2 (bottom). (b-left) Time evolution of viscous and pressure drag for cases amr1 (top) and amr2 (bottom): with adaptive without adaptive Lagrangian mapping; (b-right) Time evolution of viscous and pressure drag for cases amr1 (top) and Lagrangian mapping; amr2 (bottom) with adaptive Lagrangian mapping (solid line). The dashed line is from a computation on a uniform grid at the same resolution utilizing the forcing scheme in [23].

Table 1: Cd , Cl , and separation point for the ﬂow around a sphere at Re = 300.

Case Cd Cl θs

EF 0.645 0.066 114◦

LFF 0.682 0.064 111◦

LFMF 0.682 0.063 112◦

LFC 1.167 -0.032 -

LFMC 0.689 0.050 109◦

Ref. [4] 0.655 0.065 111◦

Ref. [7] 0.656 0.069 -

Ref. [20] 0.671 -

Two diﬀerent Lagrangian grids were used: Grid 1, composed of about 2 × 105 elements (triangles); Grid 2, composed of about 2 × 104 elements. The former is twice as ﬁne as the coarsest Eulerian cell on the cylindrical coordinate grid, and will produce smooth forcing without adaptive Lagrangian mapping. The latter is an order of magnitude coarser and needs to be adaptively reﬁned to maintain accuracy. The following computations with diﬀerent forcing strategies were conduced: a) Eulerian forcing as in [26] (EF); b) Lagrangian forcing on Grid 1, no-mapping as in [23] (LFF); c) Lagrangian forcing + adaptive Lagrangian mapping on Grid 1 (LFMF); d) Lagrangian forcing on Grid 2, no-mapping as in [23] (LFC); e) Lagrangian forcing + adaptive Lagrangian mapping on Grid 2 (LFMC). 2 2 The resulting drag, Cd = D/0.5ρU∞ πR2 , and lift, Cl = L/0.5ρU∞ πR2 , coeﬃcients (D and L are the drag and lift forces on the sphere, and ρ the ﬂuid density), and the separation angle, θ s , are given in Table 1. Reference results from the computations of Constantinescu & Squires [4], Tomboulides [20], and Johnson & Patel [7] are also included. For all forcing schemes with the exception of case LFC (Lagrangian forcing on the coarse surface grid without the adaptive reﬁnement), the predicted Cd , Cl and θ s are within 5% of each other and the reference data. This indicates that: i) the ﬁne surface grid (Grid 1) is suﬃcient to allow the Lagrangian type forcing proposed in [23] (LFF) to match the Eulerian type of forcing in [26] (EF); ii) using Grid 1 the results with and without the adaptive strategy (LFF & 12

(a)

(b)

Figure 8: Flow around a sphere at Re = 300. (a) Pressure coeﬃcient; (b) Skin-friction coeﬃcient. Both are averaged in time and along the azimuthal direction. EF, LFF, LFMF, LFC, LFMC, • body-ﬁtted results by Lee [10].

(a)

(b)

Figure 9: Flow around a sphere at Re = 300. Velocity magnitude averaged in time and along the azimuthal direction. (a) case LFC; (b) case LFMC The scale ranges from 0 (blue) to 1.125 (red) and is normalized with U∞ .

LFMF) are the same, which is a consistency check for the proposed formulation; iii) when the coarse surface grid is used (Grid 2) the mapping methodology (LFMC) maintains the accuracy of the forcing. We should note that the adaptive strategy was set up to produce an equivalent Lagrangian grid which is twice ﬁner than the local Eulerian one. For the case of surface Grid 2, for example, this requirement introduced between 6 and 36 markers per element depending on the radial location. 2 2 Fig. 8, shows the time and azimuthally averaged pressure, C p = 2(p − p∞ )/ρU∞ , and skin-friction, C f = 2τw /ρU∞ coeﬃcients, (where p∞ is the freestream pressure and τw the wall-stress). The results from the computations by Lee [10], where a body ﬁtted code was utilized, are also included for comparison. The agreement with the reference solution is excellent and the results from the diﬀerent cases are consistant with the force coeﬃcients reported above. Comparing cases LFC and LFMC one can better quantify the impact the quality of the surface Lagrangian grid has on the solution. When such a grid does not meet the requirement to be everywhere ﬁner than the underlying Eulerian grid (case LFC), then the accuracy of the solution degrades very fast, as it can be seen in the distribution of C p , and C f . This is primarily due to the excessive transpiration that occurs through the sphere for case LFC. This is clearly demonstrated in Figure 9, where the average velocity magnitude for cases LFC and LFMC is shown. These results clearly demonstrate that the adaptive strategy restores the accuracy of the force distribution and properly enforces the non-slip conditions independently of the quality of the surface grid. 3.3. Flow around an Eppler 387 airfoil The ﬂow around an Eppler 387 airfoil at Reynolds number, Re = cU∞ /ν = 30, 000 (where c, is the chord length and, U∞ , the free-stream velocity) is considered. Spedding & McArthur [19] reported detailed experiments in this Reynolds number regime for a wide range of angles of attack. Here we report large-eddy simulations (LES), of a case where the angle of attack is set to α = 10◦ . At this value of α they reported a notable jump of the lift coeﬃcient, without a signiﬁcant change in drag coeﬃcient, due to the complex dynamics of the boundary and shear layers at the leading edge of the airfoil. The latter makes this case vary challenging, especially for immersed boundary methods. 13

Table 2: Cd and Cl on the Eppler 387 airfoil at Re = 30, 000 and α = 10◦

Cd Cl

Experiment 0.0860 1.1616

LF 16.14% 18.99%

LFM1 4.87% 11.80%

(a)

LFM2 −2.52% 0.67%

LFM3 −1.90% 0.08%

(b)

Figure 10: Distributions of pressure (a) and skin-friction (b) coeﬃcients over the Eppler 387 airfoil at Re = 30, 000 and angle of attack α = 10◦ . Averages in time and along the spanwise direction.

The airfoil was always represented by 87, 000 triangles, and the following simulations were carried out: i) LF, Lagrangian forcing with no mapping, resulting in a Lagrangian grid approximately 10 times coarser than the Eulerian one, ΔL /ΔE ∼ 10; LFM1, Lagrangian forcing + mapping requiring the Lagrangian and Eulerian grids to have comparable resolution, ΔL /ΔE ∼ 1 (2.9 million markers introduced); LFM2, Lagrangian forcing + mapping, ΔL /ΔE ∼ 0.66, (6.4 million markers introduced); LFM3 Lagrangian forcing + mapping, ΔL /ΔE ∼ 0.5 (12.8 million markers introcuded). The Cartesian Eulerian grid was the same for each simulation and it was composed of 518 × 202 × 2002 nodes along the cross-stream (x), spanwise (y) and streamwise (z) directions, respectively. The grid step was equivalent to 0.001c along the streamwise and cross-stream directions in the vicinity of the airfoil, and it was stretched away from the body. Along the spanwise direction a grid step of 0.01c was adopted. The computational domain extended 10c upstream of the center of the chord and 20c downstream, 20c along the cross-stream direction and 2c along the span of the airfoil. A Dirichlet condition was enforced at the inﬂow, a convective condition at the outﬂow [14], free-stream at the top and bottom boundaries and periodic conditions along the spanwise direction. In Table 2 we report the relative errors for drag and lift coeﬃcients, in comparison to the experimental values. As expected, case LF, where no mapping is used, results in large errors in the prediction of Cd and Cl . In case LFM1 the mapping improves the results but the resolution of the Lagrangian grid is still marginal. Cases LFM2 and LFM3 are both within 2.5% of the measurements, demonstrating convergence of the numerical solution with number of Lagrangian markers. These trends are conﬁrmed in the detailed distributions of the time-spanwise-averaged pressure and skin-friction coeﬃcients shown in Figure 10 . Further details are reported in Figure 11, where the time-spanwise-averaged proﬁles of the tangential velocity are plotted along the normal direction on the suction side of the airfoil (see Figure 12 for the actual locations). In all cases a small recirculation bubble is predicted in the vicinity of the leading edge. For cases LF and LFM1, however, it is smaller than the one predicted by LFM2 and LFM3. The coarser Lagrangian grids also result in a boundary layer thickening faster than the one predicted by LFM2 and LFM3, probably a result of transpiration through the airfoil. 4. Summary In this work we propose the use of a marker paving algorithm based on isoparametric mapping, which facilitates local reﬁnement on the Lagrangian mesh used to deﬁne an immersed body in a ﬂuid ﬂow. The reﬁnement level is tied to the local Eulerian mesh in the vicinity of the body, which may be changing in space and time. Although the 14

(a)

(b)

(c)

(d)

(e)

(f)

Figure 11: Proﬁles of tangential velocity component along directions normal to the suction side of the Eppler 387 airfoil (see Fig. 12) at Re = 30, 000 and angle of attack α = 10◦ . Averages in time and along the spanwise direction.

(a)

(b)

Figure 12: Spanwise vorticity ﬁeld: instantaneous (a) and time-spanwise-averaged (b). Scale ranging from -25 (blue) to 25 (red). Solution by LFM3. Also the normal directions relative to the proﬁles in Fig. 11 are represented.

15

proposed formulation is discussed in the context of triangular surface elements, the method can be extended directly to other polygonal elements with their corresponding polynomial interpolation functions. Overall the method is very eﬃcient, as little additional information needs to be kept in memory, besides the nodes and connectivity information of the underlying triangulation, albeit the small expense of recreating the markers per surface element when needed. We should also note that the method is very robust and guaranties the continuity of the forcing independently of the quality and the Lagrangian or Eulerian grids. However, in situations with low quality Lagrangian grids, for example, (i.e. a mesh deﬁned by skewed triangles) the computational eﬃciency is reduced. This is because in both mapping schemes we present we use the same number of points in each local direction, which is selected based on the largest Lagrangian triangle edge and smallest Eulerian grid size, generating markers also in places that may not be needed. In a similar manner an excess number of markers can be generated in situations where large Lagrangian elements on the immersed body go through highly anisotropic Eulerian grids. In all cases we have examined, however, the increase in the computational cost is negligible. Anisotropic mapping and local sub-triangulation will be helpful to overcome such issues in the future. The results for a series of canonical problems demonstrate the accuracy and robustness of the method by comparing to results in the literature, as well as, to Eulerian IB schemes, which are generally not very sensitive to the selected surface grid. We should note that the numerical characteristics of the underlying IB method are not modiﬁed when marker paving is used. Acknowledgments: This research was supported by ONR Grant N000141110455, monitored by Dr K.-H. Kim, and the U.S. Department of Energy, Nuclear Energy University Program (DOE-NEUP) grant. AP is partially supported by EU-H2020 HOLISHIP - GA 689074 grant. Computational time is provided through the DoD High Performance Computing Modernization Program, and an ALCC allocation from the U.S. Department of Energy. References [1] Akselvoll, K., and Moin, P. An Eﬃcient Method for Temporal Integration of the Navier-Stokes Equations in Conﬁned Axisymmetric Geometries. Journal of Computational Physics 125, 2 (May 1996), 454–463. [2] Balaras, E. Modeling complex boundaries using an external force ﬁeld on ﬁxed Cartesian grids in large-eddy simulations. Computers and Fluids 33, 3 (2004), 375–404. [3] Bathe, K. J. Finite Element Procedures. Klaus-Jurgen Bathe, 2007. [4] Constantinescu, G., and Squires, K. LES and DES investigations of turbulent ﬂow over a sphere at Re=10,000. Flow Turbulence And Combustion 70, 1-4 (2003), 267–298. [5] Fadlun, E., Verzicco, R., Orlandi, P., and Mohd-Yusof, J. Combined immersed-boundary ﬁnite-diﬀerence methods for three-dimensional complex ﬂow simulations. Journal of Computational Physics 161, 1 (2000), 35–60. [6] Hughes, T. J. R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications, 2000. [7] Johnson, T., and Patel, V. Flow past a sphere up to a Reynolds number of 300. Journal Of Fluid Mechanics 378 (1999), 19–70. [8] Kempe, T., and Frohlich, J. An improved immersed boundary method with direct forcing for the simulation of particle laden ﬂows. Journal of Computational Physics 231, 9 (2012), 3663 – 3684. [9] Kim, J., Kim, D., and Choi, H. An immersed-boundary ﬁnite-volume method for simulations of ﬂow in complex geometries. Journal of Computational Physics 171, 1 (2001), 132–150. [10] Lee, S. A numerical study of the unsteady wake behind a sphere in a uniform ﬂow at moderate Reynolds numbers. Computers & Fluids 29, 6 (2000), 639–667. [11] Liu, W. K., Jun, S., and Zhang, Y. F. Reproducing kernel particle methods. International Journal For Numerical Methods In Fluids 20, 8-9 (Apr. 1995), 1081–1106. [12] Mei, R. Flow due to an oscillating sphere and an expression for unsteady drag on the sphere at ﬁnite reynolds number. Journal of Fluid Mechanics 270 (7 1994), 133–174. [13] Mittal, R., and Iaccarino, G. Immersed boundary methods. Annual Review Of Fluid Mechanics 37 (2005), 239–261. [14] Orlanski, I. A simple boundary condition for unbounded hyperbolic ﬂows. Journal of Computational Physics 21, 3 (1976), 251–269. [15] Peskin, C. S. Flow patterns around heart valves - numerical method. Journal of Computational Physics 10, 2 (1972), 252–&. [16] Pinelli, A., Naqavi, I. Z., Piomelli, U., and Favier, J. Immersed-boundary methods for general ﬁnite-diﬀerence and ﬁnite-volume NavierStokes solvers. Journal of Computational Physics 229, 24 (2010), 9073–9091. [17] Roma, A. M., Peskin, C. S., and Berger, M. J. An adaptive version of the immersed boundary method. Journal of Computational Physics 153, 2 (1999), 509–534. [18] Seo, J. H., and Mittal, R. A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations. Journal of Computational Physics 230, 19 (2011), 7347 – 7363. [19] Spedding, G., and McArthur, J. Span eﬃciencies of wings at low Reynolds numbers. Journal of Aircraft 47, 1 (2010), 120–128. [20] Tomboulides, A. G. Direct and large-eddy simulation of wake ﬂows: ﬂow past a sphere. PhD thesis, Princeton University, Princeton, 1993. [21] Uhlmann, M. An immersed boundary method with direct forcing for the simulation of particulate ﬂows. Journal of Computational Physics 209, 2 (2005), 448–476.

16

[22] Van Kan, J. A second-order accurate pressure-correction scheme for viscous incompressible ﬂow. SIAM Journal on Scientiﬁc and Statistical Computing 7, 3 (1986), 870–891. [23] Vanella, M., and Balaras, E. A moving-least-squares reconstruction for embedded-boundary formulations. Journal of Computational Physics 228, 18 (June 2009), 6617–6628. [24] Vanella, M., Posa, A., and Balaras, E. Adaptive mesh reﬁnement for immersed boundary methods. ASME J. Fluids Eng. 136, 4 (2014), 040909. [25] Vanella, M., Rabenold, P., and Balaras, E. A direct-forcing embedded-boundary method with adaptive mesh reﬁnement for ﬂuid-structure interaction problems. J. Comput. Phys. 229, 18 (2010), 6427–6449. [26] Yang, J., and Balaras, E. An embedded-boundary formulation for large-eddy simulation of turbulent ﬂows interacting with moving boundaries. Journal of Computational Physics 215, 1 (2006), 12–40. [27] Yang, J., and Stern, F. A non-iterative direct forcing immersed boundary method for strongly-coupled ﬂuid–solid interactions. Journal of Computational Physics 295 (2015), 779 – 804.

17

Copyright © 2021 COEK.INFO. All rights reserved.