Gridless spectral algorithm for Stokes flow with moving boundaries

Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Gridless spectral algorithm for Stokes flow with moving boundaries S.Z. Husain, J.M. Floryan * Department of Mechanical and Materials Engineering, The University of Western Ontario, London, Ontario, Canada N6A 5B9

a r t i c l e

i n f o

Article history: Received 27 December 2007 Received in revised form 20 July 2008 Accepted 23 July 2008 Available online 5 August 2008 Keywords: Gridless methods Immersed boundary conditions Moving boundary problems Spectral accuracy

a b s t r a c t A gridless, spectrally-accurate algorithm for the Stokes flow with moving boundaries is presented. The algorithm uses fixed computational domain with boundaries of the flow domain moving inside the computational domain. The spatial discretization is based on the Fourier expansions in the streamwise direction and the Chebyshev expansions in the transverse direction. Temporal discretization uses one- and two-steps implicit formulations. The boundary conditions on the moving boundaries are imposed using the immersed boundary conditions concept. Numerical tests confirm the spectral accuracy in space and theoretically-predicted accuracy in time. Different variants of the solution procedure are presented and their relative advantages are discussed. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction The first step in the most common approach used for simulations of flows in complex geometries involves numerical modeling of geometry of the flow domain. A model of the domain is constructed using grid generator and flow boundary conditions are imposed on the edges of this domain, i.e., the flow domain overlaps with the computational domain. This process, which is frequently referred to as pre-processing, could consume majority of time required for preparation of the simulation. The flow boundary conditions are enforced exactly and one is concerned with various errors associated with spatial and temporal discretization of the field equations. Complications occur when geometry of the flow domain changes as a function of time, e.g., when one deals with the moving boundary problems. In the present work, we are interested in finding a spectrally-accurate and computationally efficient solution to such problem while maintaining a sharp resolution of the location of the moving boundary. Techniques developed for handling of the moving boundary problems have been reviewed in [1] with updates given in [2]. Such techniques can be classified as Eulerian, Lagrangian and mixed. Eulerian algorithms, which are of interest in the present work, rely on a coordinate system that is either stationary in a laboratory frame of reference or moves in a prescribed manner (i.e., describable by the Galileo transformation). These algorithms can be divided for convenience into fixed grid methods, adaptive grid methods and various mapping methods.

* Corresponding author. Fax: +1 519 661 3020. E-mail address: mfl[email protected] (J.M. Floryan). 0045-7825/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2008.07.007

The adaptive grid methods adjust the grid at each time step so that one of the grid lines follows the moving boundary. The adjustment is done numerically through appropriate mappings. The process is conceptually identical regardless of whether one uses numerical grid generation combined with finite-difference discretization of the field equations or a grid generation scheme for the finite-element discretization of the same equations. These methods have high computational overhead as the grid needs to be reconstructed at each time step followed by solution of the field equations. In the case of the problem discussed in [3] and solved using finite volume discretization, the cost of grid construction reached the level of about 75% of the total cost of the computation. The boundary conditions are imposed exactly and thus the overall error includes contributions from the grid generation and the spatial and temporal discretizations of the field equations. Numerous problems occur when one is interested in high accuracy as this requires both high accuracy in grid generation and high accuracy discretization. One can reduce the cost of computations and improve accuracy by using analytical transformations that map the irregular flow domain into a regular computational domain at the cost of increased complexity of the field equations (see [4] for an example). Such mappings are available for a limited class of geometries [1]. In the fixed grid methods the boundary of the flow domain is allowed to move through a fixed grid and this eliminates the high computational cost of grid re-adjustment. Various methods for modeling of movements of boundaries based either on the surface or volume tracking procedures are discussed in [1,5]. In the former approach, a set of points is introduced to mark location of the boundary which is represented as a set of interpolated curves. These points are moved in a prescribed manner during the solution

246

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

process [5,6] and provide ability for precise identification of boundary location. In the latter approach, the information about the boundary location is not stored but the boundary is reconstructed whenever necessary on the basis of the presence of certain quantity of a convenient marker within computational cell, e.g., MAC – Marker and Cell [7], VOF – Volume of Fluid [2,8] and Level Set [9,10] methods. The reconstruction process leads to a diffused location of the boundary but the overall accuracy can be made consistent with the low accuracy discretization schemes used for the field equations. Lagrangian methods are characterized by a coordinate system that moves with the fluid. Each computational cell always contains the same fluid and its tracking requires solution of an initial value problem. These methods are well suited for moving boundary problems as they permit simple delineation of moving boundaries. The two main problems involve mesh tangling and loss of numerical accuracy associated with highly distorted meshes. Possible remedies are discussed in [1]. Mixed Lagrangian–Eulerian methods rely on the combination of concepts described above [1]. The available methods have low order accuracy and thus one is interested in the development of techniques that could overcome this limitation. The overall error is affected by the accuracy of representation of boundary location and by the accuracy associated with spatial discretization of the field equations. Spectral methods provide lowest error for the field equations but are limited to solution domains with regular geometries. The variability and complexity of geometry of the flow domain represents the main challenge for implementation of these methods, however, the use of the concept of immersed boundary conditions provides a way to combine both, i.e., the fixed, regular solution domain with a complex, time-dependent flow domain [11]. The immersed boundary conditions concept has been proposed in [12] in the context of simulation of cardiac dynamics. Boundary conditions at the edges of the physical domain are replaced by constraints imposed inside the computational domain but these constraints have been constructed using physical arguments through introduction of forcing functions that reproduce the effect of the boundary [13,14]. In this sense, the method is analogous to the fixed grid methods discussed above. The first spectrally-accurate implementation for fixed boundary problems has been given in [11] with the constraints constructed in a formal manner by imposing conditions for elimination of certain terms in the spectral expansions representing flow boundary conditions. This implementation is limited to geometries that can be represented by Fourier expansions but results in a gridless algorithm as all possible variations of boundary geometries are described in terms of the Fourier coefficients only. The programming effort associated with modeling of changes of geometry has been essentially removed as the only information required for specifying the new geometry is reduced to a set of Fourier coefficients provided as input to the code. This implementation has been extended to hydrodynamics stability problems [15] and applied successfully in the analysis of the effect of surface roughness [16]. The present work is focused on the development of spectrallyaccurate algorithms suitable for solution of time-dependent flow problems, including moving boundary problems. Such algorithms are of interest in the development of active flow control and management techniques using micro-electro-mechanical devices. As our interest is in flows through micro/nano channels, we focus our attention on the small Reynolds number limit, i.e., Stokes flow. Simpler problems, i.e., unsteady conduction problems in a corrugated channels have been considered earlier in [17,18] without encountering numerical instabilities. The model problem for the present work is described in Section 2. Spatial and temporal discretizations are discussed in Section 3. Results of numerical tests are discussed in Section 4. Section 5 provides a short summary of the main conclusions.

2. Model problem The model problem consists of an unsteady Stokes flow in a channel bounded by corrugated walls with geometry of the corrugation changing as a function of time (see Fig. 1). These geometries as well as motions of the walls are described by the following relations

yL ðx; tÞ ¼ 1 þ

n¼þ1 X

ðnÞ

HL ðtÞeinax ;

ð2:1aÞ

ðnÞ

ð2:1bÞ

n¼1

yU ðx; tÞ ¼ 1 þ

n¼þ1 X

HU ðtÞeinax ;

n¼1 ðnÞ

ðnÞ

ðnÞ

ðnÞ

where HL ¼ HL and HU ¼ HU are known and stars denote complex conjugates. The channel geometry is periodic in x with wavelength k = 2p/a and extends to ±1 in the x-direction. The dimensionless field equation has the form

o 2  r W ¼ r2 ðr2 WÞ; ot

ð2:2aÞ

where W denotes the stream function, which is defined as

uðx; y; tÞ ¼ u0 ðyÞ þ u1 ðx; y; tÞ ¼ dW0 =dy þ oW=oy ¼ oWT =oy; ð2:2bÞ vðx; y; tÞ ¼ v1 ðx; y; tÞ ¼ oW=ox ¼ oWT =ox; 2

2

2

2

ð2:2cÞ

2

and r = o /ox + o /oy . Here, u (x, y, t) and v (x, y, t) denote the total velocities in x- and y-directions, respectively, u1 (x, y, t) and v1 (x, y, t) denote velocity modifications due to the presence of boundary motion, the Poiseuille flow u0 (y) = 1  y2 is taken as the reference flow, W0 = y3/3 + y + 2/3 denotes the stream function of this flow and WT stand for the stream function of the complete flow (total stream function). The problem formulation needs to be supplemented with suitable initial conditions, which are taken to be in the form

uðx; y; 0Þ ¼ ui ðx; yÞ; vðx; y; 0Þ ¼ vi ðx; yÞ;

ð2:3aÞ ð2:3bÞ

yL ðx; 0Þ ¼ yLi ðxÞ;

ð2:3cÞ

yU ðx; 0Þ ¼ yUi ðxÞ;

ð2:3dÞ

where ui (x, y), vi (x, y), yLi (x), yUi(x) are considered to be known, and boundary conditions at the solid walls, which are taken to be in the form

u0 ðyL ðx; tÞÞ þ u1 ðx; yL ðx; tÞ; tÞ ¼ uL ðx; tÞ ¼ 0;

ð2:4aÞ

u0 ðyU ðx; tÞÞ þ u1 ðx; yU ðx; tÞ; tÞ ¼ uU ðx; tÞ ¼ 0; n¼þ1 X dHðnÞ L v1 ðx; yL ðx; tÞ; tÞ ¼ vL ðx; tÞ ¼ dyL =dt ¼ einax ; dt n¼1

ð2:4bÞ

v1 ðx; yU ðx; tÞ; tÞ ¼ vU ðx; tÞ ¼ dyU =dt ¼

n¼þ1 X n¼1

ð2:4cÞ

ðnÞ

dHU inax e : dt

ð2:4dÞ

y 1

yU(x,t)

YU x

0

YL

yL(x,t)

-1

Fig. 1. Sketch of the instantaneous form of the flow domain.

247

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

3. Discretization method We are interested in the determination of solution of (2.2)–(2.4) with spectral accuracy. The main difficulty associated with the implementation of spectral discretization arises due to the irregularity and time-dependence of the solution domain. In order to overcome this problem, we select fixed rectangular computational domain extending over one period in the x-direction and extending sufficiently far in the y-direction so that the flow domain remains always immersed inside the computational domain during the time interval of interest. If we denote the locations of extremities in the shape of the walls as YU and YL, the y-extent of the computational domain is set as (1  YL, 1 + YU) without loss of generality. The spatial discretization is based on the use of Fourier expansions in the x-direction due to periodicity of the geometry, and on expansions in terms of Chebyshev polynomials in the y-direction. We shall use standard definition of the Chebyshev polynomials and thus the y-extent of the computational domain needs to be mapped onto (1, 1) (see Fig. 1) before calculations can proceed. The required mapping has the form

^ ¼ ½y  ð1 þ Y U ÞC þ 1; y

ð3:1Þ

^ 2 h1; 1i, C = 2/(2 + YU + YL) is a constant and the governwhere y ing equation transforms into

"

2

2

#

"

4

2

4

#

o o o o o o þ C2 2 W ¼ þ 2C2 2 2 þ C4 4 W: ^ ^ ^ ot ox2 oy ox4 ox oy oy

ð3:2Þ

n¼þ1 X

n¼þ1 X

ð3:3aÞ

ðnÞ

AU ðtÞeinax ;

ð3:3bÞ

n¼1 ð0Þ

ð0Þ

ðnÞ

ðnÞ

where AL ðtÞ ¼ 1 þ C½2  Y U þ HL ðtÞ, AL ðtÞ ¼ CHL ðtÞ for n – 0, ð0Þ ð0Þ ðnÞ ðnÞ AU ðtÞ ¼ 1 þ C½Y U þ HU ðtÞ, AU ðtÞ ¼ CHU ðtÞ for n – 0. The boundary conditions at the transformed boundaries become

^L ðx; tÞÞ þ u1 ðx; y ^L ðx; tÞ; tÞ ¼ 0; u0 ðy ^U ðx; tÞÞ þ u1 ðx; y ^U ðx; tÞ; tÞ ¼ 0; u0 ðy

ð3:4aÞ ð3:4bÞ

^L ðx; tÞ; tÞ ¼ vL ðx; tÞ ¼ C1 dy ^L =dt ¼ C1 v1 ðx; y

n¼þ1 X n¼1

1

^U ðx; tÞ; tÞ ¼ vU ðx; tÞ ¼ C dy ^U =dt ¼ C v1 ðx; y

1

ðnÞ

dAL inax e ; ð3:4cÞ dt

n¼þ1 X n¼1

ðnÞ

dAU inax e : dt ð3:4dÞ

The solution can be represented in the form of Fourier expansion

Wðx; y^; tÞ ¼

n¼þ1 X

UðnÞ ðy^; tÞeinax 

n¼1

n¼þN XM

UðnÞ ðy^; tÞeinax ;

ð3:5Þ

n¼N M H

^; tÞ ¼ UðnÞ ðy ^; tÞ and star denotes complex conjugate. where UðnÞ ðy Substitution of (3.5) into the field Eq. (3.2) and separation of Fourier components lead to an uncoupled system of partial differential equations for U(n), n 2 h0, NMi, of the type

i h i oh 2 2 C D  ðnaÞ2 UðnÞ ¼ C4 D4  2C2 ðnaÞ2 D2 þ ðnaÞ4 UðnÞ ; ot

ð3:6Þ

^. Two types of temporal discretizations have been where D ¼ d=dy used. The two-step implicit method results in the following relations

h

   i C D þ 2n2 a2 C2  1:5C2 Dt 1 D2 þ n4 a4 þ 1:5n2 a2 Dt1 UðnÞ sþ1 h i h i ðnÞ 2 2 2 2 2 2 ðnÞ 1 1 ¼ 2Dt C D  ðnaÞ Us þ 0:5Dt C D  ðnaÞ Us1 ; 4

4

n 2 h0; NM i;

k¼1 X

ðnÞ

^Þ  Gk;sþ1 T k ðy

k¼0

k¼N XT

ðnÞ

^Þ; Gk;sþ1 T k ðy

ð3:8Þ

k¼0 ðnÞ

where Tk denotes the Chebyshev polynomial of kth order and Gk;sþ1 denotes the unknown coefficients of the expansion. Substitution of (3.8) into (3.7) gives

h







C4 D4 þ 2n2 a2 C2  1:5C2 Dt 1 D2 þ n4 a4 þ 1:5n2 a2 Dt 1 k¼N XT

i

h i k¼N XT ðnÞ ðnÞ Gk;sþ1 T k ¼ 2Dt 1 C2 D2  ðnaÞ2 Gk;s T k

k¼0

k¼0

h

2

þ 0:5Dt 1 C2 D  ðnaÞ

2

i k¼N XT

ðnÞ Gk;s1 T k :

ð3:9Þ

k¼0

Galerkin procedure is used to develop algebraic equations for the ðnÞ ^Þ and of (3.9) by T j ðy unknowns Gk;sþ1 , i.e., we multiply both sides pffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ^ ¼ 1= 1  y ^2 to get integrate with the weight function x

h





C4 hT j ; D4 T k i þ 2n2 a2 C2  1:5C2 Dt1 hT j ; D2 T k i

k¼0

ðnÞ AL ðtÞeinax ;

n¼1

^U ðx; tÞ ¼ y

^ UðnÞ sþ1 ðyÞ ¼

k¼N XT

^Þ plane are given as Locations of the walls in the ðx; y

^L ðx; tÞ ¼ y

while similar relations resulting from the one-step, self-starting implicit method are shown in the Appendix A. In the above the subscript s labels the time step and Dt stands for the (constant) length of the time step. Relation (3.7) has the form of an inhomogeneous orðnÞ dinary differential equation for Usþ1 . The following discussion will be carried out in the context of the two-step method, while the relevant relations for the one-step method can be readily deduced. ðnÞ The unknown function Usþ1 can be represented in terms of expansions based on Chebyshev polynomials in the form

ð3:7Þ

   ðnÞ þ n4 a4 þ 1:5n2 a2 Dt1 hT j ; T k i Gk;sþ1 k¼N i XT h 2 ¼ C hT j ; D2 T k i  ðnaÞ2 hT j ; T k i k¼0

h i ðnÞ ðnÞ j 2 h0; NT  4i;  2Dt 1 Gk;s þ 0:5Dt1 Gk;s1

ð3:10Þ

^Þ; g k ðy ^Þi ¼ where the inner product is defined as hfj ðy R1 ^ ^ ^ ^ ^ f ð y Þg ð y Þ x ð y Þd y . Eq. (3.10) leads to N  3 decoupled algebraic T k 1 j equations for each Fourier mode. Four additional equations required in order to close the system need to be derived from the boundary conditions and these conditions provide coupling between different Fourier modes. The boundary conditions are to be enforced along the lines ^U ðx; s þ 1Þ that are inside the computational do^L ðx; s þ 1Þ and y y ^ 2 h1; 1i. This main while the solution domain remains fixed at y can be accomplished using the concept of immersed boundary conditions. We begin description of the implementation of boundary conditions by noting that at time s + 1 one needs to evaluate velocity components ul (x, s + 1)  u (x, f(x,s + 1),s + 1) and vl (x, s + 1)  v (x, f(x, s + 1), s + 1) along an arbitrary time dependent line l that at ^ ¼ f ðx; s þ 1Þ, such that f is a peritime t = s + 1 occupies position y odic function with period k = 2p/a and jf (x, s + 1)j 6 1. The function f (x, s + 1) can be expressed without loss of generality as

f ðx; s þ 1Þ ¼

n¼1 X n¼1

ðnÞ

P sþ1 einax 

n¼þN XA

ðnÞ

Psþ1 einax

ð3:11Þ

n¼NA

where in the calculations the shape of the wall is approximated by the finite number of terms NA. The velocity components ul(x,s + 1) and vl(x,s + 1) are periodic in x with the same period k and thus can be expressed in terms of Fourier series as

ul ðx; s þ 1Þ  uðx; f ðx; s þ 1Þ; s þ 1Þ ¼

n¼þN Xh

ðnÞ

ð3:12aÞ

ðnÞ

ð3:12bÞ

U sþ1 einax ;

n¼N h

vl ðx; s þ 1Þ  vðx; f ðx; s þ 1Þ; s þ 1Þ ¼

n¼þN Xh n¼N h

V sþ1 einax :

248

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

It will become obvious from the follow up discussion that the summation extends to Nh = NTNA + NM. When f (x, s + 1) overlaps with the position of the wall, ul and vl are known and the coefficients in the above expansions can be determined. Since the flow representation used in the computations is limited to NM + 1 modes (see Eq. (3.5)), only on the first (NM + 1) terms in (3.12) can be accounted for. We shall now discuss the method used to enforce these conditions. Velocity components ul and vl can be evaluated along the wall at time s+1 using the discretized form of the solution, i.e.,

ul ðx; s þ 1Þ ¼ u0 ðf ðx; s þ 1ÞÞ þ C

n¼þN XM

DUðnÞ ðf ðx; s þ 1Þ; s þ 1Þeinax

 oWT ^ dy ^ sþ1 oy ^L y ^U ðx; s þ 1Þ; s þ 1Þ  WT ðx; y ^L ðx; s þ 1Þ; s þ 1Þ: ¼ WT ðx; y

n¼N M

¼ u0 ðf ðx; s þ 1ÞÞ þ C

n¼þN XM k¼N XT

and are enforced for jnj 6 Nm. Relations obtained for jnj > Nm can be used as a measure of error in the enforcement of flow boundary conditions. Similar boundary conditions can be written for the upper wall by replacing subscript L with subscript U. The above problem specification is incomplete due to the shortage of boundary conditions for the mode zero in Eq. (3.18b). While these conditions can be selected arbitrarily, we shall focus our attention on the fixed mass flow rate case, i.e., we shall require that the net mass flow rate in the x-direction is maintained constant during motion of the walls. Volume flux Q along the channel can be evaluated by integrating the x-velocity component across the channel, i.e.

Qðx; s þ 1Þ ¼

ðnÞ

Gk;sþ1 DT k ðf ðx; s þ 1ÞÞeinax ;

n¼NM k¼0

Z

^U y

ð3:13aÞ vl ðx; s þ 1Þ ¼ 

n¼þN XM

ð3:19Þ The volume flux represents an x-periodic function that can be written in the form of a Fourier expansion

inaUðnÞ ðf ðx; s þ 1Þ; s þ 1Þeinax

n¼NM

¼

n¼þN XM k¼N XT

ðnÞ

inaGk;sþ1 T k ðf ðx; s þ 1ÞÞeinax :

ð3:13bÞ

n¼NM k¼0

Chebyshev polynomials and their derivatives evaluated at the wall, i.e., Tk (f(x, s + 1)) and DTk (f(x, s + 1)), are periodic functions of x and thus can be expressed in terms of Fourier expansion as follows:

T k ðf ðx; s þ 1ÞÞ ¼

m¼þ1 X

ðmÞ

wk;sþ1 eimax ;

ð3:14aÞ

m¼1

DT k ðf ðx; s þ 1ÞÞ ¼

m¼þ1 X

ðmÞ

dk;sþ1 eimax ;

ð3:14bÞ

m¼1 ðmÞ

ðmÞ

where the method for evaluation of coefficients wk;sþ1 and dk;sþ1 is explained in Appendix B. Substitution of (3.14) into (3.13) gives

ul ðx; s þ 1Þ ¼

n¼þ1 X

F ðnÞ ðs þ 1Þeinax þ C

n¼1

n¼þ1 X m¼þN X M k¼N XT

ðmÞ

ðnmÞ

imaGk;sþ1 wk;sþ1 einax ;

ð3:15bÞ

n¼1 m¼NM k¼0

where

u0 ðf ðx; s þ 1ÞÞ ¼

n¼þ1 X

F ðnÞ ðs þ 1Þeinax :

ð3:16Þ

n¼1

ðnÞ

m¼þN X M k¼N XT

ðmÞ

ðnmÞ

Gk;sþ1 dk;sþ1 ;

¼

m¼þN XM

k X

where the zero term, i.e., Q(0), represents the net mass flux along the channel. Evaluation of this flux requires knowledge of WT at both walls. The evaluation process is very similar for both walls and thus we shall limit description to the lower wall only. Values of WT at the lower wall can be evaluated by considering velocities of material points located on the line l overlapping with this wall. These velocities are known and can be expressed as

ul ðx; s þ 1Þ ¼ 0; k¼þN ðnÞ XA ^L ðx; s þ 1Þ dy dy dAL vl ðx; s þ 1Þ ¼ ¼ C1 ^ dt dy dt n¼N ;n–0 A

ð3:21aÞ

!

einax : sþ1

The Fourier expansion on the right hand side of (3.21b) does not contain mode zero as the mean position of the wall is assumed to be independent of time. The rate of change of WT along the wall can be written as

dWT oWT ^L ðx; s þ 1ÞÞ ¼ ^L ðx; s þ 1ÞÞ ðx; y ðx; y dx ox ^ oWT dy ^L ðx; s þ 1ÞÞ L : þC ðx; y ^ oy dx

ð3:17aÞ

WT ðx; y^L ðx; s þ 1ÞÞ ¼ 

m¼NM k¼0 ðnÞ V sþ1

ð3:20Þ

ð3:22Þ

The first and second terms of the right hand side are replaced with (3.21b) and (3.21a), respectively, and the resulting expression is integrated in x resulting in

Comparison of (3.12) with (3.15) gives

U sþ1 ¼ F ðnÞ ðs þ 1Þ þ C

Q ðnÞ ðs þ 1Þeinax ;

n¼N M

ðnmÞ

n¼1 m¼NM k¼0

n¼þ1 X m¼þN X M k¼N XT

n¼þN XM

ð3:21bÞ ðmÞ

Gk;sþ1 dk;sþ1 einax ; ð3:15aÞ

vl ðx; s þ 1Þ ¼ 

Qðx; s þ 1Þ ¼

ðnmÞ ðmÞ Gk;sþ1 dk;sþ1 :

¼ NT ima

n¼þN XA

ðnÞ

ðinaCÞ

n¼NA ;n–0

ð3:17bÞ

1

dAL dt

! einax þ C L ðs þ 1Þ; sþ1

ð3:23Þ

m¼N M k¼0

The reader may note that V sþ1 is not independent but results from ðnÞ specification of U sþ1 (see Appendix C). Eq. (3.17) can be used to express boundary conditions along the lower and upper walls, i.e., for ^U ðx; s þ 1Þ. In the case of our model problem the ^L ðx; s þ 1Þ and y at y boundary conditions along the lower wall take the following form

where CL is an arbitrary function of time. An analogous expression written for the upper wall introduces another arbitrary function of time, i.e., CU. Substitution of (3.23) and an equivalent relation for the upper wall into Eq. (3.19) and extraction of mode zero results in

m¼þN X M k¼N XT

C U ðs þ 1Þ ¼ Q ð0Þ ðs þ 1Þ þ C L ðs þ 1Þ:

ð0Þ

ðmÞ

ðnmÞ

ðnÞ

Gk;sþ1 ðdL Þk;sþ1 ¼ C1 F L ðs þ 1Þ;

m¼N M k¼0 m¼þN X M k¼N XT m¼N M k¼0

ðmÞ

ðnmÞ

imaGk;sþ1 ðwL Þk;sþ1 ¼ C1

ðnÞ dAL

dt

jnj P 0

ð3:18aÞ

! ;

jnj P 1;

sþ1

ð3:18bÞ

(0)

ð3:24Þ

The value of Q (s + 1) is assumed in this analysis to be known and independent of time, and equal to the flow rate of the reference flow, i.e., Q ð0Þ ¼ 43. One of the functions CU and CL can be selected arbitrarily and the other one follows from (3.24). In the description given below

249

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

the latter one has been selected arbitrarily by introducing condition WT = 0 at a conveniently selected point x0 at the lower wall resulting in n¼þN XA

C L ðs þ 1Þ ¼ C1

ðnÞ

dAL dt

ðinaÞ1

n¼N A ;n–0

! einax0 :

ð3:25Þ

sþ1

The computationally useful form of (3.24) and (3.25) requires introduction of W as this is the unknown in our problem formulation (see Eq. (2.2); W = WT  W0). Values of W0 evaluated along the lower wall represent a known function of t and x, periodic in x, that can be expressed as

W0 ðy^L ðx; s þ 1ÞÞ ¼

n¼þN XM

ðnÞ

ðHL Þsþ1 einax ;

ð3:26Þ

n¼N M

with a similar expression for the upper wall. Values of W along the lower wall can be computed from (3.19), (3.23)–(3.26) as

Wðx; y^L ðx; s þ 1Þ; s þ 1Þ ¼

n¼þN XM

ðnÞ

ðinaCÞ

dAL dt

1

n¼N M ;n–0

þ

n¼þN XM

einax 

ðnÞ

ðinaCÞ

dAL dt

1

n¼N M ;n–0

! !

n¼þN XM

ðnÞ

ðHL Þsþ1 einax

n¼N M

sþ1

einax0 ;

ð3:27aÞ

sþ1

and for the upper wall as

Wðx; y^U ðx; s þ 1Þ; s þ 1Þ ¼

n¼þN XM

ðnÞ

ðinaCÞ

dAU dt

1

n¼N M ;n–0 n¼þN XM

þ

ðinaCÞ1

ðnÞ dAL

n¼N M ;n–0

! inax

e !



n¼þN XM

4. Testing of the algorithm ðnÞ inax U Þsþ1 e

ðH

n¼N M

sþ1

einax0 þ Q ð0Þ :

dt

ð3:27bÞ

sþ1

Stream function W at the lower wall can be expressed in terms of the Fourier expansion (3.5) and the Chebyshev expansion (3.8) in the form

Wðx; y^L ðx; s þ 1Þ; s þ 1Þ ¼

n¼þN XM

Eqs. (3.10), (3.18) and (3.30) form a complete set of algebraic ðnÞ equations for the unknown coefficients Gk;sþ1 , k = 0, . . . , NT, n = 0, . . . , NM. A solution of this system of equations moves calculations forward by one time step. Various methods of solutions will be discussed in the next section. Once the stream function has been determined, the velocity components can be computed from the definition of W (see Eq. (2.2)) while evaluation of pressure is discussed in Appendix D. The reader may be concerned whether the complete discretized system would always be solvable regardless of the locations of the moving boundaries. The potential solvability problem does not occur due to the properties of Chebyshev polynomials. The projection equations for each Fourier mode involve coefficients of all Chebyshev polynomials and form a set of coupled algebraic equations of type (3.10). The projection equations however do not provide any coupling among the Fourier modes and this coupling is provided only through the discretized boundary conditions of type (3.18) and (3.30). The complete system may thus become undetermined only if the discretized boundary conditions fail to provide the necessary couplings. This might occur if the locations of the moving boundaries coincide with zeros of the first four Chebyshev polynomials and their derivatives at both walls. This is not possible because (i) two consecutive Chebyshev polynomials do not share the same roots and (ii) a Chebyshev polynomial and its derivative cannot have overlapping roots. The solvability of the system cannot however, be guaranteed if other basis functions are used. An alternative formulation frequently found in the literature [19] uses a constant pressure gradient constraint rather than the fixed mass flux constraint and the corresponding boundary conditions are discussed in Appendix E.

inax ^ UðnÞ ; L ðyL ðx; s þ 1ÞÞe

ð3:28Þ

We shall discuss the performance of the algorithm in the context of two convenient test problems involving movements of boundaries, i.e., movements of the lower wall corresponding (i) to a traveling elastic wave and (ii) to a standing elastic wave. 4.1. Traveling elastic wave 4.1.1. Test problem Consider an elastic wave traveling along the lower wall with the upper wall being flat. The shape of the resulting channel can be described as

n¼N M

yL ðx; tÞ ¼ 1 þ

where

n¼þN XM

ðnÞ

HL einaðxctÞ ;

yU ðx; tÞ ¼ 1

ð4:1aÞ

n¼Nm ;n–0

^ UðnÞ L ðyðx; s þ 1ÞÞ ¼

m¼þN X M k¼N XT

ðmÞ

ðnmÞ

Gk;sþ1 ðwk;sþ1 ÞL :

ð3:29Þ

m¼NM k¼0

Substitution of (3.28) and (3.29) and similar expressions for the upper wall into (3.27) results in the form of the closing conditions useful for numerical implementation, i.e. m¼þN X M k¼N XT

ðmÞ

ðmÞ

Gk;sþ1 ðwk;sþ1 ÞL

m¼N M k¼0 ð0Þ

¼ HL þ

n¼þN XM

ðnÞ

ðinaCÞ1

n¼NM ;n–0 m¼þN X M k¼N XT m¼NM k¼0

dAL dt

  ðmÞ ðmÞ ð0Þ Gk;sþ1 wk;sþ1 ¼ HU þ U

and, in the simplest case of a sinusoidal wave, the shape of the lower wall becomes

yL ðx; tÞ ¼ 1 þ S cos½aðx  ctÞ ¼ 1 þ ð0:5SeiaðxctÞ þ CCÞ;

where c denotes the phase speed, a denotes the wave number and S stands for the amplitude of the wave. Variations of the location of the lower wall as a function of time are illustrated for this particular case in Fig. 2. Use of the Galileo transformation

X ¼ x  ct

! einax0 ;

ð3:30aÞ

sþ1 n¼þN XM n¼N M ;n–0

 einax0 þ Q ð0Þ :

ð0Þ

ðinaCÞ1

dAL dt

!

ð3:30bÞ

ð4:2Þ

transforms the time-dependent moving boundary problem into a time-independent fixed boundary problem. The full problem in the moving frame of reference (X, y) takes the form

"

sþ1

ð4:1bÞ

# " # o3 o4 o2 o4 c þ Wþ þ2 2 þ W ¼ 0; oX 3 oXoy2 oX 4 oX oy2 oy4 o3

with boundary conditions in the form

ð4:3Þ

250

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

t=0, T

T/8

T/4

3T/8 T/2

5T/8

3T/4

7T/8

n 11 12 13 14 15

yL

(n)

Real{DΦ (y)}

top of the corrugation

bottom of the corrugation

4.0E-06

-0.8

2S

-1

2.0E-06

0.0E+00

-1.2

S 0

2

4

-1.05

6

S -1

-0.95

x

vðyL ðXÞÞ ¼ c

vðX; þ1Þ ¼ 0; n¼þN XM

uðX; yL ðXÞÞ ¼ 0;

ðnÞ

inaHL einax :

ð4:4Þ

n¼N M ;n–0

Elastic wave propagation problem represents a good test problem as it can be solved as a steady fixed boundary problem in the moving reference frame and an unsteady moving boundary problem in the fixed reference frame, and both solutions can be compared. The former one, i.e., (4.3) and (4.4), provides testing opportunity for the analysis of accuracy of spatial discretization and, especially, accuracy of enforcement of boundary conditions, and the latter one provides opportunity for testing of accuracy of temporal discretization. The former problem is solved using the steady version of the method discussed in the previous section. 4.1.2. Solution in the moving frame of reference Problem (4.3) and (4.4) represents steady fixed boundary problem. As the first step, we wish to demonstrate the spectral accuracy of the spatial discretization used in the present work. While this discretization is fairly standard, one needs to pay special attention to the zone around the moving wall where the concept of immersed boundary conditions is used to impose flow boundary conditions due to possible problems associated with the formation of boundary layers in the distribution of modal functions. The Chebyshev expansions (3.8) with coefficients calculated using Galerkin procedure (3.10) are guaranteed to be spectrally-accurate with the increasing number of terms NT. In most cases sixty Chebyshev polynomials provide machine accuracy. When a ? 1 (corrugation with shorter wavelength) and higher Fourier modes begin to play a role, one needs to increase the number of Chebyshev polynomials in order to resolve wall boundary layers. These layers become extremely thin for larger values of a and for higher Fourier modes (see Fig. 3). Modal functions change very rapidly inside these layers while they are nearly zero in the rest of the domain. Typically one needs to use NT  80 for a = 20 and NT  160 for a = 50 in order to provide the required resolution. The second aspect of spectral accuracy involves convergence of the truncated Fourier series (3.5) describing x-variations of the unknown. In all tests dealing with this issue the number of Chebyshev polynomials NT was kept sufficiently large in order to reduce the associated error to machine accuracy. Chebyshev norm defined as

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z 1

ðnÞ

DU ¼ C2 ^ ðy ^; tÞDUðnÞ ðy ^; tÞx ^Þdy ^; DUðnÞ ðy ^ x 1

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi

^ ¼ 1= 1  y ^2 x

ð4:5Þ

-0.85

-0.8

had been adopted as a measure of the ‘‘magnitude” of the derivative of the modal function U(n). Results displayed in Fig. 4 demonstrate that this norm decreases as a function of the mode number n with the rate of decrease very rapidly reaching the (asymptotic) exponential form. Error in the enforcement of flow boundary conditions is of special interest. For convenience, we adopt the L1 norm for u- and vvelocity components evaluated at the lower wall as quantitative measures of this error. These norms are defined as

kuer ðXÞk1 ¼ kver ðXÞk1 ¼

sup

juer ðXÞj;

06X62p=a

sup 06X62p=a

jver ðXÞj;

ð4:6aÞ

where

uer ðXÞ ¼ uðX; yL ðXÞÞ;

ver ðXÞ ¼ vðX; yL ðXÞÞ  vL ðXÞ;

ð4:6bÞ

and their distributions as functions of S and a are illustrated in Figs. 5 and 6. The reader may note in judging these results that S = 2 corresponds to a situation when the upper most point of the lower wall comes into contact with the upper wall. The available results suggest that the error is at machine accuracy level if a and S are below certain critical values. Once these values are reached, the error begins to increase rapidly in a fairly universal manner. The critical values of a and S can be increased by increasing the number of Fourier modes NM used in the calculation, but the qualitative character of the error increment remains unchanged. The location of the maximum error overlaps with the location of the lowest point in the lower wall, as illustrated by distribution of error along the lower boundary displayed in Fig. 7.

(n)

uðX; þ1Þ ¼ 0;

-0.9

Fig. 3. Distribution of the real part of DU(n) as a function of y for the higher modes (n > 10) in the vicinity of the lower wall for the model problem (4.3) and (4.4) for a = 5 and S = 0.05. NM = 15 Fourier modes and NT = 80 Chebyshev polynomials used in the computations.

||DΦ ||ω

Fig. 2. Shape of the lower wall deformed by elastic traveling wave described by Eq. (4.1b) with the amplitude S = 0.2, the wave number a = 1.0 and the phase speed c = p at t = 0, T/8, T/4, 3T/8, T/2, 5T/8, 3T/4 and 7T/8, where T denotes one time period.

y

10

-2

10

-4

10

-6

10

-8

0.25 0.2 S=0.05

10

-10

10

-12

0.1 0

2

4

6

8

10

0.15 12

14

16

18

20

Fourier mode number, n Fig. 4. Variations of the Chebyshev norm of DU(n) (see Eq. (4.5)) as a function of the Fourier mode number n for the model problem (4.3) and (4.4) for different wave amplitudes S for a = 1.0. NM = 20 Fourier modes and NT = 60 Chebyshev polynomials used in the computations.

251

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

A

B 30

||uer||∞

10

10-8

10

||ver||∞

10

-8

5

-10

10

10

30

-10

α=1

5 α=1

-12

10

10 -3

10 -2

10 -1

10

10 0

-12

10-3

10-2

10-1

S

S

100

Fig. 5. Variations of the kuer(X)k1 (A) and kver(X)k1 (B) norms (see Eq. (4.6a)) as a function of the wave amplitude S for selected values of the wave number a for the model problem (4.3) and (4.4) with the phase speed c = p. The dashed and solid lines correspond to results obtained with the NM = 10, 15 Fourier modes, respectively. The reader may note that S = 2 corresponds to peak of the wave reaching the top wall. NT = 60 Chebyshev polynomials were used in the calculations.

A

B 10

-8

10

-8

10

0.15 0.1

S=0.2

-10

||ver||∞

||uer||∞

0.15

0.05

S=0.2 10

0.1

-10

0.05

0.025 10

0.025

-12

10

-1

10

0

10

α

1

10

-12

10

-1

10

0

10

α

1

Fig. 6. Variations of the kuer(X)k1 (A) and kver(X)k1 (B) norms (see Eq. (4.6a)) as a function of the wave number a for selected values of the wave amplitude S for the model problem (4.3) and (4.4) for the phase speed c = p. Other parameters as in Fig. 5.

when several waves pass through the computational box. Variations of the velocity components in a few test points displayed in Fig. 9 demonstrate the expected periodic variations in time and the expected phase differences associated with different locations in the x-direction. Distribution of shear stress at the moving boundary displayed in Fig. 10 demonstrates that the expected periodicity and phase shifts have been reproduced. Results of the temporal grid convergence studies displayed in Fig. 11 demonstrate that the one-step and two-steps algorithms reproduce the expected first- and second-order accuracy in time, respectively.

5E-09

ver

uer

0

-5E-09

0

2

x

4

6

Fig. 7. Distribution of error uer and ver (see Eq. (4.6b)) in the enforcement of flow boundary conditions at the lower wall for the model problem (4.3) and (4.4) with the phase speed c = p, the amplitude S = 0.2 and the wave number a = 1. The presented results were obtained with NM = 15 Fourier modes and NT = 60 Chebyshev polynomials.

4.2. Standing elastic wave 4.2.1. Test problem Consider standing elastic wave at the lower wall with the upper wall being flat. The shape of the resulting channel can be described as

yL ðx; tÞ ¼ 1 þ 4.1.3. Solution in the fixed frame of reference The problem discussed in the previous section becomes a moving boundary problem when expressed in a fixed reference frame. We use solution obtained in the moving frame of reference as an initial condition and track evolution of the flow as a function of time. The time history of the error of enforcement of boundary conditions is illustrated in Fig. 8. It can be seen that the maximum error is similar to that found in the case of solution obtained in the moving frame of reference, location of this error follows location of the maximum opening of the slot as it moves in the x-direction, and the magnitude of this error remains approximately constant

n¼þN XM

ðnÞ

HL ðtÞeinax ;

yU ðx; tÞ ¼ 1

ð4:7Þ

n¼NM ;n–0

and, in the simplest case of a sinusoidal wave, the shape of the lower wall becomes

yL ðx; tÞ ¼ 1 þ S cosðxtÞ cosðaxÞ ¼ 1 þ ð0:5S cosðxtÞeiax þ CCÞ;

ð4:8Þ

where a denotes the wave number, S stands for the amplitude and x denotes the frequency of the wave. Variations of the location of the lower wall as a function of time are illustrated for this particular case in Fig. 12.

A

6.5E-09

0.0E+00

2.5T

2.25T

t=2T, 3T

B

6.0E-09

ver

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

uer

252

0.0E+00

t=2T, 3T

2.75T -6.5E-09

0

2.25T

2.5T

2.75T

2

4

x

-6.0E-09

6

0

2

4

x

6

Fig. 8. Distribution of error in the enforcement of flow boundary condition uer (A) and ver (B) for the u- and v-velocity components (see Eq. (4.6b)), respectively, at the lower wall at t = 2T, 2.25T, 2.5T, 2.75T and 3T, where T stand for one time period, for the model problem (4.3) and (4.4) with the amplitude S = 0.2, the wave number a = 1.0 and the phase speed c = p. The presented results were obtained through a direct solution of the moving boundary problem in the fixed coordinates system with NM = 15 Fourier modes and NT = 60 Chebyshev polynomials. Solution corresponding to the fixed boundary problem in the moving frame of reference (4.3) and (4.4) was used as the initial condition.

A

λ/4

x=0

λ/2

B

3λ/4

0.7

λ/4

x=0

λ/2

3λ/4

1.05 0.35

v

u

0.7

0.35

0

0

-0.35

0

2

4

t

6

-0.7

0

2

t

4

6

Fig. 9. Variations of the u-velocity (A) and v-velocity (B) components for three time periods at points (x, y) = (0, 0.6), (k/4, 0.6), (k/2, 0.6), (3k/4, 0.6), k = 2p/a, for the model problem (4.3) and (4.4) solved directly as a moving boundary problem in the fixed reference frame. Solution corresponding to the fixed boundary problem in the moving frame of reference (4.3) and (4.4) was used as the initial condition. NM = 15 Fourier modes and NT = 60 Chebyshev polynomials were used in the computations. Other conditions as in Fig. 2.

5T/8

3T/4

For convenience, we use the L1 norm as a quantitative measure of error associated with the enforcement of velocity boundary conditions at the lower wall, i.e.

7T/8

1.8

τL

kuer ðx; tÞk1 ¼ kver ðx; tÞk1 ¼

1.5

sup juer ðx; tÞj; 06x62p=a

sup jver ðx; tÞj;

ð4:11aÞ

06x62p=a

where t-2T=0, T 1.2

0

T/8 2

T/4

uer ðx; tÞ ¼ uðx; yL ðx; tÞ; tÞ;

3T/8 T/2

x

4

Fig. 10. Distribution of the wall shear stress at different time levels after two time periods for the model problem (4.3) and (4.4). Other parameters are as in Figs. 2 and 9.

The full test problem has the form

"

# " # 2 o o2 o2 o4 o2 o4 þ 2 W¼ þ þ W; ot ox oy ox4 ox2 oy2 oy4

ð4:9Þ

with boundary conditions in the form

uðx; þ1; tÞ ¼ 0; uðx; yL ðxÞ; tÞ ¼ 0;

vðx; þ1; tÞ ¼ 0; vðx; yL ðx; tÞ; tÞ ¼

n¼þN XM n¼NM ;n–0

ver ðx; tÞ

¼ vðx; yL ðx; tÞ; tÞ  vL ðx; tÞÞ:

6

d  ðnÞ  inax HL ðtÞ e : dt ð4:10Þ

ð4:11bÞ

4.2.2. Solution of the test problem Fig. 13 illustrates distribution of error in the enforcement of boundary conditions at the lower wall. It can be seen that the magnitude of the error changes periodically in time with frequency equal to double the frequency of the standing wave. Variations in the magnitude of the error can be correlated with variations of the shape of the wall, i.e., the maximum of error occurs at times when the channel opening is the largest. The magnitude of error can be reduced by increasing the number of Fourier modes used in the calculations but the qualitative behaviour of error with time remains unchanged. The spatial distribution of error is illustrated in Fig. 14 after 3 and 3.5 cycles of the wave motion. It can be seen that the maximum of error shifts to location corresponding to the maximum opening of the channel. The Fourier spectra of the error shown in

253

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

A

0.003

0.002

Δt

B

2

0.001

0

0.003

0.002

Δt

2

0.001

0

0.0006

8.0E-04

1.2E-04

9E-05

Two-step method 0.0004

4.0E-04

0.0002

4.0E-05 2.0E-04

6E-05

verror

Two-step method

3E-05

One-step method

One-step method 0

0

0.015

verror

8.0E-05

uerror

uerror

6.0E-04

0.03

0.045

0

0.015

0.03

Δt

Δt

0

0.045

Fig. 11. Variations of error in the u- (Fig. 10A) and v- (Fig. 10B) velocity components as a function of the step size Dt used in the temporal discretization in the case of the model problem (4.3) and (4.4) with the amplitude S = 0.1, the wave number a = 1.0 and the phase speed c = p. The error is defined as the maximum of the absolute value of the difference between the results obtained through the direct solution of the moving boundary problem in the fixed frame of reference and the corresponding fixed boundary problem in the moving frame of reference at a time corresponding to t = 1.0. NM = 10 Fourier modes and NT = 60 Chebyshev polynomials were used in the computations. Solution of the model problem (4.3) and (4.4) in the moving frame of reference is taken as the initial condition for the direct solution method.

t=0, T

-0.8

3T/8, 5T/8 S

yL

T/4, 3T/4 -1

T/8, 7T/8 -1.2

T/2 0

2

x

4

6

Fig. 12. Shape of the lower wall modified by elastic standing wave with the wave number a = 1.0, the amplitude S = 0.2 and the frequency x = p at times t = 0, T/8, T/4, 3T/8, T/2, 5T/8, 3T/4, 7T/8 and T, where T denotes one time period.

Fig. 15 demonstrate that the NM number of Fourier modes have been eliminated according to the construction of the immersed boundary conditions algorithm described in Section 3. The largest error is associated with the first Fourier mode omitted in the enforcement of boundary conditions. The error associated with the higher omitted Fourier modes rapidly decays as the mode number increases. Fig. 15 also displays results of tests carried out in order to check if the method produces any spurious spatial oscillations.

A

NM=9 10

12

-6

12

||ver||∞

||uer||∞

B

NM=9 10-6

Three cases were considered, i.e., in case (A) the wave was represented by the principal Fourier mode and the calculations had been carried out with NM = 9 Fourier modes, in case (B) the wave was represented by the second Fourier mode (the principal mode has the wave number a = 0.5), and in case (C) the wave was represented by the third Fourier mode (the principal mode has the wave number a = 1/3). In order to have fully equivalent representations, the number of Fourier modes used in cases (B) and (C) were NM = 18 and NM = 27, respectively. The selected representations admitted subharmonics of the 1/2 type in case (B) and 1/3 type in case (C). The Fourier spectra shown in Fig. 16 demonstrate the equivalency of the results in all three cases. No subharmonics had been produced during the solution process and the modes expected to produce zero contributions in cases (B) and (C) behaved as expected. Results of temporal grid convergence studies are displayed in Fig. 16 and demonstrate that the one-step and two-steps algorithms reproduce the expected first- and second-order accuracy in time, respectively. Variations of velocity components at a test point displayed in Fig. 17 demonstrate the expected periodic variations in time. Variations of the shear stress at the moving wall are displayed in Fig. 18 and also follow the expected periodic variations in time. The instantaneous streamlines associated with a standing wave with a more complex form are displayed in Fig. 19 in order to illustrate the algorithm’s ability to deal with more complex motions.

10-9

15 10 -9

15 10-12

0

2

t

4

6

10

-12

0

2

t

4

6

Fig. 13. Variations of the kuer (x, t)k1 (A) and kver(x, t)k1 (B) norms (see Eq. (4.11a)) as a function of time over three time periods for the moving boundary problem defined by (4.8)–(4.10) for a = 1, S = 0.2 and x = p (shape of the wall is illustrated in Fig. 12). NT = 60 Chebyshev polynomials were used and the solid, dotted and dashed lines corresponds to results obtained with NM = 9, 12, 15 Fourier modes, respectively. The initial conditions used correspond to shape of the wall given by Eq. (4.8) at time t = 0 and zero flow modifications. It can be seen that the initial transient dies out and the time-periodic flow response begins to dominate after 1–2 periods.

A 1.5E-09

B 1.5E-09

ver

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

uer

254

0.0E+00

-1.5E-09

0.0E+00

0

2

4

6

-1.5E-09

0

2

4

x

6

x

Fig. 14. Distribution of the error uer (A) and ver (B) in the u- and v-velocity components (see Eq. (4.11b)) at the lower wall after three (solid line) and three and half time periods (dashed line) for the model problem (4.8)–(4.10) evaluated with NM = 15 Fourier modes. All other conditions are as in Fig. 13.

A

B at the initial state

after 3 time periods

(n)

after 3.5 time periods

3.0E-07

0

10

20

after 3 time periods

after 3.5 time periods

3.0E-07

C

B

Case A 1.0E-11

at the initial state

|V |er

|U(n)|er

6.0E-07

6.0E-07

1.0E-11

40

Fourier mode number, n

0

C

B

Case A 30

10

20

30

40

Fourier mode number, n

Fig. 15. Fourier spectra of error distributions uer and ver (see Eq. (4.11b)) in the u- (A) and v- (B) velocity components at the lower wall for the model problem (4.8)–(4.10). Three different forms of Fourier expansions were considered, i.e., case A: a = 1.0, NM = 9; case B: a = 0.5, NM = 18 and case C: a = 1/3, NM = 27. NT = 60 Chebyshev polynomials used in all cases. All other conditions as in Fig. 13.

A

0.04

Δt

0.02

0

B

0.04

Δt

0.02

0

0.4446

-0.00637

0.4446

-0.00637

Two-step method

v

u

0.444

Two-step method

0.44457

-0.0063

-0.00638

v

u

-0.00623

-0.00639 0.4434

One-step method

One-step method

0.44454

-0.0064 0

0.001

Δt

2

0.002

0.003

0

0.001

0.002

Δt

-0.00616 0.003

2

Fig. 16. Variations of the u- (A) and v- (B) velocity components as a function of the step size Dt used in the temporal discretization at a test point (x, y)=(p/2, 0.75) in the case of the model problem (4.8)–(4.10) with the amplitude S = 0.1, the wave number a = 1.0 and the frequency x = p at time t = 1.0. The initial conditions correspond to the solution of the fixed boundary problem for the wall shape given by Eq. (4.8) at t = 0. NM = 10 Fourier modes and NT = 60 Chebyshev polynomials were used in the computations.

4.3. Efficiency of the algorithm The proposed algorithm requires solution of a system of linear algebraic equations resulting from the discretization of the field equations and boundary conditions at each time step. The resulting matrix structure is displayed in Fig. 20 and underscores the fact that algebraic equations resulting from discretization of each modal equation of type (3.7) are uncoupled, i.e., the corresponding coefficients form blocks in the upper triangular form. The only coupling between Fourier modes is provided through the boundary

conditions, i.e., the corresponding entries form horizontal lines in the matrix (four lines per block). Direct solution of the matrix equation leads to a solution method that we shall refer to as the ‘direct algorithm’. The efficiency of this algorithm is very good as the only matrix entries that have to be recomputed at each time step correspond to the boundary conditions. In addition, the shape of the moving wall can be easily adjusted by changing the magnitude of the relevant Fourier coefficients. In this sense, the proposed algorithm is superior with respects to algorithms based on dynamic grid adjustments and/or mapping methods as the cost of

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

0.9

0.08

u

v

0.6

0

0.3

-0.08

0

2

4

6

t Fig. 17. The evolution of the u- (solid line) and v- (dashed line) velocity components at a test point (x, y) = (k/4, 0.6) during the first four time periods for the model problem (4.8)–(4.10) evaluated with NM = 15 Fourier modes. All other conditions are as in Fig. 13.

7T/8

3

3T/8

3T/4

T/4

τL

2.5

2 1.5 1

t-2T=0,T

T/8 0

2

5T/8 4

T/2

255

with many small matrices rather then one very large matrix as well as they open possibility for parallelization of the computations. This issue becomes significant in the case of three-dimensional problems and large number of Fourier modes and Chebyshev polynomials. The rates of convergence of all versions of the decoupled algorithm are generally very good; they decrease with an increase of the amplitude S and the wave number a characterizing the moving boundary. Their efficiency is illustrated in Fig. 21 in the context of the test problem discussed in Section 4.1 solved in the moving reference frame. It can be seen that the iterative method (version 1) requires less time than the direct method as long as the amplitude S is less than 0.23. Over-relaxation decreases the efficiency of the method for S < 0.23 but improves efficiency for higher amplitudes S. The relative efficiency of different versions is illustrated in Figs. 22 and 23 in the context of the test problem discussed in Section 4.2. Fig. 21 displays variations in the number of iterations required for convergence as a function of simulation time. The maxima in the number of iterations correspond to instants of time when the channel opening passes through the maximum and minima correspond to times when the channel becomes flat. Fig. 23 displays analogous information but expressed in terms of the computing time required for convergence. The local maxima correspond to instants of time when the channel opening passes through the maximum and minima correspond to times when the channel becomes flat. It can be seen that version 1 of the decoupled algorithm is most efficient as inclusion of additional coupling between Fourier modes increase number of iterations and time required for convergence.

6

x Fig. 18. Distribution of the wall shear stress at different time levels during the third time period. The results shown were determined using NM = 15 Fourier modes. All other conditions are as in Fig. 13.

re-gridding and matrix construction have been eliminated, and comparable to algorithms based on various forms of interface tracking while delivering much higher accuracy. The matrix of coefficients can be very large when a large number of Fourier modes are required and this motivates the search for a method of solution that avoids construction as well as inversion of the complete matrix. Structure of the matrix (see Fig. 20) suggests the use of an iterative solution algorithm based on the decoupling of Fourier modes. We have tested four versions of such algorithms which we shall refer to as the ‘decoupled algorithms’. In the first version, the unknowns corresponding to a Fourier mode of interest in equation (3.18) and (3.30) at the current time step are expressed in terms of the remaining Fourier modes using their values from the previous time step (or from the previous iteration). The solution process begins with mode 0, proceeds to the next mode using the most recent information available and continues until the last mode NM is reached, and then it is repeated until a convergence criterion is satisfied. In this way, the inversion of the complete matrix of size (NT + 1)  (2NM + 1) is replaced by a repetitive solution of system of (NT+1) equations for each Fourier mode. Versions 2–4 retain different couplings between the main modes as this might accelerate convergence. In version 2, we solve for each mode simultaneously with mode 0, e.g., we solve a system of equations for modes 0 and 1, then for modes 0 and 2, and then for modes 0 and 3, and so on. In version 3, we solve for each mode simultaneously with mode 1, i.e., mode 0 and 1, then 1 and 2, followed by 1 and 3, and so on. In the forth version, we solve for modes 0, 1 and 2 together and then separately for each of the remaining modes. It is necessary to point out that the decoupled algorithms reduce memory requirements as one needs to work

5. Conclusions A gridless algorithm for unsteady flow problems described by the biharmonic operator (Stokes flow) has been presented and tested. The algorithm uses a fixed computational domain with the flow domain completely immersed inside the computational domain. The flow boundary conditions are imposed using the concept of immersed boundary conditions. Two versions of boundary condition had been presented, i.e., corresponding to the fixed mass flux constraint and to the fixed pressure gradient constraint. The algorithm is suitable for the flow domains that have form of a channel whose geometry is periodic in the streamwise direction. The opening of the channel may change as a function of time in an arbitrary manner but resulting in an instantaneous geometry that is Fourier transformable. The algorithm uses Fourier expansions for the spatial discretization in the flow direction, Chebyshev expansions for the discretization in the transverse direction and the first- and second-order implicit temporal discretizations. Various tests confirm the spectral accuracy of the spatial discretization and the theoretically-predicted accuracy of the temporal discretization. The algorithm is very efficient as the part of the coefficient matrix corresponding to the field equations needs to be computed only once and the only changes required at each time step are limited to entries corresponding to the boundary conditions. In addition, the algorithm is very flexible as far as adaptation to different geometries is concerned as the only changes required are limited to changes in the Fourier coefficients describing motions of the walls. The numerical error associated with treatment of flow boundary conditions is well controlled during the simulation process. No numerical instabilities have been observed. Direct and decoupled versions of the algorithm had been explored. It had been found that while the direct algorithm is more efficient in situations when large-amplitude motions of the boundaries are of interest, a simple mode-decoupled algorithm is significantly faster when the interest is limited to small amplitude motions.

256

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

Fig. 19. Instantaneous streamlines in a channel bounded by yL (x, t) = 1+(0.05ieix  0.0167ie3ix + 0. 0125ie4ix + CC) cos (xt) and yU (x, t) = 1, for x = p using NM = 25 and NT = 80 evaluated after two time periods. Results are shown at times t  2T = 0, T/8, T/4, 3T/8, T/2, 5T/8, 3T/4, 7T/8 and T, where T stands for one time period. The initial conditions used correspond to shape of the walls at time t = 0 i.e. yL (x, 0) = 1 + (0.05ieix  0.0167ie3ix + 0. 0125ie4ix + CC) and yU (x, 0) = 1, and zero flow modifications.

Time required (in sec)

102

Direct Method 10

1

1.90 1.75 10

0

1.50 1.25

10-1

Relaxation factor = 1 .0

0.1

0.2

0.3

S

Fig. 20. Structure of the coefficient matrix for the test problem (4.8)–(4.10) obtained with NM = 10 and NT = 60. Non-zero entries are marked in black.

Fig. 21. Performance of the first version of the decoupled algorithm with different values of the over-relaxation factor for the model problem (4.3) and (4.4) as a function of the amplitude S of the traveling wave with the wave number a = 1.0 and the phase speed c = p. Time required by the direct method is given for reference. All calculations have been carried out with NM = 15 Fourier modes and NT = 60 Chebyshev polynomials.

257

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

Coefficients in (3.13b) can be evaluated with the help of the recur^Þ ¼ 2T k ðy ^Þ þ 2y ^DT k ðy ^Þ  DT k1 ðy ^Þ that leads to rence relation DT kþ1 ðy a relation

600

No.of iterations

500 400

1 X

ðmÞ

dkþ1;sþ1 ¼ 2

300

ðnÞ

ðmnÞ

ðB:3Þ

Method 3 Method 4

Method 2

200

whose evaluation begins at k = 0 giving the initial terms in the form ðmÞ

d0;sþ1 ¼ 0 for jmj P 0;

100

ðmÞ

P 1; 1

2

3

4

ð0Þ

d1;sþ1 ¼ 1;

Fig. 22. Number of iterations required by different versions of the decoupled algorithm as a function of the simulation time (and thus position of the wall). Simulations have been carried out over two periods of motion of the standing wave described in Section 4.2 with the wave amplitude S = 0.15, the wave number a = 1.0 and the frequency x = p using zero initial condition, NM = 15 Fourier modes and NT = 60 Chebyshev polynomials.

ðmÞ

d2;sþ1 ¼ 4Psþ1 for jmj P 0:

1 X

ðmÞ

rkþ1;sþ1 ¼ 2

ðmnÞ

ðnÞ

ðmÞ

ðmÞ

Psþ1 dk;sþ1  r k1;sþ1 þ 4dk;sþ1 for k > 3;

ðB:5Þ

n¼1

whose evaluation begins at k = 0 giving the initial terms in the form ðmÞ

ðmÞ

ðmÞ

Method 3

ðB:4Þ

Coefficients required in Appendix D can be evaluated with the ^Þ ¼ 4DT k ðy ^Þ þ 2y ^ D 2 T k ðy ^Þ help of the recurrence relation D2 T kþ1 ðy ^Þ that leads to a relation D2 T k1 ðy

r0;sþ1 ¼ r 1;sþ1 0 for jmj P 0; Method 2

ðmÞ

d1;sþ1 ¼ 0 for jmj

ðmÞ r k;sþ1

t

Required time (in sec)

ðmÞ

n¼1

Iterative Method 1

0

ðmÞ

Psþ1 dk;sþ1  dk1;sþ1 þ 2wk;sþ1 for k > 2;

ð0Þ

r 2;sþ1 ¼ 4;

ðmÞ

r2;sþ1 ¼ 0 for jmj P 1;

ðmÞ

r3;sþ1 ¼ 24Psþ1 for jmj P 0:

Method 4

ðB:6Þ

101

Appendix C ð0Þ

Demonstration that V sþ1 in Eq. (3.17b) cannot be independently specified. We can write the following relation along line l overlapping ^ ¼ f ðx; s þ 1Þ with the moving boundary y

100

Iterative Method 1 1

2

3

WT ðx þ k; f ðx þ k; s þ 1Þ; s þ 1Þ  WT ðx; f ðx; s þ 1Þ; s þ 1Þ

4

t

¼

Fig. 23. Variations of time (in sec) required to achieve convergence using different versions of the decoupled algorithm as a function of the simulation time (and thus position of the wall). All the other conditions as in Fig. 22.

x

Appendix A

Z

xþk

ðC:1Þ

 oWT ðx; f ðx; s þ 1Þ; s þ 1Þ oWT ðx; f ðx; s þ 1Þ; s þ 1Þ df þC dx ^ ox oy dx ðC:2Þ

and expressed in terms of velocity components along the line l as

Z

xþk

x

Temporal discretization of (3.6) using a one-step implicit method results in

½C4 D4 þ ð2ðnaÞ2 C2  0:5C2 Dt1 ÞD2 þ ððnaÞ4 þ 0:5ðnaÞ2 Dt 1 ÞUsþ1 ðnÞ

2

dWT ðx; f ðx; s þ 1Þ; s þ 1Þ dx ¼ 0; dx

¼0

This work has been carried out with support of SHARCNET and NSERC of Canada. SHARCNET of Canada provided computing resources.

2

xþk

which can be re-arranged as

x

Acknowledgement

1

Z

2

ðnÞ

¼ 0:5Dt ½C D  ðnaÞ Us ;

 df dx ¼ 0: vl ðx; s þ 1Þ þ ul ðx; s þ 1Þ dx

ðC:3Þ

Substitution of (3.11) and (3.12) into (C.3) and integration of the resulting expression gives ð0Þ

V sþ1 ¼ ia

n¼þN XA

ðnÞ

ðnÞ

nU sþ1 Psþ1 ;

ðC:4Þ

n¼NA

n 2 h0; NM i:

ð0Þ

and demonstrates that the mean value of vl (x, s + 1), i.e., V sþ1 , cannot be independently specified as it depends on the specification of ul (x, s + 1) as well as on the wall geometry expressed by f (x, s + 1).

Appendix B ðmÞ

ðmÞ

Evaluation of coefficients wk;sþ1 and dk;sþ1 required in Eq. (3.14) ðmÞ and r k;sþ1 required in Appendix D. Coefficients in (3.13a) can be evaluated with the help of the ^Þ ¼ 2y ^T k ðy ^Þ  T k1 ðy ^Þ that leads to recurrence relation T kþ1 ðy ðmÞ

wkþ1;sþ1 ¼ 2

n¼þ1 X

ðnÞ

ðmnÞ

ðmÞ

Psþ1 wk;sþ1  wk1;sþ1 for k > 1;

ðB:1Þ

n¼1

whose evaluation begins at k = 0 giving the initial terms in the form ð0Þ w0;sþ1

¼ 1;

ðmÞ w0;sþ1

¼ 0 for jmj P 1;

ðmÞ w1;sþ1

¼

ðmÞ Psþ1

for jmj P 0: ðB:2Þ

Appendix D Evaluation of pressure field from the known stream function field. The total pressure field is given by

pðx; y; tÞ ¼ p0 ðxÞ þ p1 ðx; y; tÞ;

ðD:1Þ

where p0(x) = 2x is the pressure field associated with the reference Poiseuille flow and p1 (x, y, t) corresponds to the modification in the pressure field due to boundary motions. The dimensionless field equations have the form

258

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259

ou1 ov1 þC ¼ 0; ^ ox oy ou1 op1 ou21 ou2  C2 21 ¼ 0; þ  ^ ot ox ox2 oy ov1 op1 ov21 ov2  C2 12 ¼ 0: þC  ^ ox2 ^ ot oy oy

ðD:2aÞ

h ðD:2bÞ ðD:2cÞ

The velocity components u1 and v1 are known and can be expressed as

^; s þ 1Þ ¼ C u1 ðx; y

n¼þN XM

ih

i

C3 D2  1:5CDt 1 Uð0Þ ðy^U ðx; s þ 1Þ; s þ 1Þ  Uð0Þ ðy^L ðx; s þ 1Þ; s þ 1Þ ^U ðx; s þ 1Þ  y ^L ðx; s þ 1Þ ¼ Ap ðs þ 1Þ½y h i ð0Þ 1 ^U ðx; sÞ; sÞ  Uð0Þ ðy ^L ðx; sÞ; sÞ  2Dt C U ðy h i ^U ðx; s  1Þ; s  1Þ  Uð0Þ ðy ^L ðx; s  1Þ; s  1Þ ; þ 0:5Dt 1 C Uð0Þ ðy

ðE:2Þ ^; s þ 1Þeinax ; D U ðy ðnÞ

n¼N M

^; s þ 1Þ ¼ ia v1 ðx; y

Integration of the above equation across the channel gives

n¼þN XM

ðD:3Þ ^; s þ 1Þeinax : nUðnÞ ðy

n¼NM

The corresponding pressure field can be expressed as

^; s þ 1Þ ¼ Ap ðs þ 1Þx þ p1 ðx; y

n¼þN XM

^; s þ 1Þeinax ; pðnÞ ðy

ðD:4Þ

which is the second condition required for the modal function U(0). The pressure gradient Ap is assumed to be known at each time step. A computationally useful form of Eq. (E.2) is obtained by expressing U(0) and D2U(0) in terms of Chebyshev expansions which are evaluated at the boundaries using the Immersed Boundary Conditions concept described in Section 3. If the boundary overlaps with the ^ ¼ f ðx; s þ 1Þ the reline l overlapping with the moving boundary y quired relations take the forms

n¼NM

where Ap(s + 1) is a function of time and appears in the expression due to the fact that an x-periodic, time-dependent velocity field may be associated with a pressure field that has a time-dependent, linear in x component. The unknowns Ap and p(n) need to be evaluated from the momentum equations. We begin with the evaluation of Ap. Substitution of Eqs. (D.3),(D.4) into (D.2b) and use of the twostep time discretization for the time derivative result in

Uð0Þ ðf ðx; s þ 1Þ; s þ 1Þ ¼

h i ð0Þ Ap ðs þ 1Þ ¼ C3 D3  1:5CDt1 D Usþ1 þ 2Dt 1 CDUð0Þ s

D2 T K ðf ðx; s þ 1ÞÞ ¼

ð0Þ

 0:5Dt 1 CDUs1 ;

k¼N XT

ð0Þ

Gk;sþ1

k¼N XT

ðsÞ

wk;sþ1 eisax ;

ðE:3aÞ

s¼Nh

k¼0

D2 Uð0Þ ðf ðx; s þ 1Þ; s þ 1Þ ¼

s¼þN Xh

ð0Þ

Gk;sþ1

k¼0

s¼þN Xh

ðsÞ

r k;sþ1 eisax ;

ðE:3bÞ

s¼Nh

where s¼þN Xh

ðsÞ

r k;sþ1 eisax :

ðE:3Þ

s¼N h

ðD:5aÞ

h i ^; s þ 1Þ ¼ CðinaÞ1 C2 D2  ðnaÞ2  1:5CDt1 DUðnÞ pðnÞ ðy sþ1

ðsÞ

Method for evaluation of coefficients rk;sþ1 is explained in Appendix B. Substitution of (E.3) written for the upper and lower walls into Eq. (E.2) results in the final form of the condition being sought.

ðnÞ

1 þ 2CDt 1 DUðnÞ s  0:5CDt DUs1 for jnj–0; ðD:5bÞ

^. Similar substitutions applied to Eq. (D.2c) lead to where D ¼ d=dy

h i ^; s þ 1Þ ¼ inaC1 C2 D2  ðnaÞ2  1:5CDt1 UðnÞ DpðnÞ ðy sþ1 ðnÞ

1  2inaDt 1 UðnÞ s þ 0:5inaDt Us1 for jnj P 0:

ðD:6Þ In the case of mode zero the right hand side of (D.6) becomes zero which results in

^; s þ 1Þ ¼ const: pð0Þ ðy

ðD:7Þ

The final expression for the complete pressure field has the form

^; s þ 1Þ ¼ ½Ap ðs þ 1Þ  2x þ pðx; y

n¼þN XM

^; s þ 1Þeinax pðnÞ ðy

n¼N M ;n–0

þ const;

ðD:8Þ

where the constant can be selected conveniently. Appendix E Specification of the closing boundary conditions corresponding to the constant pressure gradient constraint. Eq. (3.30a) represents boundary condition that can be selected arbitrarily. Eq. (3.30b) needs to be replaced with condition expressing the desired pressure gradient. Eq. (D.5a) can be rewritten as

h

i

ð0Þ 1 C3 D3  1:5CDt 1 D Uð0Þ sþ1 ¼ Ap ðs þ 1Þ  2Dt CDUs ð0Þ

þ 0:5Dt1 CDUs1 :

ðE:1Þ

References [1] J.M. Floryan, H. Rasmussen, Numerical analysis of viscous flows with free surfaces, Appl. Mech. Rev. 42 (1989) 323–341. [2] R. Scardovelli, S. Zaleski, Direct numerical simulation of free surface and interfacial flow, Annu. Rev. Fluid Mech. 31 (1999) 567–603. [3] I. Inculet, J.M. Floryan, R. Haywood, Dynamics of water droplets break-up in electric fields, IEEE Trans. Ind. Appl. 28 (1992) 1203–1209. [4] M. Hamed, J.M. Floryan, Numerical simulation of unsteady nonisothermal capillary interfaces, J. Comput. Phys. 145 (1998) 110–140. [5] J.M. Hyman, Numerical methods for tracking of interfaces, Physica 12D (1984) 396–407. [6] J. Glimm, J.W. Grove, X.L. Li, K.M. Shyue, Y. Zeng, Q. Zhang, Three-dimensional front tracking, SIAM J. Sci. Comput. 19 (1998) 703–727. [7] F.H. Harlow, J.E. Welch, Numerical study of large amplitude free surface motions, Phys. Fluids 9 (1966) 842–851. [8] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225. [9] S.J. Osher, J.A. Sethian, Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [10] J.A. Sethian, P. Smereka, Level set methods for fluid interfaces, Annu. Rev. Fluid Mech. 35 (2003) 341–372. [11] J. Szumbarski, J.M. Floryan, A direct spectral method for determination of flows over corrugated boundaries, J. Comput. Phys. 153 (1999) 378–402. [12] C.S. Peskin, The fluid dynamics of heart valves: experimental, theoretical and computational methods, Annu. Rev. Fluid Mech. 14 (1982) 235–259. [13] R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261. [14] C.S. Peskin, The immersed boundary method, Acta Numer. (2002) 479– 517. [15] J.M. Floryan, Centrifugal instability of Couette flow over a wavy wall, Phys. Fluids 14 (2002) 312–322. [16] J.M. Floryan, Three-dimensional instabilities of Laminar flow in a rough channel and the concept of hydraulically smooth wall, Eur. J. Mech. B/Fluids 26 (2007) 305–329. [17] S.Z. Husain, J.M. Floryan, Immersed boundary conditions method for unsteady flow problems described by the Laplace operator, Expert Systems in Fluid Dynamics Research Laboratory Report ESFD-1/2006, Department of

S.Z. Husain, J.M. Floryan / Comput. Methods Appl. Mech. Engrg. 198 (2008) 245–259 Mechanical and Materials Engineering, The University of Western Ontario, London, Ontario, N6A 5B9, Canada, 2006 (To appear in the Int. J. Numer. Meth. Fluids). [18] S.Z. Husain, J.M. Floryan, Implicit spectrally-accurate method for moving boundary problems using immersed boundary conditions concept, Expert

259

Systems in Fluid Dynamics Research Laboratory Report ESFD-2/2006, Department of Mechanical and Materials Engineering, The University of Western Ontario, London, Ontario, N6A 5B9, Canada, 2006. [19] J.M. Floryan, Stability of wall-bounded shear layers in the presence of simulated distributed roughness, J. Fluid Mech. 335 (1997) 29–55.