Implicit heat flux correction-based immersed boundary-finite volume method for thermal flows with Neumann boundary conditions

Accepted Manuscript Implicit heat flux correction-based immersed boundary-finite volume method for thermal flows with Neumann boundary conditions

Tongqing Guo, Ennan Shen, Zhiliang Lu, Yan Wang, Lu Dong

PII: DOI: Reference:

S0021-9991(19)30133-0 https://doi.org/10.1016/j.jcp.2019.02.016 YJCPH 8515

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

24 January 2018 25 February 2019 25 February 2019

Please cite this article in press as: T. Guo et al., Implicit heat flux correction-based immersed boundary-finite volume method for thermal flows with Neumann boundary conditions, J. Comput. Phys. (2019), https://doi.org/10.1016/j.jcp.2019.02.016

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

Highlights • Within the scopes of immersed boundary method (IBM) and finite volume method (FVM), an implicit heat flux correction-based IB-FVM for thermal flows with Neumann boundary conditions is proposed by constructing hybrid thin-plate splines (TPS) interpolation/delta function distribution. In the derivative interpolation process, the cosine delta function is replaced by the much more accurate TPS, since the use of cosine delta function results in a less accurate solution. • Different from conventional IBM, the cell face centres are defined as the Eulerian points in the framework of FVM. As compared with existing implicit IBM for Neumann conditions, the present method avoids the introduction of auxiliary layers of Lagrangian points as well as the approximate conversion from the Neumann to Dirichlet condition, and thus is suitable for an arbitrary geometry. • The heat flux of a solid cell is corrected in such a way that the same amount of heat flux as that flowing into the fluid domain across the boundary is supplemented into the solid domain.

Implicit heat flux correction-based immersed boundary-finite volume method for thermal flows with Neumann boundary conditions Tongqing Guoa, Ennan Shena, Zhiliang Lua,*, Yan Wanga, Lu Dongb a

Key Laboratory of Unsteady Aerodynamics and Flow Control, Ministry of Industry and Information Technology, Nan-

jing University of Aeronautics and Astronautics, Yudao Street 29, Nanjing, Jiangsu 210016, China b

State Key Laboratory of Disaster Prevention and Mitigation of Explosion and Impact, Army Engineering University of

PLA, Nanjing, Jiangsu 210007, China

Abstract In the frameworks of immersed boundary method (IBM) and finite volume method (FVM), an implicit heat flux correction-based IB-FVM is proposed for thermal flows with Neumann boundary conditions. With the use of a fractional-step technique, the preconditioned Navier-Stokes (N-S) equations are solved by the FVM to obtain the intermediate solution in the prediction step and the heat flux is corrected by enforcing the Neumann condition in the correction step. Different from existing IBMs, the cell face centres are defined as the Eulerian points due to the heat flux computation at each face in the FVM. The Neumann condition is implemented in such a way that the interpolated temperature gradient is equal to the specified boundary value at the same point when the corrected gradient field is interpolated from the face centres to the Lagrangian points. To achieve an implicit algorithm, the temperature derivative corrections at the Lagrangian points are set as unknowns and a system of algebraic equations is established by constructing hybrid thin-plate splines (TPS) interpolation/delta function distribution. In the derivative interpolation process, the much more accurate TPS is introduced because the use of cosine delta function yields a less accurate solution. After the distribution process, the heat flux correction of a fluid cell is evaluated by using the

solved temperature derivative corrections at the face centres, but that of a solid cell is calculated by using their additive inverses to supplement the same amount of heat flux into the solid domain as that flowing into the fluid domain across the boundary. Finally, the heat flux of a cell is corrected by adding the correction to the intermediate value and the corrected heat flux is utilized to solve the NS equations in the prediction step. As compared with the available implicit IBMs for Neumann conditions, the present method avoids the introduction of auxiliary layers of Lagrangian points as well as the approximate conversion from the Neumann to Dirichlet condition and thus is suitable for an arbitrary geometry. The proposed method is verified by simulating several benchmark thermal flows with Neumann conditions, the natural convection in an annulus and the steady or unsteady forced convection over a stationary or oscillating cylinder. All the computed results agree well with the literature data. Key words: Immersed boundary method; Finite volume method; Thermal flow; Neumann boundary condition; preconditioned Navier-Stokes equations; Fractional-step technique.

1. Introduction Due to the decoupling between the solution of governing equations and the implementation of boundary conditions, the immersed boundary method (IBM), which can simulate flows with complex and moving boundaries on simple Cartesian meshes, is receiving more and more attention. Since the IBM was originally proposed by Peskin [1] in the 1970s, a variety of IBMs have been developed and improved to study different kinds of flow problems with desirable accuracy and efficiency. The key issue in the IBM is the computation of restoring force or, more generally, the im-

plementation of boundary conditions. Earlier approaches mainly include the penalty method [1,2], the direct forcing scheme [3-6] and the momentum exchange scheme [7-10]. The conventional IBM suffers from the problem of flow penetration into the solid body, since the non-slip condition is not enforced in the solution process. To remove this drawback, Kim et al. [4] introduced a mass source/sink into the computation so that the non-slip condition was well kept. Wu and Shu [11] proposed an implicit velocity correction-based IBM, in which the intermediate velocity was corrected implicitly to satisfy the no-slip condition. Similar studies on the boundary condition-enforced IBM can also be found in [12-15]. Most of IBMs were initially proposed for incompressible flows with no-slip conditions. By extending their ideas, IBMs for thermal flows with isothermal conditions have been developed rapidly, such as multi-direct heat source scheme [17], implicit temperature correctionbased IBM [18], as well as other methods and applications [19-22]. All the above mentioned IBMs are only applicable for flows with Dirichlet conditions. There is another important type of boundary condition named as Neumann condition. As compared with Dirichlet condition, the implementation of Neumann condition is much more challenging and only several explicit [23-25] and implicit [26] IBMs have been developed until now. As for the IBM, there are two basic ways to impose the Neumann condition. One way is to convert the Neumann to Dirichlet condition and then IBMs for Dirichlet conditions can be employed. Zhang et al. [23] proposed an explicit direct heating scheme. In their scheme, an auxiliary layer was defined outside the immersed boundary and the Neumann condition was approximated by a finite difference scheme. Consequently, the Neumann condition was converted to the wall temperature condition. Wang et al. [26] present-

ed an implicit IBM by introducing two auxiliary layers inside and outside the immersed boundary to convert the Neumann condition to the Dirichlet condition on the inner auxiliary layer. This kind of method may not be suitable for geometries with sharp edges, such as an airfoil with a sharp trailing edge, since it is difficult to define such auxiliary layers for those configurations. The other way is to make a heat flux correction by imposing the Neumann condition. Ren et al. [24] established an explicit heat flux correction-based IBM by introducing a heat source/sink into the energy equation. In their work, the heat source/sink was determined from the difference between the specified heat flux and the interpolated value from the intermediate gradient field. The heat flux correction-based IBM avoids the introduction of auxiliary layers and the approximate conversion between two types of boundary conditions. So far, only explicit IBM for heat flux correction has been developed. Recently, Wang et al. [26] pointed out that most of the one-step explicit correctionbased IBMs do not accurately satisfy the Neumann condition. Moreover, all flow variable correction-based IBMs use the cell centres as the Eulerian points, which cannot be directly applied to the heat flux correction in the framework of FVM due to the requirement of heat flux computation at each cell face. The above problems motivate the present work. The Neumann thermal boundary condition directly contributes to the heat flux. Therefore, introducing a heat source/sink into the energy equation is equivalent to making a heat flux correction. In this work, a fractional-step technique will be used to decouple the solution of the N-S equations and the implementation of boundary conditions. In the prediction step, the N-S equations will be solved by the FVM to obtain the intermediate solution which neglects the effect of immersed boundary. In the correction step, the intermediate variables

will be corrected by enforcing all the physical boundary conditions. The no-slip and isothermal conditions will be imposed by employing existing implicit boundary conditionenforced IBMs for Dirichlet conditions. For the implementation of Neumann condition, all the face centres will be defined as the Eulerian points and an implicit heat flux correctionbased IB-FVM will be proposed via an interpolation process and a distribution process. After the temperature derivative correction at the face centres are solved implicitly, the heat flux correction for a mesh cell will be studied in detail. Finally, the present method will be verified by simulating several well-documented thermal flows with Neumann conditions.

2. Implicit heat flux correction-based IB-FVM 2.1 Governing equations In the Cartesian coordinate system, by introducing a source term QB into the compressible N-S equations to reflect the effect of immersed boundary, the governing equations describing a 2-D laminar thermal flow can be written as the following conservative form

∂W ∂fc1 ∂fc 2 ∂fv1 ∂fv 2 + + = + + QG + QB ∂t ∂x ∂y ∂x ∂y

(1)

together with the boundary conditions given by

−k

v = vB

on boundary Γ

(2a)

T = TB

on boundary Γ1

(2b)

∂T = q B on boundary Γ2 ∂n

(2c)

In Eq. (1), W is the vector of the conservative variables, fc1 and fc 2 are the convective flux vectors, fv1 and fv 2 are the viscous flux vectors. They are respectively defined

as

ª ρ º « ρu » « » W= « ρv » « » ¬ρ E¼

ª ρu º « 2 » ρu + p » « f c1 = « ρ uv » « » ¬ ρ uE + up ¼

ª0 º «τ » « xx » f v1 = «τ yx » « » «uτ + vτ + k ∂T » xx xy ∂x ¼» ¬«

ª ρv º « ρ vu » « » fc 2 = « ρv2 + p » « » ¬ ρ vE + vp ¼

ª0 º «τ » « xy » » f v 2 = «τ yy « » ∂T » « «uτ yx + vτ yy + k ∂y » ¬ ¼

(3)

For natural convection, the buoyancy force is the driven force and is usually calculated by using the Boussinesq approximation when the incompressible equations are solved. In the case of the compressible N-S equations, the buoyancy source term QG in Eq. (1) can be directly computed from the density variation QG = ª¬0, 0, − ( ρ − ρ ∞ ) g , − ( ρ − ρ ∞ ) gv º¼

T

(4)

In the above expressions, ρ , v , p , T and E are the fluid density, velocity vector with components u and v , pressure, temperature and total energy per unit mass, respectively. k is the thermal conductivity coefficient. τ ij represent the components of viscous stress tensor. n is the unit normal vector of the boundary. The subscripts B and ∞ refer to the immersed boundary and the free stream, respectively. q B stands for the prescribed heat flux on the boundary. g denotes the gravitational acceleration directed toward the negative y-direction. In the velocity correction-based IBM for no-slip conditions, Shu et al. [12] indicated that adding a force term in the momentum equations is equivalent to making a correction in the

velocity field. In their method, a fractional-step technique was introduced to decouple the solution of governing equations and the implementation of boundary conditions. Similarly, Eq. (1) is solved by the following prediction-correction scheme. Step 1 (Prediction step): Solve the N-S equations without the immersed boundary-related source term

∂W * ∂fc1 ∂fc 2 ∂fv1 ∂fv 2 + + = + + QG ∂t ∂x ∂y ∂x ∂y

(5)

and obtain the intermediate (predicted) solution W * . Step 2 (Correction step): Correct the intermediate solution W * by solving ∂W = QB ∂t

(6)

and get the correction δW . As a result, W n +1 = W * + δW , where W n +1 represents the physical (corrected) solution at the time level (n + 1) . In this work, Eq. (5) will be solved in the prediction step by using the FVM and the preconditioning method for low-speed flows. In the correction step, the velocity, temperature and heat flux will be corrected by enforcing the no-slip condition, the isothermal condition and the Neumann boundary condition, respectively. Both the no-slip condition and the isothermal condition will be performed by employing existing implicit IBMs for Dirichlet conditions. The Neumann thermal boundary condition will be implemented by proposing an implicit heat flux correction-based IB-FVM. 2.2 Prediction of intermediate solution by FVM It should be stated that, in the prediction step, all flow variables refer to the intermediate

variables marked with a superscript “*” in Eq. (5). For simplicity, the symbol of asterisk is omitted in this section. To simulate low-speed thermal flows accurately and efficiently, the preconditioning method is employed, which consists of the multiplication of the time-derivative item by a preconditioning matrix ī p . The preconditioned integral form of Eq. (5), formulated in primitive variables Wp = [ p, u, v, T ] , can be written as T

īp

∂ W p dΩ + ³ (Fc − Fv )dS = ³ QG dΩ ∂t Ω³ ∂Ω Ω

(7)

where Fc and Fv denote the integral forms of the convective flux vector and the viscous flux vector, respectively. To maintain the same form as that before preconditioning, the primitive variables in Eq. (7) are transformed to the conservative variables as follows ī p M −1

where M =

∂ WdΩ + ³ (Fc − Fv )dS = ³ QG dΩ ∂t Ω³ ∂Ω Ω

(8)

∂W . ∂W p

The preconditioning matrix proposed by Weiss and Smith [27] is adopted, which was found to perform very well for various flow cases. It is defined as 0 0 ρT · § Θ ¨ ¸ ρ 0 ρT u ¸ ¨ Θu īp = ¨ Θv 0 ρ ρT v ¸ ¨ ¸ ¨ ΘH − 1 ρu ρv ρc ¸ p © ¹

(9)

where H and c p are respectively the total enthalpy and the specific heat at constant

pressure,

ρT =

∂ρ ∂T

p = const

M p2 1 with a the local sound speed. ε p = . and Θ = ε pa2 1 + (γ − 1)M p2

γ is the ratio of specific heats and M p2 is determined from 2 § § (ν Δx ) ·¸ ,1¸· M p2 = Min ¨ Max ¨ M a2 , KM a2∞ , ¨ ¨ a 2 ¸¹ ¸ © © ¹

K = 1̚3

(10)

where M a and M a∞ are the local and free-stream Mach numbers, respectively, Δx is the characteristic length of the control volume and ν Δx denotes the diffusion velocity. The cell-centered FVM with central scheme is applied to discretize Eq. (8) on a Cartesian mesh, and the following system of ordinary differential equations in time can be obtained d ( Ω IWI ) = − Mī p−1 RI (W ) dt

(11)

where Ω is the volume, R represents the residual, and the index I denotes the particular control volume. For the steady flow simulation, Eq. (11) is advanced in time by using the 5-stage RungeKutta explicit scheme with a local time step to obtain a steady-state solution. In the case of unsteady computation, the dual time-stepping approach [37] is employed. Note that, the temporal accuracy will be damaged if the physical time derivative term is directly multiplied by the preconditioning matrix like Eq. (8). With the use of the dual timestepping, the preconditioning matrix will be exerted on the pseudo-time derivative term. At first, preconditioning is removed from Eq. (11) and the physical time derivative is approximated by an implicit backward difference formula of second-order accuracy 3Ω nI +1W In +1 − 4Ω nIW In + Ω nI −1W In −1 + RI (W n +1 ) = 0 2 Δt

(12)

where Δt denotes the global physical time step and the superscript n represents the time level.

A time-stepping methodology is used to solve Eq. (12), which can be written as

(

)

d Ω nI +1W I ª 3Ω n +1W I − 4Ω nIW In + Ω nI −1W In −1 º = −« I + RI W » = − R I W 2 Δt dt ¬ ¼

( )

( )

(13)

where W is the approximation to W n+1 , t denotes a pseudo-time variable and R stands for the unsteady residual. As mentioned above, the preconditioning matrix is exerted on the pseudo-time derivative term in Eq. (13) to obtain the following system similar to Eq. (11)

(

)

d Ω nI +1W I = − Mī p−1 R I W  dt

( )

(14)

The stationary solution of Eq. (14) corresponds to the flow variables at the new time level, i.e., W = W n+1 . Since R I = 0 at steady state in pseudo time, Eq. (12) is fulfilled. Hence, at each physical time step of the dual time-stepping, a steady problem governed by Eq. (14) is solved in the pseudo time by employing the same explicit time-marching scheme used in the steady flow simulation. 2.3 Corrections of velocity, temperature and heat flux by IBM 2.3.1 Velocity correction under no-slip condition As sketched in Fig. 1, the physical space and the immersed boundary are discretized by an Eulerian mesh and a set of Lagrangian points, respectively.

Eulerian mesh

Lagrangian point

Immersed boundary

Fig. 1. Immersed boundary illustration.

The velocity correction is made by applying the implicit correction-based IBM of Wu and Shu [11]. Suppose that the velocity correction is δv , then the corrected velocity is v = v * + δv . To satisfy the no-slip condition, it is required that when the corrected velocity

v is interpolated from the Eulerian points to the Lagrangian points, the interpolated velocity should be equal to the prescribed boundary value v B at the same position. In a diffuse interface IBM, the interpolation process from the Eulerian points to the Lagrangian points is conventionally done by using the following discrete delta function N

v B (X Bl ) = ¦ v (x j )D (x j − X Bl )h 2 j =1

(l = 1,2,", M )

l l 1 §¨ x j − X B ·¸ §¨ y j − YB ·¸ D (x j − X ) = 2 δ ¨ δ h © h ¸¹ ¨© h ¸¹ l B

(15)

(16)

where X B and x are respectively the coordinate vectors of a Lagrangian point and an Eulerian point, M is the number of Lagrangian points, N denotes the number of Eulerian points within the support of the delta function kernel at all Lagrangian points, and h is the mesh spacing around the immersed boundary. In Eq. (15), the corrected velocity v at an Eulerian point is defined as

v (x j ) = v * (x j ) + δv (x j )

(17)

In this work, a 4-point cosine function proposed by Peskin [28] is used, which reads

­1 ° (1 + cos(π r / 2 )) δ (r ) = ® 4 °0 ¯

r ≤2 r >2

(18)

On the other hand, the velocity correction δv (x j ) in Eq. (17) is evaluated via a distribution process, that is, the velocity corrections at the Lagrangian points are distributed back to the Eulerian points by using the following discrete delta function M

δv (x j ) = ¦ δv B (X Bi )D (x j − X Bi )Δs i i =1

( j = 1,2,", N )

(19)

where Δs i is the arc length of the i-th immersed boundary element. Substituting Eqs. (17) and (19) into Eq. (15) leads to N

N

M

j =1

j =1 i =1

v B (X Bl ) = ¦ v * (x j )D (x j − X Bl )h 2 + ¦¦ δvB (X Bi )D (x j − X Bi )Δs i D (x j − X Bl )h 2

(20)

By setting the term δvBi Δs i as an unknown, Eq. (20) can be further written as the following matrix form

AX v = B

(21)

where

[

X v = įv1B Δs1 , δvB2 Δs 2 ,", δvBM Δs M § D11 ¨ D 2¨ A = h ¨ 21 # ¨¨ © DM 1

D12 D22 # DM 2

]

T

D1N ·§ D11 ¸¨ ! D2 N ¸¨ D21 % # ¸¨ # ¸¨ " DMN ¸¹¨© DN 1 "

(22a) D12 D22 # DN 2

D1M · ¸ ! D2 M ¸ = h 2 Di Dd % # ¸ ¸ " DNM ¸¹ "

(22b)

§ v1B · § D11 ¨ 2¸ ¨ ¨ v B ¸ 2 ¨ D21 B=¨ ¸−h ¨ # # ¨ ¸ ¨¨ ¨vM ¸ © DM 1 © B ¹

D1N ·§ v1* · § v1B · § v1* · ¨ ¸ ¨ ¸ ¨ ¸ ¸ ! D2 N ¸¨ v2* ¸ ¨ vB2 ¸ 2 ¨ v2* ¸ = − h Di ¨ ¸ % # ¸¨ # ¸ ¨ # ¸ # ¨ ¸ ¸¸¨ * ¸ ¨ M ¸ * ¸ ¨ ¸ ¨ ¸ ¨ " DMN ¹© v N ¹ © v B ¹ © vN ¹

D12 D22 # DM 2

"

(22c)

where Di and Dd are the transformation matrices for the interpolation and distribution processes, respectively. Actually, Dd is the transpose of Di . After the unknown vector X v is solved from Eq. (21), the velocity field on the Eulerian mesh can be corrected by using Eqs. (19) and (17). 2.3.2 Temperature correction under Dirichlet thermal boundary condition The Dirichlet thermal boundary condition specifies a constant wall temperature. Both the no-slip condition and the constant wall temperature condition belong to the Dirichlet type, hence, the temperature correction under Dirichlet boundary condition is intrinsically consistent with the previous velocity correction. For the temperature correction, only the final equation system is given below. More details can be found in [15, 18]. AX T = B

(23)

where

[

X T = įTB1Δs1 , δTB2 Δs 2 ,", δTBM Δs M

[

B = TB1 , TB2 ,", TBM

]

T

[

]

T

− h 2 Di T1* , T2* ,", TN*

(24a)

]

T

(24b)

and the matrix A is the same as that in Eq. (21). The transformation matrix Di in Eq. (24) is for the interpolation process and has the same form as that in Eq. (22). Solve Eq. (23) and then distribute the temperature corrections at the Lagrangian points denoted by X T to the Eulerian points by using the delta function distribution

M

δT (x j ) = ¦ δTB (X Bi )D (x j − X Bi )Δs i i =1

( j = 1,2,", N )

(25)

As a result, the temperature filed on the Eulerian mesh is corrected from T (x j ) = T * (x j ) + δT (x j )

(26)

2.3.3 Heat flux correction under Neumann thermal boundary condition The Neumann thermal boundary condition prescribes a constant temperature gradient normal to the boundary and directly contributes to the heat flux. In this study, it is to be implemented strictly in the frameworks of IBM and FVM. While the FVM is applied to the energy equation, the heat flux for a mesh cell, included in the viscous flux vector Fv in Eq. (8), can be discretized as follows NF § ∂T ∂T · FD = ¦ ¨¨ k nx + k n y ¸ ΔS i ∂x ∂y ¸¹ i i =1 ©

(27)

where n x and n y are the components of the outward unit normal vector of a cell face,

ΔS represents the face area and N F denotes the face number of a cell. Eq. (27) indicates that the heat flux FD can be corrected once the temperature derivative corrections at the face centres are solved. Therefore, the key issue is to determine the temperature derivative correction by imposing the Neumann boundary condition. Different from existing IBMs using the cell centres as the Eulerian points, as for the present heat flux correction, the face centres should be defined as the Eulerian points, since the heat fluxes crossing cell faces as well as the temperature derivatives in Eq. (27) are evaluated at the face centres in the FVM. To satisfy the Neumann thermal boundary condition, it is required that when the cor-

rected temperature derivatives are interpolated from the Eulerian points to the Lagrangian points, the interpolated gradient should be equal to the given boundary value at the same position. Following the idea of previous implicit IBM for Dirichlet conditions, the temperature derivative corrections at the Lagrangian points are treated as unknowns and will be solved implicitly via interpolation and distribution processes. (1) Introduction of thin-plate splines (TPS) scheme As shown in Fig. 2, there are two sets of Eulerian meshes when the face centres are used as the Eulerian points. They are respectively composed of the face centres in the I- and Kdirection, which are denoted by rectangles and triangles in Fig. 2.

Lagrangian point

K

I

Fig. 2. Two sets of Eulerian meshes composed of face centres.

Applying the delta function interpolation to these two Eulerian meshes in Fig. 2 individually leads to

or

ªN º U B ( X Bl ) = « ¦ U ( x j ) D ( x j − X Bl ) h 2 » ¬ j =1 ¼I

(28)

ªN º U B ( X Bl ) = « ¦ U ( x j ) D ( x j − X Bl ) h 2 » ¬ j =1 ¼K

(29)

ªN º ªN º 2U B ( X Bl ) = « ¦ U ( x j ) D ( x j − X Bl ) h 2 » + « ¦ U ( x j ) D ( x j − X Bl ) h 2 » ¬ j =1 ¼ I ¬ j =1 ¼K

(30)

In the above equations, U stands for a scalar variable and the subscripts I and K represent the Eulerian meshes composed of the face centres in the I- and K-direction, respectively. Implicit immersed boundary algorithm can be then constructed according to Eq. (28) or (29) or (30). As will be shown in Section 3.1, when the cosine delta function is applied to the temperature derivative interpolation process, the convergent temperature distribution is less accurate. In this work, the more accurate TPS interpolation scheme is introduced to form hybrid TPS interpolation/delta function distribution. TPS provides a means to characterize an irregular surface by using functions that minimize an energy functional [38]. For a 2-D problem TPS reduces to the so-called infinite-plate splines (IPS) [35, 39]. The interpolated function of TPS is differentiable everywhere and the grid is not restricted to a rectangular array. In earlier studies on interpolation methods for fluid-structure interactions [39, 40], TPS turned out to be one of the most robust and accurate methods. Hence, TPS is very popular for the data transfer at the fluid-structure interface [39-41]. In Appendix A, the interpolation formula of TPS is derived from the wellestablished radial basis function (RBF) interpolation with thin-plate splines. (2) Temperature derivative interpolation by TPS For practicality and simplicity, the same face-centred Eulerian points that are used in the delta function interpolation are meanwhile used as known points for the TPS interpolation. So the number of Eulerian points involved in the TPS interpolation is still N . According to the TPS interpolation formula defined in Eq. (A.15), the temperature deriv-

ative, taking ∂T ∂x as an example, can be interpolated from the Eulerian points to the Lagrangian ones TxB = K TPS Tx

(31)

where

B x

ª ∂T T =« ¬« ∂x

∂T , X 1B ∂x

∂T ,!, ∂x X B2

º » » X BM ¼

ª ∂T Tx = « ¬« ∂x

∂T , x1 ∂x

∂T ,! , ∂x x2

º » » xN ¼

T

(32)

T

(33)

and K TPS is the transformation matrix for the TPS interpolation from the Eulerian to Lagrangian points, which is determined by the TPS function and the coordinates of known and interpolated points. Please refer to Appendix A for more details about K TPS . Note that, the dimension of K TPS is M × N , which is the same as that of the transformation matrix for delta function interpolation, such as Di in Eq. (22). Hence, Di can be directly replaced by K TPS to implement TPS interpolation. (3) Temperature derivative correction by hybrid TPS/delta function At first, the prescribed temperature gradient normal to the boundary is expressed as ∂T ∂n

= X Bl

∂T ∂x

X Bl

nx ( X Bl ) +

∂T ∂y

X Bl

n y ( X Bl )

for

( l = 1, 2," , M )

(34)

In the interpolation processes, ∂T ∂x is interpolated from the Eulerian points to the Lagrangian ones by using the following TPS interpolation

∂T ∂x

∂T j =1 ∂x N

X Bl

=¦

KljTPS xj

(35)

where ∂T ∂x

xj

§ ∂T · =¨ ¸ © ∂x ¹

*

xj

§ ∂T · +δ¨ ¸ © ∂x ¹ x j

(36)

and K ljTPS denote the corresponding elements of the matrix K TPS in Eq. (31). In the distribution process, the derivative corrections at the Lagrangian points are still distributed back to the Eulerian points by using the following delta function M § ∂T · § ∂T · = δ ¨ ¸ n x (X Bi )D (x j − X Bi )Δs i ¸ ¦ © ∂x ¹ x j i =1 © ∂n ¹ X Bi

δ¨

(37)

Substituting Eqs. (36) and (37) into Eq. (35) yields

∂T ∂x

X Bl

*

§ ∂T · = ¦¨ ¸ j =1 © ∂x ¹ N

xj

N M § ∂T · i i i TPS KljTPS + ¦¦ δ ¨ ¸ nx ( X B ) D ( x j − X B ) Δs Klj ∂ n i © ¹ j =1 i =1 XB

(38)

A similar expression as Eq. (38) can be derived for ∂T ∂y . Substituting Eq. (38) and the corresponding expression for ∂T ∂y into Eq. (34) leads to

∂T ∂n

*

*

X Bl

N § ∂T · = nx ( X Bl ) ¦ ¨ ¸ j =1 © ∂x ¹

xj

N § ∂T · KljTPS + ny ( X Bl ) ¦ ¨ ¸ j =1 © ∂y ¹

KljTPS xj

§ ∂T · i i i TPS + nx ( X ) ¦¦ δ ¨ ¸ nx ( X B ) D ( x j − X B ) Δs Klj n ∂ i © ¹ XB j =1 i =1 l B

N

M

(39)

N M § ∂T · i i i TPS + ny ( X Bl ) ¦¦ δ ¨ ¸ ny ( X B ) D ( x j − X B ) Δs Klj © ∂n ¹ X Bi j =1 i =1

By setting the term δ (∂T ∂n ) X i Δs i as an unknown, Eq. (39) can be further written as the B

following matrix form AX Tn = B

where

(40)

ª § ∂T · º § ∂T · § ∂T · X Tn = «δ ¨ ¸ Δs1 , δ ¨ ¸ Δs 2 ,", δ ¨ ¸ Δs M » © ∂n ¹ X BM © ∂n ¹ X B2 «¬ © ∂n ¹ X 1B »¼

(

) + diag ( n ( X ) , n ( X ) ,", n ( X ) ) K

T

(41a)

( ) D diag ( n ( X ) , n ( X ) ,", n ( X ) )

A = diag nx ( X 1B ) , nx ( X B2 ) ,", nx ( X BM ) K TPS Dd diag nx ( X 1B ) , nx ( X B2 ) ,", nx ( X BM ) y

ª ∂T B=« «¬ ∂n

1 B

∂T , X 1B ∂n

y

2 B

∂T ," , ∂n X B2

(

y

º » X BM » ¼

M B

TPS

y

1 B

y

2 B

M B

y

(41b)

T

)

ª ∂T * · «§¨ ¸ «© ∂x ¹ ¬

)

* ª § · T ∂ « «©¨ ∂y ¹¸ ¬

− diag nx ( X B1 ) , nx ( X B2 ) ," , nx ( X BM ) K TPS

(

d

− diag n y ( X 1B ) , n y ( X B2 ) ," , n y ( X BM ) K TPS

*

x1

§ ∂T · ,¨ ¸ © ∂x ¹

*

x2 *

x1

§ ∂T · ,¨ ¸ © ∂y ¹

§ ∂T · ," , ¨ ¸ © ∂x ¹

*

x2

§ ∂T · ," , ¨ ¸ © ∂y ¹

º » » xN ¼

T

º » » xN ¼

(41c) T

In the above expressions, the matrices K TPS and Dd are for the interpolation process and the distribution process, respectively. Dd is of the same form as that in Eq. (22). The linear system of Eq. (40) can be solved by a direct method, Gaussian method or LU decomposition. Here, the Gauss-Jordan elimination is used. Once the temperature derivative corrections at the Lagrangian points, i.e. X Tn , are solved from Eq. (40), the corrections at the Eulerian points can be evaluated via Eq. (37). (4) Heat flux correction The heat flux of a cell is to be corrected from the solved temperature derivative corrections at the Eulerian points. A constant ∂T ∂n on the boundary means that the heat flux always flows into or out of the solid domain across the boundary. In that case, the temperature in the solid domain will increase or decrease continually until a convergent state is reached and the temperature filed in the fluid domain will be seriously affected in return.

To eliminate this effect, the same amount of heat flux as that flowing into the fluid domain should be added to the solid domain. In other words, the heat flux of a solid cell should be corrected by using the additive inverses of the temperature derivative corrections in Eq. (37). Base on Eq. (27), the heat flux correction of a cell is then calculated from NF

ª

i =1

¬

∂T

§ ∂T ·

º

©

¼i

(δ FD ) F = ¦ « k ⋅ δ ¨§ ¸· nx + k ⋅ δ ¨ ¸ ny » ∂x ∂y ©

NF

ª

i =1

¬

¹

∂T © ∂x

(δ FD )S = −¦ « k ⋅ δ §¨

¹

§ ∂T · ¸ nx + k ⋅ δ ¨ ¹ © ∂y

ΔS i

· º ¸ n y » ΔSi ¹ ¼i

(42)

(43)

where the subscripts F and S refer to the fluid and solid cells, respectively. A mesh cell is defined as a fluid cell if its centre locates in the fluid domain. Otherwise, it is a solid cell. The heat flux of a mesh cell is finally corrected as

FD = FD* + δFD

(44)

where FD* is the intermediate heat flux calculated without considering the Neumann condition in the prediction step denoted by Eq. (5). In the present IBM for Neumann conditions, the temperature derivative correction and the heat flux correction are conducted in the correction step denoted by Eq. (6), and the corrected heat flux in Eq. (44) is used during solving Eq. (5) in the prediction step. 2.4 Solution procedure Hereto, an implicit heat flux correction-based IB-FVM has been proposed for thermal flows with Neumann conditions. Its main solution procedure can be outlined below: (1) Calculate the heat flux by using Eq. (44) and solve Eq. (5) to obtain the intermediate flow variables W * ;

(2) Solve Eq. (21) to get the velocity corrections at the Lagrangian points and correct the velocity field on the Eulerian mesh by using Eqs. (19) and (17); (3) Solve Eq. (23) to get the temperature corrections at the Lagrangian points and correct the temperature field on the Eulerian mesh by using Eqs. (25) and (26); (4) Solve Eq. (40) to get the temperature derivative corrections at the Lagrangian points and evaluate the heat flux corrections on the Eulerian mesh via Eqs. (37), (42) and (43); (5) Update the vector of the conservative variables W to the next time step (steady computation) or the next pseudo-time step (unsteady computation) according to the corrected velocity and temperature; (6) Repeat steps (1) to (5) until final solutions are obtained.

3. Numerical examples and discussions In this section, the proposed method will be applied to simulate several well documented thermal flows with Neumann boundary conditions, including the natural convection in a concentric horizontal annulus and the forced convection flows over stationary and oscillating circular cylinders. 3.1 Natural convection in a concentric horizontal annulus with a constant heat flux As illustrated in Fig. 3, the surface of the inner cylinder with radius R I is maintained at a constant heat flux qB = − k

∂T ∂T with = −1 in its outward normal direction and the ∂n ∂n

outer cylinder with radius RO = 2 RI is kept at a constant temperature T∞ . On the inner and outer cylinder surfaces, the no-slip condition is imposed. The flow behavior of this

problem is characterized by the Prandtl number Pr and the Rayleigh number Ra , which are respectively defined as Pr = c p

Ra =

(45)

μ k

c p ρ 2 gβL3 qB L k kμ

(46)

where L = RO − RI is the gap width of the annulus, μ is the dynamic viscosity coefficient and β is the thermal expansion coefficient.

Fig. 3. Schematic diagram for natural convection in an annulus.

A fixed Pr of 0.7 and two Ra numbers of 5700 and 50000 are considered here. As sketched in Fig. 3, the physical space is set as (5L×5L) and the annulus is placed at the centre of the square domain. At first, in order to examine the accuracy of the proposed method, four uniform Cartesian meshes are used in the simulation at Ra = 5700. The mesh spacings from fine to coarse are respectively h1=0.005, h2=0.01, h3=0.02 and h4=0.04, which have been divided by the gap width of the annulus L. To calculate norms of temperature errors, local estimates of the ex-

act solution are required. For the current second-order central scheme, estimates of the exact solution are obtained by extrapolating the solutions using the two finest meshes and the second-order Richardson extrapolation method [42]

T −T Texact = T1 − 22 1 r −1

(47)

where T1 and T2 are the numerical solutions on meshes h1 and h2, and the mesh refinement factor r = h2 h1 = 2 . The L2-norm and L’-norm are then calculated from NE

L2 -normk =

¦ (T

k ,i

− Texact ,i

i =1

)

2

(48)

NE

L∞ -normk = max Tk ,i − Texact ,i

(49)

where k indicates the mesh level and N E is the number of points included in the norm calculation. Here, all cell centres of the whole Cartesian mesh h4 are used. Table 1 lists the norms of temperature errors and the convergence rates. The norms versus the mesh spacing are plotted in Fig. 4 in a log-log scale. Because the exact solution is estimated from the solutions on meshes h1 and h2 by using the second-order Richardson extrapolation, the convergence rates at the mesh level h1 are exactly second order. At mesh levels h3 and h2, the convergence rates of L2-norm are a bit lower than second order. Table 1 Norms of temperature errors and convergence rates at Ra = 5700. Mesh spacing h h4=0.040 h3=0.020 h2=0.010 h1=0.005

L2-norm Norm 1.2992E-02 3.5631E-03 9.3030E-04 2.3257E-04

Order 1.87 1.94 2.00

L’-norm Norm 3.6746E-02 1.1771E-02 3.5832E-03 8.9581E-04

Order 1.64 1.72 2.00

Fig. 4. L2-norm and L’-norm of temperature errors versus mesh spacing at Ra = 5700.

Fig. 5 further compares the computed local temperature distributions on the inner cylinder surface at Ra = 5700 and 50000 with the results of Yoo [29] and those of Wang et al. [26], where the present results are calculated by using a mesh spacing of 0.02L. In Fig. 5, θ is the angular displacement from the uppermost point of the annulus and the dimensionless temperature T * is defined as

T* =

T − T∞ qB L k

(50)

As can be seen from Fig. 5, the present results at both Rayleigh numbers agree well with the reference data.

(a) Ra = 5700

(b) Ra =50000

Fig. 5. Comparison of local temperature distribution on inner cylinder surface.

Fig. 6 shows the isotherms and the streamlines at Ra = 5700 and 50000. In both cases, Ra is so high that the convection dominates the thermal flow. The flow is symmetric about

the vertical plane through the center of cylinders and a pair of crescent-shaped eddies are formed in the enclosure. Meanwhile, the strong convection induces a plume on the upper part of the annulus and the plume becomes stronger with the increase of Ra . The same patterns of streamline and isotherm have been reported by previous studies [24-26, 29]. 0.1

0

0.2

0

0.4

0.3

0.3 0.1

0.2

0.5

0

0.2

0.4

0.7

0.5

0

0.4

0.3 0.5 0.1

0.1

0.6

0.1

0.3

0.5

Ra = 5700

0.2

Ra = 5700 0

Pr = 0.7

0.4

Pr = 0.7

0.2

0 0.4

0.3

0.1

0.3 0.4

0.2

0.1

0.2 0 0.1

0

0

0.2

0.10.05 0.25

0 0.1 5

0.

5 0.40.45

0 .3

0.35 0.25

0.2

0

0.2

0.3

0.05 0.1

5

3 0.

2 0.

0.1 0. 05

15

0

Ra = 50000

0.15

Ra = 50000 0.3 0

0.15

Pr = 0.7

Pr = 0.7

0.1

0.0 5

2 0. 1 0.

0.05

0

0.25 0.25 0.2 0.15

0

Fig. 6. Isotherms (left) and streamlines (right) for natural convection in an annulus.

A numerical test is also conducted to further clarify the heat flux correction approach for a solid cell. At Ra = 5700, (δ FD ) S in Eq. (43) and its additive inverse − (δ FD ) S are used to correct the heat flux of a solid cell, respectively. Fig. 7 displays the temperature contours for both cases at the same numerical iteration number. In the case of − (δ FD ) S , the heat flux of a solid cell is corrected in the same way as a fluid cell and the temperature in the solid domain decreases continually over iteration, since the condition of ∂T ∂n = −1 always causes the heat flux to flow from the solid domain into the fluid domain. It can be seen from Fig. 7a that the dimensionless temperature for − (δ FD ) S is negative in the whole solid domain and gradually becomes positive around the inner cylinder. As compared with the temperature contour for (δ FD ) S in Fig. 7b, the temperature field in the fluid domain is seriously affected by the low temperature in the solid domain. Therefore, the heat flux correction for a solid cell should be conducted in such a way that the same amount of heat flux as that flowing into the fluid domain should be supplemented into the solid domain, which is formulated by Eq. (43).

(a) − (δ FD ) S

T* T

T* T

0.04 0.03 0.02 0.01 0 -0.01 -0.02 -0.03 -0.04 -0.05 -0.06 -0.07

0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

(b) (δ FD )S

Fig. 7. Temperature contour comparison between solid cell heat flux correction using − (δ FD ) S and

(δ FD )S . At the end of this section, the numerical results of delta function-based IBM are compared with those of TPS-based IBM. In Fig. 8, delta function interpolation/distribution (1) and (2) denote the results from Eqs. (28) and (30), respectively. These two distributions are nearly overlapped, which indicates that the solutions are consistent when the face centres in the I- or K- or both directions are used in the interpolation process. It also can be observed from Fig. 8 that the results of delta function-based IBM are less accurate as compared with the reference data while the results of TPS-based IBM show a good agreement. Because the only difference is the interpolation scheme used in the interpolation process of temperature derivative, current cosine delta function should be replaced by a more accurate interpolation method , such as TPS used here.

Fig. 8. Local temperature distribution comparison between delta function- and TPS-based IBMs.

3.2 Steady forced convection over a cylinder with a constant heat flux The forced convection over a stationary cylinder is steady in the low-Reynolds-number range. Its thermal flow patterns depend on the Prandtl number Pr and the Reynolds

number Re defined as

Re =

ρv∞ D μ

(51)

where v ∞ is the free-stream velocity and D denotes the cylinder diameter. On the surface of the cylinder, the uniform heat flux condition q B = − k

∂T ∂T with = −1 ∂n ∂n

in the outward normal direction and the no-slip condition are imposed. Steady cases are simulated at Pr = 0.7 and Re = 10, 20, 40, 45. The cylinder center is placed at (0, 0) and the physical space is set as [-14D, 20D]×[-15D, 15D]. A non-uniform Cartesian mesh is used for the whole rectangular domain and the subdomain of [ƺ0.6D, 0.6D]×[ƺ0.6D, 0.6D] in the vicinity of cylinder is discretized by a uniform fine mesh with a spacing of 0.01D. The cylinder is represented by 200 Lagrangian points. The heat transfer rate on the cylinder surface can be measured by the local Nusselt number Nu and the average Nusselt number Nu , which are respectively defined as

Nu ( X ) =

Nu =

qB D k [T ( X ) − T∞ ]

1 Nu( X )ds πD ³Γ

(52)

(53)

Fig. 9 compares the computed local Nusselt number distributions on the cylinder surface with the results of Bharti et al. [30] and those of Ren et al. [24], where θ is the angular displacement from the leading edge of cylinder. As compared with the results of Ren et al. [24], the present local Nusselt number distributions agree better with those of Bharti et al. [30]. The computed average Nusselt numbers are listed in Table 2, which also match well

with the reference data.

(a) Re=10

(b) Re=20

(c) Re=45 Fig. 9. Comparison of local Nusselt number distribution on cylinder surface at Re = 10, 20 and 45. Table 2 Comparison of average Nusselt number. Source Present results Bharti et al. [30] Ahmad et al. [31] Dhiman et al. [32] Ren et al. [24] Wang et al. [26]

Re = 10 2.0466 2.0400 2.0410 2.1463 2.0265 2.01

Re = 20 2.7858 2.7788 2.6620 2.8630 2.7413 2.69

Re = 40 3.7856 3.7755 3.4720 3.7930 3.7407 3.68

Re = 45 3.9915 3.9727 -

Fig. 10 shows the isotherms and streamlines for Re=10, 20 and 40, where the dimensionless temperature T* equals the reciprocal of local Nu. A pair of symmetric eddies are formed behind the cylinder and the recirculation length increases with the increase of Re. At each

Re, the maximum clustering of isotherms appears in the front surface of the cylinder, which indicates that the temperature gradient or the heat transfer rate there is higher than that in other regions. With the increase of Re, the clustering of isotherms increases and the temperature around the cylinder decreases, which manifests that the convection is enhanced. Those observations are in accordance with previous studies [24, 26, 30]. 0.0

Re=10

5

1.5

Re=10

1

5 0 .1

0.15

0.0

5

0.25 0.3 5

0

0.0

0. 15

0.55

0.2 5

0.45 0.55 5 0.6

0.45

0.25 0.15 0.05

0.35

0.5

-0.5 -1

0.15

5

-1.5 -1

Re=20

1.5

0.05

0

1

2

3

4

0

1

2

3

4

0

1

2

3

4

Re=20

1 0.15 0.0

5

0.5

5

5 0 .5

5 5 0.1 0.0

0.35

0

5

0 .2

0 .1

5

0.25

0 .3

0. 45

-0.5 0.15

-1 0.05

-1.5 -1

1.5

Re=40

0.5

5 0.0

0.15

0.45

0

0.25

0.1

0.05

5 35 0.

5 0 .2 5

0.0

Re=40

1

-0.5 0.05

-1

0.05

-1.5 -1

Fig. 10. Isotherms (left) and streamlines (right) for steady forced convection over a cylinder.

Table 3 lists the computed drag coefficients and the recirculation lengths at Re = 20, 40 as well as some reference data. The drag coefficient is defined as

CD =

FD 0.5ρv∞2 D

(54)

with FD the drag force. Benefiting from the boundary condition-enforced IBM, the force

F acting on the cylinder can be calculated by utilizing the momentum correction F =−

δ (ρv ) δt

(55)

The comparison in Table 3 demonstrates that both the drag coefficients and the recirculation lengths are in excellent agreement with the reference data. Table 3 Comparison of drag coefficient Cd and recirculation length Lw. Source Present results Dennis et al. [33] He et al. [34] Wu et al. [11] Ren et al. [24]

Cd 2.077 2.05 2.152 2.072 2.126

Re=20 Lw 0.957 0.94 0.921 0.92 0.913

Cd 1.556 1.52 1.499 1.554 1.568

Re=40 Lw 2.343 2.35 2.245 2.3 2.331

3.3 Unsteady convection over a stationary/oscillating cylinder with a constant heat flux In this section, Re is increased to 200 to investigate the unsteady thermal flows over a stationary/oscillating cylinder under the same heat flux condition as that in Section 3.2. At Re=200, the flow is characterized by the Karman vortex street shedding from the cylinder surface and the temperature field presents a similar shedding pattern due to the velocity field [23]. For moving boundary problems, only the position of moving boundary is tracked in the IBM. So it adds little extra computational cost. In this work, an oscillating cylinder in the cross-wise direction is considered and its dimensionless equation of motion is written as

y ( t ) = A sin ( 2π ft

)

(56)

Here, f is the dimensionless oscillating frequency defined as f = St = ( fD ) / v∞ . A and y are the dimensionless oscillating amplitude and displacement, respectively, which are

both non-dimensionalized by D . According to the test cases used in [23] and [36],

Re = 200 and A = 0.15 are selected, and f is set as 0.2 which is very close to the vortex shedding frequency from a stationary cylinder at Re = 200 . Based on the Cartesian mesh used in Section 3.2, for these two unsteady flow computations, the sub-domain discretized by the uniform fine mesh with a spacing of 0.01D is extended to a region of 4D×8D around the cylinder. The dual time-stepping approach is ap-

(

plied to both cases and the physical time steps are both set as Δt = 1 (180 f ) = D 180v∞ f

)

according to comparative computations. Fig. 11 shows the time history of surface-averaged Nu for the thermal flow over a stationary cylinder at Re = 200. It is clear that the surface-averaged Nu varies over time periodically. Fig. 12 compares the computed time-averaged local Nu distribution on the stationary cylinder surface with the result of Zhang et. al [23] and that of Hu et. al [36]. Fig. 13 shows a comparison of the time-averaged local T* distribution on the oscillating cylinder surface. As can be seen from Fig. 12 and Fig. 13, the present results agree well with the literature data in both cases, which indicates that the proposed method can efficiently handle the unsteady thermal flows with Neumann conditions, including moving boundary problems.

Fig. 11. Time history of surface-averaged Nu for stationary cylinder at Re = 200.

Fig. 12. Comparison of time-averaged local Nu distribution on stationary cylinder surface at Re = 200.

Fig. 13. Comparison of time-averaged local T* distribution on oscillating cylinder surface at Re = 200.

4. Conclusions In this study, an implicit heat flux correction-based IB-FVM has been proposed for ther-

mal flows with Neumann boundary conditions. With respect to the heat flux correction, the face centres should be defined as the Eulerian points within the scope of FVM. To satisfy the Neumann condition accurately, the temperature derivative corrections at the Lagrangian points are solved implicitly by constructing hybrid TPS interpolation/delta function distribution. TPS is introduced to the derivative interpolation process to improve the interpolation accuracy significantly because the use of cosine delta function leads to a less accurate solution. The heat flux of a fluid cell is then corrected by using the solved temperature derivative corrections at the cell faces, while that of a solid cell should be corrected by using their additive inverses to supplement the same amount of heat flux into the solid domain as that flowing into the fluid domain. As compared with existing implicit IBMs, the Neumann condition is treated strictly within the scope of IBM and the introduction of auxiliary layer as well as the conversion from the Neumann to Dirichlet condition is not needed. Therefore, the present method is suitable for an arbitrary geometry. The proposed method has been well verified by simulating the natural convection in an annulus and the steady or unsteady forced convection over a stationary or oscillating cylinder with Neumann boundary conditions. All the computed results show a good agreement with the reference data in the literature. In this study, only the cosine delta function is compared with TPS. Tornberg and Engquist [43] have established that alternative regularized delta functions provide better accuracy than the cosine delta function, such as the piecewise linear and piecewise cubic functions. Therefore, other much more accurate functions need to be further tested. In the future work, the present method will be extended to 3-D simulations, and will be

improved to study thermal compressible flows.

Acknowledgements This work was supported by the National Natural Science Foundation of China (Grant No. 11872212).

Appendix A. Derivation of TPS interpolation scheme The interpolation formula of TPS is to be established from the radial basis function (RBF) interpolation by using thin-plate splines as basis functions. In RBF, the interpolation function f is defined as a sum of basis functions N

(

)

f ( x ) = ¦ ai ⋅ ij x - xiC +ȥ ( x ) i=1

(A.1)

where x = [ x, y,z ] is the coordinate vector of a point, the superscript C denotes the scattered centers with known values, N represents the number of centers, ij is a given basis function with respect to the Euclidean distance < , and ȥ = b0 +b1 x +b2 y +b3 z is a linear

( i = 1, 2," , N )

polynomial. All unknown coefficients, including ai

and b0 , b1 , b2 , b3 , are de-

termined by the interpolation conditions

f ( xiC ) = diC for ( i = 1, 2," , N )

(A.2)

together with the side conditions N

N

¦a = ¦a x i

i =1

i

i =1

C i

N

N

i =1

i =1

= ¦ ai yiC = ¦ ai ziC = 0

(A.3)

C In Eq. (A.2), di is the known value at the i-th center. Eqs. (A.2) and (A.3) can be ar-

ranged into the following matrix form Ka E = d CE

(A.4)

where the extended vectors a E , d CE and the matrix K are respectively defined as

a E = [ a1 , a2 ," , aN , b0 , b1 , b2 , b3 ]

T

d CE = ª¬ d1C , d 2C ," , d NC , 0, 0, 0, 0 º¼

T

(A.5) (A.6)

ª K11 « « K 21 « # « K K = « N1 « 1 « C « x1 « yC « 1 «¬ z1C

(

C C with K ij = ϕ xi − x j

K1N

1 x1C

y1C

K 22 " K 2 N # # # K N 2 " K NN

1 x2C # # 1 xNC

y2C # yNC

K12

"

1 x2C y2C

" " "

1 xNC yNC

0 0 0

0 0 0

0 0 0

z2C

"

z NC

0

0

0

z1C º » z2C » # » » z NC » 0» » 0» 0» » 0 »¼

(A.7)

).

The unknown vector a E is then solved from Eq. (A.4) a E = K −1d CE

(A.8)

Hereto, the interpolation function f ( x ) in Eq. (A.1) is determined. Next, function valA A A A ues at M interpolated points xi = ª¬ xi , yi , zi º¼ are evaluated from

N

(

)

f ( xiA ) = ¦ a j ⋅ ij xiA - x Cj +ȥ ( xiA ) = diA for j=1

( i = 1, 2," , M )

(A.9)

The matrix form of the above equations is d A = Ka E

(A.10)

where d A and K are respective defined as

d A = ª¬ d1A , d 2A ," , d MA º¼ ª K11 « K K = « 21 « # « «¬ K M 1

(

K12 " K1N K 22 " K 2 N # # # KM 2 K MN

)

A C with K ij = ϕ xi − x j .

Substituting Eq. (A.8) into Eq. (A.10) yields

1 x1A 1 x2A # # 1 xMA

T

(A.11) y1A y2A # yMA

z1A º » z2A » # » » zMA »¼

(A.12)

d A = KK −1d CE = Td CE

(A.13)

In the above equation, T has a dimension M × ( N + 4 ) and denotes the transformation matrix for RBF interpolation from N centers to M interpolated points, which is related to the specified basis function ij ( •

)

and coordinates of centers and interpolated points.

Since the last four components of d CE are all zero, Eq. (A.13) can be further written as

d A = K ϕd C

(A.14)

T

where d C = ª¬ d1C , d 2C ," , d NC º¼ and the modified transformation matrix K ϕ with a dimension of M × N is obtained by removing the last four columns of T . Now, function values at M interpolated points can be evaluated from N centers with known values d C by directly using Eq. (A.14). Finally, by choosing the thin plate splines as basis functions, the TPS interpolation scheme can be established according to Eq. (A.14), i.e. d A = K TPS d C

(A.15)

where the transformation matrix K TPS for TPS interpolation is computed by setting ϕ ( < ) = < log < 2

(A.16)

References [1] C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (1977) 220–252. [2] Z.-G. Feng, E.E. Michaelides, Robust treatment of no-slip boundary condition and velocity updating for the lattice-Boltzmann simulation of particulate flows, Comput.

Fluids 38 (2009) 370–381. [3] E.A. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof, Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations, J. Comput. Phys. 161 (2000) 35–60. [4] J. Kim, D. Kim, H. Choi, An immersed-boundary finite-volume method for simulations of flow in complex geometries, J. Comput. Phys. 171 (2001) 132–150. [5] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comput. Phys. 209 (2005) 448–476. [6] N. Zhang, Z.C. Zheng, An improved direct-forcing immersed-boundary method for finite difference applications, J. Comput. Phys. 221 (2007) 250–268. [7] A.J.C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation Part I. Theoretical foundation, J. Fluid Mech. 271 (1994) 285–310. [8] X.D. Niu, C. Shu, Y.T. Chew, Y. Peng, A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows, Phys. Lett. A 354 (2006) 173–182. [9] Y. Hu, H. Yuan, S. Shu, X. Niu, M. Li, An improved momentum exchanged-based immersed boundary–lattice Boltzmann method by using an iterative technique, Comput. Math. Appl. 68 (2014) 140–155. [10] H.-Z. Yuan, X.-D. Niu, S. Shu, M. Li, H. Yamaguchi, A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating a flexible filament in an incompressible flow, Comput. Math. Appl. 67 (2014) 1039–1056. [11] J. Wu, C. Shu, Implicit velocity correction-based immersed boundary-lattice Boltzmann

method and its applications, J. Comput. Phys. 228 (2009) 1963–1979. [12] C. Shu, N.Y. Liu, Y.T. Chew, A novel immersed boundary velocity correction-lattice Boltzmann method and its application to simulate flow past a circular cylinder, J. Comp. Phys. 226 (2007) 1607–1622. [13] J. Wu, C. Shu, An improved immersed boundary-lattice Boltzmann method for simulating three-dimensional incompressible flows, J. Comput. Phys. 229 (2010) 5022– 5042. [14] Y. Wang, C. Shu, C.J. Teo, J. Wu, An immersed boundary-lattice Boltzmann flux solver and its applications to fluid–structure interaction problems, J. Fluids Struct. 54 (2015) 440–465. [15] Y.L. Qiu, C. Shu, J. Wu, Y. Sun, L.M. Yang, T.Q. Guo, A boundary condition-enforced immersed boundary method for compressible viscous flows, Comput. Fluids 136 (2016) 104–113. [16] J.R. Pacheco, A. Pacheco-Vega, T. Rodic, R.E. Peck, Numerical simulations of heat transfer and fluid flow problems using an immersed-boundary finite-volume method on nonstaggered grids, Numer. Heat Transfer Part B: Fundam. 48 (2005) 1–24. [17] Z. Wang, J. Fan, L. Kun, K. Cen, Immersed boundary method for the simulation of flows with heat transfer, Int. J. Heat Mass Transfer 52 (2009) 4510–4518. [18] W.W. Ren, C. Shu, J. Wu, W.M. Yang, Boundary condition-enforced immersed boundary method for thermal flow problems with Dirichlet temperature condition and its applications, Comput. Fluids 57 (2012) 40–51. [19] D.L. Young, Y.J. Jan, C.L. Chiu, A novel immersed boundary procedure for flow and

heat simulations with moving boundary, Comput. Fluids 38 (2009) 1145–1159. [20] Z. Feng, E.E. Michaelides, Heat transfer in particulate flows with direct numerical simulation (DNS), Int. J. Heat Mass Transfer 52 (2009) 777–786. [21] B.S. Kim, D.S. Lee, M.Y. Ha, H.S. Yoon, A numerical study of natural convection in a square enclosure with a circular cylinder at different vertical locations, Int. J. Heat Mass Transfer 51 (2008) 1888–1906. [22] J.M. Lee, M.Y. Ha, H.S. Yoon, Natural convection in a square enclosure with a circular cylinder at different horizontal and diagonal locations, Int. J. Heat Mass Transfer 53 (2010) 5905–5919. [23] N. Zhang, Z.C. Zheng, S. Eckels, Study of heat-transfer on the surface of a circular cylinder in flow using an immersed-boundary method, Int. J. Heat Fluid Flow 29 (2008) 1558–1566. [24] W. Ren, C. Shu, W. Yang, An efficient immersed boundary method for thermal flow problems with heat flux boundary conditions, Int. J. Heat Mass Transf. 64 (2013) 694– 705. [25] Y. Hu, D. Li, S. Shu, X. Niu, Study of multiple steady solutions for the 2D natural convection in a concentric horizontal annulus with a constant heat flux wall using immersed boundary-lattice Boltzmann method, Int. J. Heat Mass Transf. 81 (2015) 591– 601. [26] Y. Wang, C. Shu, L.M. Yang, Boundary condition-enforced immersed boundary-lattice Boltzmann flux solver for thermal flows with Neumann boundary conditions, J. Comput. Phys. 306 (2016) 237–252.

[27] J. M. Weiss, W.A. Smith, Preconditioning applied to variable and constant density flows. AIAA Journal 33 (1995) 2050-2057. [28] C.S. Peskin, The immersed boundary method, Acta Numer. 11 (2002) 479–517. [29] J.S. Yoo, Dual free-convective flows in a horizontal annulus with a constant heat flux wall, Int. J. Heat Mass Transf. 46 (2003) 2499–2503. [30] R. Bharti, R.P. Chhabra, V. Eswaran, A numerical study of the steady forced convection heat transfer from an unconfined circular cylinder, Heat Mass Transf. 43 (2007) 639–648. [31] R.A. Ahmad, Z.H. Qureshi, Laminar mixed convection from a uniform heat flux horizontal cylinder in a crossflow, J. Thermophys. Heat Transfer 6 (1992) 277–287. [32] A.K. Dhiman, R.P. Chhabra, A. Sharma, V. Eswaran, Effects of Reynolds and Prandtl numbers on the heat transfer across a square cylinder in the steady flow regime, Numer. Heat Transfer 49 (2006) 717–731. [33] S.C.R. Dennis, G.Z. Chang, Numerical solution for steady flow past a circular cylinder at Reynolds number up to 100, J. Fluid Mech. 42 (1970) 471–489. [34] X. He, G.D. Doolen, Lattice Boltzmann method on curvilinear coordinates system: flow around a circular cylinder, J. Comput. Phys. 134 (1997) 306–315. [35] R.L. Harder, R.N. Desmarais, Interpolation using surface splines, Journal of Aircraft 9 (1972) 189–191. [36] Y. Hu, D. Li, S. Shu, X. Niu, An efficient immersed boundary-lattice Boltzmann method for the simulation of thermal flow problems, Commun. Comput. Phys. 20(2016) 12101257. [37] A. L. Gaitonde, A dual-time method for the solution of the 2D unsteady Navier-Stokes

equations on structured moving meshes, AIAA-95-1877-CP. [38] J. Duchon, Splines minimizing rotation-invariant semi-norms in Sobolev spaces. Constructive theory of functions of several variables, edited by W. Schempp and K. Zeller, Oberwolfach (1976) 85-100. [39] M. J. Smith, D. H. Hodges, C. E. S. Cesnik, Evaluation of computational algorithms suitable for fluid-structure interactions, Journal of Aircraft 37 (2000) 282–294. [40] A. H. Van Zuijlen, A. de Boer, H. Bijl, Higher-order time integration through smooth mesh deformation for 3D fluid-structure interaction simulations, Journal of Computational Physics 224 (2007) 414-430. [41] S. Keye, D. Mavriplis, Summary of case 5 from sixth drag prediction workshop: coupled aerostructural simulation, Journal of Aircraft 55 (2018) 1380-1387. [42] C. J. Roy, Grid convergence error analysis for mixed-order numerical schemes, AIAA journal 41 (2003) 595-604. [43] A. K. Tornberg, B. Engquist, Numerical approximations of singular source terms in differential equations, J. Comput. Phys. 200 (2004) 462–488.