An introduction to the derivation of surface balance equations without the excruciating pain

International Journal of Heat and Mass Transfer 115 (2017) 992–999

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Technical Note

An introduction to the derivation of surface balance equations without the excruciating pain Alexandre Martin a,⇑, Huaibao Zhang b, Kaveh A. Tagavi a a b

Department of Mechanical Engineering, University of Kentucky, Lexington, KY 40506-0503, USA School of Physics, Sun Yat-Sen University, Guangzhou 510006, China

a r t i c l e

i n f o

Article history: Received 7 April 2017 Received in revised form 10 July 2017 Accepted 17 July 2017 Available online 04 August 2017 Keywords: Boundary conditions Porous flow Surface balance equations Coupled flow

a b s t r a c t Analyzing complex fluid flow problems that involve multiple coupled domains, each with their respective set of governing equations, is not a trivial undertaking. Even more complicated is the elaborate and tedious task of specifying the interface and boundary conditions between various domains. This paper provides an elegant, straightforward and universal method that considers the nature of those shared boundaries and derives the appropriate conditions at the interface, irrespective of the governing equations being solved. As a first example, a well-known interface condition is derived using this method. For a second example, the set of boundary conditions necessary to solve a baseline aerothermodynamics coupled plain/porous flow problem is derived. Finally, the method is applied to two more flow configurations, one consisting of an impermeable adiabatic wall and the other an ablating surface. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Arguably, one of the more complex and painful [1] tasks in developing a numerical simulation code, such as those in the field of Computational Fluid Dynamics (CFD), is constructing the boundary condition equations. The task becomes even more complicated when multiple codes that solve different sets of governing equations need to be coupled via physical boundary conditions at their interfaces. The traditional approach for coupling different codes is to develop the so-called ‘‘balance equations” at the interface, which are then used to solve for various parameters at that particular interface. These balance equations are often determined either by one or a combination of the following methodologies: (a) determination of empirical relations followed by phenomenological derivation of transport equations; (b) through consideration of certain physical phenomenon at the interface. A common and oft-used example of balance condition of type (a) was developed in the late 1960s by Beavers and Joseph [2]. That boundary condition, obtained experimentally, governs the behav-

ior of the velocity at the interface between a plain flow parallel to a porous medium,1 as illustrated in Fig. 1. The boundary condition at interface B is:


@u

aBJ ¼ pffiffiffiffi ðuB  QÞ: @y
B K

ð1Þ

The solution of this partial differential equation determines the ‘‘slip” velocity u at the interface, which is a function of the permeability K, the superficial velocity Q in the porous media, and the Beavers-Joseph parameter aBJ . This equation was later given a ‘‘statistical treatment” by Saffman [4], where it was shown that the velocity Q could be neglected to obtain:


@u

aBJ ¼ pffiffiffiffi uB : @y
B K

ð2Þ

A few years later [5], it was demonstrated that this equation could be generalized by relating the velocity at the wall to the shear stress:

  @u @ v aBJ þ ¼ pffiffiffiffi uB : @y @x B K

ð3Þ

A common example of type (b) boundary condition derivation is when a CFD code is coupled with a Material Response (MR) solver for studying aerothermodynamics problems, such as the surface balance equations presented by Milos and Rasky [6].

⇑ Corresponding author. E-mail addresses: [email protected] (A. Martin), [email protected]. edu.cn (H. Zhang), [email protected] (K.A. Tagavi). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.07.078 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

1 The curious reader is invited to learn about the history of the Beavers-Joseph equation by reading the well-documented paper of Nield [3].

993

A. Martin et al. / International Journal of Heat and Mass Transfer 115 (2017) 992–999

surface coupling have been based on a methodical analysis of all physical phenomena, which is then presented as surface balance equations [6]. Again, such an approach has the drawback of not being universal, varying from case to case depending on the equations that are solved in the respective domains. Moreover, certain phenomena are sometimes neglected without proper justification (e.g., the shear stress term in the energy equation [12]), or an entire conservation equation is neglected (e.g. the momentum equation reduced to zero pressure gradient [6]). However, instead of trying to document every possible pair of physical phenomenon that could occur at an interface, the present paper highlights an elegant and straightforward method to model surface equations: let the governing equations do the work. Nearly all governing equations are composed of three types of terms: the time dependent storage term, the flux term, and the source term. Both the time-dependent and the source terms are locally evaluated. By contrast, the flux terms are based on the transport of an intensive property through the enclosed boundaries of a volume. This principle has been defined by Reynolds, in his classic Reynolds Transport Theorem [14]:

Z

@Q dV þ @t

V

Fig. 1. Velocity profile of a plain flow over a porous material (image taken from Beavers and Joseph [2] reproduced with permission).

In this paper a systematic method of writing complex boundary conditions is presented. In order to demonstrate the feasibility of the method, first, a simple but well-established flow is used. The method is then applied to the much more complex problem of an aerothermodynamic flow coupled to a porous flow to demonstrate the versatility of the method.

As mentioned above, the modeling of surface phenomena between adjoining domains of different governing equations has proven to be challenging from the computational point of view. This happens when two independent solvers are specifically designed to solve different sets of equations where each is applicable only to their respective domain. Coupling two different sets of governing equations at the interface is not a trivial matter. The difficulty comes from two sources: (a) mathematically – the order of each set of differential equations results in a complex system to solve, and (b) numerically – the implementation of the resulting interface condition may require advanced cell mapping and reconstruction techniques, especially if the coupling scheme is used in a multi-processor environment. The previous works on this topic have focused on two main types of applications. The first type consists of incompressible, non-reacting flows [7–9,4,10] which results in a simplified system of equations where the mass and energy equations need not be explicitly solved. For this type of applications, the problem is simplified to determining values for pressure, normal velocity, and both tangential velocities at the interface. In most cases, the condition at the interface is obtained for very specific flows and cannot be extended to other flow configurations. A well-known example of such condition is the Beavers-Joseph equation, which, for instance, is not applicable to flows where the normal velocity component at the interface cannot be neglected. The second type of application addresses compressible, reacting flows: for example, the interaction of an atmospheric flow over the porous heat shield of a re-entry vehicle [11–13]. Previous works on

Z

Q ðu  nÞdA ¼ A

ð4Þ

S dV V

where u is the velocity vector and n the unit vector normal to surface dA. In this equation, the first term represents the timedependent storage term, the second, the flux term, and the third, the source term. In layman’s term: the change in time of quantity Q in volume V is given by the net amount of the quantity that enters through the boundaries A, plus any source of Q that may exist within the control volume. Eq. (4) can be re-written by defining a flux matrix F that includes the product Q and u, as well as other phenomena (e.g., pressure contributions in the momentum and energy equations, transport phenomena and grid convection):

Z

@Q dV þ @t

V

2. Flux-balancing at the boundary

I

I

Z

F  ndA ¼ A

ð5Þ

S dV V

If two domains – domain r and domain s – are connected via surface AB , the quantity F leaving domain r through area AB is exactly the same as that entering domain s. This concept, applied to Eq. (5), simply means that independent of the storage and source terms in each respective domain, the flux term on each side of the interface must be equal to each other. This principle can be expressed mathematically by:

Z

AB

 Z  F  n dA ¼ F  n dA r

AB

ð6Þ

s

From this analysis, it can be concluded that when imposing a boundary condition at the interface of two domains, Eq. (6) must be satisfied. This basic principle is further explained in the following sections, using two classic problems as examples. It should be noted that this derivation does not include surface source terms, some of which could take place at the interface. In such a case, Eq. (6) is rewritten as:

Z

AB

 Z  Z F  n dA ¼ F  n dA þ SA dA r

AB

s

ð7Þ

AB

where SA is the rate of the intrinsic quantity being generated at the control surface AB . 3. Example 1: coupling of Darcy’s law and Stokes equations To illustrate the concept of ‘‘flux-balancing”, as represented in Eq. (6) above, and to verify and validate its feasibility, the method is first applied to a classical case with the physical configuration illustrated in Fig. 1, which was used to derive the Beavers-Joseph

994

A. Martin et al. / International Journal of Heat and Mass Transfer 115 (2017) 992–999

equation. For this type of configuration, the porous flow is solved using Darcy’s law, and the plain flow using the Stokes equations, as detailed below.

where n is the normal vector of surface A bounding volume V. In this formulation, Eq. (14) therefore describes the momentum fluxes that leave (or enter) a porous volume V through surface A.

3.1. Darcy’s law

3.2. Stokes equations

A porous medium is usually defined as a fixed-in-space solid matrix with connected void spaces through which a fluid can flow [15]. The law governing the flow thorough a porous medium was obtained experimentally by Henry Darcy in 1856, while observing the filtration system of the aqueducts of the city of Dijon [16]. Darcy’s law is given by

For the proposed example, the plain flow is governed by the Stokes equations in 2D, assuming a Newtonian fluid. In matrix form, the steady-state form of these equations become:

@p l ¼  vs; @x K

ð8Þ

where l is the viscosity of the fluid and K is the permeability of the porous medium. The superficial velocity, v s , is calculated by averaging the velocity over the entire cross section of the porous medium. The superficial velocity is therefore defined as v s ¼ /u, where u is the fluid velocity averaged over the fluid area of the cross section, and / is the porosity. Although the rigorous way to derive Darcy’s law is to average the Navier–Stokes equation over the porous volume [17], it is also possible to obtain it using scale analysis. Assuming that the flow is Newtonian, steady-state, one-dimensional, and incompressible, the momentum equation in the x-direction, reduces to

0¼

@p @ þ @x @y



l

 @u : @y

ð9Þ

The solution of this equation with a no-slip velocity boundary condition is a Poiseuille flow, and takes the form of:

u¼

@p ðhy  y2 Þ ; @x 2l

ð10Þ

where h is the height (in y-direction) of the uniform cross section of the flow enclosure. An expression for the superficial velocity in the x-direction can be obtained by integrating over the cross section:

vs

2

@p h ¼ : @x 12l

ð11Þ

By comparing Eqs. (8) and (11), it can be seen that a permeability 2

can be derived as K ¼ h =12. Going back to Eq. (9), it can therefore be deduced that the characteristic length of the derivative is pro@ portional to the inverse of the square root of K, that is @y  p1ffiffiKffi or @ a p ffiffi ffi more specifically @y ¼  , where a is a proportionality constant. K

Utilizing the same flux form as for the plain flow, the 2D Darcy’s laws can be written as:



r

p

0

0

p

 

 lK u

 lK v

!

¼ 0:

ð12Þ

Using the result from the above scale analysis and setting @ @ ¼ @y ¼ paffiffiKffis , this equation can be formulated as: @x



r

p

0

0

p



 r  as l

pffiffivffi 2 puffiffiKffi uþ K vpþu ffiffiffi 2 pvffiffiffi K K

! ¼ 0:

ð13Þ

After integrating over volume V, and applying the divergence theorem, the following is obtained:

I "

p

0

0

p



 n  as l

pffiffivffi 2 puffiffiKffi uþ K vpþu ffiffiffi 2 pvffiffiffi K K

!

#

 n dA ¼ 0:

ð14Þ



r

p

0

0

p



2 @u @x

 lr 

@v @x

@u @y

þ @u @y

þ @@xv

! ¼ 0:

2 @@yv

ð15Þ

Similarly, this equation can be expressed as a momentum flux by integrating over volume V, while applying the divergence theorem:

I "

p 0 0



p

nl

2 @u @x

@v @x

@u @y

þ @u @y

þ @@xv

2 @@yv

!

#  n dA ¼ 0;

ð16Þ

3.3. Flux-balancing: parallel flow over a permeable wall It can be shown that the Beavers-Joseph boundary condition can be derived analytically by coupling the plain flow and porous flow at the interface. For the interface between the plain flow and the porous medium, the fluxes on both sides must be balanced, i.e. no fluxes are destroyed or created at the surface and that the fluxes only pass through the interface. Because the two adjacent domains share a contact interface of area AB and a normal vector nB , an equation at the interface is obtained by ‘‘balancing”, or equating, the fluxes from Eqs. (16) and (14):

"

¼

p 0



 0 p "  p 0

2l @u @x

l @u þ l @@xv @y

l @u þ l @@xv @y 2l @@yv

!#  nB f

as plffiffiKffi ðu þ v Þ  l 0 p as pffiffiKffi ðv þ uÞ 2as plffiffiKffi v 2as plffiffiKffi u

!#

ð17Þ  nB

p

Since the plain flow is overlaid on top of the porous medium, as illustrated in Fig. 1, then nB ¼ ð0; 1Þ, and the coupling equation becomes:

  @u @ v 1 þ ¼ as pffiffiffiffi ðu þ v Þjp @y @x f K     @v l p þ 2l ¼ p þ 2as pffiffiffiffi v @y f K p

ð18aÞ ð18bÞ

where the subscript f and p denote the fluid and porous medium sides, respectively. It can be seen that Eq. (18a) is the same as Eq. (3). If it is assumed that the velocity in the y-direction is zero, the above set of equations reduces to:

@uf as ¼ pffiffiffiffi up @y K pf ¼ pp ð19bÞ

ð19aÞ

Eq. (19a) is exactly the Beavers-Joseph boundary condition as expressed in Eq. (2). The present method also yields the continuity of pressure at the interface, Eq. (19b), which is implicitly assumed in most boundary condition treatments. This derivation, however, clearly shows that this equation is only valid for very specific types of flow (Stokes, in the plain flow domain, Darcy, in the porous media domain), for a very specific geometrical configuration (stratified flow over a porous medium), with a very specific assumption (no velocity in the y-direction). It is clear that this interface condition would not be appropriate to couple, for example, the Navier–Stokes

995

A. Martin et al. / International Journal of Heat and Mass Transfer 115 (2017) 992–999

equations with the Forchheimer equation for porous flow [18] since it is missing some of the key terms, such as the advection term. 4. Example 2: Coupling of Darcy-Brinkman and Navier-Stokes equations Another common flow situation is that of a porous material interacting with a reacting compressible flow. One typical example is when a spacecraft enters an atmosphere while traveling at hypersonic speeds [19,20,12]. Under such extreme conditions, the spacecraft’s thermal protection system (TPS) experiences high levels of aerodynamic heating from the flow around it. One type of TPS that is widely used is a charring ablative TPS, where a carbon porous structure is infused with a resin that pyrolyzes at high temperatures. As a result, relatively cooler pyrolysis gases are introduced into the flow field while leaving behind a carbon matrix, forming a porous medium [21]. Traditionally, a CFD code and a Material Response (MR) code are used to model the flow field and the TPS side, respectively. To fully and accurately analyze an atmospheric entry, however, a fully-coupled approach that can allow complete solutions for both the material and the fluid is desired. Such full scale analysis on an entire trajectory has not yet been performed. Both set of governing equations in the CFD and MR solvers can be written in the following form

@Q þ $  ðF adv  F diff Þ ¼ S; @t

and e and h are respectively the total energy and total enthalpy of the gas per unit mass. 4.2. Material response governing equations Material Response codes solve for gaseous mass, solid mass, momentum, and energy conservation equations. Thus, the terms in Eq. (20) for a porous medium are given by:

0 B B B B B B B B B B B Q ¼B B B B B B B B B B @

ð20Þ

where Q is a vector of the conserved variables, F adv and F diff are advective and diffusive fluxes, respectively, and S is a source vector. Integrating this last equation over a control volume results in Eq. (4).

F diff

4.1. Fluid dynamics governing equations For a reacting compressible flow in thermal equilibrium, the terms of Eq. (20) are:

0

qg

1

0

qg1 u

qg1 v

qg1 w 1

1 C B B . C .. .. .. C B B . C . . . C B B . C C B C B C B q u q v q w Bq C g ngs g ngs C B gngs B gngs C C B C F adv ¼ B 2 Q ¼B C q u þ p q v u q wu B qg u C g g C B g C B C B Bq v C 2 B qg uv qg v þ p qg wv C B g C C B C B B @ qg w A A @ qg uw qg v w qg w2 þ p C qg e qg uh qg v h qg wh 1 0 J1 1 0 x_ 1 C B .. C B B . C . C B B .. C C B C B C B J C B ngs C B C Bx _ C B ngs C B s s C B sxx xy xz C B S ¼ F diff ¼ B C B 0 C C B syx s s C B yy yz C B B 0 C C B C B szy szz C B szx C B C B 0 A @ ngs C B X A @ su  q_00  ðJi hi Þ 0

where s is the viscous tensor, J1 ; . . . ; Jngs are mass diffusion for each

species, and q_00 is the heat conduction vector. The subscript i, from 1 to ngs (total number of gas species), represents each gas species, _ i the mass production of species i. In these expreswhich makes x sions, qg is overall gas density, u ¼ ðu; v ; wÞ is the bulk velocity,

/qgngs

qs1 .. .

qsnss /qg u /qg v /qg w /qg e þ qs es

/qg1 u

/qg1 v

/qg1 w

/qg v h

/qg wh

1

C C B .. .. .. C C B . . . C C B C C B C C B /q u / q v / q w g ngs g ngs g ngs C C B C C B C C B 0 0 0 C C B C C B .. .. .. C C B . . . C C F adv ¼ B C C B 0 0 0 C C B C C B C C B /q u2 þ p / q v u / q wu g g C C B g C C B C B /qg uv /qg v 2 þ p /qg wv C C C B C C B A @ /qg uw /qg v w /qg w2 þ p A

/qg uh 1 0 x_ g1 B . C B . C 1 0 B . C J1 C B C Bx C B B _ gngs C .. C B C B . C B _ C B x C B B s1 C C B Jngs B . C C B C CS¼B ¼B B .. C 0 C B C B C B Bx _ snss C C B 0 C B C B C B ngs C B B Dx C A @ _ 00 X C B q  ðJi hi Þ B Dy C C B i¼1 C B @ Dz A SD

where es and qs is the total energy and total density in the solid phase, respectively. In the Darcy-Brinkman formulation, the porous medium viscous forces are treated as a source term in the momentum equations, denoted by Dx ; Dy , and Dz in Eq. (22), and may be calculated by solving the following linear equation:

K xx

B @ K yx K zx

i¼1

.. .

0

ð22Þ

0

ð21Þ

1

/qg1

K xy K yy K zy

K xz

10

Dx

1

0

u

1

B C CB C K yz A@ Dy A ¼ /l@ v A; K zz Dz w

ð23Þ

The three-by-three matrix on the left hand side is an anisotropic tensor of solid permeability. If the material is orthotropic, Eq. (23) can be simplified:

0

K xx

0

B @ 0 0

K yy

Dx ¼ 

/l u; K xx

0

0

10

Dx

1

0

u

1

B C CB C 0 A@ Dy A ¼ /l@ v A; K zz Dz w Dy ¼ 

/l v; K yy

Dz ¼ 

ð24Þ /l w: K zz

ð25Þ

The diffusive effect of porous media is modeled as a source term in the energy equation, which is given as:

SD ¼ Dx u þ Dy v þ Dz w

ð26Þ

As was done in the previous example with Eq. (13), a dimensional analysis moves back the Darcy forces and thermodynamic work from the source term to the diffusive fluxes:

996

A. Martin et al. / International Journal of Heat and Mass Transfer 115 (2017) 992–999

0

F adv

    @Y @Y _ 00 Y i  qg D i _ 00 Y i  qg D i ¼ m m @x f @x p

1

/qg1 w

C B .. .. .. C B . . . C B C B B /q u /qgngs v /qgngs w C g ngs C B C B C B 0 0 0 C B C B .. .. .. C B . . . ¼B C C B 0 0 0 C B C B B /q u2 þ p /qg v u /qg wu C C B g C B 2 C B /qg uv / q v þ p / q w v g g C B C B 2 @ /qg uw /qg v w /qg w þ p A /qg v h

/qg uh

0

F diff

/qg1 v

/qg1 u

which does not make any assumptions about the individual contributions of mass transport phenomenon. But if the mass diffusion term is neglected on the porous side, Eq. (5a) of Milos and Rasky [6] is obtained:



qg D

/qg wh 1

J1 .. .

ð27Þ

C B C B C B C B C B Jngs C B C B B0 0 0C C B .. .. C B .. C B. . . C B ¼B 0 0C C B0 C B C B C B C B D C B F C B C B C B ngs C B X A @ 00 WF  q_  ðJi hi Þ

@Y i @x



_ 00 ðY i;f  Y i;p Þ: ¼m

However, it might not always be proper to neglect mass diffusion in the porous medium flow [22,23]. It should be noted that Eq. (5b) of Milos and Rasky [6] can be obtained by including a surface source term for surface chemical reactions.

h

K xx

K zz

2u K xx

f

h

B v B þ puffiffiffiffiffi DF ¼ as /lB pffiffiffiffiffi K yy @ K xx pwffiffiffiffiffi þ puffiffiffiffiffi



i

qg u2 þ p  sxx ¼ /qg u2 þ p  as /l pffiffiffiffiffiffiffi

where

2 puffiffiffiffiffi K xx

ð34Þ

f

4.3.2. Momentum balance Most Material Response (MR) codes do not solve the momentum equation explicitly, but rather implicitly include Darcy’s law in energy and mass conservation. Thus, a Surface Momentum Balance (SMoB) equation is usually not necessary to close the system of equations. However, when coupling to a CFD code, a momentum boundary condition needs to be applied. Using the approach adopted in this work, the momentum fluxes at the interface are therefore:

i¼1

0

puffiffiffiffiffi K yy

þ pvffiffiffiffiffi K xx

2 pvffiffiffiffiffi K yy

pwffiffiffiffiffi K yy

þ pvffiffiffiffiffi K zz

ð33Þ

puffiffiffiffiffi K zz

þ

i

ð35aÞ p

v

!#

u K yy p    h i w u qg uw  szx ¼ /qg wu  as /l pffiffiffiffiffiffiffi þ pffiffiffiffiffiffiffi f K zz p K xx

qg uv  syx ¼ /qg v u  as /l pffiffiffiffiffiffiffi þ pffiffiffiffiffiffiffi

pwffiffiffiffiffi 1 ð28Þ

K zz

ð35bÞ

K xx

f

K

xx C pvffiffiffiffiffi þ pwffiffiffiffiffi C K yy C K zz A 2 pwffiffiffiffiffi

"



ð35cÞ

and WF ¼ uDF .

This set of equations is much more detailed than the continuity of pressure condition [24–26] that is used in most coupling approaches:

4.3. Flux-balancing: aerothermodynamic flow over a porous material

pf ¼ pp

Using the technique demonstrated in Section 2, the balance of fluxes yields:

½F adv  n  F diff  nf ¼ ½F adv  n  F diff  np

ð29Þ

Without making any assumptions on the direction of the flow, the interface between the porous flow and the plain flow can be arbitrarily oriented in the x-direction. The normal vector of that surface is therefore n ¼ ð1; 0; 0Þ. 4.3.1. Mass balance Applying Eq. (29) to the conservation of mass for species i, using the terms of Eqs. (21) and (27), the interface condition is therefore:

h

i

h

qgi u þ Jx;i ¼ /qgi u þ Jx;i f

i

p

ð30Þ

The overall mass flow rate, obtained by summing over all species, is conserved at the interface, that is:

h i h i _ 00 ¼ qg u ¼ /qg u m f

p

ð31Þ

and the mass diffusion term for the species on the fluid side is calculated by Fick’s law as

J x;i ¼ qg D

@Y i @x

ð32Þ

Combining Eqs. (30) with (32) yields the Surface Mass Balance (SMB) equation:

ð36Þ

or the slightly more accurate total pressure condition [27,11]:

h

i

h

i

qg u2 þ p ¼ /qg u2 þ p f

ð37Þ

p

4.3.3. Energy balance Applying the technique to the energy terms yields the Surface Energy Balance (SEB) equation:

h

qg uh  ðsxx u þ sxy v þ sxz wÞ þ q_ 00x þ 

l 

X

¼ /qg uh  as / pffiffiffiffiffiffiffi u þ v þ w K xx 2

2

2



J x;i hi þ

q_ 00x

i f

þ

X

 ð38Þ

J x;i hi p

In this equation, it is assumed that puffiffiffiffiffi þ pvffiffiffiffiffi þ pwffiffiffiffiffi  0, per mass K xx

K yy

K zz

conservation of an incompressible fluid. As with the SMB, simplifying using the total conservation of mass (i.e., Eq. (31)) and neglecting the work performed by shear on both sides of the surface reduces this equation to Eq. (6) of Milos and Rasky [6]:

" kT

@T T @T v X @Y q g hi D i  kv  @x @x @x i

# _ 00 ðhp  hf Þ ¼ q_ 00x þ m

ð39Þ

f

It is to be noted that Milos and Rasky [6] include radiative fluxes and emission, both omitted in the present derivation. They also express conduction and thermal diffusion explicitly using Fourier’s and Fick’s law, respectively.

997

A. Martin et al. / International Journal of Heat and Mass Transfer 115 (2017) 992–999

5. Example 3: Navier–Stokes equations with an impermeable, adiabatic wall The procedure outlined in Section 2 can also be used to define an adiabatic, impermeable wall boundary condition for the Navier–Stokes equations. For the sake of simplicity, a single–species incompressible flow is considered. Starting from Eq. (7), domain r fluxes are thus represented by the Navier–Stokes equations, and domain s fluxes are null. For this example, the surface source term SA is needed for the momentum equations. This term, defined in the following as R00 , represents the reaction forces per unit area on the surface (for instance, in aerodynamic applications, the drag, lift, and lateral forces). The surface conditions are obtained by applying the fluxes of Eq. (21) to a surface oriented in the x-direction. Therefore, for the mass balance results in:

h

i qg u ¼ 0

Considering that the gas density cannot be 0, this leads to a zero velocity in the x-direction at the surface. The momentum flux balance then becomes:

h h

i

qg u2 þ p  sxx ¼ R00x

ð41aÞ

qg uv  syx ¼ R00y

ð41bÞ

qg uw  szx ¼ R00z

ð41cÞ

f

i

f

i

f

and the energy flux becomes:

h

F grid ¼ Qug

0

f

½p  sxx f ¼ R00x

 syx f ¼ R00y

ð43bÞ

½szx f ¼ R00z

ð43cÞ

6. Example 4: Coupling of Navier-Stokes equations and a surface ablating material In the field of aerothermodynamics, it is often necessary to account for the interaction of an ablating surface with a high enthalpy flow. This is especially true when studying thermal protection systems, whether it is from the fluid perspective [28,29] or from the material response perspective [12,11]. For this type of problems, the surface of the material is recessing, and mass, energy and momentum are transferred into the flow field. Therefore, domain r is a fluid flow, and domain s, a recessing solid wall. The fact that the surface is recessing introduces an additional flux which was not included in Example 3. This flux has taken various names in the literature, but is often described as a grid advection flux. In terms of control volume analysis, this flux is obtained

C C C C qgngs wg C C qg wug C C C qg wv g C C C qg wwg A .. .

C B .. .. .. C B . . . C B C B B q u qgngs v qgngs w C C B gngs C B ¼ B q u2 þ p C q u v q uw g g C B g C B 2 C B qg v u q v þ p q v w g g C B B A @ qg wu qg wv qg w2 þ p C 0

F diff

ð43aÞ

ð44Þ

.. .

qg1 wg 1

qgngs v g qg v ug qg vv g qg v wg qg eug qg ev g qg ewg 0 qg1 u qg 1 v qg 1 w 1

ð42Þ

This set of equations represents the condition at the wall when solving the Navier–Stokes equations. It should be noted that no assumptions have been made about the velocity at the wall, which makes these equations valid whether or not a no-slip condition is enforced. By inserting the mass conservation equation (qg u ¼ 0) into both the momentum and energy equations above, and by enforcing a no-slip velocity (u ¼ 0), a more familiar set of equations is obtained

 

00  @T q_ x f  k ¼0 @x f

F grid

qg1 v g

qg1 ug

B . B . B . B B q ug B gngs ¼B B qg uug B B q uv g B g B @ qg uwg

F adv

i

qg uh  ðsxx u þ sxy v þ sxz wÞ þ q_ 00x ¼ 0

ð45Þ

The fluxes from the fluid domain are therefore composed of the grid advection flux, and both advective and diffusive fluxes of Eq. (21):

ð40Þ

f

h

by implicitly solving the Geometric Conservation Law [30], and is often used in CFD codes in order to modify a computational grid while calculating the solution. More details of this technique are given in Ref. [31], and applications to control volume problems are presented in Ref. [32]. From a physical point of view, the grid movement only affects conserved variables, and are therefore described similarly to the advective fluxes. More details on the procedure necessary to obtain this result can be found in Ref. [12]. Using the nomenclature previously defined, and defining ug as the wall (or grid) velocity vector, the grid convection flux F grid can be expressed as:

B B B B B B B B ¼B B B B B B B @

qg hv

qg hu

J1 ...

sxx syx szx

Jngs

sxy syy szy

su  q_ 00 

sxz syz szz

ngs X ðJi hi Þ

1

ð46Þ

qg hw

C C C C C C C C C C C C C C C A

i¼1

From the material side, the surface fluxes can simply be obtained by setting the gas velocity to 0 in Eq. (22):

0

F grid

qs1 ug B B B B B q ug B sngs ¼B B B B B B @ qs es ug

1 0 B .. C B . C C B B 0 C C B C B ¼B 0 C C B B 0 C C B C B @ 0 A 0

F diff

q_ 00

qs1 v g .. .

qsngs v g 0 0 0

qs es v g

1

0 1 0 C B .. C C C B C B.C C B0C C qsngs wg C B C C C F adv ¼ B B0C C B C C B0C C B C C B C C @0A A 0 qs es wg

qs1 wg

ð47Þ

998

A. Martin et al. / International Journal of Heat and Mass Transfer 115 (2017) 992–999

These fluxes can thus be used in Eq. (7) to obtain the condition at the surface. Similar to the previous examples, a surface oriented in the x-direction, and a surface reaction force R00 are used. For the mass balance, the following surface condition is obtained:

h

i



qgi ðu  ug Þ þ Jx;i ¼ qi ug



f

ð48Þ

s

For the momentum balance, the equation becomes:

h h h

i

qg uðu  ug Þ þ p  sxx ¼ R00x

ð49aÞ

qg v ðu  ug Þ  syx ¼ R00y

ð49bÞ

qg wðu  ug Þ  szx ¼ R00z

ð49cÞ

f

i

f

i

f

And for the energy:

h

qg ðuh  ug eÞ  ðsxx u þ sxy v þ sxz wÞ þ q_ 00x þ

¼ qs ug es þ q_ 00x



X

J x;i hi

i f

ð50Þ

s

To relate these equations to practical applications, ug could either be imposed for an inert receding wall, or calculated using reaction kinetics for an ablating wall. It should be noted that these equations are related to boundary conditions often found in the relevant literature [29,28]. For instance, by neglecting the gas velocity at the wall, these equations reduce to:

h

qgi ug þ J x;i

i f

 ¼ qi ug s

ð51Þ

½p  sxx f ¼ R00x

 syx f ¼ R00y

ð52bÞ

½szx f ¼

ð52cÞ

h

ð52aÞ

R00z

qg ug e þ q_ 00x þ

X

J x;i hi

i f

 ¼ qs ug es þ q_ 00x s

ð53Þ

7. Conclusion It is demonstrated that complex boundary conditions need not be derived one example at a time, a process that is normally ‘‘painful” and time consuming. Rather, a systematic and universal method of matching governing equations at the boundaries and interfaces can provide a rigorous approach that can consistently result in the correct set of boundary and interface conditions. This method is successfully tested against four well-known sets of flow conditions, namely Beavers-Joseph flow, an aerothermodynamic interface flow, an adiabatic impermeable wall, and an ablating surface. It is expected that when applied to other flow configurations, independent of the complexity of the flows, the underlying method presented here will result in an accurate and useful set of boundary and interface conditions. Conflict of interest The authors declared that there is no conflict of interest. Acknowledgments Financial support for this work was provided by NASA Award NNX15AD73G and NNX13AN04A. The authors are also grateful to G.V. Candler and S.C.C. Bailey for insightful discussions and comments.

References [1] J.R. Shewchuk, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain, Technical Report CMU-CS-TR-94-125, Carnegie Mellon University, 1994. [2] G.S. Beavers, D.D. Joseph, Boundary conditions at a naturally permeable wall, J. Fluid Mech. 30 (01) (1967) 197–207, http://dx.doi.org/10.1017/ S0022112067001375. [3] D.A. Nield, The Beavers–Joseph boundary condition and related matters: a historical and critical note, Transp. Porous Media 78 (3) (2009) 537–540, http://dx.doi.org/10.1007/s11242-009-9344-y. [4] P.G. Saffman, On the boundary condition at the surface of a porous medium, Stud. Appl. Math. 50 (2) (1971) 93–101, http://dx.doi.org/10.1002/ sapm197150293. [5] I.P. Jones, Low Reynolds number flow past a porous spherical shell, Math. Proc. Camb. Philos. Soc. 73 (1) (1973) 231–238, http://dx.doi.org/10.1017/ S0305004100047642. [6] F.S. Milos, D. Rasky, Review of numerical procedures for computational surface thermochemistry, J. Thermophys. Heat Transf. 8 (1) (1994) 24–34, http://dx. doi.org/10.2514/3.497. [7] J. Ochoa-Tapia, S. Whitaker, Momentum transfer at the boundary between a porous medium and a homogeneous fluid – I. Theoretical development, Int. J. Heat Mass Transf. 38 (14) (1995) 2635–2646, http://dx.doi.org/10.1016/00179310(94)00346-W. [8] J. Ochoa-Tapia, S. Whitaker, Momentum transfer at the boundary between a porous medium and a homogeneous fluid – II. Comparison with experiment, Int. J. Heat Mass Transf. 38 (14) (1995) 2647–2655, http://dx.doi.org/10.1016/ 0017-9310(94)00347-X. [9] M. Cieszko, J. Kubik, Derivation of matching conditions at the contact surface between fluid-saturated porous solid and bulk fluid, Transp. Porous Media 34 (1–3) (1999) 319–336, http://dx.doi.org/10.1023/A:1006590215455. [10] M.L. Bars, M.G. Worster, Interfacial conditions between a pure fluid and a porous medium: implications for binary alloy solidification, J. Fluid Mech. 550 (2006) 149–173, http://dx.doi.org/10.1017/S0022112005007998. [11] A. Martin, I.D. Boyd, Modeling of heat transfer attenuation by ablative gases during the Stardust re-entry, J. Thermophys. Heat Transf. 29 (3) (2015) 450– 466, http://dx.doi.org/10.2514/1.T4202. [12] A. Martin, I.D. Boyd, Strongly coupled computation of material response and nonequilibrium flow for hypersonic ablation, J. Spacecraft Rock. 52 (1) (2015) 89–104, http://dx.doi.org/10.2514/1.A32847. [13] D. Kuntz, B. Hassan, D. Potter, Predictions of ablating hypersonic vehicles using an iterative coupled fluid/thermal approach, J. Thermophys. Heat Transf. 15 (2) (2001) 129–139, http://dx.doi.org/10.2514/2.6594. [14] O. Reynolds, The Sub-mechanics of the Universe, vol. III, Cambridge University Press, 1903, pp. 12–13. [15] D.A. Nield, Modelling fluid flow and heat transfer in a saturated porous medium, J. Appl. Math. Decis. Sci. 4 (2) (2000) 165–173, http://dx.doi.org/ 10.1155/S1173912600000122. [16] H. Darcy, Les Fontaines Publiques de la Ville de Dijon, Dalmont, Paris, France, 1856, p. 570. [17] S. Whitaker, Flow in porous media I: A theoretical derivation of Darcy’s law, Transp. Porous Media 1 (1) (1986) 3–25, http://dx.doi.org/10.1007/ BF01036523. [18] P.H. Forchheimer, Wasserbewegung durch Boden, Z. Vereines Deutscher Ingenieure 45 (1901) 1782–1788. [19] A. Martin, I.D. Boyd, Non-Darcian behavior of pyrolysis gas in a thermal protection system, J. Thermophys. Heat Transf. 24 (1) (2010) 60–68, http://dx. doi.org/10.2514/1.44103. [20] A. Martin, L.C. Scalabrin, I.D. Boyd, High performance modeling of an atmospheric re-entry vehicles, J. Phys.: Conf. Ser. 341 (1). http://dx.doi.org/ 10.1088/1742-6596/341/1/012002 (Article 012002). [21] W.H. Bowman, R.M. Lawrence, Ablative materials for high-temperature thermal protection of space vehicles, J. Chem. Educ. 48 (10) (1971) 690–691, http://dx.doi.org/10.1021/ed048p690. [22] A. Martin, I.D. Boyd, I. Cozmuta, M.J. Wright, Chemistry model for ablating carbon-phenolic material during atmospheric re-entry, in: 48th AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper 2010-1175, Orlando, FL, 2010. http://dx.doi.org/10.2514/6.2010-1175. [23] J.R. Lachaud, N.N. Mansour, A. Ceballos, D. Pejakovic´, L. Zhang, J. Marschall, Validation of a volume-averaged fiber-scale model for the oxidation of a carbon-fiber preform, in: 42nd AIAA Thermophysics Conference, AIAA Paper 2011-3640, Honolulu, HI, 1–10, 2011. http://dx.doi.org/10.2514/6.2011-3640. [24] T. Suzuki, K. Fujita, T. Sakai, Numerical analysis of graphite ablation in nitrogen flow, in: 46th AIAA Aerospace Sciences Meeting and Exhibit, AIAA-2008-1217, Reno, NV, 2008. http://dx.doi.org/10.2514/6.2008-1217. [25] H. Kihara, N. Hirata, K. Abe, A study of thermal response and flow field coupling simulation around Hayabusa capsule loaded with light-weight ablator, Open J. Fluid Dyn. 3 (2A) (2013) 100–107, http://dx.doi.org/10.4236/ ojfd.2013.32A016. [26] O. Joshi, P. Leyland, Stability analysis of a partitioned fluid–structure thermal coupling algorithm, J. Thermophys. Heat Transf. 28 (2014) 59–67, http://dx. doi.org/10.2514/1.T4032. [27] P.A. Gnoffo, C.O. Johnston, R.A. Thompson, Implementation of radiation, ablation, and free energy minimization modules for coupled simulations of hypersonic flow, in: 47th AIAA Aerospace Sciences Meeting and Exhibit, AIAA

A. Martin et al. / International Journal of Heat and Mass Transfer 115 (2017) 992–999 Paper 2009-1399, Orlando, FL, 12, 2009. http://dx.doi.org/10.2514/6.20091399. [28] L. Trevino, G.V. Candler, Numerical simulation of surface patterns on sublimating ablative materials, in: 53rd AIAA Aerospace Sciences Meeting, AIAA Paper 2015-1452, Kissimmee, FL, 2015. http://dx.doi.org/10.2514/6. 2015-1452. [29] C.H. Mortensen, X. Zhong, Real-gas and surface-ablation effects on hypersonic boundary-layer instability over a blunt cone, AIAA J. 54 (3) (2016) 980–998, http://dx.doi.org/10.2514/1.j054404.

999

[30] P.D. Thomas, C.K. Lombard, Geometric conservation law and its application to flow computations on moving grids, AIAA J. 17 (10) (1979) 1030–1037, http:// dx.doi.org/10.2514/3.61273. [31] J.-Y. Trepanier, M. Reggio, H. Zhang, R. Camarero, A finite-volume method for the Euler equations on arbitrary Lagrangian-Eulerian grids, Comput. Fluids 20 (4) (1991) 399–409, http://dx.doi.org/10.1016/0045-7930(91)90081-R. [32] H. Zhang, M. Reggio, J.-Y. Trépanier, R. Camarero, Discrete form of the GCL for moving meshes and its implementation in CFD schemes, Comput. Fluids 22 (1) (1993) 9–23, http://dx.doi.org/10.1016/0045-7930(93)90003-R.