A coupling approach to modeling heat transfer during a full transient flight cycle

International Journal of Heat and Mass Transfer 110 (2017) 587–605

Contents lists available at ScienceDirect

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

A coupling approach to modeling heat transfer during a full transient flight cycle M.-P. Errera a,⇑, M. Lazareff a, J.-D. Garaud a, T. Soubrié b, C. Douta c, T. Federici d a

ONERA, The French Aerospace Lab, 29 Avenue de la Division Leclerc, 92322 Châtillon Cedex, France Andheo, 92322 Châtillon Cedex, France c Safran Aircraft Engines, Site de Villaroche, Rond-Point René Ravaud, 77550 Moissy Cramayel, France d Safran Tech, Pôle M&S, Rue des Jeunes Bois-Châteaufort, 78114 Magny-Les-Hameaux, France b

a r t i c l e Article history: Received 7 October Received in revised Accepted 14 March Available online 27

i n f o 2016 form 3 February 2017 2017 March 2017

Keywords: Fluid-structure interaction Conjugate heat transfer Thermal coupling Transient Numerical stability Flight cycle

a b s t r a c t The purpose of the present study is to describe novel numerical coupling schemes to analyze the transient temperature field in a solid via a conjugate heat transfer procedure. Emphasis is put on the interfacial treatment based on two complementary treatments: Dirichlet-Robin and Neumann-Robin transmission conditions. The numerical methods are first presented on the basis of a stability analysis in an aerothermal model problem. Stability conditions are expressed and the mathematical expression of the most relevant coupling parameters are provided for the first time. Furthermore, an overview of all the coefficients that can be used in a transient thermally-coupled procedure are given and a unified approach for steady and unsteady ramps is proposed. Then, these interfacial schemes are applied to the problem of convective heat transfer over, and transient conduction heat transfer within, a flat plate. A comparative study with realistic operating conditions is carried out, at low and large Biot numbers. It is shown that certain choices of coupling coefficients, even if physically reasonable, may result in nonconverging algorithms. This confirms that a model problem provides insight to the behavior of complicated heat transfer cases and constitutes an invaluable aid for generating efficient interfacial schemes. Indeed, the numerical computations demonstrate the efficiency of the numerical schemes based on the main theoretical results. The trends predicted by the model problem are recovered and excellent convergence properties are observed in all cases. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction In recent years many studies have been devoted to analyze the behavior of various conjugate heat transfer (CHT) problems, and in general, they are most often limited to steady cases, i.e. when a fluid-solid steady state is sought [1–6]. The simulation of the transient heat load in solid structures is much less common, and the literature on transient CHT has, so far, not been very extensive. It is beginning to be employed in turbomachinery applications to account for the time-dependent thermal response of structures to ambient conditions, for instance in the remarkable pioneering studies at the University of Surrey [7–9]. This approach is essential nowadays to describe accurately the unsteady heat load which could lead to substantial gains in engine performance and component reliability. Note that unsteady CHT is not restricted to turbine

⇑ Corresponding author at: ONERA/DAAA/NFLU, Châtillon, France. E-mail address: [email protected] (M.-P. Errera). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.03.048 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

blade cooling applications. It can also be found in modeling heating, cooling and ventilating flows in building simulations [10–12]. Unsteady CHT calculations are rather rare in practice because they imply significant computational costs. If in a dynamic coupling, the exchange of information between the fluid and the solid is carried out with a period equal to the fluid time step, transients are calculated and this leads to a very expensive solution, especially over a long period of time. As we are interested here only in the transient temperature analysis in a solid domain, it is clearly not a viable approach to use a transient fluid solution. Accordingly, a two-way coupling of a dynamic thermal modeling in the solid and a sequence of steady states in the fluid is generally considered [7–9,13] to analyze transient conduction. This approach is widely justified by the significant differences between time constants in the two media. It is interesting to note that this quasi-steady approach is also used in other investigations in order to improve the computational time [14–16]. The key point in the coupling implementation is the choice of relevant interfacial conditions. These conditions have a direct

588

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

Nomenclature a Bi BiD C D
D g h K L q T U y+ z

thermal diffusivity [m2 s
1] Biot number mesh Biot number heat capacity [W s K
1] Fourier number normalized Fourier number temporal amplification factor heat transfer coefficient [W m
2 K
1] thermal conductance [W m
2 K
1] solid thickness [m] heat flux [W m
2] temperature [K] velocity [m s
1] non-dimensional wall distance temporal amplification factor

Greek symbols a coupling coefficient [W m
2 K
1] e convergence tolerance j spatial amplification factor

impact on numerical properties of the coupling methodology. However, most of the interfacial numerical treatments in CHT problems are based on a Dirichlet transmission condition (temperature imposed on the fluid side) supposedly due to stability reasons. This was recommended by Giles [1] in a pioneering stability analysis. Yet, it should be kept in mind that three particular conditions were used in his analysis (time-marching procedure in both media, coupling between two unsteady media, update at every time step) and thus it is reasonable to assume that the conclusions of this investigation cannot be extended to all CHT cases. There is ample evidence that there are different levels of dependency between the fluid and the solid and many types of CHT problems can be identified. For instance, as noted by Verstraete et al. [17], at large Biot number, thermal gradients may become important, and thus it seems ‘‘natural” to impose a Neumann condition. This is confirmed by CHT test cases [18] and by stability analysis [17,19]. Another recent paper points in this direction. A coupling methodology for weakly transient conjugate heat transfer problems based on two complementary interfacial conditions was proposed by Gimenez et al. [20] and this study showed the interest in using various interfacial conditions. In the same vein, the coupling of a transient solid solution with a sequence of steady states is a very specific issue (neither of the above three particular conditions are met) that needs to be treated as such. Here, we consider this specific problem. In CHT analysis, adaptive coupling coefficients have been highlighted and expressed for the first time [6]. They are derived from a 1D prototype coupled model in which a numerical transition can be identified. This fundamental result was obtained through a normal mode stability analysis based on the theory of GodunovRyabenkii [21–25]. The performance of these coupling coefficients in the framework of a Dirichlet-Robin procedure was tested recently [19,26] and the relevance of the predictive model was fully confirmed. Indeed, the so-called ‘‘optimal coefficients” provided oscillation-free solutions and always the smallest convergence rate. Moreover, stability criteria were proposed, for selecting the most relevant interfacial conditions to implement. However, the previous work was devoted to steady CHT solutions only, and it is not possible to employ this interface model

k

m q

Dy

thermal conductivity [W m
1 K
1] unit normal vector density [Kg m
3] grid size [m]

Subscripts f fluid domain s solid domain 0+ solid side at the interface 0
fluid side at the interface Superscripts n temporal index min minimum max maximum opt optimal (^) unknown values

directly. It is then essential to carry out a new stability analysis. The primary goal of this paper is to perform an extensive analysis specially adapted to problems in a full transient flight cycle. In effect, steady and unsteady CHT have very little in common as already outlined in [27,28] where the main numerical characteristics of these two approaches were summarized and compared. Real-world systems are almost always characterized by multiple interacting physical phenomena and the interactions occur on a wide range of spatial and temporal scales. Therefore, even if the numerical solution of the individual physics is stable and robust, there is no assurance that the coupled system will be also stable and robust. It is common to observe unexpected convergence behaviors (slow convergence, instabilities, oscillations and even bifurcations). This is all the more true in unsteady CHT where this issue is more difficult and sensitive and requires a careful definition of the nature of the thermal coupling implementation at the interface since a relevant fluid-solid solution is required at each time coupling. The goals of this paper are to propose a numerical methodology, accurate and rapidly converging, to analyze the transient temperature field in a solid via a CHT method and to provide optimal local coefficients at the fluid-solid interface when a transient solid solution is coupled to a sequence of fluid steady states. Two complementary interfacial treatments, namely Dirichlet-Robin (D-R) and Neumann-Robin (N-R) conditions, will be considered. This paper is structured as follows. The problem specification and the algorithm is presented first (Section 2). Then, two interface treatments are analyzed, their similarities and differences are identified, and emphasis is put on the stability bounds and optimal coefficients (Section 3). On this theoretical basis, a unified treatment to deal with steady and unsteady ramps of a full transient cycle is described (Section 4). Then, a test case is presented (Section 5). A comparative study of convective heat transfer over, and transient conduction within, a flat plate is considered at various operating conditions. The impact on stability and convergence is then described at large (Section 6) and low (Section 7) Biot numbers. At the end of this paper, heteroclite ramps are considered (Section 8).

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

2. Problem specification and solution approach 2.1. Example of a transient cycle The aim of this study is to define relevant interfacial conditions of a two-way coupling procedure for practical CHT applications during a full transient flight cycle. This cycle is above all characterised by a long period of time, that is, the entire duration of a flight (i.e. several hours). This is needed since high requirements about the life-span of the engine elements are required in particular heat load characteristics during all the stages of a flight. These stages are expressed by ramps and these ramps are defined from ramp points, namely specific instants representing changes in operating conditions. Fig. 1 gives an example of a transient cycle where a typical quantity U (engine speed, temperature, flow rate, . . .) is plotted over time. One can see that this cycle covers a broad range of usual operating conditions from stand-still, idle, engine acceleration to maximum take-off (MTO) and cruise conditions. The time intervals, also simply called ‘‘ramps” are preestablished instants from the time position of ramp points where environment conditions are available. In these ramps, linear distributions of the environment parameters are assumed. It is only at ramp points that a fluid-structure coupling may take place and then these points can also be regarded as coupling times. Outside these specific and limited instants, flow conditions are not known and must be interpolated. In contrast, there are no restrictions on the number of solid time steps. In other words, each ramp can be subdivided into as many temporal divisions as are required to analyze the transient metal heat conduction in solid regions. 2.2. Time scale disparities A transient coupling is challenging due to the great time-scale disparities of the physical models between the fluid and the solid. The fluid flow requires usually a much smaller temporal resolution than the structure. Accordingly, a description of the transients in the fluid would be almost impractical over a period comparable to the thermal response time of the solid. It is then necessary to develop simplified coupling strategies to minimize the number of CFD runs. 2.3. A quasi-dynamic strong approach An alternative and computationally cheaper procedure is to account for the time scale disparities and consider the flow solution as a sequence of steady states. In most cases, this assumption is valid since the fluid and the solid operate on different time constants. We therefore may assume that the influence of unsteadiness in the fluid domain is negligible and thus the flow solution is considered as a sequence of steady states and the solid as a pro-

cess allowed to evolve over time and space. Such an approach is referred to as a ‘‘quasi-dynamic” procedure but we are aware that the terminology is far from standard. The fluid-solid interaction exists at both a physical level and an algorithmic level. In a loose coupling of physical models, the physical components of the global system are coupled by a weak interaction. The equilibrium of the fluid-solid interface is not enforced. This treatment is suitable in some engineering applications such as steady CHT problems or aeroelasticity. However, in the transient aerothermal coupling as defined in Fig. 1, very few coupling instants are available and thus, for accuracy reasons, it is essential to impose the equilibrium at the interface through a strong fluidsolid interaction. Strong coupling can be achieved with monolithic procedures. In this approach, the interactions between the fluid and the solid are solved simultaneously, and thus can be very time consuming in a transient problem. Furthermore, this approach often leads to poorly-conditioned systems. Alternatively, in a partitioned method, separate solvers are used for the fluid and the solid. As a result, appropriate methods at the interface must be investigated to take into account characteristic time discrepancies and to define relevant coupling strategies. As a strong-coupling algorithm must be used, a sub-iteration process is set up at every coupling time to reach a converged solution. If this process converges, then a ‘strongly coupled’ algorithm is obtained. This iteration procedure, briefly described in the next section, can be regarded as a means to reduce the computational time of a fully transient coupling and as an efficient procedure based on the time step dictated by the dynamics of the solid domain. Note that convergence can be accelerated if the influence of the flow on the structure is included in the structure [29] and in the same way, the influence of the structure is included in the flow solver [30].

2.4. Numerical algorithm The fluid-solid interaction is achieved by partitioning the problem into fluid and solid parts solved separately with boundary conditions (BC) calculated by the other part [31,32,35], even when a black box is used [36], which is the case in the present study. From a software standpoint, the most practical means of coupling two individual physics together is based on some form of iterative methods whose variants are known under different names as the Picard iteration, successive substitution or fixed-point iteration. The algorithm is rather simple to implement and this method is robust, although it is subject to stability requirements [33,34]. The most common and practical method of iteration schemes is Gauss-Seidel-like coupling iterations associated with inner iterations repeated until convergence. This sequential treatment that

U MTO cruise

idle

0

idle

4000

6000

8000 Time (sec)

10000

12000

Fig. 1. Example of a transient flight cycle - Solid lines: unsteady ramps – Dotted lines: steady ramps.

589

Fig. 2. Transient fluid-solid algorithm.

590

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

can be seen as a Conventional Serial Staggered (CSS) procedure [37]. Fig. 2 shows this sequence of iteration steps. The various steps in Fig. 2 can be defined as follows : ①: Transient calculation in the solid ②: Exchange from the solid to the fluid ③: Steady fluid computation ④: Exchange from the fluid to the solid Convergence Test – Comparison between solid & fluid states at the interface. Convergence? ⑤ No: Go back to tc (step ①). One more iteration between tc and tc + Dtc in the solid until convergence ⑤0 Yes: first transient calculation in the solid from tc + Dtc. The transient temperature in the solid is calculated in step ① using interpolated fluid thermal states as BC at every temporal solid increment. In this paper, only a linear temporal interpolation is implemented. Other choices are, however, possible. For instance a three-point interpolation based on the two thermal states at the ramp points and on the current temporal point. Another possibility is to employ a linear interpolated relationship between the heat flux and the temperature. These interpolation procedures will be implemented in the near future. After the calculation of the temperature field in the solid, the temperatures are compared to the previous ones at the coupled interfaces. If the local maximum temperature difference does not reach the criterion for convergence, new boundary temperatures are prescribed to the fluid and a CFD calculation is carried out again until the sequence converges. Then, the solid calculation moves on to the next time step. For further details see [13,38]. Note that a Gauss-Seidel procedure is inherently sequential in execution order of fluid and solid solver. This sequential algorithm is not necessarily a major problem for parallel computing as the fluid solver is much more costly and can be parallelized separately. 2.5. The key point: The fluid-solid interfacial treatment Coupling two individual physics codes together to achieve multiphysics applications may lead to unexpected numerical behaviors such as slow convergence problems, oscillations or instabilities. Even in the case where the numerical techniques of the individual physics are stable, robust, and fast-converging, there is no assurance that the coupled system will be stable and fast. This should cause us to carefully define the nature of the fluid-solid mechanisms. The fixed-point iteration procedure in Fig. 2 may converge slowly in fluid-structure interaction problems [39] or not converge at all and thus, some acceleration techniques have been considered [40]. The most critical and salient point in the implementation of a partitioned procedure is the data exchange. If a fluid-solid coupling is prone to instability or oscillations, we have to look at the root causes of such behavior and the first question to be examined is the nature of the conditions at the interface. These conditions are generally enforced to ensure stability or to increase the convergence speed by using relaxation coefficients. Unfortunately, these coefficients are most often set in a rather arbitrary and global manner. Our emphasis in this study is to examine the stability behavior of local adaptive coefficients at the fluid-solid interface taking into account the most important quantities involved in a thermal numerical coupling. These quantities are the physical and numerical parameters involved in a fluid-solid thermal interaction: the thermal properties, the characteristic dimensions, and the time step. These features are already included in a 1D model problem and it’s one of the reasons why a simplified model can give insight into the potential instabilities of 3D computations. The other reason is that one may reasonably assume that the modes that may be unstable are those whose variation is the direction normal to

the shared boundary. As a result, the analysis carried out in a simplified model should remain valid for real 2D/3D CHT calculations. Unfortunately, we will show that the results obtained in this steady case cannot be generalized. Thus, a new stability analysis has to be carried out and the resulting models must then be tested. 3. Fluid-solid interfacial conditions 3.1. Robin-Robin conditions The potential stability and convergence problems can be clarified and alleviated by using efficient coupling procedures in the general 2-coefficient framework of Robin-Robin conditions. Robin conditions exhibit many interesting features. They are based on coupling coefficients that accurately take into account the interactions across a shared boundary. These coefficients mimic the boundary of each sub-system and can be regarded as local interface stiffness. The fluid-solid conditions will now be presented, followed by a thorough and detailed study of these conditions applied to a model problem representing a transient coupling. 3.2. Fluid-solid interfacial equations at coupling points Below and throughout the text, the order of the sequences of the CSS algorithm is the one defined in Fig. 2 (the solid is performed before the fluid). However, the convergence analysis would not be changed in any significant way by assuming that the fluid problem goes first. On the solid side, a Robin condition at the shared boundary has the following form, at the nth iteration at the coupling time in the temporal interval t c þ Dtc

h i ^ns þ an
1 q T^ ns f

t c þDt c

h i ¼ qn
1 þ an
1 T n
1 f f f

ð1Þ

t c þDt c

where q is the heat flux, T the temperature and ðaf ; as Þ are coupling coefficients, with (^) denoting the sought values. qs ¼ q0þ ¼
K s @T s =@ m, is the normal heat flux where m is the inward normal to the solid domain (See Appendix A). Similarly, qf ¼ q0
¼
K f @T f =@ m where m points to the outward normal direction of the fluid domain (See Appendix A). With these convention signs, the Robin condition on the fluid side may be expressed as

h i ^nf
ans T^ nf q

tc þDt c

  ¼ qns
ans T ns tc þDtc

ð2Þ

Conditions (1) and (2) act to couple the two models and domains at their shared interface. At convergence, to make the overall model well-posed, we require heat flux normal to the interface and temperature to be continuous across the interface independently of the coefficients, provided af þ as –0 (or a1 þ a1s –0) f

for well-posedness. These coefficients are two adjustable parameters that control how energy is communicated across the common interfacial surface. They are key factors of any numerical CHT procedure. Note that in most commercial codes, Robin thermal conditions written in the form of (1) and (2) are not directly available. It is straightforward, however, to express them in terms of (h, T⁄), where h is the heat transfer coefficient and T⁄ is a reference temperature [41]. 3.3. Solid interfacial equations at ‘‘non-coupling” points At intermediate solid time steps, there is no adjacent fluid state at the same time. The BC on the solid are then interpolated between two known fluid steps

½qns þ anf T ns t

c
¼ interp: of BC from the fluid states

ð3Þ

591

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

These interpolated BC must be viewed as ‘‘pseudo-couplings” that modify the solid state at each temporal increment without modifying the stability of the coupling process.

4.1. Reminder of the optimal coefficient for steady CHT For the sake of brevity, details of the 1D stability model are not presented; the interested reader is referred to [6,19] for a more complete discussion. We will just insist here on the major differences due to the transient nature of the CHT problem addressed in this paper. Firstly, it is useful to point out the main result in steady CHT analysis. When a Dirichlet-Robin procedure is implemented (i.e., as ¼ 1), it has been shown that the temporal amplification factor of the model problem is a non-monotonic function that reaches its absolute minimum, always less than unity, at a single coupling parameter referred to as the ‘‘optimal” coefficient. Its expression is given by

aðoptÞ f

ð4Þ


f is a normalized Fourier number that will be defined in where D Section 4.5 and K f denotes the thermal conductance of the 1st fluid cell at the interface (K f ¼ 2kf =Dyf in a finite-volume approach). This expression, efficient in a steady context, is neither pertinent nor useful to the unsteady issues presented in this paper. Eq. (4) is just provided for comparison purposes. 4.2. Normal mode analysis A stability analysis is carried out (see Appendices A and B) by considering a coupling prototype problem composed of two subdomains (j P 0 and j 6 0) with a common boundary ðj ¼ 0Þ. Then, normal mode solutions are expressed in the following form

(

T nj

¼

zn
1 jþj ; j P 0 z

n

j

j
;

j60

ð7Þ

This solution avoids an increasing exponential solution as j ! þ1.

4. Optimal procedures for steady and unsteady CHT

Kf
f Þ ¼ ð1
D 2

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 z
1 z
1 js ðDs ; zÞ ¼ 1 þ
1
1þ 2Ds z 2Ds z

ð5Þ

4.4.2. Fluid spatial amplification factor In the fluid domain, the spatial amplification factor has the following form (See B-6)

j
¼ jf ¼

h 1
2kf =Dyf

!
2

 
2 h ¼ 1
Kf

ð8Þ

Consequently, the condition jjf j j > 1 avoids an increasing exponential solution as j !
1, if the following stability condition is verified

h <1 2K f

ð9Þ

This simple expression is worth considering. In the stability analysis performed in [6], the trivial solution j
¼ 1 was obtained. This time, the fluid domain provides a stability condition jj
j > 1 verified only if h < 2K f , which can be assumed in the vast majority of cases. However, this condition should not be overlooked. It means that diffusion in the fluid near-wall region must prevail over convection. As a result, when the convective coefficient is expected to be high (impinging jets for example), the mesh grid in the nearwall region must be refined enough such that condition (9) holds. It should also be noted that this condition is independent of the coupling process and must be verified independently in the fluid domain. We do not know whether this condition has been expressed before. 4.4.3. Temporal amplification factor The expression of the temporal amplification factor is given by (see details in Appendix B)

z ¼ gðz; af ; as Þ ¼

ðaf þ as ÞK s ðh
af ÞðK s
as Þ js ðz; Ds Þ þ ðh þ as ÞðK s þ af Þ ðh þ as ÞðK s þ af Þ ð10Þ

where j and z are two complex functions. j is the ‘‘spatial amplifi-

cation factor” (jþj in the sub-domain j P 0 and j
j in the subdomainj 6 0). z is the ‘‘temporal amplification factor”, which means that from one time step n
1 to the next n, each mode increases in amplitude by a ratio jzj.

where K s ðK s ¼ ks =Dys Þ denotes the thermal conductance of the 1st solid cell (see Appendix A). Ds represents the mesh Fourier number in the dynamic solid domain defined by

4.3. Stability in the sense of Godunov-Ryabenkii (G-R)

where as is the solid thermal diffusivity, Dt is the time step and Dys is the solid grid size. The mesh Fourier number characterizes heat conduction in a transient medium. It can be thought of as the dimensionless measure of time used in transient conduction discrete problems characterized by the space interval Dys . At this point, the following three important results are recalled (see [6] for more details):

It can be shown (see Appendix B) that the temporal amplification factor can be written as follows

z ¼ gðz; af ; as Þ

ð6Þ

The 1D discrete model problem will be stable in the sense of Godunov-Ryabenkii provided there are no solutions to z ¼ gðzÞ, for all z, with jzj P 1. 4.4. Stability conditions for unsteady CHT All stability conditions, expressed in Appendices A and B, are discussed here. 4.4.1. Solid spatial amplification factor In the solid domain, the spatial amplification factor js is the solution to a quadratic equation, and the root that fulfils the condition jjs j j < 1 (see Appendix B), is given by

Ds ¼ as Dt=Dy2s

ð11Þ

(a) the maximum of jgj for jzj P 1 (max jgj), also referred to as jzjP1

the temporal amplification factor, is reached on the unit circle jzj ¼ 1 (maximum principle). (b) on the unit circle, a transition occurs when

jgðz ¼ 1; af ; as Þj ¼ jgðz ¼
1; af ; as Þj

ð12Þ

and at this transition, the minimum of max jgj is reached, with jzjP1

max jgj ¼ max jgj < 1 always verified. jzjP1

jzj¼1

592

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

(c) as a result, at this transition, the G-R stability criterion gðzÞ–z for jzj P 1 is automatically fulfilled. The specific coefficients for this transition to occur are referred to as the ‘‘optimal coefficients” and their mathematical expressions will be determined in § 5.2 ðas ¼ 1Þ and § 5.3 ðas ¼ 0Þ in the transient CHT analysis. Note that for Dirichlet-Robin conditions in steady CHT analysis, the transition equation (12) always holds. However, we will see in this present article (Section 5.2) that this is not always true. The coefficients ðaf ; as Þ for which max jgj is located at z ¼
1 correspond to a fast process, prone to instability. On the contrary, the coefficients for which max jgj is located at z ¼ þ1 can generate a slow convergence but the process is always stable. The intersection combines the advantages of the first two options, without suffering their disadvantages: the coupling process is fast and is always stable. It is the reason why the conditions where this intersection occurs (transition) play a fundamental role in CHT problems. 4.5. A normalized Fourier number A last point to conclude this section. For simplicity, we define the function representing a normalized Fourier number in terms of Ds

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 1 1 Ds
s ¼ 1 þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi D 1þ
1¼ Ds Ds 1 þ Ds þ 1 þ 2Ds

ð13Þ


s is a monotonic increasing function that plays an essential role in D transient heat transfer. It is defined in the domain ½0; 1½ and varies in the range ½0; 1½. Introducing this normalized Fourier number simplifies the writing of various expressions in the stability analysis. For instance,
s can also be employed this number has been used in Eq. (4) and D to simplify the expression of the complex function js at the specific (and strategic) points, as we can see from (7)



js ðDs ; z ¼ þ1Þ ¼ 1 js ðDs ; z ¼
1Þ ¼ D
s

ð14Þ

5.2. Dirichlet-Robin (D-R) transmission procedure 5.2.1. Conditions This transmission scheme is widely used in the literature. Temperature is imposed on the fluid side and a Robin condition on the solid side. These conditions are obtained by introducing the coefficients



af P 0 as ¼ 1

ð15Þ

into (1) and (2). The following interfacial conditions are thus obtained :

(

^ns þ an
1 q T^ ns ¼ qn
1 þ an
1 T n
1 s f f f T^ nf ¼ T ns

ð16Þ

Note that the temperature on the right-hand side of the first equation in (16) is located in the solid domain. Thus, only the fluid heat flux is sent to the solid and has to be interpolated. 5.2.2. Stability analysis Substituting (15) into (10), the amplification factor of this specific transmission procedure can be written simply as

gðz; af Þ ¼

K s  js ðz; Ds Þ þ af
h K s þ af

ð17Þ

Using (14) and (17), the stability condition jgj < 1 becomes


1 <


s þ af
h K sD < þ1 K s þ af

ð18Þ

The inequality on the right-hand side always holds. Thus, the stability condition is readily obtained

af >

h Ks
sÞ
ð1 þ D 2 2

ð19Þ

By introducing the ‘‘mesh Biot number” defined by

BiD ¼

h Ks

ð20Þ

Then, the lower stability bound can be expressed as 5. Mathematical study of two complementary interfacial transmission schemes 5.1. Two transmission conditions considered The general formulation of a 2-coefficient method covers a very large family of schemes. It is preferable to initially limit our investigation to two opposite transmission conditions in the fluid domain, while considering a general Robin condition in the solid domain. This approach has already been used in [20] for weakly transient CHT problems and as result, the present paper can be regarded as a continuation of the process. These two complementary approaches will be used at the fluid-solid boundary: (a) Dirichlet-Robin conditions: solid temperature imposed on the fluid side. (b) Neumann-Robin conditions: solid heat flux imposed on the fluid side. The condition as ¼ 0 results in pure flux (Neumann) BC (perfectly insulating boundary) and as ¼ 1 causes pure temperature (Dirichlet) BC (perfectly conducting). As a result, it seems, in theory, that these two opposite conditions create an ideal basis for better understanding and addressing complex transient heat transfer problems.

min amin ¼ unsteady ¼ af

Ks
s Þ ½BiD
ð1 þ D 2

ð21Þ

The dimensionless number BiD is used here to characterize the heat transferred from the first solid grid cell to the surrounding fluid. Condition (21) reveals the existence of three different zones in the numerical behavior of the D-R transmission procedure in transient CHT that can be obtained by simple algebra and defined in the following way:
s  is characterized by a lower sta(1) The first region ½BiD > 1 þ D bility limit as clearly outlined by (21). In this region, when af > amin , jgj first decreases from 1 to a certain limit and f afterwards increases from this extremum to 1. This limit is the absolute minimum of the amplification factor, always located in the stable zone and provided by the optimal coefficient.
s  does not exhi
s Þ=2 < BiD < 1 þ D (2) The second region ½ð1 þ D bit any stability restriction. It is a narrow region where as previously, jgj has an absolute minimum, namely when af is equal to the optimal coefficient.
s Þ=2 is characterized by a con(3) The third region ½BiD < ð1 þ D tinually increasing function in the domain af P 0. In this case, as in the preceding one, there is no stability restriction.

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

The minimum of the amplification factor is obtained for af ¼ 0. In other words, in this third region, there is no optimal coefficient. Note that only in this region, the minimum is obtained when the derivative vanishes. In the first two regions, the minimum is the one and only intersection of two functions. The mathematical expression of the absolute minimum, which only exists in the first two regions, the ‘‘optimal” coefficient is readily provided by (12)

aopt ¼h
f


 Ks
s Þ ¼ K s BiD
1 þ Ds ð1 þ D 2 2

ð22Þ

It can therefore be seen that the optimal coefficient is defined by the difference between the heat transfer coefficient (at steady state) and a combination of the solid conductance (during the solid transients) and the Fourier number. At very large Fourier number
s  1), the expression (22) simply reduces to the difference (D between the fluid and solid thermal conductances (h
K s ). It is the only result that could have been ‘‘intuitively” found. Let e0 P 0 be given. Then, the best choice for the coupling coefficient, adapted for all conditions, can be simply expressed as

af ¼ maxðe0 ; K s  ½BiD
ð1 þ D
s Þ=2Þ

ð23Þ

All these results are summarized in Table 1, where the domain and range of jgj are provided. The conditions depend upon both the mesh Biot number and the Fourier number. In the bottom row of Table 1, the values g 01 ; g 02 ; g opt are easily obtained from (17). 5.3. Neumann-Robin (N-R) transmission procedure 5.3.1. Conditions This transmission scheme is less common in the literature. The heat flux is imposed on the fluid side and a Robin condition on the solid side. These conditions are obtained by introducing the coefficients



af P 0 as ¼ 0

ð24Þ

(

^ns þ an
1 q T^ ns ¼ qn
1 þ an
1 T n
1 s f f f ^nf ¼ qns q

593

ð25Þ

Note that, the heat flux on the right-hand side of the first equation in (25) is located in the solid domain. Thus, only the fluid temperature is sent to the solid and has to be interpolated. 5.3.2. Stability analysis Substituting (24) into (10), the amplification factor of this specific transmission procedure can be written simply as

gðzÞ ¼

af D
s þ ðh
af Þ ðK s þ af ÞBiD

ð26Þ

Now, taking into account the guidelines and remarks put forward in the previous section, it is easy to show that the stability condition jgj < 1 in the numerical behavior of N-R transmission procedure is composed of two different zones:
s  does not exhibit any stability (1) The first region ½BiD > 1
D restriction. In this region, jgj has an absolute minimum, namely when af is equal to the optimal coefficient. This limit is the absolute minimum of the amplification factor, always located in the stable zone and defined by the optimal coefficient.
s  is characterized by an upper (2) The second region ½BiD 6 1
D stability limit. max amax ¼ unsteady ¼ af

2h
s Þ
BiD ð1
D

ð27Þ

In this region, jgj first decreases from 1 to a certain limit and afterwards increases from this extremum to 1. This limit is the absolute minimum of the amplification factor provided by the optimal coefficient. As we can see, an optimal coefficient can this time be defined over the whole range af P 0 and for all mesh Biot numbers. Its expression, obtained directly from (12) is given by

aopt ¼ f

into (1) and (2). The following interfacial conditions are thus obtained: Table 1 Dirichlet-Robin procedure for transient analysis ðas ¼ 1Þ – Stability conditions.

2h
s 1
D

ð28Þ

594

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

It can be seen that the optimal coefficient corresponds to the
s Þ). quotient between 2h and a solid transient effect (term ð1
D
At very large Fourier number (Ds  1Þ, the optimal coefficient becomes infinite, which means that a Neumann-Dirichlet condition is appropriate. These results are summarized in Table 2, where the domain and range of jgj are provided. In the bottom row of Table 2, the values g 1 and g opt are easily obtained from (26).

impact on stability of a N-D procedure that can turn the coupling process into an unconditionally stable method. Instead, this cannot be achieved with the D-R procedure in which the stability limit can
s Þ is only never be removed and the influence of the term ð1 þ D slightly stabilizing. This indicates that the stability behavior of CHT conditions is not intuitive and obvious. Understanding the reasons for differences in performance is crucial and implies a detailed analysis taking into account the main physical and numerical parameters such as the local Biot and Fourier numbers.

5.4. Comparative summary 6. A unified approach for steady and unsteady ramps The interfacial schemes presented in the preceding section, as well as the stability bounds and the optimal coefficients have been obtained from a normal mode analysis on the basis of the Godunov-Ryabenkii theory. Even if this theory is well known, the results derived from this methodology in the context of transient CHT problems are completely new, presented for the first time in this paper. An analysis of these results is convincing in this respect. For instance, the optimal coefficient for steady CHT (Eq. (4)) has nothing in common with the one for a transient analysis (Eq. (22)). Indeed, it is most instructive to notice that in the former expression, no ‘‘solid” terms are present while in the latter, fluid and solid terms are directly competiting. More generally, the stability bounds are different and there is a zone (3rd row in Table 1) where no optimal coefficient can be defined. A similar argument can be made about the Neumann-Robin condition. This clearly shows that the same boundary condition exhibits properties that will inevitably depend upon the problem sought to be addressed. It is rather clear from Table 1 (lower stability bound at large Biot number & unconditional stability at large Biot number) and Table 2 (upper stability bound at low Biot number & unconditional stability at low Biot number) that D-R and N-R conditions are complementary conditions with converse properties. This complementary nature of the two transmission procedures allows us to find, for a given mesh Biot number, a zone of unconditional stability. This is a very positive feature which cannot be achieved by a single transmission procedure. Another major difference can also be observed from Tables 1 and 2. In the first case, the time step (i.e. the Fourier number) has only a small effect on stability and in general cannot remove the existing stability restriction (amin in Table 1). On the contrary, f there is obviously a temporal threshold (i.e. a Fourier number) in Table 2) does not exist above which the stability restriction (amax f anymore. This aspect is less important in the scope of this study. Indeed, looking at the transients implies taking a small Fourier number (Ds ¼ 0:5; 1). For higher values of Ds , in another context
s Þ can have a significant and for a different purpose, the term ð1
D

6.1. Steady and quasi-steady states It should first be emphasized that the ‘‘steady strategy” in which only a steady state is sought has nothing in common with the numerical approach proposed in this paper. For instance, the Dirichlet-Robin procedure proposed for the first time in [6] based on the ‘‘optimal coefficient” concept is radically different in several ways with the Dirichlet procedure presented here where the influence of unsteadiness in the fluid domain is negligible [41]. However, the ‘‘unsteady” stability analysis is the most general approach and thus may be simply extended to steady states by replacing BiD by Bi (a dimensionless number evaluated with the
s ! 0. As a consecharacteristic thickness L of the body) and D quence, it is an easy matter to obtain from Tables 1 and 2 the stability bounds and optimal coefficients at steady state. 6.2. Stability bounds and optimal coefficient at steady-state for Dirichlet-Robin condition At large Biot number, Table 1 shows (1st row) that there is a lower stability bound that becomes at steady state

h 2


amin steady ¼

  Ks Ks K s hL
1 ¼ ðBi
1Þ ¼ 2 2 2 ks

ð29Þ

In the same manner, the optimal coefficient becomes at steady state

aopt steady ¼ h


Ks Ks ¼ ð2Bi
1Þ 2 2

ð30Þ

Biot numbers larger than one characterize difficult problems due to the non-uniformity of temperature fields within the solid. Thermal gradients within the solid may become important, even though the material is a good conductor. It is the reason why a temperature sent to the fluid side exhibits stability limits precisely defined in the present paper.

Table 2 Neumann-Robin procedure for transient analysis ðas ¼ 0Þ – Stability conditions.

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

This stability bound is an increasing function of the thermal penetration depth as illustrated by Eq. (21) or by Eq. (29). Consequently, the greater the penetration depth, the greater the value of the stability bound amin . Put another way, the stability limit prof vided by a given penetration depth is able to stabilize phenomena with shorter penetrations but is unstable for higher values. As at steady state, the thermal penetration depth becomes the solid thickness (expressed in the conventional Biot number), the resultmin min ing stability limit amin steady ðasteady > aunsteady Þ provides a valid limit for all the lower values of penetration and as result for all transient phenomena. At low Biot numbers, there is no stability bound as shown in Table 1 (2nd and 3rd row). 6.3. Stability bounds and optimal coefficient at steady-state for Neumann-Robin condition

2h 1
Bi

ð31Þ

In the same manner, the optimal coefficient becomes at steady state

aopt steady ¼ 2h

they are complementary and provide the coefficients to stabilize any of the above interface treatments. At low or moderate Biot numbers, it is recommended to adopt a D-R procedure. At large Biot numbers, a N-R procedure is preferable. If a D-R procedure is adopted, the general numerical conditions are provided by Table 1. To avoid any risk of instability at large Biot numbers, it is recommended to use, for steady and unsteady ramps, coupling coefficients such that af > amin steady and ideally

af ¼ aopt steady . However, it should be recalled that this condition is well suited for low or moderate Biot numbers in which no stability problems are theoretically detected. If a N-R procedure is adopted, the general numerical conditions are provided by Table 2. To avoid any risk of instability at low or moderate Biot numbers, it is recommended to use, for steady and unsteady ramps, coupling coefficients such that af < amax unsteady and ideally af ¼ aopt steady . However, it should be recalled that this condi-

At low Biot number, Table 2 shows (2nd row) that there is an upper stability bound that becomes at steady state

amax steady ¼

595

ð32Þ

At moderate or low Biot numbers (smaller than 1), problems are generally thermally simple, due to relatively uniform temperature fields inside the body although this temperature may be changing as heat passes into the solid from the surface. It is therefore quite understandable that a heat flux sent to the fluid side exhibits stability limits. The same reasoning applies in this case. This stability bound is also an increasing function of the thermal penetration depth as illustrated by Eq. (27) or by Eq. (31). Consequently, the greater the penetration depth, the greater the value of the stability bound amax . But this time, we arrive at exactly the opposite conclusion. f The stability limit provided by a given penetration depth is able to stabilize phenomena with larger penetrations but is unstable for lower values. As a result, the upper stability bound, for all penetrations depths, is given by the smallest penetration Dys . As at steady state, the penetration depth becomes the solid thickness, max max the stability limit amax unsteady ðaunsteady < asteady Þ provides a valid limit for all greater values of penetration and thus for all transient and opt steady phenomena. Interestingly, when Dys ! 0, amax unsteady ! asteady . This means that the optimal value is very close to the upper stability bound based on the smallest scale. At large Biot numbers, there is no stability bound as shown in Table 2 (1st row). It is worth noting that Verstraete [17,18] has developed another numerical approach, limited to steady CHT, and has drawn similar conclusions with the same stability bounds. 6.4. Strategy for a unified approach in a full transient cycle In a transient CHT procedure, the thermal penetration depth in the solid body is not known. However, the stability bounds and optimal coefficients, just established in Sections 6.2 and 6.3, take into account all penetration depths. Consequently, they allow us to define in a simple way stable conditions that adapt easily to all ramps in a full transient cycle. First and foremost, it is necessary to determine an average Biot number in order to consider a single method (either D-R or N-R) at a given fluid-solid interface. Tables 1 and 2 have clearly shown that

tion is well suited for large Biot numbers in which no stability problems are theoretically detected. This approach allows us to define in a simple way stable conditions for all penetration depths in transient CHT analysis. In theory, this approach is valid, i.e. stable, for all kinds of ramps in a transient flight cycle and thus adapts easily to all operating scenarios. Nevertheless, we should bear in mind that strong transient conditions may justify the use of the unsteady coefficients, precisely defined in this paper. Otherwise the unified treatment could be slightly less favourable, in terms of CPU time. 7. Test case 7.1. Numerical tools The computer code in the solid, called Z-set [42], is a threedimensional finite-element solver. It is an object-oriented code for material and structural mechanics simulations, with many non-linear solution capabilities. In the present study, only the thermal solver has been used. Assuming that there are no heat sources, the temperature field within the solid body is modeled according to the balance equation

qs C s @ t T
r  ðKrTÞ ¼ 0

ð33Þ

where T = T(t, x, y, z) is the unknown unsteady temperature field within the body, qs the density and Cs the solid heat capacity. K is a thermal conductivity matrix. For this work, K is treated as a scalar constant, e.g. K = ks. All CFD computations have been performed using the commercial code Fluent. This code is based on a finite-volume approach and the k
x SST turbulence model is employed. It is a robust two-equation eddy-viscosity turbulence model that combines the k
x turbulence model and k
e turbulence model such that the k
x is used in the inner region of the boundary layer in which it is more accurate and switches to the k
e in the free shear flow [43–45]. No wall functions have been employed. Indeed, a refined fluid mesh in the boundary layer has been used such that the criterion of yþ in the first cell is yþ  1. The coupling between these two computer codes have been carried out through the coupling library OpenPalm developed jointly by Cerfacs and ONERA [46]. 7.2. Geometry and operating conditions In an attempt to demonstrate the efficacy of the numerical strategy proposed in this paper, a simple CHT test case with operating conditions covering all the main elements highlighted, is presented in Fig. 3.

596

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

y

U T

80 m.s-1 1200 K = =

L=350 mm L=3 mm

x

Time-dependent temperature T=f(t) Fig. 3. Geometry & operating conditions.

7.3. Fluid mesh and convective coefficient

7.4. Unsteady conditions

Sections 5 and 6 have shown that the numerical interfacial methods are directly dependent on the Biot number. As the interest here is on implementing adaptive local schemes, it is fundamental to cover a broad range of different Biot numbers and ideally within the same flowfield. The mesh illustrated in Fig. 4 has been adapted because it produces a fairly high exchange coefficient at the leading edge. Fig. 5 shows the heat transfer coefficient h along the flat plate obtained with the grid of Fig. 4 and using the reference temperature T1. It is seen that it is sufficiently refined at the leading edge to obtain a coefficient that varies from 2800 to 70. This interval (from 40 to 1!) is large enough to obtain Biot numbers in the same proportions. Thus, all the conditions, identified in the last column in Tables 1 and 2, will be represented within the same flow field. The aim of this preliminary study focused on CFD only, is not to extract a precise estimation of the convective coefficient h to eventually implement it during the CHT process. We just want to generate a well-suited mesh, especially at the leading edge and at the interface. The values represented in Fig. 5 will not be used in the CHT computation. As described earlier, the convective coefficient will be estimated during the coupled simulations. The solid mesh is coincident at the interface. This precaution was taken to avoid any spatial grid-to-grid interpolation error and to focus the numerical study on stability issues. In the ydirection, the solid contains 11 mesh-points uniformly distributed. Fig. 6 shows a part of the fluid-solid mesh which is extremely

Two ramps are considered as shown in Fig. 7. The first ramp could be regarded as a sudden acceleration during 9 s, followed by cruise conditions.

refined at the interface ðDyf ¼ 2  10
6 mÞ to capture heat transfers without the need of wall functions. This is done by choosing the first grid point away from the wall such that the nondimensional distance yþ is close to unity. This is a recommended option in the case of the k-x SST turbulence model.

7.5. Characteristic numbers Two cases are considered in this study, as shown in Table 3. 7.6. Convergence criterium The convergence criterion (error tolerance or stopping criterion) used in this paper is based on the infinity norm of the absolute interfacial temperature, defined by

DT s ¼ maxjT kþ1
T ki j i i¼1;N s

ð34Þ

where k is the coupling iteration and i represents the spatial index along the solid interface (Ns faces). We impose

DT s < e with e ¼ 10
3

ð35Þ

7.7. ‘‘Fluid” stability condition The stability analysis has highlighted the ‘‘fluid” stability condition (9) derived from the ‘‘spatial” amplification factor. This very particular threshold to be respected is not connected to the coupled temporal amplification factor and thus can be verified independently. All the values required by Eq. (9) are provided in Table 4. The maximum value of the heat transfer coefficient in Table 4 was estimated from a CFD test (see Fig. 5). In a CHT computation, this coefficient may vary but not enough to be significantly

Fig. 4. Fluid mesh at the leading edge.

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

597

Fig. 5. Heat transfer coefficient, h, along the flat plate.

Fig. 6. Fluid-solid mesh at the interface.

8. Comparative study of interfacial treatments at large Biot number 8.1. Dirichlet-Robin procedure

Fig. 7. Cycle composed of unsteady (9 s) and steady ramps.

affected. It is then clear that the mesh grid in the near-wall region is refined enough such that condition (9) holds. As stated above, this requirement is independent of the ‘‘coupled” stability conditions provided in Tables 1 and 2, that we are going to examine in the next sections.

This study, at large Biot number, employs the values indicated in the first row of Table 3 (case I) and the corresponding numerical properties are summarized in the first two rows of Table 1. Fig. 8 shows the maximum of temperature error (see Eq. (34)) of the coupling using the Dirichlet-Robin procedure with a coupling coefficient equal the heat transfer coefficient (af ¼ h). It can be seen that large oscillations arise, leading to a nonconvergent (and also non-divergent) algorithm. This is a striking example because it shows that the most widely used condition (Dirichlet on the fluid side) associated with the most widely used coefficient (af ¼ h) may lead to serious numerical problems. This exemplifies that the choice of certain coupling coefficients, which seem to be physically reasonable and used routinely in industry for design purposes, may result in non-converging behavior. This is fully explained by the model problem. Indeed, if we look at the 1st row of Table 1, af  h is inside the stable zone. This would seem

598

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

Table 3 Physical constants and Biot numbers. ks (W/m/K) 0.15 13.6

qC p (J/m3/K)

Bi

BiD 4.6 – 0.2 0.05 – 0.0022

Case

4

1.5  10 1.36  106

46 – 2 0.5 – 0.022

Large Biot number (Case I) Low Biot number (Case II)

Table 4 Stability condition in the fluid domain. hmax (W m
2 K
1) 3

2.8  10

kf (W m
1 K
1)
2

2.5  10

Dyf (m)
6

2  10

K f ¼ 2kf =Dyf 4

2.5  10

hmax =2K f 0.056

Fig. 8. Large oscillations for the interface temperature.

Fig. 9. Convergence behavior for various coefficients. Case I – Large Biot numberDirichlet-Robin condition.

to be a paradox. Actually, this Table also shows an instability at large Biot number and the corresponding lower stability bound is amin  2h. Thus, the interval, between this lower stability bound f

Other CHT computations have been carried out and are summarized in Fig. 10. The qualitative and quantitative prediction provided by the 1D model (1st row of Table 1) is in perfect

and af  h is very narrow which increasingly limits the ability of the D-R procedure to take into account regions where the Biot number is large. As a result, for stability reasons, a Dirichlet condition on the fluid side is clearly not appropriate even if it can always be stabilized as shown in the 1st row of Table 1. Indeed, it can be seen why this might work in some cases and also why it might be very limited and possibly hazardous in others since h is not always clearly defined. It is, however, necessary to extend the analysis to determine the exact numerical properties of the D-R scheme. A detailed analysis has been performed and illustrated in Fig. 9 which shows the behavior of various coupling coefficient parameters. Of specific interest are the small or large oscillations at af  h and also at af  1:4h. However, these oscillations do not lead to an instability. This is really problematic in industrial calculations. Naturally, an oscillation of an amplitude of 10–100 degrees will never be acceptable, but it can be dangerous to accept an oscillation of just one degree, as depicted in Fig. 9. It is much better to employ the optimal coefficient that leads to a fast and monotone convergence process as illustrated by the figure. Another important point must be specified. The convective coefficient h varies significantly near the leading edge and within a single cell; it can vary by a maximum factor of 2, 4 between the upstream and down stream faces. Thus, this variation (h becomes h ) has been considered in the calculation of h.

Fig. 10. Convergence vs coupling coefficient. Case I – Large Biot number- DirichletRobin condition.

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

599

agreement with the numerical CHT computations. The theoretical  ¼ h2 and perfectly confirmed by stability bound is actually amin f the computations. 8.2. Neumann-Robin procedure The mathematical model predicts (see Table 2 – row 1): – an unconditional stable behavior 8af P 0

 2h – an optimal coefficient aopt f

– a ‘‘frozen” convergence for af ! 0

Constant values of the coupling coefficients (af ¼ 10 and af ¼ 100 in Fig. 11) were initially tested showing better convergence for the larger of these two values. However, further increase in this coefficient (af ¼ 1000) results in a deterioration of the convergence because large and constant values slow down the convergence rate at the trailing edge. Contrary to what one might expect, Fig. 11 indicates that af ¼ h seems to be the preferable choice for a

convergence rate used in this paper ðe ¼ 10
3 Þ (better that the optimal coefficient !). Note that another study with coupling coefficients close to h, has shown an optimal coefficient of aopt  1:3h. f The question arises as to whether the value of the coupling coefficient depends upon the adopted convergence rate

ðe ¼ 10 Þ. The only way to find out this information is to investigate all the way to the smallest value of e to provide a precise answer to this question !. By way of example, Fig. 12 shows CHT computations with a convergence rate e as small as possible (the smallest that can be obtained). Note that, in general, convergence performances of coupling coefficients do not remain uniform at every instant of the iterative procedure. By and large we can say that the coefficients have rather good convergence properties, at least in the initial moments. However, as the solution evolves, these performances may be impaired more or less dramatically. For instance, the sudden break of the slope for af ¼ 1:4h in Fig. 9 illustrates this tendency and shows that this coefficient exhibits excellent convergence properties during the first five iterations but suddenly a slight oscillatory behavior is observed. An acceptable level of convergence cannot be achieved. The same phenomenon occurs in
3

Fig. 12. Convergence behavior for various coefficients with e ¼ 10
10 . Case I – Large Biot number- Neumann-Robin condition.

Fig. 12 where a kink can be observed at the 7th iteration for af ¼ 1:3h but in less dramatic fashion since a convergence solution is obtained. In both cases, remarkable convergence properties were observed at the initial instants, even with higher performances than the optimal coefficient. These phenomena are due to a less than ideal energy transfer between the fluid and solid media that ultimately appears as the solution evolves. In the case of Fig. 9, the fluid-solid system is under-relaxed (af < aopt in D-R procef dure) which generates oscillations quickly. At lower coefficients such as af < amin , the system becomes unstable. In the case of f in N-R procedure) Fig. 12, the system is over-relaxed (af < aopt f and energy transfer is unnecessarily slightly frozen. At the lower limit af ¼ 0, the system would remain static.

 2h proWhat is remarkable is that the optimal coefficient aopt f

vided by the model problem theory, is found to be the fastest coupling coefficient beyond a certain threshold of convergence. It is worth noting that very small convergence rates ðe ¼ 10
10 Þ can be reached. 8.3. Comparative study A Dirichlet condition is not satisfactory for large Biot numbers. An instability zone is identified theoretically and confirmed by the CHT computations. However, it can be argued that an ‘‘optimal” treatment exists, but it could be dangerous to employ it directly as it involves a heat transfer coefficient depending on a reference temperature that can rarely be correctly estimated. Anyway, employing coupling coefficients clearly on the ‘‘stable branch” ðaf  aopt f Þ, a monotonic converging process can always be

Fig. 11. Convergence behavior for various coefficients. Case I – Large Biot numberNeumann-Robin condition.

obtained. To sum up, a Neumann condition is ideal for large Biot numbers as shown in the theoretical study (first row of Table 2), which does not predict any instability, and is amply supported by the present CHT test cases. The optimal coefficient provides, perhaps surprisingly, the fastest convergence rate. Indeed, as already mentioned the process is unconditionally stable and does not present any pathological behavior (oscillations, instabilities, slow convergence). Furthermore, it is interesting to note that a simple ‘‘Neu mann-Dirichlet” condition, with no relaxation at all, is not the best-performing procedure, even if it is always stable.

600

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

In summary, when the coupling coefficients are employed on the basis on the 1D model, no numerical problems have been identified. 9. Comparative study of interface treatments at low Biot number 9.1. Dirichlet-Robin procedure This study, at low Biot number, employs the values indicated in the 2nd row of Table 3 (Case II) and the corresponding numerical properties are summarized in the third row of Table 1. Other systematic studies have been performed and the resulting maximum error is shown in Fig. 13 where the outcomes of 7 different computations (7 different constant values of the coupling coefficient af , from 0 to 105) are presented. From this convergence plot, it can be seen that: – a stable process is found for any af P 0 – the fastest convergence is obtained for af ¼ 0 (8 coupling iterations needed) – the number of coupling iterations tends to increase as af increases

Fig. 14. Zoomed view on convergence behavior for various coefficients. Case II – Low Biot number Dirichlet-Robin condition.

These numerical results are in perfect agreement with the results obtained from the 1D model summarized in the 3rd row of Table 1. Clearly, the amplification factor rises as af rises until the limit jgj ¼ 1. Notice that this row is the only one that does not exhibit any optimal coupling coefficient in theory. In practice, if we look a little more closely, by implementing two small values (af ¼ 100 & af ¼ 300), a moderate improvement from the previous results is obtained (see zoomed view of Fig. 14). These trends are not predicted by the theoretical model and despite their small impact, it is useful to report them. This implies that rather than imposing directly the heat flux on the solid side, with no relaxation, a small amount of relaxation has a positive impact on the convergence rate. A summary of all these test cases with a Dirichlet-Robin condition is provided in Fig. 15 which gives the number of iterations needed to converge (with a maximum absolute errore ¼ 10
3 ) against af . As predicted by the 1D model, this function varies con-

Fig. 15. Convergence vs coupling coefficient. Case II – Low Biot number DirichletRobin condition.

tinuously, with no unstable zone. A very slight disagreement is noted around af  300, i.e. af  hmax =10,where the bestperforming process is obtained (gain of just one iteration) whereas the 1D model predicted it at af ¼ 0. 9.2. Neumann-Robin procedure Unlike the preceding case, the mathematical model predicts (see Table 2 – row 2) a numerical behavior far from being monotonic: – a ‘‘frozen” convergence for af ¼ 0

– an increasing convergence rate until af ¼ aopt f

Fig. 13. Convergence behavior for various coefficients. Case II – Low Biot number Dirichlet-Robin condition.

– the convergence process slows down from aopt to amax f f – an unstable zone for af > amax f

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

601

We are now aiming to accurately assess the Neumann-Robin condition through many numerical tests, to see if we can confirm the results predicted by the 1D model theory, quantitatively and qualitatively. First, small and constant values of af were used, because according to the model, the optimal value is around 2h, and h varies greatly from hmax  2800 (leading edge) to hmin  70 (trailing edge). Fig. 16 illustrates these tests with 7 outcomes. For af ¼ 0, a very slow convergence is observed. After the initialization of the coupling (independent of the coefficient) everything seems to be effectively frozen. Fig. 16 also shows: – an increasing convergence rate from af ¼ 0 – an opposite trend around af ¼ 150 – for higher coefficients, this trend is being accentuated (convergence in 58 iterations for af ¼ 200) This demonstrates that we cannot increase the coefficient globally, regardless of the local fluid-solid interaction. The coefficient for which the convergence slows down is very low, and close to 2hmin . This zone is detrimental to the performance of the entire interface. This means that it would be very costly to stabilize the coupling process by taking a constant coefficient, which might be adjusted to the specific leading edge. After that, the analysis has been refined further to find the optimal coefficient of the CHT computation. Fig. 17 shows a zoomed view of all these tests (the outcomes of 10 computations are represented !). The optimal coefficient of 1:4h is acceptable (in theory, aopt  2h). f

Fig. 17. Zoomed view on convergence behavior around af ¼ h. Case II – Low Biot number- Neumann-Robin condition.

 2h, until an unstaOther tests have been carried out from aopt f ble solution is obtained. In quantitative terms, the CHT calculations show that aopt  1:4h and amax  2:6h, whereas the model gives f f

aopt ¼ 2h and amax ¼ 2:6h. f f A summary of all these test cases with a Neumann-Robin condition is provided in Fig. 18 which gives the number of iterations needed to converge (with a maximum absolute error e ¼ 10
3 ) in terms of af . As predicted by the model problem, this function is a non monotone function which varies continuously, with an unstable zone for af > amax in perfect agreement with the theoretical f model (2nd row of Table 2).

Fig. 18. Convergence vs coupling coefficient. Case II – Low Biot number- NeumannRobin condition.

9.3. Comparative study A Neumann condition is not recommended for low Biot numbers. An instability zone is identified theoretically (Table 2) and confirmed by the CHT computations. As for large Biot numbers associated to the Dirichlet condition, an ‘‘optimal” coefficient exists, and the model shows that the zone between the optimal and the unstable coefficients ½2h; 2:6h is extremely narrow. Nevertheless, if a Neumann condition is still desired, it would be best to choose a coefficient largely on the ‘‘stable branch” ðaf  aopt f Þ,

Fig. 16. Convergence behavior for various coefficients. Case II – Low Biot number Neumann-Robin condition.

employing coupling coefficients clearly on this branch. If so, there is always a risk attached to slowing down the coupling process. Without any doubt, a Dirichlet condition is ideal for low Biot numbers as shown in the theoretical study (3rd row of Table 1),

602

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

which does not predict any instability, and amply supported by the CHT test cases in which the process is unconditionally stable. In summary, when the coupling coefficients are employed on the basis of the 1D model, no numerical problems were identified. 10. CHT computation with composite ramps 10.1. Objective Our goal here is to perform a reference test case combining composite ramps and to employ the numerical interface treatments presented in this paper. The geometry and operating conditions are displayed in Fig. 3. Two ‘‘probes” will be considered: – the temperature at the leading edge calculated on the upper side (F-S interface) – the temperature at the trailing edge on the upper side The input temperature level, f(t), imposed on the lower side (input signal) is given in Table 5. Fig. 19. Solid temperature profile at the F-S interface at low Biot number. In black: the imposed temperature profile f(t) at the lower side.

10.2. Temporal coupling instants The temporal coupling instants are those indicated in Table 5, namely 12 instants from 0.5 s to 12 s (ramp points). Each coupling period has been divided into 10 ‘‘solid” time steps. In this way, the maximum mesh Fourier number based on the solid grid size is around unity (see Eq. (11)). 10.3. Results The CHT computation presented in Fig. 19 for a full transient cycle is performed at small Biot numbers. As reported in this paper, a Dirichlet-Robin condition is well-suited and no relaxation parameter is needed (3rd row of Table 1). At small Biot numbers, the heat conduction inside the body is much faster than the heat convection away from its surface due to the small solid resistance (high conductance) or to the small distance the heat has to diffuse into the solid and temperature gradients are negligible inside of it. This trend is observed in Fig. 19 where the solid temperature profile is represented at the interface. The other CHT computation presented in Fig. 20 for the same full transient cycle is carried out at large Biot numbers. We know that a Neumann-Robin condition is then well-suited and that this time, relaxation parameters are not needed but recommended (1st row of Table 2). Thus, we have employed the optimal coefficient whose expression is given by (32). Note that the solid conductivities and heat capacities have been selected so as to obtain the same thermal diffusivity (see Table 3) at low and large Biot numbers. As a result the same solid time steps have been used (10 solid increments during a coupling period). At large Biot numbers, thermal gradients within the solid may become important, even though the material is a good conductor and more complicated heat transfer equations for ‘‘transient heat conduction” are required to describe the time-varying and nonspatially-uniform temperature field within the material body. It is exactly what is observed in Fig. 20 where the temperature along the upper face of the plate (fluid-solid interface) diverges dramatically from the boundary condition imposed on the lower part,

Fig. 20. Solid temperature profile at the interface at large Biot number. In black: the imposed temperature profile f(t) at the lower side.

especially at the leading edge where the thermal gradients are more significant. 11. Summary and conclusions Two new numerical interfacial models have been presented for transient CHT models. They have been derived from a model problem and two interface conditions have been analyzed. These fundamentally different transmission procedures do not necessarily compete and can be combined to work together for a more efficient overall algorithm. When the first method is unconditionally stable, the second one exhibits an upper stability bound. When the second

Table 5 Temperature values on the lower side T (K). Time (s) f(t)

0. 1000

0.5 1120

1.0 1120

1.5 1190

2.0 1190

2.5 1120

5.9 1120

6.0 1190

8.0 1190

8.1 1120

10.0 1120

11.0 1000

12.0 1000

603

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

method is unconditionally stable, the first one exhibits a lower stability bound. The numerical criteria establishing the nature and character of the most relevant interfacial treatment have been identified and mathematically expressed. The numerical studies presented in this paper illustrate and confirm the main theoretical results derived from the mathematical normal mode analysis. The stability bounds, the optimal coefficient and the general numerical behavior are largely supported by the CHT test cases. The trends predicted by the prototype problem are recovered and excellent convergence properties are observed in all cases. On the basis of this coupling prototype problem, an efficient and unified coupling methodology has been defined. This shows unequivocally that a coupling prototype greatly improves the prospect for fast convergent algorithms and provides a useful framework to address stability issues in realistic CHT problems. It appears that the root cause of convergence problems in CHT problems are undoubtedly due to the interfacial treatments. Energy transfer at the interface depends on many numerical and physical local factors. It is therefore fundamental to employ adaptive coefficients to optimize energy transmission, reflecting local conditions and stability considerations. Otherwise, many numerical problems may arise, such as a low convergence rate, large oscillations or even instabilities. This stresses that the interfacial treatment is the key point in a numerical CHT procedure. Of course, it should be kept in mind that a fixed-point iteration procedure exhibits only a linear convergence behavior. As a result, there is no doubt that this can be improved and acceleration techniques could be considered. Emphasis must be put, however, on the nature of the treatment at the interface. In this regard, the use of adaptive coupling coefficients always have a positive impact on the convergence process in a CHT computation, especially when a configuration exhibits many different features that cannot be taken into account efficiently even by ‘‘static” coupling parameters or by procedures routinely used in industry. Acknowledgments The authors would like to acknowledge the Direction Générale de l’Aviation Civile (DGAC) for their financial support. We also like to thank the ‘‘first reviewer” for his valuable comments and recommendations. Finally, thanks to Andrew Mayne for his thorough rereading of this article. Appendix A The 1D fluid-solid model is composed of the fluid domain Xf (steady state approach) and the solid domain Xs (unsteady approach). This discrete model is presented in Fig. A1. The stability analysis in this section is performed using general Robin-type interface conditions on either side of the fluid-solid interface, from which any other particular condition may be obtained.

A.1. Robin conditions on the solid side Let us consider that the fluid and solid domains are thermally coupled through a conventional staggered (CSS) algorithm. In the first step, in the temporal interval ½t n ; t nþ1 , with the sign convention of Fig. A1 (y denotes the direction normal to the interface, the unit vector, m, points towards the inward normal of the solid domain), a Robin condition is applied on the solid side ðj ¼ 0þÞ n nþ1 n n n qnþ1 0þ þ af T 0þ ¼ q0
þ af T 0


ðA-1Þ

A.2. Robin conditions on the fluid side In the same manner and with the same unit normal (this time, m is an outer-pointing normal of the fluid domain), a Robin condition is applied on the fluid side ðj ¼ 0
Þ nþ1 nþ1 nþ1 nþ1 qnþ1 T 0
¼ qnþ1 T 0þ 0

as 0þ
as

ðA-2Þ

where af and as are positive coupling coefficients. Note that even though only a fluid steady-state solution is sought, it is important to consider the fluid discretization to determine the spatial characteristic equation, as will be seen in Appendix B. The solid heat flux at the interface, q0þ , may be expressed as follows

q0þ ¼
K s ðT 1
T 0þ Þ

ðA-3Þ

where K s is the solid thermal conductance (K s ¼ ks =Dys ). The fluid heat flux at the interface, q0
, has a similar expression

q0
¼


kf ðT 0

T
1=2 Þ ¼
K f ðT 0

T
1=2 Þ Dyf =2

ðA-4Þ

The fluid domain runs until steady-state conditions are reached. After that, it is usual to express the local convective heat transfer as follows

q0
¼ hðT 
T 0
Þ where h is the heat transfer coefficient and T temperature.

ðA-5Þ 

is a reference

A.3. Combined interface condition Note that, at each coupling, the fluid solution is performed until steady state is achieved, via a finite-volume method, as indicated in the zone j 6 0 of Fig. A1. At convergence in the fluid side, the convective coefficient h is computed from the known heat flux q0
and a reference temperature T⁄, as expressed by (A-5). Combining (A-1) and (A-3), one obtains nþ1 n ðK s þ anf ÞT nþ1 þ anf T n0
0þ ¼ þq0
þ K s T 1

Substituting the definition of q0
(A-5) into (A-6), yields

Fig. A1. Discrete model.

ðA-6Þ

604

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605

T nþ1 0þ ¼

1  þ ðanf
hÞT n0
þ hT  ½þK s T nþ1 1 ðK s þ anf Þ

ðA-7Þ

In the same manner, combining (A-2) and (A-5), one obtains nþ1 nþ1 nþ1 ðh þ anþ1 ÞT nþ1 T 0þ þ hT s 0
¼
q0þ þ as



ðA-8Þ

Substituting the definition of q0þ (A-3) into (A-8), the interfacial temperature becomes

T nþ1 0
¼

1 nþ1 ½ðanþ1
K s ÞT nþ1 0þ þ K s T 1  þ constant ðh þ anþ1 Þ s s

ðA-9Þ

Substituting the definition of T nþ1 0þ (A-7) into (A-9), the above equation becomes

T nþ1 0
¼

ðaf þ as ÞK s ðh
af ÞðK s
as Þ n T nþ1 þ T þ constant ðh þ as ÞðK s þ af Þ 1 ðh þ as ÞðK s þ af Þ 0
ðA-10Þ

The ‘‘constant” is a function at the external BC ðj ¼
JÞ that plays no role in the stability analysis. In Eq. (A-10), the temporal index is omitted for the coupling coefficients af and as .

B.3. Amplification factor in the ‘‘steady” fluid domain In the fluid domain, the solution is considered only at steadystate. From (A-4), we obtain at time n+1:

! kf kf 
h T nþ1 T nþ1
hT 0
¼ Dyf =2 Dyf =2
1=2

Inserting the ‘fluid’ normal mode solutions (B-1), one obtains

j
¼ jf ¼ 

f

This stability analysis is similar to the standard Fourier stability method except that the Fourier analysis ignores boundary conditions and, as these may affect the stability, the theory of Godunov & Ryabenkii (G-R) is preferable. We introduce the normal mode solution by considering eigensolutions of the form

T nj ¼

zn
1 jsj ; j P 0 zn j

j f;

ðB-1Þ

j60

where z is the ‘‘temporal amplification factor” and amplification factor”.

j is the ‘‘spatial

In the solid domain, the solution is advanced in time using the following discrete equation (centrally discretized in space and treated with a forward implicit scheme)

B.4. Temporal amplification factor of the coupled problem Inserting the eigensolutions (B-1) into the single interface condition (A-10), the following amplification factor is readily obtained

z ¼ gðz; af ; as Þ ðaf þ as ÞK s ðh
af ÞðK s
as Þ js ðz; Ds Þ þ ðh þ as ÞðK s þ af Þ ðh þ as ÞðK s þ af Þ

n As z ¼ T nþ1 0
=T 0
, from (A-10), we obtain:

ðaf þ as ÞK s T nþ1 1 z¼ ðh þ as ÞðK s þ af Þ T n0


ðB-2Þ

z
1 Ds z



js þ 1 ¼ 0 j > 0

ðB-4Þ

There is one acceptable solution to the Eq. (B-4) given by the choice of the root with the minus sign. This topic has already been discussed in [6].

ðh
af ÞðK s
as Þ ðh þ as ÞðK s þ af Þ

ðB-8Þ

ðB-9Þ

Note that the solid amplification factor, js , is just apparently dependent on the fluid temperature T n0
. Indeed, combining (A-1) and (A-2), it is possible to extract T n0
n n nþ1 n n n n T n0
¼ ½qnþ1 0þ
q0þ þ af T 0þ þ as T 0þ =ðaf þ as Þ

ðB-10Þ

And finally js can be expressed with physical solid data only, as follows

References

where as is the solid thermal diffusivity,Dt is the time step and Dxs is the solid grid size. Inserting the ‘solid’ eigensolutions (B-1) into the temporal discrete solution (B-2) leads to a quadratic expression in the solid domain

þ

T nþ1 1 T n0


Ds ¼ as Dt=Dx2s

ðB-3Þ

!

and thus, comparing (B-7) and (B-8), we see that the spatial amplification factor in the solid domain takes the simple form

js ¼



ðB-7Þ

B.5. Amplification factor in the ‘‘transient” solid domain

This central implicit scheme is used from the interior point j = 1 and not at the interface j = 0+ (A-3). At this boundary, the general Robin interface condition (A-1) is implemented. The expression of the Fourier number Ds in Eq. (B-2) is given by

j2s
2 þ

f

and the corresponding stability condition is: jjf j > 1. This ‘fluid’ stability condition is discussed in Section 4.4 and in Section 7.7.

js ¼

B.2. Characteristic equation in the ‘‘transient” solid domain

nþ1 T nþ1
T nj ¼ Ds ðT nþ1 þ T nþ1 j>0 j jþ1
2T j j
1 Þ

ðB-6Þ

2

and the corresponding stability condition is jgj < 1. This stability condition of the coupled problem is discussed in Sections 4.3 and 4.4 and in Section 5.

B.1. Eigensolutions

(

1

1
2k =hDy

¼ Appendix B

ðB-5Þ

T nþ1 T nþ1 1 1 ¼ ðanf þ ans Þ nþ1 n n n T 0
q0þ
qn0þ þ anf T nþ1 0þ þ as T 0þ

ðB-11Þ

[1] M.B. Giles, Stability analysis of numerical interface conditions in fluidstructure thermal analysis, Int. J. Numer. Meth. Fluids 25 (1997) 421–436. [2] B. Roe, R. Jaiman, A. Haselbacher, P.H. Geubelle, Combined interface boundary method for coupled thermal simulations, Int. J. Numer. Meth. Fluids 57 (2008) 329–354. [3] 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. [4] V. Kazemi-Kamyab, 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. Meth. Fluids 74 (2014) 113–133. [5] J. Lindström, J. Nordström, A stable and high-order conjugate heat transfer problem, J. Comput. Phys. 229 (2010) 5440–5456. [6] M.-P. Errera, S. Chemin, Optimal solutions of numerical interface conditions in fluid-structure thermal analysis, J. Comput. Phys. 245 (2013) 431–455.

M.-P. Errera et al. / International Journal of Heat and Mass Transfer 110 (2017) 587–605 [7] Z. Sun, J. Chew, N. Hills, K. Volkov, C. Barnes, Efficient finite element analysis/computational fluid dynamics thermal coupling for engineering applications, J. Turbomach. 132 (2010). [8] V. Ganine, U. Javiya, N. Hills, J. Chew, Coupled fluid-structure transient thermal analysis of a gas turbine internal air system with multiple cavities, J. Eng. Gas Turbines Power 134 (2012). [9] Z. Sun, J. Chew, N. Hills, L. Lewis, C. Mabilat, Coupled aerothermomechanical simulation for a turbine disk through a full transient cycle, J. Turbomach. 134 (2012). [10] Z. Zhai, Q. Chen, Impact of determination of convective heat transfer coefficient on the coupled energy and CFD simulation for buildings, Proc. Build. Simul. Confer. 3 (2003) 1467–1474. [11] I. Beausoleil-Morrison, The adaptive coupling of computational fluid dynamics with whole-building thermal simulation, in: 7th International IBPSA Conference, Rio de Janeiro, Brazil, 2001. [12] Q. Chen, Z. Zhai, L. Wang, Computer modeling of multiscale fluid flow and heat and mass transfer in engineered spaces, Chem. Eng. Sci. 62 (2007) 3580–3588. [13] M.-P. Errera, B. Baqué, A quasi-dynamic procedure for coupled thermal simulations, Int. J. Numer. Meth. Fluids 72 (2013) 1183–1206. [14] D.J. Culler, J.J. Mc Namara, Studies on fluid-thermal-structural coupling for aerothermoelasticity in hypersonic flow, AIAA J. 48 (8) (2010) 1721–1738. [15] L. He, M.L.G. Oldfield, Unsteady conjugate heat transfer modeling, J. Turbomach. 133 (3) (2011) 031022. [16] C. Shen, F.-X. Sun, X.-L. Xia, Analysis on transient conjugate heat transfer in gap-cavity-gap structure heated by high speed airflow, Int. J. Heat Mass Transf. 67 (2013) 1030–1038. [17] T. Verstraete, S. Scholl, Stability analysis of partitioned methods for predicting conjugate heat transfer, Int. J. Heat Mass Transf. 101 (2016) 852–869. [18] T. Verstraete, Multidisciplinary Turbomachinery Component Optimization Considering Performance, Stress, and Internal Heat Transfer PhD thesis, Universiteit Gent, 2008. [19] M.-P. Errera, F. Duchaine, Comparative study of coupling coefficients in Dirichlet-Robin procedure for fluid-structure aerothermal simulations, J. Comput. Phys. 312 (2016) 218–234. [20] G. Gimenez, M.P. Errera, D. Baillis, Y. Smith, F. Pardo, A coupling numerical methodology for weakly transient conjugate heat transfer problems, Int. J. Heat Mass Transf. 97 (2016) 975–989. [21] S.K. Godunov, V.S. Ryabenkii, The Theory of Difference Schemes, An introduction, Amsterdam, North-Holland, 1964. [22] H-O. Kreiss, Stability theory for difference approximations of mixed initial boundary value problems, I., Math. Comp. 22 (1968) 703–714. [23] S. Osher, Systems of difference equations with general homogeneous boundary conditions, Trans. Am. Math. Soc. 137 (1969) 177–201. [24] B. Gustaffson, H.-O. Kreiss, A. Sundström, Stability theory of difference approximations for mixed initial boundary value problems, Math. Comput. 26 (1972) 649–686. [25] B. Gustafsson, The Godunov-Ryabenkii Condition: the Beginning of a New Stability Theory Technical report 1999-014, Department of Information Technology (1999), 1999. [26] R. El Khoury, M. Errera, K. El Khoury, M. Nemer, Efficiency of steady state coupling schemes for the treatment of Fluid-Structure thermal interactions, Int. J. Therm. Sci. 115 (2017) 225–235.

605

[27] M.-P. Errera, G. Turpin, Temporal multiscale strategies for conjugate heat transfer problems, J. Coupled Syst. Multiscale Dyn. 1 (2013) 89–98. [28] M.-P. Errera, Efficient coupling procedures in steady and unsteady thermal analysis, in: Proceedings of the 15th International Heat Transfer Conference, IHTC-15, August 10-15 (2014) Kyoto, Japan, 2014. [29] B. Hübner, D. Dinkler, A simultaneous solution procedure for strong interactions of generalized Newtonian fluids and viscoelastic solids at large strains, Int. J. Numer. Meth. Eng. 64 (2005) 920–939. [30] W.G. Dettmer, D. Peric´, A fully implicit computational strategy for strongly coupled fluid-solid interaction, Arch. Comput. Methods Eng. 14–3 (2007) 205– 247. [31] C.A. Felippa, K.C. Park, Staggered transient analysis procedures for coupled dynamic systems formulation, Comput. Methods Appl. Mech. Eng. 24 (1980) 61–111. [32] S. Piperno, C. Farhat, B. Larrouturou, Partitioned procedures for the transient solution of coupled aeroelastic problems, Comput. Methods Appl. Mech. Eng. 124 (1995) 79–112. [33] H.G. Matthies, J. Steindorf, Partitioned strong coupling algorithms for fluidstructure interaction, Comput. Struct. 81 (8–11) (2003) 805–812. [34] J.M. Ortega, W.C. Rheinbolt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970. [35] K.C. Park, A. Felippa, Partitioned analysis of coupled systems, in: T. Belytschko, H.T.J.R. (Eds.), Computational Methods for Transient Analysis, North-Holland Pub. Co., 1983, pp. 157–219. [36] J. Degroote, Partitioned simulation of fluid-structure interaction, Arch. Comput. Methods Eng. 20–3 (2013) 185–239. [37] C. Farhat, M. Lesoinne, On the accuracy, stability, and performance of the solution of three-dimensional nonlinear transient aeroelastic problems by partitioned procedures, AIAA/ASME/ASCE/AHC/Struct., Struct. Dyn. Mater. Conf. 1 (1996) 629–641. [38] B. Baqué, M.-P. Errera, A. Roos, F. Feyel, Simulation of conjugate heat transfer via a temporal multiscale approach, J. Multiscale Comput. Eng. 11 (4) (2013) 333–345. [39] P. Causin, J.-F. Gerbeau, F. Nobile, Added-mass effect in the design of partitioned algorithms for fluid-structure problems, Comput. Methods Appl. Mech. Eng. 194 (2005) 4506–4527. [40] U. Küttler, W. Wall, Fixed-point fluid-structure interaction solvers with dynamic relaxation, Comput. Mech. 43 (1) (2008) 61–72. [41] M.-P. Errera, R.R. El Khoury, A numerical efficient approach for coupled aerothermal simulations in fluid-structure systems, in: 51th International Conference on Applied Aerodynamics, 2016, 4-6 April 2016, Strasboug France. [42] http://www.zset-software.com. [43] F.R. Menter, Zonal two-equation k-x turbulence model for aerodynamic flows, AIAA Paper 1993–2906, 1993. [44] F.R. Menter, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA J. 32 (8) (1994) 269–289. [45] F.R. Menter, M. Kuntz, R. Langtry, Ten years of industrial experience with the SST turbulence model, Turbulence, Heat and Mass Transfer 4 (2003) Begell House Inc. [46] http://www.cerfacs.fr/globc/PALM_WEB/.