Conjugate heat transfer using modified interface conditions
for the Navier-Stokes equations
Jan Nordstr¨om†∗
Department of Mathematics, Link¨oping University, SE-581 83 Link¨oping, Sweden
Jens Berg†
Department of Information Technology, Uppsala University, SE-751 05, Uppsala, Sweden
Abstract
This paper evaluates the use of the compressible Navier-Stokes equations, with pre-scribed zero velocities, as a model for heat transfer in solids. In particular in con-nection with conjugate heat transfer problems.
We derive estimates, and show how to choose and scale the coefficients of the energy part in the Navier-Stokes equations, such that the difference between the energy equation and the heat equation is minimal.
A rigorous analysis of the physical interface conditions for the conjugate heat transfer problem is performed and energy estimates are derived in non-standard L2 -equivalent norms. The numerical schemes are proven energy stable with the physical interface conditions and the stability of the schemes are independent of the order of accuracy.
We have performed computations of a conjugate heat transfer problem in two different ways. One where, as traditionally, the heat transfer in the solid is governed by the heat equation. The other where the heat transfer in the solid is governed by the Navier-Stokes equations. The simulations are compared and the numerical results corroborate the theoretical results.
1. Introduction
Heat transfer is an important factor in many fluid dynamics applications. Flows are often confined within some material with heat transfer properties. Whenever there is a temperature difference between the fluid and the confining solid, heat will be transfered and change the flow properties in a non-trivial way. This interaction and heat exchage is referred to as the conjugate heat transfer problem [1, 2, 3, 4]. Examples of application areas include cooling of turbine blades and nuclear reactors, atmospheric reentry of spacecrafts and gas propulsion micro thrusters for precise satellite navigation.
The standard method of computing conjugate heat transfer problems is to solve the Navier-Stokes equations in the fluid and couple it with the heat equation in the
∗Corresponding author: Jan Nordstr¨om, E-mail: jan.nordstrom@liu.se
†Parts of this work were completed while the authors visited the Centre for Turbulence Research at Stanford University
solid. This intuitive method is appealing due to the simplicity of the heat equation and has been considered before in e.g. [1, 2, 3, 4]. There are, however, some drawbacks that we shall address in this paper:
• The fluid and solid domain requires different physics solvers.
• Additional codes might be required to handle the communication between the solvers.
• Conjugate heat transfer cannot be easily implemented into existing multi-block Navier-Stokes codes.
As an attempt to circumvent these drawbacks, we will investigate how, and if, the compressible Navier-Stokes equations themselves can model heat transfer in solids. We will show how to scale and choose the coefficients of the energy part of the Navier-Stokes equations, such that it is as similar to the heat equation as possible. Numerical simulations of heat transfer is solids are performed to show the similari-ties, and differences, of the temperature distributions obtained by the Navier-Stokes equations and the heat equation. We will not overwrite, or strongly force, the veloc-ities in the Navier-Stokes equations to zero in each time integration stage since that would ruin the stability of the numerical method that we use. Instead, the velocities will be enforced weakly at the boundaries and interfaces only.
Conjugate heat transfer problems have been computed using a variety of meth-ods. For the stationary problem, methods include the finite volume method [5], the finite element method [6] and the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) [7]. For the unsteady problem, overlapping grids have been used [3]. The interface conditions have been imposed either strongly, weakly or by a mixture of both.
The method we use in this paper is a finite difference method on Summation-By-Parts (SBP) form with the Simultaneous Approximation Term (SAT) for imposing the boundary and interface conditions weakly. The SBP-SAT method has been used for a variety of problems and has proven to be robust [8, 9, 10, 11, 12, 13, 14, 15, 16, 17].
The SBP finite difference operators were originally constructed by Kreiss and Sc-hearer [18] with the purpose of constructing an energy stable finite difference method [19]. Together with the weak imposition of boundary [20] and interface [21] condi-tions, the SBP-SAT provides a method for constructing energy stable schemes for linearly well-posed initial-boundary value problems [22]. There are SBP operators based on diagonal norms for the first [23] and second [24, 25] derivative accurate of order 2, 3, 4 and 5, and the stability analysis we will present is independent of the order of accuracy.
In the previous literature, a mathematical investigation of the interface condi-tions in terms of well-posedness of the continuous equacondi-tions and stability of the resulting numerical scheme has not been performed to our knowledge. In [1], stabil-ity is proven for a second order accurate finite difference method by analysing the simplified case of two coupled heat equations. We shall in this paper hence focus on the mathematical treatment of the fluid-solid interface rather than computing physically relevant scenarios.
The paper consists of two major parts. Section 2 and 3 are devoted to the Navier-Stokes equations and the heat equation on a single domain to verify that the Navier-Stokes equations are indeed able to accurately model heat transfer in solids. Sections 4 and 5 are devoted to the conjugate heat transfer problem. Since there are no analytical solutions available, or easily constructed solutions, we shall numerically solve the conjugate heat transfer problem in two ways:
1. By coupling the Navier-Stokes equations across an interface and let the solu-tions in the different subdomains represent the fluid and solid, respectively. 2. By, as conventionally, coupling the Navier-Stokes equations for the fluid with
the heat equation for the solid.
By having the two solutions to the conjugate heat transfer problem, we can study the effects of using the Navier-Stokes equations throughout the entire computational domain, either in the fluid or the solid.
2. Similarity of the Navier-Stokes equations and the heat equation The two-dimensional compressible Navier-Stokes equations in conservative form are
qt+ Fx+ Gy = 0 (1)
where q = [ρ, ρu, ρv, e]T are the density, x- and y-directional momentum and energy, respectively. The energy is given by
e = p
γ − 1+ 1 2ρ u
2+ v2 , (2)
where p is the pressure and γ is the ratio of specific heat. Furthermore, we have
F = FI− FV, G = GI− GV. (3)
where the superscript I denotes the inviscid part of the flux and V the viscous part. The components of the flux vectors are given by
FI = [ρu, p + ρu2, ρuv, u(p + e)]T, GI = [ρv, ρuv, p + ρv2, v(p + e)]T, FV = [0, τxx, τxy, uτxx+ vτxy+ κTx]T,
GV = [0, τxy, τyy, uτyx+ vτyy + κTy]T,
(4)
where we have the temperature T and the thermal conductivity coefficient κ. The stress tensor is given by
τxx = 2µ ∂u ∂x+λ ∂u ∂x + ∂v ∂y , τyy = 2µ ∂v ∂y+λ ∂u ∂x + ∂v ∂y , τxy = τyx = µ ∂u ∂y + ∂v ∂x , (5) where µ and λ are the dynamic and second viscosity, respectively. To close the system we need to include an equation of state, for example the ideal gas law
Here R = cP− cV is the specific gas constant and cp and cv the specific heat capacity
for constant pressure and volume, respectively.
If we let u = v = 0, all the convective terms and viscous stresses are zero and by using (2) and (6), equation (1) reduces to
ρt = 0 px = 0 py = 0 Tt = 1 cVρ (κTx)x+ (κTy)y . (7)
The last equation is similar, but not identical, to the variable coefficient heat equa-tion.
For ease of comparison with the heat equation we transform to non-dimensional form as follows: u = u ∗ c∗ ∞ , v = v ∗ c∗ ∞ , ρ = ρ ∗ ρ∗ ∞ , T = T ∗ T∗ ∞ , (8) p = p ∗ ρ∗ ∞(c∗∞)2 , e = e ∗ ρ∗ ∞(c∗∞)2 , λ = λ ∗ µ∗ ∞ , µ = µ ∗ µ∗ ∞ , (9) cP = c∗P c∗P ∞, cV = c∗V c∗P ∞, R = R∗ c∗P ∞, κ = κ∗ κ∗ ∞ , (10) x = x ∗ L∗ ∞ , y = y ∗ L∗ ∞ , t = c ∗ ∞ L∗ ∞ t∗, (11)
where the ∗-superscript denotes a dimensional variable and the ∞-subscript the reference value. L∗∞ is a characteristic length scale and c∗∞ is the reference speed of
sound. The equation of state (6) becomes in non-dimensional form
γp = ρT. (12)
By using (8)-(11), the last equation in (7) becomes
Tt= 1 P ec 1 cVρ (κTx)x+ (κTy)y (13) where P ec= c∗∞L∗∞ α∗ ∞ , α∗∞ = κ∗∞ ρ∗ ∞c∗P ∞ (14) are the P´eclet number based on the reference speed of sound and the thermal diffu-sivity respectively.
The heat equation, non-dimensionalized using (8)-(11), can be written as
˜ Tt= 1 P ec 1 csρs κsT˜x x+ κsT˜y y (15)
where P ec is defined in (14) and cs, ρs, κs are the specific heat capacity, density
and thermal conductivity of the solid, respectively. In this case, all coefficients are constant but rewritten in a form which resembles (13).
In order to compare (13) and (15), we define β = P ecρcV, βs = P ecρscs and rewrite (13) and (15) as βTt= (κTx)x+ (κTy)y, (16) β ˜Tt= β βs κsT˜x x +κsT˜y y . (17)
Note that βs is constant for the solid. Since β > 0 and ∂β∂t = 0 we can estimate the
difference T − ˜T in the β-norm defined by ||T − ˜T ||2β = Z Ω T − ˜T 2 βdΩ (18)
where Ω is the computational domain. By subtracting (17) from (16), multiplying with T − ˜T and integrating over Ω we obtain
1 2 d dt||T − ˜T || 2 β = − Z Ω κ∇T · ∇T + κs βs β∇ ˜T · ∇ ˜T + I ∂Ω T − ˜T κ∇T −κs βs β∇ ˜T · nds + Z Ω κs βs T − ˜T∇β · ∇ ˜T dΩ + Z Ω κ +κs βs β ∇T · ∇ ˜T dΩ. (19)
In order to obtain as similar temperature distributions from the heat equation and Navier-Stokes equation as possible, the right-hand-side of (19) has to be less than or equal to zero. Note that we specify the same boundary data for T and ˜T , in which case the boundary integral is zero. By further assuming that ∇β = 0 we can rewrite (19) as the quadratic form
1 2 d dt||T − ˜T || 2 β = − Z Ω ∇T ∇ ˜T T κ −1 2 κ + κs βs β −1 2 κ +κs βs β κs βs β ∇T ∇ ˜T . (20) By computing the eigenvalues of the symmetric matrix in the quadratic form and requiring that they be non-negative, we can conclude that we need κ − κs
βsβ = 0.
Thus, if the relations
κ β − κs βs = 0, ∇β = 0 (21) hold, then d dt||T − ˜T || 2 β ≤ 0 (22)
and the Navier-Stokes equations and the heat equation produces the exact same solution for the temperature if given identical initial data. However, the relations (21) hold only approximately and not in general due to compressibility and possible non-zero velocities in the Navier-Stokes equations. We will numerically investigate the deviations in later sections.
3. The discrete problem on a single domain
The Navier-Stokes equations and the heat equation are discretized using a finite difference method on Summation-By-Parts (SBP) form with weak boundary condi-tions using the Simultaneous Approximation Term (SAT). In the basic formulation, the first derivative is approximated by
ux ≈ Dv = P−1Qv, (23)
where v is the discrete grid function approximating u. The matrix P is symmetric, positive definite and defines a discrete norm by ||v||2 = vTP v. In this paper, we consider diagonal norms only. The matrix Q is almost skew-symmetric and satisfies the SBP property Q + QT = diag[−1, 0, . . . , 0, 1]. There are SBP operators based on
diagonal norms with 2nd, 3rd, 4th and 5th order accuracy, and the stability analysis does not depend on the order of the operators [23, 26]. The second derivative is approximated here by using the first derivative twice, i.e.
uxx ≈ D2v = (P−1Q)2v. (24)
This yields a non-compact second derivative operator which simplifies the program-ming. Compact formulations of the second derivative operator are derived in [24, 25]. In order to extend the operators to higher dimensions, it is convenient to intro-duce the Kronecker product. For arbitrary matrices A ∈ Rm×n and B ∈ Rp×q, the
Kronecker product is defined as
A ⊗ B = a1,1B . . . a1,mB .. . . .. ... an,1B . . . am,nB . (25)
The Kronecker product is bilinear, associative and obeys the mixed product property
(A ⊗ B)(C ⊗ D) = (AC ⊗ BD) (26)
if the usual matrix products are defined. For inversion and transposing we have
(A ⊗ B)−1,T = A−1,T ⊗ B−1,T (27)
if the usual matrix inverse is defined. The Kronecker product is not commutative in general, but for square matrices A and B there is a permutation matrix R such that
A ⊗ B = RT(B ⊗ A)R. (28)
Let Px,y, Qx,y and Dx,ydenote the difference operators in the coordinate direction
indicated by the subscript. The extension to multiple dimensions is done by using the Kronecker product as follows:
¯ Px = Px⊗ Iy, Q¯x = Qx⊗ Iy, ¯ Py = Ix⊗ Py, Q¯y = Ix⊗ Qy, ¯ Dx = Dx⊗ Iy, D¯y = Ix⊗ Dy, ¯ D2 x = D2x⊗ Iy, D¯y2 = Ix⊗ Dy2. (29)
Due to the mixed product property (26), the operators commute in different co-ordinate directions and hence differentiation can be performed in each coco-ordinate direction independently. The norm is defined by
||T||2 = TTP T¯ (30)
where ¯P = ¯PxP¯y = Px⊗ Py.
3.1. The heat equation
The main theoretical tool for analyzing well-posedness of the partial differential equations, or stability of the resulting numerical schemes, is the energy method in continuous or discrete form.
We multiply (15) with T and integrate over the spatial domain. By using the Green-Gauss formula for higher dimensional integration by parts we obtain
||T ||2
t + 2 ˜αs ||Tx||2+ ||Ty||2 = 2˜αs
I
∂Ω
(T Tx, T Ty) · nds (31)
where ˜αs = P e1ccκsρss and n is the outward pointing normal vector. The norm is the
standard L2-norm defined by ||T ||2 =R
ΩT
2dΩ. The terms T T
x and T Ty are
indefi-nite everywhere along the boundary and hence needs to be supplied with boundary conditions to obtain a bounded energy growth. We shall use the Dirichlet boundary condition
T (x, y, t) = g(x, y, t), x, y ∈ ∂Ω. (32)
By assuming that g ≡ 0 we obtain the energy estimate ||T ||2
t + 2 ˜αs ||Tx||2+ ||Ty||2 = 0 (33)
and a well-posed problem [27, 28, 29, 22].
Let the spatial domain of interest be the unit square, 0 ≤ x, y, ≤ 1, discretized with M + 1 and N + 1 equidistant grid points, respectively. A semi-discretization of (15) using the SBP-SAT technique can be written as
Tt= ˜αs D¯2xT + ¯D 2 yT + σ1P¯−1E¯x0,y(T − gx0,y) + σ2P¯ −1¯ ExM,y(T − gxM,y) + σ3P¯−1E¯x,y0(T − gx,y0) + σ4P¯ −1¯ Ex,yN(T − gx,yN) (34)
where ¯Ex0,y = Ex0,y ⊗ Iy, Ex0,y = diag[1, 0, . . . , 0], is an operator which selects the
west boundary elements as indicated by the subscript. gx0,y = gx0,y(t) is the
time-dependent boundary data at the west boundary as indicated by the subscript. The other operators and boundary data are defined analogously. The parameters σ1,...,4
has to be determined such that the scheme (34) is stable in the ¯P -norm. The choice of penalty parameters are given in
Proposition 3.1. The scheme (34) is stable by choosing
σ1 ≤ − ˜ αs 4p1 , σ2 ≤ − ˜ αs 4p2 , σ3 ≤ − ˜ αs 4p3 , σ4 ≤ − ˜ αs 4p4 , (35)
where p1,...,4 are small enough such that Px(1, 1) − p1 ≥ 0, Px(M, M ) − p2 ≥ 0,
Proof. To simplify the algebra we consider only the north boundary. All other boundaries are treated in the same way independently of each other. We apply the discrete energy method by multiplying with TTP and using the SBP properties of¯ the operators. By assuming that gx,yN = 0 we obtain
||T||2t = 2 ˜αsTTP¯xE¯x,yND¯yT + 2σ4T T ¯ PxE¯x,yNT − 2 ˜αs || ¯DxT|| 2 + || ¯DyT||2 (36) and the penalty parameter σ4 has to be chosen such that ||T||2t ≤ 0. We rewrite
|| ¯DyT||2 = p4( ¯DyT)TP¯xE¯x,yN( ¯DyT) + ( ¯DyT) ˜P ( ¯DyT) (37)
where ˜P = ¯P − p4P¯xE¯x,yN. By using (37) we can write (36) as the quadratic form
||T||2t + 2 ˜αs || ¯DxT||2P˜ + || ¯DyT||2 = (Rξ)T (Px⊗ M4) Rξ (38)
where R is a permutation matrix, ξ = [T, ¯DyT]T, Px = diag[ ¯PxE¯x,yN, ¯PxE¯x,yN] and
M4 = 2σ4 α˜s ˜ αs −2 ˜αsp4 . (39)
The eigenvalues λ1,2 of M4 are non-postive if (35) holds.
Note that σ1,2 and σ3,4 are proportional to ∆x1 and ∆y1 , respectively. Alternative
formulations which avoids the dependence of mesh size can be found in [30, 31, 32]. 3.2. The Navier-Stokes equations
The compressible Navier-Stokes equations in two space dimensions requires three boundary conditions at a solid wall [22]. Since we are aiming for modelling heat transfer in a solid using (1), we must let both the tangential and normal velocities be zero. As the the third boundary condition we specify the temperature along the boundary [22, 10, 33].
In order to analyze (1) in terms of well-posedness and stability, we use the non-dimensional, linear, constant coefficient and symmetric formulation
wt+ (A1w)x+ (A2w)y = (A11wx+ A12wy)x+ (A21wx+ A22wy)y (40) where = M a
Re, Re is the Reynolds number and M a is the Mach number. The
coefficient matrices can be found in [22, 34]. The symmetrized variables are
w = " ¯ c √ γ ¯ρρ, u, v, 1 ¯ cpγ(γ − 1)T #T , (41)
where an overbar denotes the constant state which we have linearized around. More details can be found in [22, 9, 10, 34, 35, 8].
By applying the energy method to (40) and using the boundary conditions
u = 0, v = 0, T = Tw (42)
we get the energy estimate ||w||2
t ≤ 0 where ||w||2 = R wTwdΩ. Details of the
We consider again the unit square discretized with M + 1 and N + 1 grid points, respectively. For simplicity we consider the south boundary only. An SBP-SAT discretization of (1) non-dimensionalized using (8)-(11) can be written as
qt+ ( ¯Dx⊗ I4)F + ( ¯Dy ⊗ I4)G = S, (43)
where q is the discrete grid function approximating q and F, G are the discrete flux vectors corresponding to (4). A SAT S which imposes the boundary conditions (42) at the south boundary is given by
S = (P¯y−1E¯x,y0 ⊗ I4) ¯Σ (q − 0) + εσ2( ¯P −1 y E¯x,y0 ⊗ I4) ¯H2q − 0 + εσ3( ¯Py−1E¯x,y0 ⊗ I4) ¯H3q − 0 + ε( ¯P −1 y E¯x,y0 ⊗ I4) ¯Θ ( ¯Dx⊗ I4)q − 0 + εσ4( ¯Py−1E¯x,y0 ⊗ I4) ¯H4q − gw , (44)
where ¯Hi = Ix ⊗ Iy ⊗ Hi and Hi are 4 × 4 matrices that have the only non-zero
element 1 at the (i, i) position on the diagonal. The bold fonted 0 denotes a zero vector of the same size as the solution vector to indicate that the boundary data is zero. The data vector gw can be defined as
gw=
1
γ(γ − 1)(ρTw⊗ e4), (45)
where Tw is the prescribed wall temperature and e4 = [0, 0, 0, 1]T. The penalty
matrices ¯Σ and ¯θ and penalty coefficients σ2,3,4 are given in
Proposition 3.2. The scheme (43) for the compressible Navier-Stokes equations, imposing the boundary conditions (42), given by (44) is stable using the penalty matrices ¯ Σ = Ix⊗ Iy⊗ Σ (46) with Σ = σ 2γ 1 0 √γ pγ(γ − 1) 0 0 0 0 √ γ 0 γ pγ(γ − 1) pγ(γ − 1) 0 pγ(γ − 1) γ − 1 , σ ≤ −c¯ 2, (47) ¯ Θ = Ix⊗ Iy⊗ Θ (48) where Θ = 0 0 0 0 0 0 λ+µ2 ¯ρ 0 0 0 0 0 0 0 0 0 (49)
and the penalty coefficients σ2 ≤ − µ3 p3(λ + 3µ)(µ − λ) ¯ρ , σ3 ≤ − 1 4p3 λ + 2µ ¯ ρ , σ4 ≤ − 1 4p3 γµ P r ¯ρ. (50) Proof. See [33] for a proof and further details.
With the schemes (34) and (43) being energy stable we can solve the two equa-tions simultaneously and study the difference due to compressibility effects and possible non-zero velocities in the Navier-Stokes equations.
0 2 4 6 8 10 0 0.02 0.04 0.06 0.08 0.1 Time l2ïdifference (a) l2-difference 0 2 4 6 8 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Maximum difference Time (b) l∞-difference
Figure 1: Difference of temperature in time between the Navier-Stokes equations and the heat equation with constant boundary data
3.3. Numerical comparisons of the temperature distributions
Since we consider the compressible Navier-Stokes equations, relation (21) does not hold in general since neither the thermal conductivity of the fluid nor the density is constant. To study the deviation we consider the extremely different materials air and gold at near room temperature. The equations in the solid domain will be non-dimensionalized using reference values from the fluid. The material properties were obtained from [36] and can be seen in in Table 1.
Table 1: Properties of air and gold at near room temperature in SI-units
Material Density Specific heat capacity Thermal conductivity
Air ρ∗ = 1.205 c∗P = 1005 κ∗ = 0.0257
Gold ρ∗s = 1.93 × 104 c∗
s = 130 κ
∗
s = 300.0
We will numerically solve the Navier-Stokes equations and the heat equation simultaneously. The coefficients of the Navier-Stokes equations have been chosen such that relation (21) holds initially. The solutions will deviate over time due to compressibility effects and we compute the difference as a function of time.
Two cases will be considered, both using 65 × 65 grid points. In the first case we construct the initial data for the temperature using the function
f (x, y, t) = 1
2 sin(x
2+ y2− t) + cos(x2+ y2− t) + 1 (51)
and we let
ρ(x, y, 0) = 1, u(x, y, 0) = 0, v(x, y, 0) = 0 (52) for the Navier-Stokes equations. The boundary data for the temperature is set equal to one along the entire boundary, while the velocities are set to zero for the Navier-Stokes equations. The temperature will reach the same steady-state solution and we plot the l2- and l∞-difference as a function of time in Figure 1.
As an even more severe test, we use the function (51) to create the initial- and time-dependent boundary data for the temperature. The other initial data for the
0 20 40 60 80 100 0 1 2 3 4 5x 10 ï3 l2ïdifference Time (a) l2-difference 0 20 40 60 80 100 0 0.002 0.004 0.006 0.008 0.01 0.012 Maximum difference Time (b) l∞-difference
Figure 2: Difference of temperature in time between the Navier-Stokes equations and the heat equation with time-dependent boundary data
Navier-Stokes equations are given as before by (52). The result can be seen in Figures 2.
We can see from Figures 3-4 that the temperature distribution from the Navier-Stokes equations and the heat equation are not identical in the time-dependent case. This is because the relations (21) do not hold in general. In particular, we see in Figure 5 that one of the reasons is the non-constant density. The variations are, however, very small.
A mesh refinement study shows that the difference does not converge to zero with the mesh size in either case. This is expected since we solve two similar, but different, equations. However, we must make sure that the difference remains small. Note that since we are using weak boundary conditions, the boundary values are not identical. The difference in boundary values does, however, converge to zero with the mesh size.
0 0.2 0.4 0.6 0.8 1 0 0.5 1 0.5 1 1.5 2 y Temperature x (a) Time t = 20 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0.8 1 1.2 1.4 1.6 1.8 y Temperature x (b) Time t = 100
Figure 3: Temperature from the Navier-Stokes equations with time-dependent boundary data
In the first case shown in Figure 1, the large difference in temperature at the beginning of the computation is due to the fact that the initial- and boundary data
0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.5 1 1.5 2 x 10ï3 y Difference between NavierïStokes and Heat equation
x (a) Time t = 20 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.5 1 1.5 x 10ï3 y Difference between NavierïStokes and Heat equation
x
(b) Time t = 100
Figure 4: Difference in temperature between the Navier-Stokes equations and the heat equation with time-dependent boundary data
0 0.2 0.4 0.6 0.8 1 0 0.5 1 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005 y Density x (a) t = 20 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0.9998 0.9999 1 1.0001 1.0002 1.0003 y Density x (b) t = 100 Figure 5: Density from the Navier-Stokes equations
does not match. This large initial discontinuity induces velocities which decreases with time. There is no significant difference in the computational time it takes to reach the steady-state solution.
In the second case shown in Figure 2, the initial- and boundary data match, and hence there are no significant velocities induced. The maximum deviation is approximately 0.2% while on average it is approximately 0.08%, and the difference does not increase with time.
These two examples show that the temperature distribution obtained by the Navier-Stokes equations does not differ a lot from that obtained by the heat equation. Initial velocities due to discontinuities in the initial-/boundary data decrease fast and the compressibility of the Navier-Stokes equations only has a small (zero in the steady-state case) effect on the temperature.
x y -5 0 5 10 -10 -5 0 5 10 x y 0 0.2 0.4 0.6 0.8 1 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4
Figure 6: Typical geometry of a CFD computation. Only a small part of the computational domain is solid.
4. Two different coupling techniques
In [4], a linear, symmetric model of the Navier-Stokes equations was coupled with the heat equation in 1D. The coupling was proven well-posed, stable and high-order accurate. The drawback of the coupling technique used in [4], and also in [1, 3, 2, 17], is that the discretization required two separate solvers - one for the Navier-Stokes equations and one for the heat equation. Even a third code might be needed which handles the communication between the two solvers [17]. This complicates the programming and implementation into existing multi-block Navier-Stokes solvers.
Here we aim at a different approach. Instead of coupling the Navier-Stokes equations to the heat equation, we will let the Navier-Stokes equations themselves govern the heat transfer in the solid as described in the previous sections. In this way, the programming effort is restricted to modifying the interface conditions of an already existing multi-block Navier-Stokes solver.
The complexity of solving the Navier-Stokes equations compared to the heat equation is of course significantly larger. However, a common scenario in a compu-tational fluid dynamics (CFD) computation is that only a small part of the com-putational domain is solid, for example when computing the flow field around an airfoil, see Figure 6.
The verification of the new procedure will be done by comparing it to the stan-dard method of coupling the Navier-Stokes equations to the heat equation. We shall denote the Navier-Stokes equations coupled with the heat equation by NS-HT and the coupled Navier-Stokes equations by NS-NS.
4.1. Coupling of the Navier-Stokes equations (NS-NS)
We consider the Navier-Stokes equations in the two domains Ω1 = [0, 1] × [0, 1]
and Ω2 = [0, 1] × [−1, 0] with an interface at y = 0. Denote the solution in Ω1 by
q = [ρ, ρu, ρv, e] and in Ω2 by ˜q = [ ˜ρ, ˜ρ˜u, ˜ρ˜v, ˜e].
The interface will be considered as a solid wall and hence we impose the general interface conditions for the velocities
φ1u − ψ1uy = 0, v = 0,
φ2u + ψ˜ 2u˜y = 0, ˜v = 0,
(53)
where we require
φ1ψ1 ≥ 0, φ2ψ2 ≥ 0. (54)
The conditions (53) includes the no-slip as well as slip-flow solid wall conditions. These conditions are thoroughly discussed in [33] and we will not go into details here. In this paper we use φ1 = φ2 = 1 and ψ1 = ψ2 = 0.
4.1.1. Well-posed coupling conditions for the temperature
To couple the temperature of the two equations we will use continuity of tem-perature and heat fluxes according to Fourier’s law,
T = ˜T , κ1Ty = κ2T˜y. (55)
The well-posedness of the coupled Navier-Stokes equations using the interface con-ditions (53) and (55) will be shown by the energy method.
Consider the linearized, frozen coefficient and symmetric Navier-Stokes equations
wt+ (A1w)x+ (A2w)y = ε (A11wx+ A12wy)y+ (A21wx+ A22wy)y , ˜ wt+ (B1w)˜ x+ (B2w)˜ y = ε (B11w˜x+ B12w˜y)y+ (B21w˜x+ B22w˜y)y . (56)
The energy estimates of w and ˜w will be derived in the L2-equivalent norms
||w||2 H1 = Z Ω1 wTH1wdΩ, || ˜w||2H2 = Z Ω2 ˜ wTH2wdΩ˜ (57) where H1,2 = diag[1, 1, 1, δ1,2], δ1,2 > 0 (58)
are to be determined. We apply the energy method and consider only the terms at the interface y = 0. We get by using the conditions (54) that
d dt ||w|| 2 H1 + || ˜w|| 2 H2 ≤ −2ε 1 Z 0 δ1µ¯1 ¯ ρ1¯c21(γ1− 1)P r1 T Ty− δ2µ¯2 ¯ ρ2c¯22(γ2− 1)P r2 ˜ T ˜Ty dx, (59) where the bar denotes the state around which we have linearized and the subscript 1 or 2 refer to values from the corresponding subdomain Ω1 or Ω2. We shall choose
to linearize around the same speed of sound, that is ¯c1 = ¯c2 and let γ1 = γ2. By
requiring continuity of temperature (T = ˜T ) equation (59) reduces to
d dt ||w|| 2 H1 + || ˜w|| 2 H2 ≤ −2ε ¯ c2 1(γ1− 1) 1 Z 0 T δ1κ¯1 ¯ ρ1cP1 Ty− δ2κ¯2 ¯ ρ2cP2 ˜ Ty dx. (60)
In order to obtain an energy estimate by using continuity of the heat fluxes, we need to choose the weighs
δ1 = ¯ρ1cP1, δ2 = ¯ρ2cP2 (61) since then d dt ||w|| 2 H1 + || ˜w|| 2 H2 ≤ −2ε ¯ c2 1(γ1− 1) 1 Z 0 Tκ¯1Ty − ¯κ2T˜y dx = 0. (62)
Hence the interface conditions (55) gives an energy estimate and no unbounded energy growth can occur.
Note that the physical interface conditions (55) requires an estimate in a different norm than the standard L2-norm. From a mathematical point of view, any interface
condition which give weights that define a norm will result in a well-posed coupling. 4.1.2. The discrete problem and stability
In [14], a stable and conservative multi-block coupling of the Navier-Stokes equa-tions were shown. The coupling was done by considering continuity of all quantities and of the fluxes. In our case, the velocities are uncoupled and the equations are coupled only by continuity of temperature and heat fluxes.
We consider again the formulation (56) and discretize using SBP-SAT for im-posing the interface conditions (53) and (55) weakly. We let for simplicity the subdomains be discretized by equally many uniformly distributed gridpoints which allow us to use the same difference operators in both subdomains. We stress that the subdomains can have different discretizations as can be seen from i.e. [14, 37]. This assumption simplifies the notation and avoids the use of too many subscripts.
We discretize (56) using the SBP-SAT technique as wt+ ˆDxF + ˆDyG = S,
˜
wt+ ˆDxF + ˆ˜ DyG = ˜˜ S,
(63)
where the discrete fluxes are given by
F = ˆA1w − ε ˆA11wx+ ˆA12wy , G = ˆA2w − ε ˆA21wx+ ˆA22wy , ˜ F = ˆB1w − ε˜ ˆB11w˜x+ ˆB12w˜y , ˜ G = ˆB2w − ε˜ ˆB21w˜x+ ˆB22w˜y . (64)
The hat notation denotes that the matrix has been extended to the entire system as ˆ Dx = Dx⊗ Iy⊗ I4, Dˆy = Ix⊗ Dy ⊗ I4, ˆ Aξ= Ix⊗ Iy ⊗ Aξ, Bˆξ = Ix⊗ Iy⊗ Bξ, (65)
where ξ is a generic index ranging over the indicies which occur in (64).
The SATs imposing the interface conditions (53) and (55) can be written as
S =Pˆy−1Eˆx,yNΣˆ1 w − g I + εσ 2Pˆy−1Eˆx,yN φ1Hˆ2w − ψ1Hˆ2Dˆyw − g1 + εσ3Pˆy−1Eˆx,yN ˆH3w − g2 + ε ˆPy−1Eˆx,yNΘˆ1 ˆ H3Dˆxw − ∂g2 ∂x + εσ4Pˆy−1Eˆx,yN ˆI T 1w − ˆI T 2w˜ + εσ5Pˆy−1Dˆ T yEˆx,yN ˆI T 1w − ˆI T 2w˜ + εσ6Pˆy−1Eˆx,yN ¯ κ1Iˆ1TDˆyw − ¯κ2Iˆ2TDˆyw˜ (66) and ˜ S =Pˆy−1Eˆx,y0Σˆ2 w − ˜˜ g I + ε˜σ 2Pˆy−1Eˆx,y0 φ2Hˆ2w − ψ˜ 2Hˆ2Dˆyw − ˜˜ g1 + ε˜σ3Pˆy−1Eˆx,y0 ˆH3w − ˜˜ g2 + ε ˆPy−1Eˆx,y0Θˆ2 ˆ H3Dˆxw −˜ ∂ ˜g2 ∂x + ε˜σ4Pˆy−1Eˆx,y0 ˆI T 2w − ˆ˜ I T 1w + ε˜σ5Pˆy−1Dˆ T yEˆx,y0 ˆI T 2w − ˆ˜ I T 1w + ε˜σ6Pˆy−1Eˆx,y0 ¯ κ2Iˆ2TDˆyw − ¯˜ κ1Iˆ1TDˆyw . (67)
Here ˆP = ¯P ⊗ I4, ˆEx,y0 = ¯Ex,y0⊗ I4, ˆHj = Ix⊗ Iy⊗ Hj and Hj is a 4 × 4 matrix with
the only non-zero element 1 at the (j, j)th position on the diagonal and the operators ˆ
I1,2 selects the interface elements. The penalty matrices ˆΣ1,2 = Ix⊗ Iy⊗ Σ1,2, ˆΘ1,2 =
Ix⊗ Iy ⊗ Θ1,2, and the penalty coefficients σ2,...,6 and ˜σ2,...,6 has to be determined
such that the scheme is stable.
Remember that in the energy estimates for the continuous coupling, a non-standard L2-equivalent norm was used. The same modification to the norms has to be done in the discrete case. Thus, the discrete energy estimates will be derived in the norms ||w||2 H1 = w TP ˆˆH 1w, || ˜w||2H2 = ˜w TP ˆˆH 2w,˜ (68) where, ˆ H1 = Ix⊗ Iy⊗ H1, Hˆ2 = Ix⊗ Iy⊗ H2, (69)
and the matrices H1,2 are defined in (58) with the weights given in (61). Note that
ˆ
P ˆH1,2 = ˆH1,2P .ˆ
By applying the energy method to (63) and adding up we get d dt||w|| 2 H1+ d dt|| ˜w|| 2 H2 + DI = IT (70)
where the dissipation term, DI, is given by DI = 2ε ˆ Dxw ˆ Dyw T ˆ P ˆH1 0 0 P ˆˆH1 ˆ A11 Aˆ12 ˆ A21 Aˆ22 ˆ Dxw ˆ Dyw + 2ε ˆ Dxw˜ ˆ Dyw˜ T ˆ P ˆH2 0 0 P ˆˆH2 ˆ B11 Bˆ12 ˆ B21 Bˆ22 ˆ Dxw˜ ˆ Dyw˜ . (71)
The interface terms can be split into three parts as IT = IT1+ IT2+ IT3 where IT1
are the inviscid terms, IT2 the velocity terms and IT3 the coupling terms related to
the temperature.
In [33] it is shown how to choose Σ1,2, Θ1,2, σ2,3and ˜σ2,3, with small modifications,
such that the inviscid and velocity terms are bounded. Here we shall focus on the coupling terms. With appropriate choices of Σ1,2, Θ1,2, σ2,3 and ˜σ2,3 as described in
[33] we get d dt||w|| 2 H1 + d dt|| ˜w|| 2 H2 + DI ≤ IT3, (72)
where IT3 can be written as the quadratic form
IT3 = −ε(Rξ)T(Px⊗ M ) Rξ. (73)
To obtain (73), we have used the permutation similarity property of the Kro-necker product, R is a permutation matrix and ξ = [Ti, ˜Ti, (DyT)i, (DyT)˜ i]T where
the subscript i denotes the values at the interface. Note that we do not need the specific form of R, it is sufficient to know that such a matrix exists. Furthermore, we have Px = diag[δ1Px, δ2Px, δ1Px, δ2Px], (74) and M = −2σ4 σ4+ ˜σ4 κ¯1γ1− σ5− ¯κ1σ6 κ¯2σ6+ ˜σ5 σ4+ ˜σ4 −2˜σ4 σ5+ ¯κ1σ˜6 −¯κ2γ1− ˜σ5− ¯κ2σ˜6 ¯ κ1γ1− σ5− ¯κ1σ6 σ5+ ¯κ1σ˜6 0 0 ¯ κ2σ6+ ˜σ5 −¯κ2γ1− ˜σ5− ¯κ2σ˜6 0 0 . (75) Since Px is positive definite and the Kronecker product preserves positive
definite-ness, the necessary requirement for (72) to be bounded is that the penalty coefficients are chosen such that M ≥ 0. The penalty coefficients are given in
Proposition 4.1. The coupling terms IT3 in (72) are bounded using
˜
σ4 = σ4 ≤ 0, σ5 = −¯κ1r, σ6 = γ + r, ˜σ5 = −¯κ2(γ1+ r), ˜σ6 = r, r ∈ R (76)
and hence the scheme (63) is stable.
Proof. With the choices of penalty coefficients given in Proposition 4.1, the matrix M in (75) reduces to M = 2σ4 −1 1 0 0 1 −1 0 0 0 0 0 0 0 0 0 0 (77)
4.2. Coupling of the Navier-Stokes equations with the heat equation (NS-HT) Consider again the heat equation (15) and the Navier-Stokes equations in the constant, linear, symmetric formulation (40). We will again use the interface condi-tions
T = ˜T , κTy = κsT˜y (78)
and show that they lead to a bounded energy growth. The estimates of w and ˜T will again be derived in the L2-equivalent norms
||w||2 J1 = Z Ω1 wTJ1wdΩ1, || ˜T ||2ν2 = Z Ω2 ˜ T2ν2dΩ2 (79)
where J1 = diag[1, 1, 1, ν1] and ν1,2 > 0 are to be determined.
Remember that the symmetrized variables for the Navier-Stokes equations are
w = " ¯ c √ γ ¯ρρ, u, v, 1 ¯ cpγ(γ − 1)T #T . (80)
and note that there is a scaling coefficient in the temperature component. To simplify the analysis, we rescale (15) by multiplying the equation with 1
¯
c√γ(γ−1. To apply
the energy method, we rewrite the speed of sound based P´eclet number P ecin (14)
as P ec= P r · Re M a = P r ε (81)
where P r = c∗P ∞µ∗∞/κ∗∞ is the Prandtl number and thus (15) becomes ˜ Tt ¯ cpγ(γ − 1) = εκs P r¯cpγ(γ − 1)ρscs ˜Txx+ ˜Tyy . (82)
By applying the energy method to (82) and (40), and adding the results we obtain d dt ||w||2 J1 + 1 ¯ c2γ(γ − 1)|| ˜T || 2 ν2 ≤ −2ε ¯ c2γ(γ − 1)P r 1 Z 0 ν1γµ ¯ ρ T Ty − ν2κs ρscs ˜ T ˜Ty dx. (83) If we choose ν1 = ¯ κ ¯ρ γµ, ν2 = ρscs (84)
and apply the interface conditions (78) we get
d dt ||w||2 J1 + 1 ¯ c2γ(γ − 1)|| ˜T || 2 ν2 ≤ −2ε ¯ c2γ(γ − 1)P r 1 Z 0 TκT¯ y − κsT˜y dx = 0 (85)
and hence the conditions (78) does not contribute to unbounded energy growth. Note again that the application of the physical interface conditions (78) requires the use of a non-standard norm in the energy estimates.
4.2.1. The discrete problem and stability
The discretization of the coupled system is analogous to that which is presented in [4], and extended to multiple dimensions as described before. We hence only present the numerical scheme and the choice of interface penalty coefficients such that the scheme is stable.
An SBP-SAT discretization of the Navier-Stokes equations coupled with the heat equation is given by, when only considering the interface terms,
wt+ D¯x⊗ I4 F + ¯Dy⊗ I4 G = S,
˜
Tt− ¯D2xT + ¯˜ D2yT˜
= ˜S. (86)
The penalty terms are given by S = P¯y−1E¯x,yN ⊗ ¯Σ1 w − gI + εσ2 P¯y−1E¯x,yN ⊗ I4 φ1H¯2w − ψ1H¯2 D¯y⊗ I4 w − g1 + εσ3 P¯y−1E¯x,yN ⊗ I4 ¯ H3w − g2 + ε ¯Py−1E¯x,yN ⊗ I4 ¯ Θ1 ¯ H3 D¯x⊗ I4 w − ∂g2 ∂x + ε ¯Py−1E¯x,yN ⊗ Σ4 ¯ I1Tw − ¯I2T( ˜T ⊗ e4) + ε ¯Py−1D¯TyE¯x,yN ⊗ Σ5 ¯ I1Tw − ¯I2T( ˜T ⊗ e4) + ε ¯Py−1E¯x,yN ⊗ Σ6 ¯ κ ¯I1T D¯y⊗ I4 w − κsI¯2T( ¯DyT ⊗ e˜ 4) , (87)
where Σ4,5,6= diag[0, 0, 0, σ4,5,6] and
˜ S = ετ4P¯y−1E¯x,yN ˜T − T + ετ5P¯y−1D¯ T yE¯x,yN ˜T − T + ετ6P¯y−1E¯x,yN κsD¯yT − ¯˜ κ ¯DyT . (88)
The choices of penalty parameters such that the scheme is stable is given in
Proposition 4.2. The scheme (86) for coupling the Navier-Stokes equations with the heat equation is stable with the SATs given by (87), (88) where the penalty coefficients for the coupling terms are given by
r ∈ R, σ4 = τ4 ≤ 0, σ5 = −κsr, σ6 = −1 + rP r P r , τ5 = − ¯ κ (−1 + rP r) P r , τ6 = r. (89)
Proof. We apply the energy method, using the modified discrete norms, ||w||2 J1 = w T( ¯P ⊗ J 1)w, || ˜T||2ν2 = ν2 ˜ TTP T,¯ (90)
where J1 = diag[1, 1, 1, ν1] and ν1,2 are given in (84). Using appropriate penalty
terms for the inviscid part and the velocity components of the Navier-Stokes equa-tion, see [10, 33], we obtain the energy estimate
d dt||w|| 2 J1 + d dt|| ˜T|| 2 ν2 ≤ 0 (91)
Table 2: Difference between NS-NS and NS-HT at time t = 500
Difference
N l∞ l2 Interface
32 1.1514e-03 6.8992e-04 1.1514e-03 64 2.4612e-04 1.4491e-04 2.4612e-04 128 4.3440e-05 2.5329e-05 4.3440e-05
5. Numerical results
In this section we will quantify the difference between the two methods, NS-NS and NS-HT, in two representative cases. The computational domain is Ω = Ω1∪ Ω2
where Ω1 = [0, 1] × [0, 1] and Ω2 = [0, 1] × [−1, 0]. All computations are done using
3rd-order accurate SBP operators and the time integration is done using the classical 4th-order Runge-Kutta method.
The computations are initialized with zero velocities everywhere and temperature T = 1 in both subdomains. In the x-direction we have chosen periodic boundary conditions. At y = −1 we specify T = 1.5 and at y = 1 we have T = 1. For the Navier-Stokes equations we have no-slip solid walls as described in [33] for the veloc-ities. These choices of boundary conditions renders the solution to be homogeneous in the x-direction.
Under the assumption of identically zero velocities and periodicity in the x-direction, the exact steady-state solution can be obtained as
T = − k 2(k + 1)y + 3k + 2 2(k + 1), ˜ T = − 1 2(k + 1)y + 3k + 2 2(k + 1), (92)
where k is the ratio of the steady-state thermal conductivities. We can see from (92) that the only occurring material parameter is the ratio between the thermal conductivity coefficients. Neither the density nor the thermal diffusivity has any effect on the steady-state solution. The larger the ratio of the thermal conductivities is, the stiffer the problem becomes. Hence to simplify the calculations we have chosen two materials such that k = 5 (rather than air and gold for which k ≈ 11600).
The temperature distribution at time t = 500, which is the steady-state solution, is seen in Figure 7 when using 65 grid points in each coordinate direction and subdomain. In Figure 8 we show an intersection of the absolute difference along the line −1 ≤ y ≤ 1 at x = 0.5 together with the time-evolution of the l∞- and
l2-differences. In Figure 8(b) we can see that the large initial discontinuity gives
differences in the beginning of the computation. As the velocities are damped over time, the difference decreases rapidly towards zero.
In Table 2 we list the results for different number of grid points. As we can see, the differences are very small. Even for the coarsest mesh, the relative maximum and interface differences are less than 0.1% while the relative l2-difference is about
0.05%. The differences will not decrease asymptotically to zero with the design order of accuracy of the scheme since we compare two different systems of equations.
0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x y NS−NS, time t = 500.00, Nx=64, Ny=64 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5
(a) Temperature distribution from NS-NS −11 −0.5 0 0.5 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 y Temperature t=500.00, N=64 NS−NS NS−H
(b) Temperature distribution for NS-NS and NS-HT along y at x = 0.5.
Figure 7: Temperatures at time t = 500 from NS-NS and NS-HT using 65 grid points in each coordinate direction and subdomain
Next we consider an unsteady problem. The boundary data at the south bound-ary is given by the time-dependent perturbation
f (x, t) = 1 + 0.25 ∗ sin(t) ∗ sin(πx) (93) and hence there will be no steady-state solution. In the x-direction in the solid domain, we have changed from periodic boundary conditions to solid wall boundary conditions with prescribed temperature T = 1. The no-slip solid wall boundary condition has the additional benefit of damping the induced velocities in the Navier-Stokes equation. The results can be seen in Figure (9) where we plot the l∞- and
l2-difference as a function of time. As we can see, the difference does not approach
zero but it remains bounded and the relative mean difference is still less than 0.5% while the maximum difference is less than 1.5%. Thus, despite of the rather large variation in the boundary data, NS-NS and NS-HT produces very similar solutions. The choice of (93) was done to prevent missmatching data at the corners and the introduction of additional numerical errors.
As a final example, we consider the NS-NS coupling with a flow over a slab of material for which the ratio of the thermal conductivities is k = 100. The initial condition is temperature T = 1 in the fluid domain and T = 1.5 in the solid domain. The boundary conditions are periodic in the x-direction. At the south boundary in the solid domain we let T = 1.5 and at the north boundary in the fluid domain, there is a Mach 0.5 free-stream boundary condition with T = 1, as described in [9]. Figure 10 shows a snapshot of the solution at time t = 2.5. The time-step had to be decreased with a factor of four, compared to the previous cases, to compensate for the increased stiffness of the system.
−10 −0.5 0 0.5 1 1 2 x 10−4 y Difference t=500.00, N=64
(a) Temperature difference between NS-NS and NS-NS-HT along y at x = 0.5. 0 100 200 300 400 500 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 t Max Mean
(b) l∞- and l2-difference in time
Figure 8: Temperature differences for NS-NS and NS-HT using 65 grid points in each coordinate direction and subdomain
0 20 40 60 80 100 0 0.005 0.01 0.015 0.02 t Max Mean
Figure 9: l∞- and l2-difference in time between NS-HS and NS-HT for an unsteady problem
6. Conclusion
We have evaluated the compressible Navier-Stokes equations as a model for heat transfer, and in particular conjugate heat transfer problems. For the one domain case, we have derived conditions for how to choose the coefficients and scalings of the energy equation in the Navier-Stokes equations such that it is as similar to the heat equation as possible. Computations using material data for air and gold at near room temperature showed that the l2-difference was less than 0.08% and did
not increase in time.
A thorough investigation of the interface conditions leading to well-posedness of the continuous equations were performed. The physical interface conditions required that a non-standard L2-equivalent norm was used to derive the energy estimates.
All numerical schemes were proven energy stable and the stability does not de-pend on the order of accuracy. Higher order operators can thus be used without modifications to the scheme.
We have performed numerical simulations of a conjugate heat transfer problem using two different approaches. In the steady-state computation, both the
temper-−0.5 0 0.5 1 1.5 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x y
Temperature and velocity field
1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45
(a) Temperature distribution and velocity profile
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 y Temperature intersection Fluid Solid
(b) Intersection of the temperature dis-tribution
Figure 10: Temperature and velocity profiles for a flow past a slab of material using NS-NS with 65x65 grid points
ature difference at the interface and the l∞-difference were less than 0.1% and the
l2-difference was about 0.05%. In the unsteady case, the difference does not
ap-proach zero over time, but does not increase with time. The relative l2-difference
was less than 1.5% while the relative l2-difference was less than 0.5%.
Instead of writing separate solvers for the fluid and solid parts, an existing multi-block Navier-Stokes solver can be modified. In this way there will be no additional communication overhead between the different codes and the programming is re-stricted to modifying the interface conditions and scalings in the subdomains.
Thus, using the compressible Navier-Stokes equations for conjugate heat transfer problems is a viable and efficient solution if a multi-block Navier-Stokes solver exist for which the interface conditions can be modified.
Acknowledgments
The computations were performed on resources provided by SNIC through Up-psala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under Project p2010056.
References
[1] M. B. Giles. Stability analysis of numerical interface conditions in fluid-structure thermal analysis. International Journal for Numerical Methods in Fluids, 25:421–436, August 1997.
[2] B. Roe, R. Jaiman, A. Haselbacher, and P. H. Geubelle. Combined interface boundary condition method for coupled thermal simulations. International Journal for Numerical Methods in Fluids, 57:329–354, May 2008.
[3] William D. Henshaw and Kyle K. Chand. A composite grid solver for conjugate heat transfer in fluid-structure systems. Journal of Computational Physics, 228(10):3708 – 3741, 2009.
[4] Jens Lindstr¨om and Jan Nordstr¨om. A stable and high-order accurate conjugate heat transfer problem. Journal of Computational Physics, 229(14):5440–5456, 2010.
[5] Michael Sch¨afer and Ilka Teschauer. Numerical simulation of coupled fluid– solid problems. Computer Methods in Applied Mechanics and Engineering, 190(28):3645 – 3667, 2001.
[6] Niphon Wansophark, Atipong Malatip, and Pramote Dechaumphai. Stream-line upwind finite element method for conjugate heat transfer problems. Acta Mechanica Sinica, 21:436–443, 2005.
[7] Xi Chen and Peng Han. A note on the solution of conjugate heat transfer problems using simple-like algorithms. International Journal of Heat and Fluid Flow, 21(4):463 – 467, 2000.
[8] Jan Nordstr¨om and Mark H. Carpenter. Boundary and interface conditions for high-order finite-difference methods applied to the Euler and Navier-Stokes equations. Journal of Computational Physics, 148(2):621 – 645, 1999.
[9] Magnus Sv¨ard, Mark H. Carpenter, and Jan Nordstr¨om. A stable high-order finite difference scheme for the compressible Navier-Stokes equations, far-field boundary conditions. Journal of Computational Physics, 225(1):1020–1038, 2007.
[10] Magnus Sv¨ard and Jan Nordstr¨om. A stable high-order finite difference scheme for the compressible Navier-Stokes equations: No-slip wall boundary conditions. Journal of Computational Physics, 227(10):4805 – 4824, 2008.
[11] Magnus Sv¨ard, Ken Mattsson, and Jan Nordstr¨om. Steady-State Computa-tions Using Summation-by-Parts Operators. Journal of Scientific Computing, 24(1):79–95, 2005.
[12] Ken Mattsson, Magnus Sv¨ard, Mark Carpenter, and Jan Nordstr¨om. High-order accurate computations for unsteady aerodynamics. Computers and Flu-ids, 36(3):636 – 649, 2007.
[13] X. Huan, J.E. Hicken, and D.W. Zingg. Interface and Boundary Schemes for High-Order Methods. In the 39th AIAA Fluid Dynamics Conference, AIAA Paper No. 2009-3658, San Antonio, USA, 22–25 June 2009.
[14] Jan Nordstr¨om, Jing Gong, Edwin van der Weide, and Magnus Sv¨ard. A stable and conservative high order multi-block method for the compressible Navier-Stokes equations. Journal of Computational Physics, 228(24):9020–9035, 2009. [15] Ken Mattsson, Magnus Sv¨ard, and Mohammad Shoeybi. Stable and accurate schemes for the compressible Navier-Stokes equations. Journal of Computa-tional Physics, 227:2293–2316, February 2008.
[16] Jan Nordstr¨om, Sofia Eriksson, Craig Law, and Jing Gong. Shock and vortex calculations using a very high order accurate Euler and Navier-Stokes solver. Journal of Mechanics and MEMS, 1(1):19–26, 2009.
[17] Jan Nordstr¨om, Frank Ham, Mohammad Shoeybi, Edwin van der Weide, Mag-nus Sv¨ard, Ken Mattsson, Gianluca Iaccarino, and Jing Gong. A hybrid method for unsteady inviscid fluid flow. Computers & Fluids, 38:875–882, 2009.
[18] Heinz-Otto Kreiss and Godela Scherer. Finite element and finite difference methods for hyperbolic partial differential equations. In Mathematical Aspects of Finite Elements in Partial Differential Equations, number 33 in Publ. Math. Res. Center Univ. Wisconsin, pages 195–212. Academic Press, 1974.
[19] Heinz-Otto Kreiss and Godela Scherer. On the existence of energy estimates for difference approximations for hyperbolic systems. Technical report, Uppsala University, Division of Scientific Computing, 1977.
[20] Mark H. Carpenter, David Gottlieb, and Saul Abarbanel. Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: Methodol-ogy and application to high-order compact schemes. Journal of Computational Physics, 111(2):220 – 236, 1994.
[21] Mark H. Carpenter, Jan Nordstr¨om, and David Gottlieb. A Stable and Con-servative Interface Treatment of Arbitrary Spatial Accuracy. Journal of Com-putational Physics, 148(2):341 – 365, 1999.
[22] Jan Nordstr¨om and Magnus Sv¨ard. Well-Posed Boundary Conditions for the Navier–Stokes Equations. SIAM Journal on Numerical Analysis, 43(3):1231– 1255, 2005.
[23] Bo Strand. Summation by Parts for Finite Difference Approximations for d/dx. Journal of Computational Physics, 110(1):47 – 67, 1994.
[24] Ken Mattsson and Jan Nordstr¨om. Summation by parts operators for finite difference approximations of second derivatives. Journal of Computational Physics, 199(2):503–540, 2004.
[25] Ken Mattsson. Summation by parts operators for finite difference approxima-tions of second-derivatives with variable coefficients. Technical Report 2010-023, Uppsala University, Division of Scientific Computing, 2010.
[26] Magnus Sv¨ard and Jan Nordstr¨om. On the order of accuracy for difference approximations of initial-boundary value problems. Journal of Computational Physics, 218(1):333–352, 2006.
[27] Heinz-Otto Kreiss. Stability theory for difference approximations of mixed initial boundary value problems. I. Mathematics of Computation, 22(104):pp. 703–714, 1968.
[28] Bertil Gustafsson, Heinz-Otto Kreiss, and Arne Sundstr¨om. Stability theory of difference approximations for mixed initial boundary value problems. II. Math-ematics of Computation, 26(119):pp. 649–686, 1972.
[29] Bertil Gustafsson, Heinz-Otto Kreiss, and Joseph Oliger. Time Dependent Problems and Difference Methods. Wiley Interscience, 1995.
[30] Jens Lindstr¨om and Jan Nordstr¨om. Spectral analysis of the continuous and dis-cretized heat and advection equation on single and multiple domains. Technical Report 2010-030, Department of Information Technology, Uppsala University, December 2010.
[31] Jing Gong and Jan Nordstr¨om. Interface procedures for finite difference ap-proximations of the advection-diffusion equation. Journal of Computational and Applied Mathematics, 236(5):602 – 620, 2011.
[32] Mark H. Carpenter, Jan Nordstr¨om, and David Gottlieb. Revisiting and extend-ing interface penalties for multi-domain summation-by-parts operators. Journal of Scientific Computing, 45(1-3):118–150, 2010.
[33] Jens Berg and Jan Nordstr¨om. Stable Robin solid wall boundary conditions for the Navier-Stokes equations. Journal of Computational Physics, 230(19):7519 – 7532, 2011.
[34] Saul Abarbanel and David Gottlieb. Optimal time splitting for two- and three-dimensional Navier-Stokes equations with mixed derivatives. Journal of Com-putational Physics, 41(1):1 – 33, 1981.
[35] Jan Nordstr¨om. The use of characteristic boundary conditions for the Navier-Stokes equations. Computers & Fluids, 24(5):609 – 623, 1995.
[36] The engineering toolbox. http://www.engineeringtoolbox.com, May 2011. [37] Ken Mattsson and Mark H. Carpenter. Stable and accurate interpolation
op-erators for high-order multiblock finite difference methods. SIAM Journal on Scientific Computing, 32(4):2298–2320, 2010.
