A high order time-accurate loosely-coupled solution algorithm for unsteady conjugate heat transfer problems

Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217

Contents lists available at SciVerse ScienceDirect

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

A high order time-accurate loosely-coupled solution algorithm for unsteady conjugate heat transfer problems V. Kazemi-Kamyab ⇑,1, A.H. van Zuijlen, H. Bijl Aerodynamics Section, Faculty of Aerospace Engineering, Delft University of Technology, P.O. Box 5058, 2600 GB Delft, The Netherlands

a r t i c l e

i n f o

Article history: Received 11 December 2012 Received in revised form 22 April 2013 Accepted 24 May 2013 Available online 6 June 2013 Keywords: Unsteady conjugate heat transfer High order implicit time integration Implicit–explicit (IMEX) Partitioned method Loosely-coupled

a b s t r a c t Thermal interaction of fluids and solids, or conjugate heat transfer (CHT), is encountered in many engineering applications. Noting that time-accurate computations of transient CHT problems can be computationally expensive, we consider the use of high order implicit time integration schemes which have the potential to be computationally more efficient relative to the commonly used second order implicit schemes. For thermally weak couplings, we present a loosely-coupled solution algorithm where high order implicit–explicit (IMEX) Runge–Kutta schemes are employed for time integration. The IMEX schemes consist of the explicit first-stage singly diagonally implicit Runge–Kutta (ESDIRK) schemes, for advancing the solution in time within each separate fluid and solid subdomain, and the explicit Runge–Kutta (ERK) schemes, for explicit integration of part of the coupling terms. By considering a numerical example (an unsteady conjugate natural convection in an enclosure), temporal order preservation of the coupling algorithm (without subiterating) is demonstrated. In addition, the stability of the loosely-coupled algorithm is investigated numerically for the CHT test-case; when the ratio of the thermal effusivities of the fluid and solid subdomains is much smaller than unity, using large Fourier numbers of the subdomains is possible, indicating that time-step size is restricted by accuracy rather than stability. Furthermore, the (computational) work-(temporal) precision character of several time integration schemes in solving the CHT test-case is compared over a range of accuracy requirements; for time-accurate solutions, the fourth and fifth order IMEX schemes are 1.5 times more efficient than Crank–Nicolson and 2.7 times more efficient than BDF2. The computational gain is higher for smaller tolerances. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Thermal interaction of flows and structures, also referred to as conjugate heat transfer, arises in many engineering applications. Examples are: the thermal interaction among various phases and materials in the mold region of a continuous casting process, the thermal interaction between the fluid in a U-tube pipe with the surrounding soil in a geothermal heat exchanger system, the cooling of gas turbine blades, and the cooling of electronic chips. In order to obtain a better understanding of the physics of the coupled problem and hence to increase the efficiency and/or safety of designs, numerical simulations serve as a viable tool. The monolithic and partitioned approaches are two commonly used methods for solving the thermal coupling of flows and structures. In the monolithic method, the solution in the global domain is obtained by solving the governing equations within the subdomains as well as the interface equations simultaneously [1,2]. The monolithic method requires the production of a single code ⇑ Corresponding author. 1

E-mail address: [email protected] (V. Kazemi-Kamyab). This author funded by STW Grant 10113.

0045-7825/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2013.05.021

specifically tailored for conjugate heat transfer problems [3]. In the partitioned method on the other hand, a separate physics solver is associated with each subdomain [4]. The solution in the subdomains are coupled at the interface through a set of transmission conditions, and a coupling algorithm is required for the transfer of data between the subdomains. By solving the coupled problem in a partitioned manner, one can take advantage of the already existing efficient and highly optimized separate fluid and solid codes. In engineering applications, typically an implicit time integration is preferred over an explicit one in order to circumvent time step restrictions due to probable stiffness in the problem. Stiffness in a system can, for example, arise due to the nature of the governing equations or due to the generated grid (such as clustering of nodes near an area of interest [5]). Performing time-accurate computations of transient CHT problems can be computationally demanding, in particular when low order time integration schemes are used. The obtained solution can suffer from low levels of temporal accuracy, and in order to increase the accuracy of the solution, smaller time steps must be taken. This results in an increase in the computational cost of solving the coupled problem. As a potential solution, we consider the use of high order implicit time integration schemes for advancing the coupled problem in time.

206

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217

Computational efficiency of the high order time integration schemes relative to second order backward difference scheme (BDF2) has been demonstrated in Bijl et al. [6] for fluid flow computations and in van Zuijlen and Bijl [7] for the partitioned simulation of the mechanical coupling of flows and structures. Using high order spatial discretization in combination with high order temporal discretization can also contribute to the overall efficiency of the method. In this paper however, following the method of lines approach, we will assume a spatial discretization, and focus on time integration. When the coupled problem is solved monolithically, using an implicit time integration scheme, the interface equations are resolved implicitly. However, in the partitioned method, some or all of the interface terms are treated explicitly, depending on the arrangement with which the two coupled domains are solved (parallel (Block Jacobi) or sequential (Block Gauss–Seidel)). Therefore, in the partitioned approach, the interface equations are treated in a segregated manner where one of the interface equations is applied as a boundary condition for one subdomain and the other as boundary condition for the second subdomain [4]. If at each time step (or stage of an implicit Runge–Kutta scheme), a single interface iteration (subiteration or fixed-point iteration) [4,8] is performed, the partitioned algorithm is referred to as looselycoupled, otherwise it is referred to as strongly-coupled. This paper focuses on loosely-coupled solution algorithms; see [4,9] for examples on strongly coupled algorithms, and [3,4,10–12] for examples on partitioned algorithms with explicit time integration. Loosely-coupled solution algorithms can provide an efficient way of solving time-accurate CHT problems relative to the monolithic approach when the thermal coupling between the subdomains is weak (when the ratio of the thermal effusivities pffiffiffiffiffiffiffiffiffiffi (e ¼ kqcp ) of the coupled domains is much smaller than unity) [13]. To the authors’ knowledge, loosely-coupled solution algorithms with up to second order implicit time integration schemes have been reported in the literature (see [3,9,13,14]). In this paper, a loosely-coupled solution algorithm is presented in which a family of high order implicit–explicit (IMEX) Runge–Kutta schemes are used for time integration. The IMEX schemes consist of the explicit first-stage singly diagonally implicit Runge–Kutta (ESDIRK) schemes, which are used for advancing the solution in time within each separate fluid and solid subdomain, and equal order and number of stages explicit Runge–Kutta (ERK) schemes for explicit integration of part of the coupling terms. The IMEX schemes considered here were originally developed for solving time-accurate convection–diffusion-reaction (CDR) problems [15] and later employed for the loosely coupled simulation of the mechanical coupling of flows and structures [7]. The solution obtained using a loosely-coupled algorithm, contains an additional source of error compared to the monolithic solution, denoted as the partitioning error. As a result, the temporal accuracy and stability of the coupling algorithm is influenced by the partitioning error. Therefore, in designing loosely-coupled solution algorithms, a number of issues needs to be considered. One, whether the design order of the time integration scheme is preserved without subiterating. Second, what are the stability properties of the algorithm; for practical computations, it is preferred that Dt is restricted by accuracy rather than stability. These two issues are investigated numerically by considering a CHT problem (unsteady conjugate natural convection in an enclosure). While in the loosely coupled multi-stage IMEX schemes a single interface iteration is performed at each (implicit) stage, in the second order loosely coupled Crank–Nicolson scheme [13], only one is performed per time-step. However, for the same time-step, the high order IMEX schemes generally provide temporally more accurate solutions. For the CHT test-case, the (computational) work-(temporal) precision character of the high order IMEX and

second order Crank–Nicolson and BDF2 schemes is compared over a range of accuracy requirements. We investigate whether the high order IMEX schemes can compete with the second order schemes for a reasonable portion of the work-precision spectrum, i.e. whether the additional work per time-step of the IMEX schemes is compensated by the gain in temporal accuracy. In what follows, first the equations governing conjugate heat transfer are discussed. After a brief overview of the ESDIRK and IMEX time integration schemes, the details of the loosely coupled solution algorithm are presented. Next, numerical examples are considered, in order to demonstrate the applicability of the algorithm, to investigate the temporal order preservation and stability of the algorithm, and finally its computational efficiency relative to the second order time integration schemes. 2. Governing equations In the conjugate heat transfer problem considered here, the fluid domain is modeled using the Boussinesq approximation of the Navier–Stokes system which in primitive variables is given by:

r  u ¼ 0;

ð1Þ

@u ¼ ðu  rÞu þ mr2 u  rp  bgjðT f  T ref Þ; @t

ð2Þ

@T f 1 ¼ ðu  rÞT f þ r  ðkf rT f Þ; ðqcp Þf @t

ð3Þ

where u is the velocity vector, p is the kinematic pressure, m the kinematic viscosity, T f the temperature, kf the thermal conductivity, cp;f the heat capacity, qf the density, b the compressibility factor of the fluid, g the acceleration due to gravity, and j is a vector indicating the direction in which the gravity acts. The solid domain is modeled using unsteady heat conduction:

@T s 1 ¼ r  ðks rT s Þ: ðqcp Þs @t

ð4Þ

The governing equations are accompanied by appropriate initial and boundary conditions. For a well-posed problem, the continuity of the temperature and heat flux are imposed at the common interface (I ) of the domains:

T f ðI ; tÞ ¼ T s ðI ; tÞ;

ð5Þ

qs ðI; tÞ ¼ qf ðI ; tÞ;

ð6Þ

where qm ðI ; tÞ ¼ km rT m  n with m the index of the subdomain and n the outward normal of the interface. To identify the governing parameters in the conjugate heat transfer problem described by (1)–(6), the equations are nondimensionalized using appropriate non-dimensional quantities (see [13]). Based on the dimensionless form of the equations, it is observed that the Prandtl number (Pr ¼ am ), the Rayleigh number f   gbðT H T C ÞL3ref , the ratio of the thermal conductivities (kks ), and Ra ¼ maf f   as the ratio of thermal diffusivities of the domains a are the govf

erning parameters of the problem. 3. Model problem In this section, the description of a one dimensional model problem is presented which will be used to discuss the details of the loosely-coupled solution algorithm. The model problem has been commonly used in the literature (for example [3,4,10]) to analyze stability of numerical algorithms for thermal coupling of domains. The model problem consists of thermal coupling of two domains X1 ¼ ½L1 ; 0 and X2 ¼ ½0; L2 , with their common

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217

interface I at x ¼ 0. The governing equation within each sub-domain is transient linear heat conduction. Since there is no flow normal to the interface in the CHT problem of Section 2, convection can be neglected from the 1-D problem [3]. For simplicity, each sub-domain has constant material properties. The initial boundary value problem is given by:

@T m ðx; tÞ @q ¼  m ; m ¼ 1; 2; @t @x @T m ðx; tÞ qm ðx; tÞ ¼ km ; m ¼ 1; 2; @x T 1 ðx ¼ L1 ; tÞ ¼ T 1;lbc ;

ðqcp Þm

T 2 ðx ¼ L2 ; tÞ ¼ T 2;rbc ; T 2 ðI; tÞ ¼ T 1 ðI; tÞ; q1 ðI; tÞ ¼ q2 ðI ; tÞ;

ð9Þ

4.1. Imposing interface boundary conditions

ð10Þ ð11Þ ð12Þ ð13Þ ð14Þ

where T m and qm respectively represent the temperature and heat flux fields within each domain, with m the index of the subdomain. The non-interface boundaries are constant Dirichlet conditions. An analytical solution to the model problem can be obtained in terms of Fourier series using the method of separation of variables (see [16, chapter 6]). However, a simpler form of the analytical solution can be obtained if the two coupled domains are considered as semi-infinite. The non-interface boundary conditions (x ! 1) are given by:

T 1 ðx ! 1; tÞ ¼ T 1;lbc ;

ð15Þ

T 2 ðx ! 1; tÞ ¼ T 2;rbc :

ð16Þ

The solution to this coupled problem is given by Duchaine [12] and Carslaw and Jaeger [17]:

  jxj p ffiffiffiffiffiffiffiffiffiffi ; ðT  T Þerfc 2;rbc 1;lbc 1 þ ee21 4a1 t   1 x T 2 ¼ T 2;rbc  e2 ðT 2;rbc  T 1;lbc Þerfc pffiffiffiffiffiffiffiffiffiffi ; 1 þ e1 4a2 t

ð17Þ ð18Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where em ¼ km qm cp;m is the thermal effusivity within each subdomain. The interface temperature remains constant with time and is given by:

T I ;exact ¼

T 1;lbc þ ee21 T 2;rbc 1 þ ee21

;

ð19Þ

while the heat flux at the interface decreases with time and is given by:

qI;exact ¼

T 2;rbc  T 1;lbc e2 pffiffiffiffiffiffi : 1 þ ee21 pt

ð20Þ

It is noted that the ratio of the thermal effusivities of the coupled domains,

r¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffi k2 q2 cp;2 k2 a1 e2 ¼ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; e1 k1 q1 cp;1 k1 a2

4. Semi-discretized form of the coupled problem

ð8Þ

ð7Þ

T 1 ðx; t ¼ 0Þ ¼ T 1;lbc ;

e2 e1

r ! 1, the

The method of lines [18] is used to solve the PDEs, in which the spatial operators are discretized first in order to obtain the semidiscrete form of the equations. Since the global domain X is split into X1 and X2 , prior to space discretization, interface boundary conditions must be imposed on the two subdomains.

T 2 ðx; t ¼ 0Þ ¼ T 2;rbc ;

T 1 ¼ T 1;lbc þ

when r  1, the thermal coupling is weak, and as strength of the coupling increases.

207

ð21Þ

is one of the governing parameters of the coupled problem. Assuming that r  1, the variation of the temperature field in X1 from the imposed initial condition is small, and changes mainly occur in X2 . As this ratio approaches unity, both the subdomains experience temperature variations. Therefore, it appears that r can be used as a measure of the strength of the thermal coupling:

In this paper, the Dirichlet–Neumann transmission conditions are considered for coupling the two domains at the interface. For partitioned algorithms, the stability of the coupling algorithm depends on the correct assignment of the interface conditions. Using a two-dimensional variant of the model problem in Section 3, Henshaw and Chand [4] analyzed the stability and the rate of convergence of Dirichlet–Neumann interface iterations, using the h scheme for time integration and without any particular assumption on the spatial discretization. They arrived at the following estimate for the rate of convergence of the iterations for the smooth components of the solution:

U

k2 k1

rffiffiffiffiffi a1 ;

ð22Þ

a2

where km and am (with m ¼ 1; 2) represent the thermal conductivity and thermal diffusivity of each subdomain. The estimate was obtained by imposing the Dirichlet (temperature) condition on X2 with the interface temperature of X1 being prescribed as its value, and the Neumann (the heat flux) condition on X1 with the interface heat flux of X2 being prescribed as its value. Performing interface iterations based on the assigned interface boundary conditions is stable if U < 1, otherwise the two interface boundary conditions must be interchanged. While the estimate (22) has been obtained for a stronglycoupled solution algorithm, it can also be used as a criterion for imposing interface boundary conditions for loosely-coupled solution algorithms; that is, for stability of a loosely-coupled partitioned algorithm, satisfying (22) is necessary but might not be sufficient. It is interesting to note that the expression in (22) is equal to the ratio of thermal effusivities of the coupled domains, i.e. U ¼ r. Therefore, the domain with the higher e is assigned the Dirichlet condition and the one with the lower e the heat flux condition. Here, it is assumed that based on the material properties of the two subdomains, the imposed interface boundary conditions satisfy the criterion (22), i.e. the temperature condition is assigned to X2 and the heat flux condition to X1 . 4.2. Space-discretization and semi-discrete form With the interface boundary conditions imposed, the two domains are now discretized in space using the cell-centered finite volume (CFV) method (see Fig. 1). Applying the volume integral form of the governing Eq. (7) to each control volume results in:

ðDxqcp Þm

d T m;j ¼ qm jjþ1=2 j1=2 dt

m ¼ 1; 2;

ð23Þ

where for simplicity cells of equal size are used within each subdomain. Substituting (8) into (23) for qm and using central difference to approximate the temperature gradients at the cell faces, the semi-discrete form of Eq. (7) for the interior cells, after some rearrangements, is given by:

208

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217

Fig. 1. Discretization of X1 and X2 subdomains using CVF.

ðDxqcp Þm

d km 2km km T m;j ¼ T m;jþ1  T m;j þ T m;j1 dt Dxm Dxm Dxm ¼ F m;j

m ¼ 1; 2;

two systems is through the interface temperature T I and heat flux qI .

ð24Þ

where F m;j is the cell-residual obtained as a result of discretizing the governing equation in space. In X2 , the discretization of the cell adjacent to the interface, T 2;1 , is given by (noting that one-sided differencing is used to approximate the temperature gradient in the interface heat flux  2  q2;I ¼ ðk2 @T Þ ): @x I

ðDxqcp Þ2

d k2 k2;I T 2;1 ¼ ðT 2;2  T 2;1 Þ  ðT 2;1  T 2;I Þ dt Dx2 Dx2;b ¼ F 2;1 ;

ð25Þ

where T 2;I is the interface temperature on the side of X2 and Dx Dx2;b ¼ 22;1 is the distance from the cell center of T 2;1 to the interface. k2;I is the interface conductivity on the side of X2 and unless otherwise mentioned is equal to constant value of k2 prescribed to the complete field (X2 ). Rewriting (25), we obtain:

  d k2 k2 k2;I T 2;1 þ f21 ¼ F 2;1 ; ðDxqcp Þ2 T 2;1 ¼ T 2;2  þ dt Dx2 Dx2 Dx2;b

ð26Þ

k

where f21 ¼ Dx2;I T 2;I is the contribution from the Dirichlet boundary 2;b (interface) condition. The value of the interface temperature is obtained by prescribing the temperature of the interface node in X1 as its value T 2;I ¼ T 1;I . However, it is noted that X1 does not have a node at the interface. Therefore, in order to solve the coupled system an approximation to the interface temperature is required. Here, linear extrapolation from cells in X1 close to the interface is used to approximate the interface temperature:

T 1;I ¼ lT 1;2 þ nT 1;1 ;

ð27Þ

where l ¼ 1=2 and n ¼ 3=2 for a uniform grid in X1 . In X1 , the discretization of the cell adjacent to the interface, T 1;1 , is:

ðDxqcp Þ1

d k1 T 1;1 ¼ q1;I  ðT 1;1  T 1;2 Þ ¼ F 1;1 ; dt D x1

ð28Þ

where q1;I is the interface heat flux on the side of X1 and based on the Dirichlet-Neumann formulation, the interface heat flux of X2 is prescribed as its value:

q1;I ¼ q2;I ¼ 

k2;I ðT 2;1  T 2;I Þ: Dx2;b

ð29Þ

The semi-discrete form of the coupled problem can be expressed by the following two coupled ODE systems:

M1

d T 1 ¼ F 1 ðT 1 ; qI ; tÞ; dt

ð30Þ

M2

d T 2 ¼ F 2 ðT 2 ; T I ; tÞ; dt

ð31Þ

where T 1 and T 2 are vectors containing the discrete solution in the FV cells, M1 and M2 are diagonal matrices containing the product of the cells’ volumetric heat capacities (ðqcp Þm ) with the cell volumes, and F 1 and F 2 are the residual vectors. The coupling between the

5. Time integration In this paper, for each subdomain, the ESDIRK schemes are considered for time integration which can be made of arbitrary high order while retaining the L-stability property. While BDF1 and BDF2 are L-stable, the Crank–Nicolson scheme is only A-stable. It is possible to construct higher order BDF methods, but they are only LðaÞ stable and their stability region drops as the order of the scheme is increased (see [19]). For an ODE system of the form M dT ¼ F ðT; tÞ, the solution at each stage of the ESDIRK scheme can dt be written as:

MT ðkÞ ¼ MT n þ Dt

k X aIki F ðtn þ ci Dt; T ðiÞ Þ i¼1

k X ¼ MT n þ Dt aIki F ðiÞ ;

ð32Þ

i¼1

P where aIki are the stage weights and ci ¼ j aIij are the quadrature ðiÞ n nodes of the scheme t ¼ t þ ci Dt. While the stages can be of low order, high order solution at the next time level can be achieved by the weighted sum of the stage residuals such that the lower order errors cancel out:

T nþ1 ¼ T n þ M1 Dt

s X bi F ðiÞ ;

ð33Þ

i¼1

P where bi are the scheme’s main weights with i bi ¼ 1, and s is the number of stages. In this paper, stiffly accurate ESDIRK schemes are considered where aIsi ¼ bi and thus the solution of the last stage is equal to the solution of the next time-level, T nþ1 ¼ T ðsÞ . Therefore, when a fully coupled approach is used to solve the coupled problem (performing interface iterations until convergence to the monolithic solution), computing (33) becomes unnecessary. The fact that the stages can be of lower order compared to the design order of accuracy, makes the Runge–Kutta schemes susceptible to order reduction in the presence of substantial stiffness (see [20]). The coefficients and weights are usually arranged in a Butcher tableau (see Fig. 2). For the ESDIRK schemes, as the name implies, the diagonal coefficients are equal ðaIkk ¼ cÞ. Furthermore, it is possible to incorporate the ESDIRK schemes into solvers which already include the Backward Euler, since from an implementation view point, each stage of the ESDIRK scheme resembles the BDF1 scheme with a source term. In solving time-accurate advection–diffusion-reaction problems, where separable stiff and non-stiff components are identifiable, Kennedy and Carpenter [15] present a family of high order

Fig. 2. A four stage additive RK method consisting of an (a) ESDIRK and (b) ERK scheme.

209

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217

IMEX (additive) Runge–Kutta schemes where the ESDIRK schemes are used for integrating the stiff components of the problem (to retain stability) and ERK schemes (which are computationally less expensive than the corresponding implicit schemes) are used for integrating the non-stiff components of the problem. In order to be consistent with the ESDIRK schemes, the introduced ERK schemes have the same order and number of stages as the ESDIRK schemes (they have the same weight factors bi and ci coefficients). For the time-step solution to have the designed high order accuracy, evaluating (33) is necessary. In this paper, the coefficients of the implicit ESDIRK schemes are denoted by aIki and that of the ERK schemes by aEki as shown in Fig. 2. Analogously, in solving unsteady CHT problems using the loosely-coupled solution algorithm, the coupling terms between the fluid and solid subdomains are assumed as separable components which can be integrated with the explicit ERK schemes; the fluid and solid systems are integrated with the implicit ESDIRK schemes. This approach follows from [7] where the mechanical coupling of fluids and structures is considered. Furthermore, in the present study, similar to [7], the high order IMEX schemes developed in [15] are used for time integration. This class of IMEX schemes has also been used in [5] to overcome geometry induced stiffness in a single (fluid) domain; the non-stiff portion of the domain is integrated using ERK schemes and the stiff portions using the ESDIRK schemes. By expressing the predictor–corrector h scheme, used in [13] for solving unsteady CHT problems, in terms of a two-stage IMEX scheme (its Butcher tableau is given in Fig. 3), similarities can be observed between the h scheme with h ¼ 0:5 (Crank–Nicolson) and the second stage of the high order IMEX schemes. For the second stage of the IMEX schemes, the coefficients of the corresponding ESDIRK and ERK schemes are respectively given by aI21 ¼ c; aI22 ¼ c, and aE21 ¼ 2c and those of the (predictor– corrector) Crank–Nicolson scheme are aI21 ¼ 0:5; aI22 ¼ 0:5; aE21 ¼ 1. 6. Solution algorithm

ðkÞ

ð34Þ

i¼1

where A2 represents the discretization of the diffusion term and f 2 is a vector containing the contribution of the non-interface boundary. The vector b2 contains zero entries except for the discretization of the cell next to the interface, T 2;1 , for which its entry is 1. Integrating the coupling term using the ERK schemes, is equivalent to using the following flux predictor: ðÞ

f21 ¼

k1 E X aki  aIki ðiÞ f21 ; aIkk i¼1

Fig. 3. Butcher tableau for the predictor-corrector h scheme.

ð35Þ

ð36Þ

i¼1 ðiÞ

where F 2 is the residual vector of any of the stages defined by: ðiÞ

ðiÞ

ðiÞ

ðiÞ

F 2 ¼ A2 T 2 þ f21 b2 þ f 2 :

ð37Þ

The equivalence between (34) and (36) can be readily observed by substituting (35) into (36). When the thermal conductivity and the mesh is fixed in time, k k the coefficient Dx2;I in f21 ¼ Dx2;I T 2;I is constant and therefore the 2;b 2;b flux-predictor in (35) is reduced to predicting the interface temperature (i.e. state-predictor): ðÞ

T 2;I ¼

k1 E X aki  aIki ðiÞ T 2;I ; aIkk i¼1

ð38Þ

ðkÞ

With T 2 computed, we assign: ðÞ

ðÞ

q1;I ¼ q2;I ¼ 

k2;I ðkÞ ðÞ ðT  T 2;I Þ Dx2;b 2;1

ð39Þ

ðkÞ

and solve for T 1 using the ESDIRK schemes: k1   X ðkÞ ðkÞ ðÞ ðkÞ ðiÞ M1 ðT 1  T n1 Þ ¼ DtaIkk A1 T 1  q1;I b1 þ f 1 þ Dt aIki F 1 ;

ð40Þ

i¼1

where A1 represents the discretization of the diffusion term and f 1 is a vector containing the contribution of the non-interface boundary. The vector b1 contains zero entries except for the discretization ðiÞ of the cell next to the interface for which its entry is 1. F 1 is the residual vector of any of the stages: ðiÞ

ðiÞ

ðiÞ

F 1 ¼ A1 T 1  q1;I b1 þ f

ðiÞ 1 :

ð41Þ

With the stage temperature field T ðkÞ m computed, the interface conditions are updated: ðkÞ

k k1     X X ðiÞ ðiÞ ðiÞ aIki A2 T 2 þ f 2 þ Dt aEki f21 b2 ; i¼1

k1   X ðkÞ ðkÞ ðÞ ðkÞ ðiÞ M2 ðT 2  T n2 Þ ¼ DtaIkk A2 T 2 þ f21 b2 þ f 2 þ Dt aIki F 2 ;

ðkÞ

ðkÞ

T 1;I ¼ lT 1;2 þ nT 1;1 ;

In this section, the details of the loosely-coupled solution algorithm is presented where the IMEX schemes (presented in Section 5) are used for time integration. We consider solving the coupled system at each stage integrating X2 first, i.e. using Block Gauss–Seidel (BGS-21). In computing the temperature field in X2 at an implicit stage of the ESDIRK schemes (the semi-discretized form of T 2 is given in k (24)), the coupling term f21 ¼ Dx2;I T 2;I is treated explicitly and inte2;b grated with the ERK schemes. In other words, the semi-discrete form of T 2 is integrated in time using the IMEX schemes:

M 2 ðT 2  T n2 Þ ¼ Dt

ðkÞ

in computing the stage solution T 2 :

ðkÞ

ð42Þ

ðkÞ

T 2;I ¼ T 1;I ; ðkÞ

q2;I ¼  ðkÞ

ð43Þ

k2;I ðkÞ ðkÞ ðT  T 2;I Þ; Dx2;b 2;1

ð44Þ

ðkÞ

q1;I ¼ q2;I

ð45Þ ðkÞ

ðkÞ

and the stage residual vectors F 1 and F 2 are computed using (41) and (37) respectively. After solving all s stages, the solution to the time-level T nþ1 is m obtained by using (33) where through cancelation of the lower order error terms, T nþ1 has the designed high order accuracy. m The solution procedure is the same when the coupled problem is solved by integrating X1 first (BGS-12). In solving for the temperature field in X1 , a prediction of the flux at the interface is required which is given by (38), with T 2;I replaced with q1;I . After solving for T 1 , the interface temperature is evaluated using (27) and is imposed as the interface condition on X2 . 6.1. Stability of the solution algorithm Because of partitioning, there is the possibility of numerical instability not inherent in the monolithic approach [3]. By applying the stability theory of Godunov and Ryabenkii to the 1-D model problem of Section 3 (with the modification that the two coupled domains are considered to be semi-infinite), [3,13] analyzed the stabilities of loosely-coupled solution algorithms where respectively BDF1 and the h scheme were used for time integration. By correct assignment of the interface boundary conditions, the solution algorithm was found to be unconditionally stable using

210

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217

BDF1 [3]. For stable computations using the predictor–corrector Crank–Nicolson scheme [13], the criterion rd2 < c needs to be ðqc DxÞ

satisfied, where r ¼ ðqcpp DxÞ2 , and d2 ¼ DDtxa22 (Fourier number of X2 ). c 1

2

is equal 1 for the vertex-based discretization and 2=3 for the cell-center discretization with linear extrapolation of the interface temperature. By substituting the definitions of r and d2 into the criterion, the stability limit can be written as,

rd2 < c )

Dt ðqkc2p Þ

1

Dx1 Dx2

< c;

ð46Þ

which imposes restriction on Dt given the material properties and grid size of the coupled domains. The stability criterion can also be expressed in the following useful form:

pffiffiffiffiffiffiffiffiffiffi

r d1 d2 < c;

ð47Þ

which directly relates the strength of the coupling r to the stability of the algorithm. Applying the same stability theory to the IMEX schemes is not straightforward, given the multiple stages (Butcher tableau coefficients). However, it is possible to identify parameters that may influence the stability of the coupling algorithm based on the fully-discretized equations. For the purpose of analysis, it is more convenient to cast the system of difference equations (for a stage) of the loosely-coupled IMEX into the form shown in Appendix A. The equations, correspond to integrating X1 first (BGS-12). It is noted that in the discretization of the interface node T 1;0 , (A.2), pffiffiffiffiffiffiffiffiffiffi the product r d1 d2 is part of the coupling term. Furthermore, for the predictor–corrector Crank–Nicolson, as shown earlier, the product also appears in the stability criterion (47). Hence, it can be expected that this product is an influential parameter in the stability of the high order IMEX schemes. This issue will be further investigated in the results section. 7. Numerical examples In this section, the loosely coupled solution algorithm presented in Section 6 is used to solve a conjugate heat transfer problem in order to demonstrate the applicability of the algorithm, to investigate the temporal order preservation and stability of the algorithm, and finally its computational efficiency relative to commonly used second order time integration schemes. Initially results are shown for a case with constant material properties. This is followed by investigating the temporal order when the thermal conductivities are temperature dependent. 7.1. Case of constant material properties 7.1.1. Problem description and specifications As a test-case, unsteady conjugate natural convection in a square enclosure is considered. The configuration of the problem, shown in Fig. 4, follows from Kaminski and Prakash [21] (where the steady-state solution of the coupled problem was sought). The solid domain occupies the rectangular region X1 ¼ ½0; L ½h; 0 and the fluid domain occupies the square region X2 ¼ ½0; L ½0; L with the interface between the two domains located at I ¼ ½0; L 0. The vertical boundaries (at x ¼ 0, and x ¼ L) of both the solid and fluid domains are insulated. Dirichlet conditions are prescribed to the top boundary of the fluid domain (Tðx; LÞ ¼ T C ) and the lower boundary of the solid domain (Tðx; hÞ ¼ T H ). The gravity acts in the x direction (j ¼ ð1; 0; 0Þ). The accuracy and stability of the loosely-coupled algorithm are investigated for the following two cases of different strength in thermal coupling:

Fig. 4. The configuration of the model problem.

Case-a

k1 ¼ 1600;

q1 ¼ 7500; cp;1 ¼ 0:5; k2 ¼ 1; q2 ¼ 1;

m2 ¼ 0:7; b ¼ 0:7 105 ; g ¼ 1; L ¼ 1; h ¼ 0:2; k1 a1 ¼ 1; T H ¼ 2; T C ¼ 1; ¼ 1600; ¼ 0:4; Pr ¼ 0:7; k2 a2

cp;2 ¼ 1; T ref

Ra ¼ 1 105 ;

ð48Þ

Case-b

k1 ¼ 80;

q1 ¼ 7:5; cp;1 ¼ 0:12; k2 ¼ 1; q2 ¼ 1; cp;2 ¼ 1;

m2 ¼ 7; b ¼ 4:9 105 ; g ¼ 1; L ¼ 1; h ¼ 0:2; T ref ¼ 1; k1 a1 T H ¼ 2; T C ¼ 1; ¼ 80; ¼ 88:9; Pr ¼ 7; k2 a2 Ra ¼ 7 104 ;

ð49Þ

where the Ra numbers are evaluated based on DT ¼ T H  T C . As pointed out earlier, the governing non-dimensional parameters of the problem are given by Pr; Ra;

a1 a2

and

k1 . k2

The ratio between the

material properties of the coupled domains in (48) and (49) correspond to those of steel-air and steel-water couplings respectively. Thermal coupling is weak (r ¼ 4:1 104  1) for Case-a, and relatively strong for Case-b (r ¼ 0:1). Based on the prescribed values of the material properties (in both Case-a and Case-b), for stability X2 takes the Dirichlet condition at the interface, while X1 the Neumann condition. 7.1.2. Solution procedure The solid domain X1 is discretized using the finite volume method with the temperature stored in the cell center. The fluid domain X2 is also discretized using the finite volume method with staggered arrangement of the variables [22]. In each of the subdomains, a uniform grid-spacing is used in both directions: for the solid domain: N 1;x ¼ 80; N 1;y ¼ 20 (with N: number of cells), and for the fluid domain N 2;x ¼ 80; N 2;y ¼ 80. Linear extrapolation from the cells in X1 close to the interface, (27), is used to approximate the interface temperature. In both the momentum and energy equations, the second order centered scheme is used for discretizing the convective terms and the second order differencing scheme is used to discretize the diffusive terms. In the flow solver, at each stage, Picard iterations are used for solving the non-linear coupled equations [22]. The flow solver is iterated to a strict tolerance (iter < 109 ; iter : maximum of the residuals of the momentum

211

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217

and energy equations) to eliminate iteration error as a contaminating variable in the study. For both cases a and b, computations were carried out to steady-state from an initial state of unit temperature imposed on the global domain (Tðx; t ¼ 0Þ ¼ 1 for x 2 X ¼ X1 [ X2 ). The purpose was to obtain information regarding the time-history of the solution during the transient phase. The time-history of the temperature at the interface and cavity centers, and x-component of the velocity at the cavity center are shown in Fig. 5. The solution reaches steady-state at t  0:3 for Case-a and t  0:2 for Case-b. The monolithic approach of Schäfer and Teschauer [23] is used for the computations (a single temperature equation is solved across both the subdomains while using a separate flow field solver in the fluid subdomain). The fourth order ESDIRK (ESDIRK4) with time step-size of Dt ¼ 5 104 was used for time integration. 7.1.3. Temporal accuracy In order to perform the accuracy analysis, computations are first performed with the monolithic solver to time level tIC from an initial state of unit temperature imposed on global domain. This is done to step over the time that is required for the heat to travel from the bottom boundary to the interface (thermal interaction between the two subdomains has been initiated at this time). The fifth order ESDIRK (ESDIRK5) with time step-size of DtIC is used for time integration. For Case-a, t IC ¼ 0:01 and DtIC ¼ 1 104 , and for Case-b, t IC ¼ 1 104 and Dt IC ¼ 2 105 . The obtained solution at t IC is used as initial condition for the problem. The coupled problem is advanced in time to tfinal with the loosely-coupled algorithm and by using various time integration schemes (for Case-a, t final ¼ 0:06, and for Case-b, tfinal ¼ 0:02). As Fig. 5 demonstrates, the solution at tfinal þ t IC is within the transient phase. For each case, the partitioned solution at t final was compared to that of the temporally exact solution (defined as the monolithic solution at t final obtained using a fine time step: ESDIRK5 is used for time integration with Dtfine ¼ 5 104 for Case-a, and Dtfine ¼ 2 105 for Case-b). To investigate how well the coupled problem is solved compared to the monolithic approach, computations are also carried out with the monolithic solver. The time integration error (t ) is defined as the difference between the temporally exact solution and the monolithic solution at tfinal obtained with a certain time step which is coarser relative to Dtfine . The difference between the temporally exact solution and the partitioned solution at t final is defined as the total error in the partitioned solution (total ) and represents the sum of the time integration error (t ) and partitioning error (), (total ¼ t þ ). For total to be of the same order as t , the partitioning error  needs to be at least of the same order as t . In Fig. 6, for both Case-a and Case-b, the total error in the global temperature field obtained with

the loosely-coupled algorithm (integrating X1 first, BGS-12) is shown as a function of time step size. Also shown, is the corresponding (time integration) error in temperature field obtained using the monolithic approach. The design orders of the IMEX schemes are clearly observable for Case-a, as the results in Fig. 6(a) demonstrate. Therefore, there is no need to perform additional interface iterations at each stage to retrieve the time integration’s design order. In retaining the design order, as Fig. 7(a) and (b) demonstrate, the updated interface conditions must be used in evaluating the stage residual functions, F ðkÞ m , otherwise order reduction occurs (in the figures, compare third and fourth order IMEX (3-4),BGS-12 with IMEX (3-4)-OR,BGS-12). For a time-accurate partitioned solution, the partitioning error () should be as low as possible, desirably below the time integration error such that the partitioning error is not the dominant source of error in the partitioned solution (hence total  t ). Otherwise the temporal accuracy of the partitioned solution is reduced. Below we investigate the possible influence that the order with which the domains are integrated, and the strength of the thermal coupling may have on accuracy. In Fig. 7, by comparing the total error in the temperature field obtained with the two different integration sequences with the corresponding (time integration) error in the monolithic solution, it is observed that the partitioned solution obtained by integrating X1 first (BGS-12) is more accurate than that of BGS-21 and for this test-case is as accurate as the monolithic solution. In solving the coupled problem using BGS-21, the partitioning error (as a result of segregated solve of the interface equations) is placed entirely in the continuity of the interface temperature equation and therefore this equation is solved approximately while the continuity of the interface heat flux is satisfied exactly. This is in opposite to solving the problem using BGS-12. In general it is not clear which approximation (having a discontinuity in the interface heat flux or interface temperature) is less disadvantageous and seems to be an issue which is problem dependent. At least for this test case, as the results in Fig. 7 demonstrate, using BGS-12 provides a solution which is more accurate than BGS-21. Furthermore, the temporal convergence character of the total error of the solution using BGS-21 is also influenced by the partitioning error; to observe the full-design orders smaller time-step sizes need to be considered. It is noted that by solving the test-case on coarser grids of the fluid and solid subdomains (hence smaller Fourier numbers DDtxa2m relative to the current mesh) full-design m orders has been observed. By comparing the results of cases a and b in Fig. 6(a) and (b), it can be inferred that as the strength of the thermal coupling increases (r ! 1), so does the partitioning error. As a result, the accuracy of the partitioned solution accuracy deteriorates (this

2

0

1.9

−5

1.8 −10

1.6

−15

1.5

−20

u

T

1.7

1.4

−25

1.3

Case−a, Cavity center Case−b, Cavity center Case−a, Interface center Case−b, Interface center

1.2 1.1 1 0

0.05

0.1

0.15 t

0.2

0.25

0.3

−30 Case−a Case−b

−35 −40 0

0.05

0.1

0.15 t

0.2

0.25

0.3

Fig. 5. Time-history of (a) the temperature at the interface and cavity centers and (b) x-component of the velocity at the center of the square enclosure.

212

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217

−2

−1

−3

CN ESDIRK3 ESDIRK4 ESDIRK5 pred−corr CN IMEX3 IMEX4 IMEX5

−2 2

O(Δ t )

−3 || || ε

−6 −7

CN ESDIRK3 ESDIRK4 ESDIRK5 pred−corr CN IMEX3 IMEX4 IMEX5

4

O(Δ t )

−8 5

−9 −3.4

−4

T total

3

O(Δ t )

O(Δ t )

−3.2

−3

−2.8 −2.6 log Δ t

−2.4

−2.2

10

−5

log

log

10

|| ε

T total

||

2

2

−4

−5 2

−6

O(Δ t )

−7 −8

3

O(Δ t )

4

O(Δ t ) 5

O(Δ t )

−9 −5

−2

−4.5

−4 log

−3.5 Δt

−3

−2.5

10

10

Fig. 6. Total error in the global temperature field (Ttotal ) as a function of Dt for (a) Case-a, (b) Case-b. In the figures, CN refers to monolithic solution obtained with Crank– Nicolson, ESDIRK 3  5 refer to the monolithic solutions obtained using the 3rd -5th order ESDIRK schemes. The solution obtained using the predictor–corrector Crank– Nicolson scheme of [13] is denoted by pred-corr CN, and the solutions obtained using the 3rd–5th order IMEX schemes are denoted by IMEX 3-5.

−2

0

−3 −4 log10 || ε Ttotal ||

−3 −4

log

10

|| ε T

total

||

2

−2

ESDIRK3 IMEX3, BGS−12 IMEX3, BGS−21 IMEX3−OR, BGS−12

−5

O(Δ t3)

−5 −6

O(Δ t4)

−7

−6

−8

−7 −8 −5

ESDIRK4 IMEX4, BGS−12 IMEX4, BGS−21 IMEX4−OR, BGS−12

2

−1

−4.5

−4

−3.5 log Δ t

−3

−2.5

−2

10

−9 −4.5

−4

−3.5

log10 Δ t

−3

−2.5

−2

Fig. 7. The total error in the global temperature field as a function of time-step size for Case-a. In the figures, the partitioned solutions obtained by integrating X1 first are indicated by BGS-12, while those obtained by integrating X2 first are indicated by BGS-21. IMEX (3-4)-OR,BGS-12 refer to the partitioned solutions obtained using IMEX (3-4) for time integration where the updated interface conditions are not used in the evaluation of the stage residual functions, and hence order reduction occurs.

holds irrespective of the sequence with which domains are integrated). Furthermore, it is observed that as the order of the time integration scheme increases, the partitioning error increases. It is noted that similar to results of Case-a obtained using BGS-21, to observe the full-design orders smaller time-step sizes need to be considered. 7.1.4. Stability By neglecting flow in the fluid domain X2 , and noting that the left and right boundaries of the problem are insulated (Fig. 4), the coupled problem reduces to the 1-D model problem of Section 3 where the governing equation within each subdomain is conduction. Following the discussion in Section 6.1, first the stability of the IMEX schemes is investigated numerically for the 1-D model problem. We investigate whether similar to the stability of Crank–Nicolson (47), a relation between r and stability of the IMEX schemes can be observed; of particular interest is a condition which will indicate the possibility of using large Fourier number of subdomains. The results of the 1-D model problem will also serve as reference for the 2-D CHT test-case. 7.1.4.1. Stability of the 1-D model problem. For the analysis, three different values of r are considered: 4:1 104 (Case-a), 0:1 (Case-b), and 1. For case of r ¼ 1, material properties of both

subdomains are set to 1. A step initial condition is imposed on the global domain: Tðx; t ¼ 0Þ ¼ 2; x 2 X1 and Tðx; t ¼ 0Þ ¼ 1; x 2 X2 . The non-interface boundary conditions are kept unchanged and set to T 1 ðy ¼ L1 ; tÞ ¼ 2 and T 2 ðy ¼ L2 ; tÞ ¼ 1. To minimize the influence of the non-interface boundary conditions on the analysis, the length of the subdomains in the y-direction (in Fig. 4) are zzextended to sufficiently large values; the thermal penetration pffiffiffiffiffiffiffiffi depth dm am t within each subdomain at the time at which instability occurs is far away from the boundaries. The number of cells in both domains were accordingly increased to keep the grid size in the y-direction unchanged compared to the CHT case. For each case, Dt is incremented until the computations become unstable. The approximate Dt and the corresponding Fourier numbers of the subdomains at which simulations become unstable are presented in Table 1(a)–(c). Computations using both the integration sequences (BGS-12, BGS-21) provided similar results. For the predictor–corrector Crank–Nicolson scheme, substituting the model parameters of Case-a and Case-b, and case of r ¼ 1 into (46) (with c ¼ 23), the expected Dt values at which instabilities initiate are respectively Dt ¼ 0:31; Dt ¼ 7:4 105 , and 8:33 105 . As the results in Table 1(a)–(c) show, for the three cases instability initiates when the criterion (46) is not satisfied. This also verifies the numerical procedure used to investigate the stability of the algorithm.

213

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217 Table 1 Approximate Dt and corresponding Fourier numbers dm at which simulations become unstable for (a) Case-a, (b) Case-b, and (c) Case r ¼ 1.

Dt (a) r ¼ 4:1 104 Pred-Corr CN

3:2 10

1:4 10

IMEX-3

5:9 102

2:5 106

3:8 106

IMEX-4

1:4 103

6:0 106

9:0 106

IMEX-5

1:3 102

5:5 105

8:3 105

(b) r ¼ 1 101 Pred-Corr CN

7:5 105

6:7 101

4:8 101

IMEX-3

6:7 104

6:0 102

4:3 100

IMEX-4

2:8 103

2:5 103

1:8 101

IMEX-5

3:1 104

2:7 102

1:9 100

(c) r ¼ 1 Pred-Corr CN IMEX-3 IMEX-4 IMEX-5

8:34 105 4

3:2 10

4

1:8 10

4

2:1 10

r

Dtref

ðd1 d2 Þref

d2

d1 1

Table 2 Influence of grid spacing on stability.

3

8:3 101 3:2 10

0

1:8 10

0

2:1 10

0

2:0 10

3

5:3 101 2:0 10

Dy1;C ¼ 2Dy1

Dy2;C ¼ 2Dy2

Dtref Dt C

ðd1 d2 Þref ðd1 d2 Þ

Dtref Dt C

ðd1 d2 Þref ðd1 d2 Þ

5:8 104

1:4 103

5:3 1013

2:1

1:1

2:1

1:1

1:3 101

2:8 103

4:5 104

2:4

1:4

1:8

0:8

1 100

1:8 104

2:1 100

2

1

2

1

Table 3 Approximate Dt and the corresponding Fourier numbers of the subdomains at which simulations become unstable for Case-b and finite length of subdomains.

Dt

r ¼ 1 101 IMEX-3

d2

d1 4

7:1 10

2

4:5 100

3

6:4 10

0

IMEX-4

2:6 10

2:3 10

1:7 101

0

IMEX-5

3:1 104

2:8 102

2:0 100

1:2 10

3

0

1:3 10

For the high order IMEX schemes, as the results in Table 1(a)–(c) demonstrate, the loosely-coupled algorithm becomes unstable for sufficiently large values of the Fourier numbers of the coupled domains. A similar observation is made for the predictor–corrector Crank–Nicolson. In addition, the IMEX schemes demonstrate better stability properties compared to the Crank–Nicolson scheme for all the considered r values; while they outperform Crank–Nicolson for case of r  1, they are marginally better for r ¼ 1. For Crank–Nicolson, the results in the tables can be explained pffiffiffiffiffiffiffiffiffiffi using the available stability criterion (r d1 d2 < 23). When r  1, pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi it follows that d1 d2 can take large values ( d1 d2  1) while still adhering to the criterion and vice versa. Analyzing the results in the tables, suggest a similar trend for the IMEX schemes. Comparing the stability results of Case-a with the corresponding ones for Case-b and case of r ¼ 1, demonstrates that when the thermal coupffiffiffiffiffiffiffiffiffiffi pling between the subdomains is weak (r  1), the product d1 d2 is large and the IMEX schemes remain stable to large Fourier numpffiffiffiffiffiffiffiffiffiffi bers. However, as r ! 1 (Case-b, r ¼ 1), d1 d2 reduces significantly, and the Fourier numbers of the domains are small. Hence, pffiffiffiffiffiffiffiffiffiffi the results imply that the product r d1 d2 , which in (A.2) was shown to be part of the coupling term in the discretization of the interface node, has an influential effect on the stability of the IMEX schemes. More importantly, r can be identified as the single parameter which provides the condition (r  1) for using large Fourier numbers (the reason for using implicit time integration). For the Crank–Nicolson scheme, based on the stability criterion (46), the Dt for stability depends linearly on the mesh spacings Dxm . Alternatively, from (47), it is observed that a change in the mesh spacing results in a change in Dt such that the product of the Fourier  2numbers of the subdomains remains constant, i.e. d1 d2 < 2=3 . A stability investigation was performed using r IMEX4, to see whether similar observations can be made for the high order IMEX schemes; the grid in one domain is increased by a factor of two while retaining the grid size fixed in the other. As the results in Table 2 show, by increasing the mesh spacing by a factor of two, the approximate Dt at which instability initiates (denoted in the table by DtC ) reduces by approximately the same factor relative to Dtref . Furthermore, the product of the Fourier numbers of the coupled domains is approximately the same as the reference ðd1 d2 Þref . The stability results in Table 1 and Table 2 were obtained assuming semi-infinite subdomains. Computations were also carried out for cases a and b to investigate the extend to which the stability results of the IMEX schemes may be influenced by the

non-interface boundary conditions placed at a finite distance from the interface (the top and bottom boundaries (at y ¼ L and y ¼ h) in Fig. 4). The same length of the subdomains Lm;y and same number of cells N m;y in the y-direction of the CHT problem were used for this 1-D test case. For Case-b, the results are shown in Table 3. A maximum difference of 6% with the results of the semi-infinite subdomains in Table 1 is observed. The difference was found to be larger for Case-a; as an example, for IMEX-4 instability initiates at a Dt approximately 3 times larger than the corresponding result in Table 1(b). The governing non-dimensional parameters of the 1-D model problem, due to lack of convection, are aa12 and kk12 . While the stability of the high order IMEX schemes is influenced by the two ratios, however, for the purpose of our analysis, it is more appropriate  qffiffiffiffi to use r ¼ kk21 aa12 . For a given strength of thermal coupling, we are interested in whether it is possible to take large Fourier numbers as desired using the loosely coupled IMEX schemes. r (an indicator of the strength of the coupling), as demonstrated in this section, does serve the purpose and provides an estimate of cases (r  1) where the high order IMEX schemes should be applied to. 7.1.4.2. Stability of the CHT problem. Returning to the CHT test-case, the problem is two dimensional. Furthermore, the fluid domain includes more complex physics, i.e. coupled-nonlinear equations. While the governing non-dimensional parameters of 1-D model problem are

a1 a2

and

k1 , k2

the CHT test-case includes Pr and Ra as

two additional parameters. Hence, for 2-D (and 3-D) computations, satisfying the stability criteria for the 1-D model problem is necessary, but may not be a sufficient one [3]. Therefore, stability of the coupling algorithm in solving the CHT problem was also investigated numerically, in a similar manner as the 1-D problem. For Case-a, numerical computations using the high order IMEX schemes with Dt values sufficient to capture the transients were performed without encountering stability issues (maximum considered Dt was 1 102 ). For Case-b, a noticeable difference with the stability results of the 1-D model problem in Table 3 was not observed. Similar observations were made when the Pr and Ra numbers were separately reduced by a factor of 10 while all the other parameters were fixed. Therefore, at least for this particular test-case, the results of the CHT computations suggest a marginal influence of Pr; Ra, and other possible contributing factors such as two dimensionality on stability. 7.1.4.3. Intermediate conclusions. The main conclusion which should be drawn from the stability investigation on the simplified

214

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217

1-D and the 2-D CHT test-cases is the possibility to use the single variable r, to identify a class of problems which the loosely coupled IMEX schemes are appropriate to apply to. This is a particularly useful result, given the absence of an analytical stability criterion. Based on the analysis, it may be concluded that when the thermal coupling between the subdomains is weak (r  1), the IMEX schemes remain stable to large Fourier numbers dm . However, as r approaches unity, the stability and accuracy of the schemes reduce. The possibility to use large Fourier numbers for thermally weak couplings, is an indication that time-step size is restricted by accuracy rather than stability, at least on reasonable grids. This is demonstrated for the CHT test-case under study: since the transients to steady state for cases a and b are completed within  t ¼ 0:3 (see Fig. 5), taking time-step sizes in the order of Dt ¼ 101 can be considered large and correspond to a time which most of the transients have passed. For Case-a (r  1), computations with the largest considered time-step size Dt ¼ 1 102 sufficient to capture the transients were stable. However, for case-b (r ¼ 0:1), the time-step sizes at which instability occur, Table 3, are adequate to capture the transient phase.

per time-step. However, for the same time-step, the high order IMEX schemes generally provide temporally more accurate solutions. For the CHT test-case considered here, we investigate whether the additional work per time step of the IMEX schemes is compensated by the gain in accuracy compared to Crank–Nicolson and BDF2 schemes. Based on the results of the accuracy and stability investigation for the IMEX schemes discussed above, only Case-a where the thermal coupling between the two domains is weak (r  1) is considered. Prior to discussing the results, as an illustration of the quality of the solution at different error levels, the variation with time of the temperature and x-component of the velocity at the center of the cavity is demonstrated in Fig. 8. The partitioned solutions are obtained using the predictor–corrector Crank–Nicolson and as reference, the data from the temporally exact solution are also included. To capture the transients with reasonable accuracy using the Crank–Nicolson scheme, a time-step size in the range of 1 103 6 Dt 6 2 103 should be considered (as Fig. 8(b) shows). The total error in the temperature field and x-component of the velocity field are plotted against the work in respectively Fig. 9(a) and (b) for several time integration schemes. The work is defined as the total number of implicit calculations during the simulation

7.1.5. Computational efficiency assessment While in the loosely coupled multi-stage IMEX schemes a single interface iteration is performed at each (implicit) stage, in the predictor–corrector Crank–Nicolson scheme, only one is performed

t

(W ¼ ðs  1Þ final , with s  1 being the number of implicit stages Dt per time-step). The temporal accuracies of the velocity and temperature fields using a time-step size in the range of

1.18 1.16

0 −4

ESDIRK5: Δ t = 5 ⋅ 10 −3

1.14 1.12

CN : Δ t = 1 ⋅ 10

−5

−3

CN : Δ t = 2 ⋅ 10

−3

CN : Δ t = 5 ⋅ 10

−2

CN : Δ t = 1 ⋅ 10

−10 u

T

1.1 1.08

−15

1.06

ESDIRK5: Δ t = 5 ⋅ 10−4 −3

CN : Δ t = 1 ⋅ 10

1.04

CN : Δ t = 2 ⋅ 10−3

−20

CN : Δ t = 5 ⋅ 10−3

1.02

CN : Δ t = 1 ⋅ 10−2

1 0

0.01

0.02

0.03 t

0.04

0.05

0.06

0

0.01

0.02

0.03 t

0.04

0.05

0.06

Fig. 8. Time-history of (a) the temperature and (b) x-component of the velocity at the center of the square enclosure. The results are obtained using the predictor–corrector Crank–Nicolson scheme. As reference, the results of the temporally exact solution obtained with ESDIRK-5 are also included.

−2 −1 −3 −2 2 u || ε total ||

−6 −7 −8

−3 −4

10

−5

log

log

10

T || ε total ||

2

−4

pred−corr CN IMEX3 IMEX4 IMEX5 Mono BDF2

−5

pred−corr CN IMEX3 IMEX4 IMEX5 Mono BDF2

−6 −7

0.5

1

1.5 2 log10 Work

2.5

3

1

1.5

2 log

10

2.5

Work

Fig. 9. Work-accuracy character of several time integration schemes for obtaining the solution to: (a) the temperature field, (b) x-component of the velocity.

215

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217

1 103 6 Dt 6 2 103 are approximately in the range of 3

2

4

3:5

1 10 6 utotal 6 1 10 and 1 10 6 Ttotal 6 1 10 respectively. Comparing Fig. 8(a) and (b), shows that for a given time step size, the error in capturing the transients of the solution is higher in the velocity field. Hence, we continue the discussion only considering results in Fig. 9(b). As Fig. 9(b) demonstrates, for reasonable accuracies of

ðkÞ

Furthermore, in solving for T 2 , an approximation to the interface conductivity, k2;I , which appears in the coefficient matrix A2 (as a result of spatial discretization of the cells next to the interface (26)) needs to be considered. One option is to make use of the available interface conductivities evaluated at the previous stages. Here, an expression similar to state predictor given by (38) is considered ðÞ for computing k2;I :

1 103 6 utotal 6 1 102 , both IMEX4 and IMEX5 reduce the computational work relative to BDF2 by approximately a factor of 1:9  2:7. IMEX-4 and IMEX-5 are as efficient as Crank–Nicolson

k2;I ¼

for utotal  102 and for an accuracy of utotal  103 are more efficient than Crank–Nicolson by a factor of 1:5. As the figure demon-

After the computation of T 2 , the interface flux is evaluated using:

strates, for high-accuracy solutions, utotal < 104 , all the three IMEX schemes are computationally more efficient than both BDF2 and Crank–Nicolson. The A-stable Crank–Nicolson scheme performs well compared to the L-stable time-integration schemes, in particular BDF2. Decline in the efficiency of the Crank–Nicolson may be expected when the solution exhibits large gradients and/or a convection dominated flow is present in the fluid domain. Furthermore, in the present study, a fixed time step was used for each simulation. For a fixed temporal accuracy however, the time step usually varies during a typical simulation [6]. Therefore, a fixed time step may not be the optimal choice. The Runge–Kutta schemes include embedded schemes which allow for the use of automatic error-based time-step controllers; for more details and examples, see [9] (a CHT problem), and [5,6,15,20] (CDR and fluid flow problems using the Runge–Kutta schemes of [15] considered in this paper).

ðÞ

k1 E X aki  aIki ðiÞ k2;I ðT 2;I Þ: aIkk i¼1

ð52Þ ðkÞ

ðÞ

1 ðÞ ðkÞ ðÞ ðk T  f21 Þ Dx2;b 2;I 2;1

ðÞ

q1;I ¼ q2;I ¼  ¼

ð53Þ

ðiÞ k1 E X aki  aIki k2;I ðkÞ ðiÞ ðT 2;1  T 2;I Þ: I D x a 2;b kk i¼1

ð54Þ

The rest of the algorithm remains unchanged. An alternative but not equivalent approach is to use the state ðÞ

ðÞ

predictor in the evaluation of f21 . By expressing f21 ðÞ f21

¼

ðÞ k2;I

ðÞ T , Dxb 2;I

it is observed that approximations to

ðÞ k2;I and

as ðÞ T 2;I

ðÞ

are required. In this approach, first T 2;I is evaluated using the state ðÞ k2;I

predictor (38). is then computed by substituting the obtained value for the interface temperature into corresponding conductivðÞ

ðÞ

ðÞ

ity function (k2;I ¼ k2;I ðT 2;I Þ). The computed value of k2;I is used both in the coupling term and in the appropriate term in the coefðkÞ

ficient matrix A2 . After solving for T 2 , the interface flux is evalu-

7.2. Case of temperature dependent thermal conductivity

ðÞ

ðÞ

ðÞ k2;I

2;b

Now the complexity of the coupled problem is increased by assuming that conductivities of the two subdomains vary with temperature according to the function:

km ðT m Þ ¼ km;0 ð1 þ km;1 ðT m  T ref Þcm;1 þ km;2 ðT m  T ref Þcm;2 Þ;

ð50Þ

where km;0 is the conductivity at the initial state (t ¼ 0), and km;1 ; km;2 ; cm;1 , and cm;2 are constant. We wish to examine the order-preservation of the loosely-coupled solution algorithm under the condition of nonlinear coupling terms. The value of the thermal conductivity at a control volume (or cell) face is evaluated by using the harmonic averaging of the cell-center conductivities (obtained with (50)) adjacent to the face kface ¼

2ki kiþ1 ki þkiþ1

[24]. For the faces located at the boundaries of the sub-

domains at which the temperature condition is imposed, the conductivity is obtained by using (50). When the domain which requires the flux condition at the interface is integrated first, BGS-12, the loosely-coupled algorithm presented in Section 6 can be used without any further considerations in addition to taking into account that k2;I ¼ k2;I ðT 2;I Þ. However, when the domain which requires the temperature at the interface is integrated first, a few issues need to be discussed. In ðkÞ

solving for T 2 , the flux predictor (35) and state (temperature) preðÞ

dictor (38) no longer lead to the same prediction of f21 in (36) due to the nonlinearity of the interface heat conductivity (k2;I ¼ k2;I ðT 2;I Þ). As pointed out in Section 6, for (36) to be equivalent to (34) where the coupling term is integrated using an ExpliðÞ

cit RK, the following prediction of f21 needs to be considered (flux predictor): ðÞ

f21 ¼

k1 E k1 E aki  aIki ðiÞ aki  aIki 1 X 1 X ðiÞ ðiÞ f21 ¼ k2;I ðT 2;I ÞT 2;I : I Dx2;b i¼1 akk Dx2;b i¼1 aIkk

ðkÞ

ðÞ

ðÞ

ated using q1;I ¼ q2;I ¼  Dx ðT 2;1  T 2;I Þ where kI

ð51Þ

ðÞ

and T 2;I are

known. The model parameters of Case-a along with the following values of the parameters in (50) are used in this test case (the conductivities are weak functions of temperature):

k1;0 ¼ 1600; k1;1 ¼ 0:0062; k1;2 ¼ 0; c1;1 ¼ 1:5; c1;2 ¼ 0; T ref ¼ 1; k2;0 ¼ 1; k2;1 ¼ 0:64; k2;2 ¼ 0:085; c2;1 ¼ 1:0; c2;2 ¼ 1:5: ð55Þ The CHT test-case is solved using the loosely-coupled algorithm noting the modifications discussed above. Both the flow and solid solvers are iterated to a specified tolerance (here 109 was used) to resolve the non-linearities present within each subdomain. The convergence of the total error in the global temperature field using the loosely-coupled algorithm with time-step size is investigated in the same manner discussed earlier. The results are shown in Fig. 10. The design order of IMEX-3 is clearly observable for BGS12. Similar to the case of linear conductivities considered earlier, the accuracy of the solution using BGS-21 is deteriorated due to the partitioning error (and smaller time-steps need to be considered to observe the full design order). As a result of the non-linear coupling terms (temperature dependent properties), the two integration sequences may also exhibit different stability behaviors, as demonstrated in [13] (for the Crank–Nicolson scheme). In solving the coupled problem using BGS-21, the approach which makes use of the state predictor, is not equivalent to integrating the complete coupling term using an ERK and therefore, there is the possibility of order reduction. However, as the results in Fig. 10 demonstrate, at least for this test-case order reduction is not observed (IMEX3,BGS-21, State Pred). Furthermore, the partitioned solutions obtained with both the alternatives (flux and state predictors) are approximately of equal accuracy relative to the monolithic solution.

216

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217

to be analyzed. Furthermore, the computational efficiency of the strongly-coupled algorithm (with ESDIRK schemes for time integration) relative to the loosely-coupled IMEX schemes needs to be investigated. These issues are the subject of future work.

−1 −2

−4

Appendix A. Parameters influencing the stability of the coupling for IMEX time schemes

−5 3

O(Δ t )

log

10

|| ε T

total

||

2

−3

ESDIRK3 IMEX3,BGS12 IMEX3,BGS 21, State Pred IMEX3,BGS 21, Flux Pred

−6 −7

3

O(Δ t )

−8 −9 −5

−4.5

−4

−3.5 log Δ t

−3

−2.5

−2

10

Fig. 10. Total error in the global temperature field as a function of Dt for case with nonlinear conductivity of subdomains.

In this section, for the one dimensional model problem described in Section 3, the parameters that may influence the stability of the loosely-coupled algorithm presented in Section 6 are identified. To simplify the analysis typically a vertex based rather than cell-centered spatial discretization is preferred (see [13] for details). By applying the IMEX schemes to the semi-discrete form of the problem, the following system of equations at each implicit stage of the loosely-coupled algorithm is obtained (assuming that X1 is integrated first (BGS-12)): ðkÞ

T 1;j ¼ T n1;j þ aIkk d1 ðT 1;j1  2T 1;j þ T 1;jþ1 ÞðkÞ 8. Conclusions and future work

þ

k1 X aIki d1 ðT 1;j1  2T 1;j þ T 1;jþ1 ÞðiÞ ;

j < 0;

ðA:1Þ

i¼1

In order to reduce the computational expense of solving timeaccurate unsteady CHT problems using a loosely-coupled partitioned approach, we considered the use of high order implicit time integration schemes which have the potential to improve the computational efficiency relative to the commonly used second order implicit schemes. In particular, we presented a loosely-coupled solution algorithm where high order implicit–explicit (IMEX) Runge–Kutta schemes are employed for time integration. The high order IMEX schemes consist of the explicit first-stage singly diagonally implicit Runge–Kutta (ESDIRK) schemes, for advancing the solution in time within each separate fluid and solid subdomain, and the explicit Runge–Kutta (ERK) schemes, for explicit integration of part of the coupling terms. To demonstrate the applicability and the temporal order preservation (without subiterating) of this loosely-coupled algorithm, an unsteady conjugate natural convection in an enclosure was considered as a numerical example. Furthermore, the stability of the algorithm was investigated numerically and compared with the results of the analytical stability analysis available in the literature for the Crank-Nicolson scheme. The analysis demonstrated that the loosely-coupled IMEX algorithm becomes unstable for sufficiently large Fourier numbers. Furthermore, when the ratio of the thermal effusivities of fluid and solid subdomains is much smaller than unity, it is possible to use large Fourier numbers, indicating that time-step size is restricted by accuracy rather than stability. In addition, the results showed better stability properties of the IMEX schemes compared to the Crank–Nicolson scheme. Furthermore, the (computational) work-(temporal) precision character of the high order IMEX and the commonly used second order implicit schemes in solving the CHT test-case was compared over a range of accuracy requirements; for time-accurate solutions, the fourth and fifth order IMEX schemes are 1:5 times more efficient than the Crank–Nicolson scheme and 2:7 times more efficient than the second order backward difference (BDF2) scheme. The computational gain is higher for smaller tolerances. As the thermal coupling between the domains becomes stronger, the stability and the accuracy of the loosely-coupled IMEX schemes reduce. Under such circumstances, it is more appropriate to use a strongly-coupled solution method where subiterations are performed to stabilize and/or increase the accuracy of the partitioned solution. When the high order ESDIRK schemes are used for time integration in a strongly-coupled solution algorithm, the stability and rate of convergence of the interface iterations need

ðkÞ

T 1;0 ¼ T n1;0  2aIkk d1 ðT 1;0  T 1;1 ÞðkÞ 

k1 X 2aIki d1 ðT 1;0  T 1;1 ÞðiÞ i¼1

k1 pffiffiffiffiffiffiffiffiffiffi X þ 2 aEki r d1 d2 ðT 2;1  T 2;0 ÞðiÞ ;

ðA:2Þ

i¼1 ðkÞ

T 2;j ¼ T n2;j þ aIkk d2 ðT 2;j1  2T 2;j þ T 2;jþ1 ÞðkÞ ; þ

k1 X aIki d2 ðT 2;j1  2T 2;j þ T 2;jþ1 ÞðiÞ

j > 0;

ðA:3Þ

i¼1 ðiÞ

ðiÞ

T 2;0 ¼ T 1;0

i ¼ 1; . . . ; k:

ðA:4Þ

References [1] T.J. Heindel, S. Ramadhyani, F.P. Incropera, Conjugate natural convection from an array of protruding heat sources, Numer. Heat Transfer Part A: Appl. 29 (1996) 1–18. [2] X. Chen, P. Han, A note on the solution of conjugate heat transfer problems using simple-like algorithms, Int. J. Heat Fluid Flow 21 (2000) 463–467. [3] M.B. Giles, Stability analysis of numerical interface conditions in fluid– structure thermal analysis, Int. J. Numer. Methods Fluids 25 (1997) 421– 436. [4] W.D. Henshaw, K.K. Chand, A composite grid solver for conjugate heat transfer in fluid–structure systems, J. Comput. Phys. 228 (2009) 3708–3741. [5] A. Kanevsky, M.H. Carpenter, D. Gottlieb, J.S. Hesthaven, Application of implicit–explicit high order Runge–Kutta methods to discontinuous-Galerkin schemes, J. Comput. Phys. 225 (2007) 1753–1781. [6] H. Bijl, M. Carpenter, V. Vatsa, C. Kennedy, Implicit time integration schemes for the unsteady compressible Navier–Stokes equations: laminar flow, J. Comput. Phys. 179 (2002) 313–329. [7] A.H. van Zuijlen, H. Bijl, Implicit and explicit higher order time integration schemes for structural dynamics and fluid–structure interaction computations, Comput. Struct. 83 (2005) 93–105. [8] J. Degroote, R. Haelterman, S. Annerel, P. Bruggeman, J. Vierendeels, Performance of partitioned procedures in fluid–structure interaction, Comput. Struct. 88 (2010) 446–457. [9] P. Birken, K.J. Quint, S. Hartmann, A. Meister, A time-adaptive fluid–structure interaction method for thermal coupling, Comput. Visual. Sci. 13 (2010) 331– 340. [10] B. Roe, A. Haselbacher, P.H. Geubelle, Stability of fluid–structure thermal simulations on moving grids, Int. J. Numer. Methods Fluids 54 (2007) 1097– 1117. [11] B. Roe, R. Jaiman, A. Haselbacher, P.H. Geubelle, Combined interface boundary condition method for coupled thermal simulations, Int. J. Numer. Methods Fluids 57 (2008) 329–354. [12] F. Duchaine, S. Mendez, F. Nicoud, A. Corpron, V. Moureau, T. Poinsot, Conjugate heat transfer with large eddy simulation for gas turbine components, C.R. Mec. 337 (2009) 550–561.

V. Kazemi-Kamyab et al. / Comput. Methods Appl. Mech. Engrg. 264 (2013) 205–217 [13] V. Kazemi-Kamyb, A.H. van Zuijlen, H. Bijl, Accuracy and stability analysis of a second order time-accurate loosely-coupled partitioned algorithm for transient conjugate heat transfer problems, Int. J. Numer. Methods Fluids (2013) (submitted for publication). [14] A.R. Crowell, B.A. Miller, J.J. McNamara, Computational modeling for conjugate heat transfer of shock–surface interactions on compliant skin panels, in: Second AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, 4–7 April 2011, Denver, Colorado. [15] C.A. Kennedy, M.H. Carpenter, Additive Runge–Kutta schemes for convection– diffusion-reaction equations, Appl. Numer. Math. 44 (2003) 139–181. [16] M.N. Özisßik, Boundary Value Problems of Heat Conduction, International Textbook Company, 1968. [17] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Oxford University Press, 2003.

217

[18] R. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems, SIAM, 2007. [19] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, third ed., Springer, 2004. [20] M.H. Carpenter, C.A. Kennedy, H. Bijl, S.A. Viken, V.N. Vatsa, Fourth-order Runge–Kutta schemes for fluid mechanics applications, J. Sci. Comput. 25 (2005) 157–194. [21] D.A. Kaminski, C. Prakash, Conjugate natural convection in a square enclosure: effect of conduction in one of the vertical walls, Int. J. Heat Mass Transfer 29 (1986) 1979–1988. [22] P. Wesseling, Principles of Computational Fluid Dynamics, Springer, 2000. [23] M. Schäfer, I. Teschauer, Numerical simulation of coupled fluid–solid problems, Comput. Methods Appl. Mech. Engrg. 190 (2001) 3645–3667. [24] P. Majumdar, Computational Methods for Heat and Mass Transfer, Taylor and Francis, 2005.