Efficient monolithic projection method with staggered time discretization for natural convection problems

International Journal of Heat and Mass Transfer 144 (2019) 118677 Contents lists available at ScienceDirect International Journal of Heat and Mass T...

Download PDF

3MB Sizes 1 Downloads 39 Views

Report

PDF Reader
Full Text

International Journal of Heat and Mass Transfer 144 (2019) 118677

Contents lists available at ScienceDirect

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

Efficient monolithic projection method with staggered time discretization for natural convection problems Xiaomin Pan, Ki-Ha Kim, Jung-Il Choi ⇑ Department of Computational Science and Engineering, Yonsei University, Seoul 03722, Republic of Korea

a r t i c l e

i n f o

Article history: Received 4 March 2019 Received in revised form 1 September 2019 Accepted 1 September 2019

Keywords: Monolithic projection method Staggered time discretization Numerical stability Temporal and spatial second-order accuracy

a b s t r a c t For a more efficient algorithm, we introduce staggered time discretization to improve the previous method (Pan et al., 2017), fully decoupled monolithic projection method with one Poisson equation (FDMPM-1P), to solve time-dependent natural convection problems. The momentum and energy equations are discretized using the Crank–Nicolson scheme at the staggered time grids, in which temperature and pressure fields are evaluated at half-integer time levels (n þ 12), while the velocity fields are evaluated at integer time levels (n þ 1). Numerical simulations of two-dimensional (2D) Rayleigh–Bénard convection show that the proposed method is more computationally efficient and stable than FDMPM-1P, while preserving the second-order spatial and temporal accuracy. Further, the proposed method provides accurate predictions of 2D Rayleigh–Bénard convection under different thermal boundary conditions for a Rayleigh number up to 1010 , three-dimensional turbulent Rayleigh–Bénard convection in the range of 1 105 —2 107 in horizontal periodic domain, and three-dimensional turbulent Rayleigh–Bénard convection in the range of 1 106 —1 107 in bounded domain. Finally, we theoretically confirmed that the proposed method is second-order in time and it is more stable than FDMPM-1P for small Ra. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction The phenomenon of natural convection has attracted considerable attention [1–10] owing to its wide range of engineering applications. To fundamentally understand and develop efficient heat transfer systems, it is essential to investigate the temperature and fluid distributions occurring in natural convection problems, which involve not only incompressible and nonlinear momentum equations but also strong coupling between incompressible flows and heat transfers. This coupling is based on the fact that the thermo-fluid flow is driven by a buoyancy force depending on the temperature distribution, meanwhile, the temperature is convected by the background fluid flow. Yee’s method [11] used staggered time discretization for nonstationary Maxwell’s equations. Inspired by this method, to obtain a more efficient method, we improved the fully decoupled monolithic projection method with one Poisson equation (FDMPM-1P) [8] based on staggered time discretization to solve the time-dependent natural convection problem. To the best of our knowledge, this might be the first attempt to use the staggered time discretization in projection methods for thermo-fluid ⇑ Corresponding author. E-mail addresses: [email protected] (X. Pan), [email protected] (K.-H. Kim), [email protected] (J.-I. Choi). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118677 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

problems. The buoyancy, nonlinear convection, and linear diffusion terms are discretized using the Crank–Nicolson scheme, while the velocity and temperature fields are evaluated at staggered time intervals, unlike [7,8]. This staggered time discretization ensures that the nonlinear convection part in the energy equation is naturally linearized and enables the immediate update of the buoyancy term once the temperature is obtained. In this study, we compute the temperature field at half-integer time levels (n þ 12) by solving a series of tri-diagonal linear systems, and then, we solve the incompressible momentum equations using a momentum solver (e.g., velocity-components decoupled projection method (VDPM) [12]) with the updated temperature field. The above-mentioned processes produce a non-iterative monolithic velocity-components decoupled projection method with staggered time discretization (MVDPM-STD). Various numerical simulations are conducted for evaluation. Numerical tests of two-dimensional (2D) Rayleigh– Bénard convection and 2D periodic forced flow are carried out to verify the numerical performance. 2D periodic flows with fixedtemperature and fixed-flux boundary conditions are investigated to examine the ability of the proposed method to deal with problems under different boundary conditions. Finally, the proposed method is adopted to solve three-dimensional (3D) turbulent Rayleigh–Bénard convection problems in both periodic and bounded domains. Theoretical analysis of the stable property and temporal accuracy is also presented.

2

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

Nomenclature

Greek symbols thermal diffusivity (m2/s)

a

Dimensional variables L characteristic length (m) g gravitational acceleration (m/s2) hc cold wall temperature (K) hot wall temperature (K) hh

thermal expansion coefficient K1 3 density kg=m kinematic viscosity m2 =s wavenumber computational domain boundary of domain X

b

q

m j

X C

Dimensionless variables D discrete divergence operators G discrete gradient operator L discrete Laplacian operator Nu local Nusselt number Ra Raleigh number ¼ gbðhh hc ÞL3 =ðmaÞ Pr Prandtl number ð¼ m=aÞ Q local heat flux t time Dt time step h grid spacing ui velocity components ði ¼ 1; 2; 3Þ p pressure h temperature xi Cartesian coordinates ði ¼ 1; 2; 3Þ

Superscripts and Subscripts it time averaged values D horizontal plane-time averaged values iA;t D iV;t volume-time averaged values D ix2 ;t vertical-time averaged values 0

fluctuation values root-mean-square values values on the vertical walls complex conjugate Fourier transform dimensional values

rms w ^

2. Proposed monolithic projection method with staggered time discretization 2.1. Governing equations Under the Oberbeck–Boussinesq approximation, we consider the natural convection of incompressible fluids in a bounded domain, X 2 Rd (d ¼ 2or3), with a sufficiently smooth boundary, C, over a finite time interval 0; T f :

rﬃﬃﬃﬃﬃﬃ @u Pr 2 r u þ FðhÞ; þ u ru ¼ rp þ @t Ra r u ¼ 0; @h 1 þ u rh ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ r2 h; @t RaPr

ð1Þ ð2Þ ð3Þ

where u; p, and h are the velocity vector, pressure, and temperature, respectively. FðhÞ denotes an external force that is proportional to the temperature variation. Here, we consider a more general hypothesis on F.

sion coefficient of the fluid, kinematic viscosity, thermal diffusivity, density, and gravitational acceleration, respectively. The variables with symbol are dimensional. Two dimensionless parameters (i.e., Prandtl number (Pr) and Raleigh number (Ra)) characterize natural convection problems. 2.2. Numerical discretization To discretize Eqs. (1)–(3) in time, we apply the Crank–Nicolson scheme with a staggered time discretization algorithm [11,13]:

unþ1 un 1 nþ1 Gunþ1 þ un Gun þ u 2 Dt rﬃﬃﬃﬃﬃﬃ 1 1 Pr nþ1 1 nþ1 Lu þ Lun þ F hnþ2 þ mbcCN2 ; ¼ Gpnþ2 þ 2 Ra nþ1

Dunþ1 ¼ cbcCN ;

0 F ðh1 Þ 6 C F 1 0

for all h1 2 H ðXÞ, where F ðh1 Þ ¼ 1

dFðhÞ j . dh h¼h1

We employ the following non-dimensionalization in Eqs. (1)– (3):

xi ¼

e xi ; L

ui ¼

ei u ; U0

t¼

et ; L=U 0

m gbðhh hc ÞL3 ; Pr ¼ ; Ra ¼ a ma

p¼

e p U 20

q

;

h¼

e h hc ; hh hc

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ where xi ; L; ui ; U 0 ¼ gbLðhh hc Þ; hh , and hc represent the Cartesian coordinates, characteristic length, velocity components, characteristic velocity, hot wall temperature, and cold wall temperature, respectively. In addition, b; m; a; q, and g denote the thermal expan-

ð5Þ

1 1 1 hnþ2 hn2 1 n nþ12 1 n þ u Gh þ Ghn2 ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Lhnþ2 þ Lhn2 þ ebcCN ; 2 Dt 2 RaPr ð6Þ 1

Hypothesis 1. We assume that F satisfies

ð4Þ

1

where G; L, and D represent the discrete gradient, Laplacian, and divergence operators, respectively. These discrete spatial operators are evaluated using the second-order central difference scheme on a staggered spatial grid in discrete domain, in which pressure and temperature are stored at the cell centers and velocities are located at the cell interfaces. Here, Dt is the time step, h is the uniform grid spacing, and the superscript n denotes the nth time level. Note that the boundary condition of velocity has been incorporated into nþ1

nþ1

mbcCN2 and cbcCN , and the boundary condition of temperature n has been incorporated into ebcCN . Because the staggered time discretization guarantees the linearization of the nonlinear convection part in the energy equation, we only need to apply the linearization technique [14,15] to the convection term in discretized momentum Eq. (4):

3

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

1 n 1 nþ1 u Gunþ1 þ un Gun ¼ N unþ1 þ unþ1 Nn þ O Dt2 ; 2 2

I are identity matrices, I d1 ¼ ð1; ; 1Þtr , and is the Kronecker |ﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄ}

ð7Þ n

where N :¼ un G and Nn :¼ Gun only contain the information at time level n and second-order temporal accuracy is maintained. According to linearization (7), by ignoring the higher-order terms, Eqs. (4)–(6) can be approximated as

nþ1

ð8Þ

;

ð9Þ

1 1 1 hnþ2 hn2 1 n nþ12 n N h þ N hn2 þ 2 Dt 1 1 1 n ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Lhnþ2 þ Lhn2 þ ebc ; 2 RaPr

ð10Þ

where the boundary conditions have been incorporated into nþ12

mbc

; cbc

nþ1

1

mation ðAn Þ

1

DtI, which enables us to rewrite Eq. (14) as

1 n12 A O O h C C B B C C B B n C C ¼ B F An B R O A A @ @ nþ1 1 cbc Dun O D DðAn Þ G 0 10 nþ1 1 I O O dh 2 C B CB n 1 CB B @ O I ðA Þ G A@ dunþ1 C A 1 O O I dpnþ2 10 0 1 n12 I O O A O O h CB B C CB B B F An C@ O I DtG C O A A @ O O I O D DtDG 0 nþ1 1 dh 2 B C B @ dunþ1 C A: 0

1

n12

0

Rh

unþ1 un 1 n nþ1 þ unþ1 Nn þ N u 2 Dt rﬃﬃﬃﬃﬃﬃ 1 1 Pr nþ1 1 nþ1 Lu þ Lun þ F hnþ2 þ mbc 2 ; ¼ Gpnþ2 þ 2 Ra Dunþ1 ¼ cbc

d

product. Further, to deal with the coupling between unþ1 and pnþ2 , we apply block lower-upper (LU) decomposition with the approxi-

n

, and ebc . Here, we have

1 nþ1 nþ1 u Gunþ1 þ un Gun þ mbcCN2 2 1 n nþ1 ¼ N unþ1 þ unþ1 Nn þ mbc 2 þ O Dt2 2

1

dpnþ2 ð15Þ

ð11Þ We can obtain the following procedure from system (15):

and nþ1

ebcCN2 ¼ ebc

nþ12

:

ð12Þ nþ1

Because both cbcCN and cbc

nþ1

depend on the discrete operator D

and the boundary condition of velocities at the ðn þ 1Þth time level, we have nþ1

cbcCN ¼ cbc

nþ1

:

ð13Þ

n1

n1

1

A h 2 dhnþ2 ¼ Rh 2 ; A du n

;nþ1

nþ12

DtDGdp du

nþ1

ð16Þ

1 ¼ R þ F dhnþ2 ; n

¼ cbc

;nþ1

¼ du

nþ1

ð17Þ

Dun Ddu;nþ1 ; nþ12

DtGdp

ð18Þ

;

ð19Þ

Similar to Appendix A in [16], we present the explicit forms of nþ12

mbc

; cbc

nþ1

where du;nþ1 is intermediate variable of dunþ1 . For simplicity, we refer to procedure (16)–(19) as monolithic velocity-components coupled projection method with staggered time discretization (MVCPM-STD).

n

, and ebc in Appendix A.

2.3. Decoupling process of variables While maintaining the second-order temporal accuracy, we can rewrite discrete formulations (8)–(10) in a coupled matrix-vector nþ12

form by introducing increments of variables, dh

nþ1

, du

, and

1

dpnþ2 :

1 10 1 n1 dhnþ2 A 2 O O C CB B h B F An G CB dunþ1 C A A@ @ 0

O 0

1

dpnþ2

D O

1 1 1 0 n1 n n1 þ ebc A h 2 hn2 Rh 2 C B 1 C B C B nþ12 n n n n1 C n C; B ¼B Rn A B r þ mbc þ F h 2 A u Gp 2 C :¼ @ A @ nþ1 n nþ1 cbc Du n cbc Du n12

rh

ð14Þ n12

1

n n 1 where A h ¼ D1t I þ 12 N 2pﬃﬃﬃﬃﬃﬃﬃ L, An ¼ D1t I þ I d1 2 N RaPr qﬃﬃﬃﬃ 1 1 n 1 1 n n 1 Pr 1 LÞ þ 12 Nn , r h 2 ¼ hDt2 12 N hn2 þ 2pﬃﬃﬃﬃﬃﬃﬃ Lhn2 , and 2 Ra RaPr qﬃﬃﬃﬃ n Pr rn ¼ uDt þ 12 Ra Lun . This results in a monolithic approach for solving nþ1 1 du dhnþ2 and . Here, O denotes the generic zero matrix, I and 1 dpnþ2

1

To obtain dhnþ2 from Eq. (16), similar to Eq. (31) in [8], we apply n12

approximate factorization to the coefficient matrix A h

in each

1

direction. Once dhnþ2 is solved, the updated buoyancy term along with the notation of Rn will be explicitly given. Further, Eqs. (17)–(19) can be solved using VDPM [12], in which the approximate block LU decomposition along with approximate factorization is introduced. Because we use the second-order central difference scheme in space, the inversions of the tri-diagonal matrices must be obtained instead of the inversions of the large n12

sparse matrices A h

and An (iterative methods are usually used 1

for obtaining inversions) for obtaining dhnþ2 and du;nþ1 . Avoiding the iterative procedures enables us significantly accelerate the calculation while maintaining the second-order temporal accuracy. This leads to the MVDPM-STD, where no iterative procedure is 1

required for obtaining dhnþ2 and du;nþ1 . It is worthy to note that compared with FDMPM-1P [8], the proposed algorithm requires no additional operations related to the buoyancy part in momentum Eq. (17); further, the procedure of the temperature field need not to be updated. This enables faster computation than FDMPM1P. Moreover, in Appendix B, we theoretically confirmed that MVDPM-STD has slightly better numerical stability than FDMPM1P, while maintaining second-order temporal accuracy.

4

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

3. Numerical simulations We perform numerical simulations of 2D and 3D thermal convection problems to validate the proposed algorithm, verify its numerical stability and computational efficiency, investigate its temporal and spatial accuracy, demonstrate its ability to deal with different boundary condition problems, and confirm its 3D capability. Here, we consider a linear form of the buoyancy term, i.e., FðhÞ ¼ ð0; hÞT for 2D cases and FðhÞ ¼ ð0; h; 0ÞT for 3D case. Five different problems were considered: 2D Rayleigh–Bénard convection; 2D periodic forced flow, which has an analytical solution; 2D periodic flows with different temperature boundary conditions (i.e., fixed temperature and fixed heat flux); 3D periodic Rayleigh– Bénard convection; 3D bounded Rayleigh–Bénard convection. As a reference for the numerical performance of the proposed MVDPM-STD, we implemented FDMPM-1P and semi-implicit projection method (SIPM) [17]. SIPM is a classical projection method that discretizes linear diffusion terms using the Crank–Nicolson scheme and buoyancy and nonlinear convection terms using the explicit Adams–Bashforth scheme. In addition, to solve the Poisson equation, we introduce a multi-grid method with a convergence criterion in terms of the relative l -residual as 1010 for 2D Rayleigh–Bénard convection, while a Fourier diagonalized Poisson solver is employed to solve the periodic flows. All the simulations were performed on an Intel Xeon E5430 (2.66 GHz) CPU. 2

3.1. Validation and stability: 2D Rayleigh–Bénard convection To validate the proposed algorithm, we considered a 2D Rayleigh–Bénard convection problem presented in [8], which models the heat transfer occurring via a fluid between two horizontal plates in an enclosure, X ¼ ½0; L ½0; H . In the current simulation, we set L ¼ H ¼ 1 with uniform grid spacing (h ¼ 1=128 for Ra ¼ 105 and 106 and h ¼ 1=256 for Ra ¼ 2 106 ) for both x1 direction and x2 -direction. The initial static state with p ¼ 0 and h ¼ 0:5ðhh þ hc Þ was heated from the bottom wall with hh ¼ 0:5 and cooled at the upper wall with hc ¼ 0:5, whereas the two side walls were thermally insulated. A no-slip velocity boundary condition was applied to all the walls of the enclosure. We stop the calculation when the numerical solution reaches a steady state: sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 1 X nþ1 1 X nþ12 n1 7 n #i;j #i;j 2 6 107 vi;j vi;j 6 10 and N1 N2 N1 N2

where v refers to u1 or u2 and # refers to p or h. Here, N1 and N2 denote the total number of grid intervals in x1 - and x2 -directions, respectively. For comparison with the FDMPM-1P proposed in [8], we considered a Boussinesq fluid (air) with Pr ¼ 0:71 and a Rayleigh number ranging from 105 to 2 106 ; all the calculations were performed until the steady state was achieved. Table 1 lists Dtmax , the Courant-Friedrichs-Lewy (CFL) number, the averaged computation time, T step , for one time step, and the relative computational cost, R, of obtaining steady-state solutions for MVDPM-STD and FDMPM-1P compared to SIPM. Here, Dt max is estimated such that for the largest Dt, the numerical solutions do not blow up. The CFL number is defined as CFL ¼ UmaxhDtmax , where U max ¼ max ðju1 j þ ju2 jÞ in a steady state. A bisection method is ðx1 ;x2 Þ2X

adopted to find the approximate value of Dt max . In addition, R is defined as the ratio of T step =Dtmax for MVDPM-STD (or FDMPM1P) to T step =Dtmax for SIPM. From Table 1, we can see that MVDPM-STD has a much larger Dtmax than SIPM for all three cases. MVDPM-STD is more stable than FDMPM-1P for small Ra (eg., Ra ¼ 105 ), while its stability is comparable to that of FDMPM-1P as Ra increases. This is consistent with the conclusion from the stability analysis in Appendix B.2. Moreover, we found that T step of MVDPM-STD is comparable to that of SIPM and lesser than that of FDMPM-1P. This is because FDMPM-1P requires additional operations for obtaining intermediate velocities as well as additional computation in the temperature-correction step [8]. Accordingly, MVDPM-STD is more efficient than FDMPM-1P and SIPM. In particular, MVDPM-STD is more than two times faster than FDMPM-1P for Ra ¼ 105 , while it is more than twenty times faster than SIPM for Ra ¼ 106 and 2 106 . Overall, the proposed MVDPM-STD outperforms FDMPM-1P and SIPM in terms of stability and computational efficiency. Table 2 compares the characteristic values of MVDPM-STD, SIPM, and FDMPM-1P with Dt listed in Table 1. Table 2 includes the maximum horizontal velocity u1max , maximum vertical velocity u2 max over the computational domain, and the averaged Nusselt R 1 @h number, Nu, on the hot wall (x2 ¼ 0). Here, Nu ¼ 0 @x j dx1 . x ¼0 2 2 From Table 2, it can be seen that the characteristic values of the proposed MVDPM-STD are in good agreement with those of SIPM and FDMPM-1P. Thus, MVDPM-STD provides accurate predictions of steady-state solutions with considerably large Dt.

Table 1 Maximum time step Dtmax and CFL number allowed for stable computation in SIPM, FDMPM-1P, and MVDPM-STD at different values of Ra with uniform grid spacing h. The computational efficiencies of the three schemes are presented. Note that Dt is non-dimensional in the setting, where the domain length is 1 and the temperature difference is 1. T step ðsÞ is the averaged computation time for one time step, and R is defined as the ratio of T step =Dtmax for MVDPM-STD (or FDMPM-1P) to T step =Dt max for SIPM. 105 ðh ¼ 1=128Þ

Ra

106 ðh ¼ 1=128Þ

2 106 ðh ¼ 1=256Þ

Dtmax

CFL

T step

R

Dtmax

CFL

T step

R

Dtmax

CFL

T step

R

SIPM

3:11 102

1:79

0:028

1:000

1:25 102

0:88

0:027

1:000

4:45 103

0:66

0:117

1:000

FDMPM 1P

3:14 101

18:11

0:041

0:143

2:78 101

19:50

0:041

0:068

1:71 101

25:56

0:177

0:039

MVDPM STD

4:77 101

27:44

0:027

0:063

3:07 101

21:50

0:029

0:044

1:77 101

26:47

0:125

0:027

Table 2 Comparison of characteristic values obtained from MVDPM-STD with benchmark solutions obtained from SIPM and FDMPM-1P [8] for different Ra with uniform grid spacing h and Dt listed in Table 1. Ra

u1 max u2 max Nu

105

106

2 106

MVDPM-STD

SIPM [8]

FDMPM-1P [8]

MVDPM-STD

SIPM [8]

FDMPM-1P [8]

MVDPM-STD

SIPM

FDMPM-1P

0:344 0:376 3:914

0:344 0:376 3:912

0:344 0:376 3:913

0:368 0:404 6:320

0:368 0:404 6:309

0:368 0:404 6:311

0:356 0:398 7:029

0:355 0:398 7:022

0:356 0:398 7:027

5

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

Fig. 1. Errors of u; p, and h from MVDPM-STD for ðaÞ h ¼ 1=2048 and ðbÞ Dt ¼ 1=4096at t ¼ 0:5.

3.2. Temporal and spatial accuracy: 2D periodic forced Rayleigh– Bénard convection A 2D periodic forced flow, which is an unsteady-state flow with an analytical solution occurring in X ¼ ½0; 1 ½0; 1 , was considered to investigate the temporal and spatial accuracy of the proposed MVDPM-STD. Without loss of completeness, all the numerical settings can be referred to the Section 4.2 in [8]. Fig. 1 displays the errors of u; p, and h obtained using the proposed MVDPM-STD at t ¼ 0:5. Fig. 1 (a) displays the results with a fixed grid spacing (h ¼ 1=2048) along with various computational time steps, and it confirms that second-order temporal accuracy is achieved for all variables. Furthermore, we present the results with a sufficiently small computational time step (Dt ¼ 1=4096) along with various uniform grid spacings in Fig. 1 (b), which shows that second-order spatial accuracy is obtained for all variables. Thus, linearizations, approximate LU decompositions, and the approximate factorization technique applied in the proposed MVDPM-STD do not degrade the second-order temporal and spatial accuracy. This is consistent with the theoretical discussion on the second-order temporal accuracy of MVDPM-STD in Appendix B.1. 3.3. Validations on different temperature boundary conditions: 2D periodic Rayleigh–Bénard convection We obtained the results of high-resolution direct numerical simulation (DNS) of 2D Rayleigh–Bénard convection for Ra up to 1010 on the basis of constant-temperature and constant-flux problems governed by Eqs. (1)–(3) with Ra replaced by Ra . Note that, as a control parameter, Ra is the measure of the imposed thermal force whose definition depends on the thermal boundary conditions. The flows are simulated in X ¼ ½0; 2 ½0; 1 . Periodic boundary conditions are considered for all variables in the horizontal direction and no-slip boundary conditions are implemented at top and bottom walls. For fixed temperature at top and bottom walls, Ra ¼ Ra with temperature boundary conditions hjx2 ¼0 ¼ 1 and hjx2 ¼1 ¼ 0; for the case with fixed vertical flux, Ra ¼

RahNuiX;t with

@h @x2

@h ðx2 ¼ 0Þ ¼ @x ðx2 ¼ 1Þ ¼ 1. Here, to evaluate 2

the two scenarios, we introduce the space-time averaged Nusselt pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ number, hNuiX;t , i.e., hNuiX;t ¼ 1 þ PrRahu2 hiX;t and hNuiX;t ¼ 1 pﬃﬃﬃﬃﬃﬃﬃﬃ for the fixed-temperature and fixed-flux cases, 1 PrRa hu2 hiX;t

respectively, where hiX;t indicates an average over the entire domain X and time [2]. All the results reported here are for Pr ¼ 1. Following the same procedure as in [2], simulations of fixed temperature and fixed flux were carried out using the proposed MVDPM-STD. Fig. 2 shows Ra hNuiX;t for two different boundary conditions obtained from MVDPM-STD. It is observed that the present results are in good agreement with those presented in [2]. We found that for smaller Ra < 106 , thermal convection in the fixedflux case is stronger than that in the fixed-temperature case, while for larger Ra > 106 , thermal convection in the two cases is comparable. For further details, we report another quantitative comparison, i.e., horizontally and temporally averaged temperature profiles for fixed-flux and fixed-temperature thermal convection with the DNS results in Fig. 3. The proposed MVDPM-STD provided excellent temperature distribution close to the hot bottom wall compared with those in [2]. In addition, the temperature profiles for the two different boundary cases become indistinguishable as Ra increases [2]. All these simulations demonstrate the ability of the proposed MVDPM-STD to deal with problems under different boundary conditions for a wide range of Ra up to Ra ¼ 1:05 1010 . 3.4. 3D Rayleigh–Bénard convection We solved the 3D Rayleigh–Bénard convection problem in both periodic and bounded domains to confirm that the proposed algorithm can reasonably solve 3D turbulent problems with different boundary conditions. 3.4.1. 3D Rayleigh–Bénard convection in periodic domain For 3D Rayleigh–Bénard convection in periodic domain, we refer the reader to Section 4.3 in [8] for governing equations, detailed numerical settings, and definitions. We present the validation of MVDPM-STD in Table 3 by comparing the Péclet number (Pe) and Nusselt number with DNS predictions in [8] and the Ra-scaling laws in [1]. The Péclet number Pe ui;rms ¼ ui;rms L2 =a is taken on the basis of root-mean-square v* ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ + u D E 2 u , where hiA;t denotes velocity ui;rms ðx2 Þ ¼ t ui ui A;t

A;t

the average over the horizontal plane and time. The Nusselt number Nuw ðtÞ is defined by the normalized wall-normal derivative

6

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

Fig. 2. Comparison of heat transport data for fixed-temperature and fixed-flux simulations using MVDPM-STD with the results in [2].

Fig. 3. Comparisons of horizontally and temporally averaged temperature profiles (from top to bottom) obtained from MVDPM-STD with the DNS results in [2]: Ra ¼ 1:03 104 ; 1:03 105 ; 1:21 106 ; 1:22 108 , and 1:05 1010 . ðaÞ Fixed-temperature and ðbÞ fixed-flux boundary conditions.

D

E

a @x@h2 jx2 ¼1 A Nuw ðtÞ ¼ aðhh hc Þ=L2

ð20Þ

at the hot wall x2 ¼ 1, and NuX ðt Þ is obtained by

NuX ðt Þ ¼

< Q >V

aðhh hc Þ=L2

;

ð21Þ

where iV denotes the average over the entire domain. The nonuniformity of heat transfer over time can be measured in terms of the standard deviation

rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ D 2 E Nuw ðt Þ Nuw and t rD ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 E ¼ NuX ðtÞ NuX

rðNuw Þ ¼

rðNuX Þ

t

where Nuw ¼ hNuw ðt Þit , NuX ¼ hNuX ðt Þit , and it denotes the average c are over time. In Table 3, Pe uc1;rms ; Pe uw 1;rms , and Pe u2;rms

defined from the center of u1;rms , the peak of u1;rms near the wall, and the center of u2;rms , respectively. We found that the results obtained by the proposed method are in excellent agreement with the FDMPM-1P predictions and the Ra-scaling laws for all Ra. In addition, consistent with [18], the standard deviations of Nusselt numbers of both FDMPM-1P and MVDPM-STD monotonically increase with an increase in Ra, indicating that the heat transfer uniformity degrades. In contrast, for each Ra; rðNuX Þ is considerably greater than rðNuw Þ, especially for larger Ra. This is because a larger Ra promotes an increase in velocity, which strengthens convective heat transfer. Then, convective heat transfer becomes more complex, which significantly increases the standard deviation of the convective Nusselt number (NuX ). Further, we confirmed that the maximum CFL number allowed for stable computation in MVDPM-STD ranges from 5:0 to 7:0, which is comparable to that of FDMPM-1P (3:9 to 6:3) and larger than that in [1] (less than 1:5). Thus, the proposed MVDPM-STD reasonably simulates 3D turbulent Rayleigh–Bénard convection with considerably large Dt.

7

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

Table 3 c Comparison of characteristic values obtained from MVDPM-STD with values from FDMPM-1P [8] and/or Ra-scaling in [1] (Pe uc1;rms ¼ 0:074Ra0:52 and Pe uw 1;rms and Pe u2;rms 0:46 with 0:27Ra ) for different Ra with different meshes. Ra

Pe uc1;rms

105

32:35

29:46

32:12

51:31

4:37

5:06

4:34

0:14

0:14

4:35

4:60

4:33

0:10

0:10

106

106:78

97:55

105:78 142:65 155:37 141:37 151:26 155:37 150:73

8:48

8:49

8:40

0:26

0:26

8:44

8:58

8:39

0:15

0:15

12:91

13:46

0:43

0:42

13:53

13:04

13:44

0:18

0:18

19:00

20:31

0:72

0:72

20:39

19:64

20:29

0:22

0:21

Pres.

Pe uw 1;rms

Ref. [1] Ref. [8]

Pres. 45:53

Pe uc2;rms

Ref. [1] Ref. [8] 53:87

45:21

Pres. 51:33

53:87

We also compare some statistically steady-state profiles obtained from the proposed MVDPM-STD with the DNS predictions [1,8] for Ra ¼ 2 107 . Fig. 4 (a) shows profiles of Peðu1;rms Þ and Peðu2;rms Þ, and Fig. 4 (b) shows profiles of the horizontally timeaveraged root-mean-square temperature hrms and horizontally time-averaged temperature hhiA;t . All the profiles in Fig. 4 are in good agreement with the previous results [1,8]. As shown in Fig. 4 (a), Pe at the walls (x2 ¼ 1 and x2 ¼ 1) is zero owing to the no-slip boundary condition of velocities at the vertical wall. However, because of the influence of the large-scale convection, we can observe an increase in Pe in the bulk region; consequently, a boundary layer where the velocity increases steeply from the noslip wall is also observed. In Fig. 4 (b), the temperature profile shows a drastic change within the thermal boundary layer that can be determined by the gradient of the mean temperature at the wall [1,19]. To study the contribution of convection and conduction to turbulent heat transfer, we report the temporally and horizontally

a@x@h 2

and conductive heat flux profiles, aðh

A;t

h hc Þ=L2

rðNuw Þ

Nuw

Ref. [1] Ref. [8] Pres. Ref. [1] Ref. [8] Pres. Ref. [8] Pres. Ref. [1] Ref. [8] Pres. Ref. [8]

5 106 233:90 225:27 233:29 312:05 325:75 308:28 312:08 325:75 311:84 13:58 2 107 460:49 463:20 451:30 618:84 616:36 605:47 551:66 616:36 551:62 20:46

averaged dimensionless convective heat flux profiles, D E

rðNuX Þ

NuX

hu2 hiA;t

aðhh hc Þ=L2 ,

, in Fig. 5 for different

Ra. Note that both the convective and conductive heat fluxes are symmetric [20] about the horizontal center plane (x2 ¼ 0); we show them only in the interval of 1 6 x2 6 0. We found that the convective heat flux is small in the near-wall region owing to the no-slip boundary condition, and it increases and becomes flat in the bulk region because of the fully developed large-scale motion. The conductive heat flux dominates in the near-wall

region owing to the large temperature gradient and reduces to zero in the bulk region. The maximum value of the convective heat flux increases and the boundary layer thickness decreases as Ra increases. To investigate the scaling of the energy spectra of the turbulent Rayleigh–Bénard convection, we illustrate the time-averaged spectra of velocity and temperature as functions of horizontal wavenumber for Ra ¼ 2 107 in Fig. 6. Throughout this paper, the total kinetic energy spectrum is given by

Ek ðjH Þ ¼

1 X D b0 ! b0 ! E ui k ; t ui k ; t x2 ;t 2! j k j¼jH

and the temperature spectrum is denoted as

Eh ðjH Þ ¼

X D 0 ! 0 ! E hb k ; t hb k ; t x2 ;t ! j k j¼jH

where u0i ¼ u uiA and h0 ¼ h hiA are fluctuations, ub0i and hb0 denote Fourier transforms of u0i and h0 , respectively, and indicates qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ! complex conjugate. Here, k ¼ ðj1 ; j2 ÞT ; jH ¼ j21 þ j22 , and j1 and

j2 are wavenumbers in the x1 - and x2 -directions, respectively. Consistent with Kerr’s work [1], for the kinetic energy spectrum (see Fig. 6 (a)), the roll-off rate is 5=3 in low-wavenumber region, while a high-wavenumber spectrum of j3 H is obtained. On the other hand, for the temperature spectrum (see Fig. 6 (b)), a roll-off rate of 7=5 is obtained for low wavenumbers, whereas for higher wavenumbers, a

Fig. 4. Comparisons of statistically steady-state solutions from MVDPM-STD with DNS data in [1,8] for Ra ¼ 2 107 . ðaÞ Time-averaged normalized root-mean-square D velocities Peðu1;rms Þ and Peðu2;rms Þ and ðbÞ time-averaged root-mean-square temperature hrms and time-averaged temperature hiA;t .

8

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

Fig. 5. Temporally and horizontally averaged dimensionless ðaÞ conductive heat flux profiles and ðbÞ convective heat flux profiles along the vertical direction for various Ra. These profiles are symmetric about x2 ¼ 0.

Fig. 6. Time-averaged horizontal spectra for Ra ¼ 2 107 . ðaÞ Total kinetic energy spectrum and ðbÞ temperature spectrum. 5=3

spectrum of jH holds. This is in good agreement with the results demonstrated in [21]. Fig. 7 shows the instantaneous horizontal planes of vertical velocity and temperature fields close to the bottom plane (x2 ¼ 0:9) and in the mid-layer plane (x2 ¼ 0) and the vertical plane (x3 ¼ 8) of the vertical velocity and temperature fields for Ra ¼ 2 107 . We set the flow and temperature contour levels in the range of 0:3 0:3 to better illustrate the large-scale motion of turbulent Rayleigh–Bénard convection. From Fig. 7, it can be clearly seen that the temperature and vertical velocity appear to vary at different attitudes: large-scale structures dominate the inside of the cell and small-scale structures become prominent in the near-wall region, and the motion of updraft and downdraft arises from the hot and cold planes, respectively. In addition, the hot plumes are taken from the bottom plane to the top plane by the ascending flows, and the cold plumes fall to the bottom plane because of the descending flows. 3.4.2. 3D Rayleigh–Bénard convection in bounded domain To validate the proposed algorithm for solving 3D turbulent flows in the bounded domain, we considered a 3D Rayleigh– Bénard convection problem presented in [22] in bounded domain

X ¼ ½0; 1 ½0; 1 ½0; 1 with Pr ¼ 0:7 for Ra ¼ 106 ; 3 106 , and 107 . The top and bottom walls are kept at constant cold (hc ¼ 0)

and hot (hh ¼ 1) temperatures, respectively. The other four vertical walls are adiabatic. All six walls had the no-slip velocity boundary condition. Similar to [22], we implemented N 1 ¼ N 2 ¼ N 3 ¼ 128; 192, and 256 for Ra ¼ 106 ; 3 106 , and 107 , respectively. The initial static state with p ¼ 0 and h ¼ 0:5ðhh þ hc Þ was implemented, and governing Eqs. (1)–(3) were considered. To evaluate the proposed MVDPM-STD, we define the following variables measured in [22]. The Péclet numbers

PeðuÞ ¼

rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Ra 2 Ra 2 u1 þ u22 þ u23 V;t and Peðu2 Þ ¼ u Pr Pr 2 V;t

are taken based on values of u and vertical velocity u2 , respectively. The following Nusselt numbers are involved: D E D E pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ @h @h NuX ¼ 1 þ PrRahu2 hiV;t ; Nuh;c ¼ 0:5 @x j þ j ; x ¼0 x ¼1 @x 2 2 2 2 V;t

V;t

Nukin ¼ 1 þ Pr hu iV;t and Nuh ¼ hh iV;t ;

where

u ¼ 12

P @ui i;j

@xj

@u

þ @xj i

2

and

h ¼

P @h 2 i

@xi

are the kinetic and

thermal energy dissipation rates, respectively. Table 4 presents comparison of the characteristic values obtained from MVDPMSTD with the results from [22]. It is confirmed that the results obtained from MVDPM-STD show good agreement with [22] for all values of Ra.

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

9

Fig. 7. Instantaneous contour plots in the horizontal planes ða; b; c; dÞ and vertical planes ðe; f Þ of vertical velocity ða; c; eÞ and temperature ðb; d; f Þ for Ra ¼ 2 107 . ða; bÞ : x2 ¼ 0:9; ðc; dÞ : x2 ¼ 0, and ðe; f Þ : x3 ¼ 8.

4. Conclusion In summary, we proposed a more efficient and simpler method, MVDPM-STD, by improving the FDMPM-1P [8] with the application of staggered time discretization to solve time-dependent natural convection problems. In the current algorithm, temperature and pressure fields are evaluated at half-integer time levels (n þ 12), and the velocity fields are evaluated at integer time levels

(n þ 1). This enables us to solve the coupled thermal transport and momentum equations using existing algorithms, and also reduce the computational cost because of no additional operations on the buoyancy part required for the intermediate solutions. Numerical tests on 2D Rayleigh–Bénard convection showed that the MVDPM-STD significantly mitigates the time-step restriction, compared with SIPM, and its stability is slightly better than that of FDMPM-1P for all Ra. However, its computational cost for each

10

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

Table 4 Comparison of characteristic values obtained from MVDPM-STD with values from [22] for different Ra with different meshes. Ra

PeðuÞ

Peðu2 Þ

NuX

Nuh;c

Nukin

Nuh

Pres.

Ref. [22]

Pres.

Ref. [22]

Pres.

Ref. [22]

Pres.

Ref. [22]

Pres.

Ref. [22]

Pres.

106

210:31

208:80

144:84

145:58

8:33

8:33

8:30

8:35

8:24

8:24

8:26

8:26

3 106

358:11

357:11

249:00

249:96

11:38

11:46

11:42

11:48

11:33

11:35

11:35

11:37

107

654:00

654:86

454:68

454:92

16:30

16:22

16:27

16:27

16:16

16:07

16:21

16:10

time step is nearly the same as that of SIPM, resulting in a computational speedup. Further, 2D periodic forced Rayleigh–Bénard convection was simulated to demonstrate the second-order temporal and spatial accuracy of MVDPM-STD. We also performed highresolution DNS of 2D Rayleigh–Bénard convection for Ra up to 1010 under fix-temperature or fix-flux boundary conditions, which confirmed the robustness of MVDPM-STD for types of thermal boundary conditions in a wide range of Ra. Finally, numerical predictions of 3D turbulent Rayleigh–Bénard convection indicated that MVDPM-STD is capable of providing reasonable 3D predictions compared with previous numerical results [1,8,22]. In addition, we theoretically confirmed that the proposed MVDPM-STD is second-order in time and has better stability than FDMPM-1P. Because the staggered time discretization in MVDPM-STD induces a non-iterative monolithic approach for the momentum-energy coupling problem, MVDPM-STD can be extended to problems involving the scalar transport and Navier–Stokes equations, such as multiphase flows [23,24] and chemotaxis-driven bioconvection problems [25]. Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgments This work was supported by National Research Foundation of Korea (NRF) grant funded by the Korean government (Ministry of Science and ICT) (No. NRF-20151009350) and KISTI (K-18-L12C08), and in part by Yonsei University (Yonsei Frontier Lab.–Young Researcher Supporting Program) of 2018. 1

Appendix A. Notations of mbcnþ2 ; cbc

nþ1

, and ebcn

Without loss of generality, in 2D domain X ¼ ½0; L1 ½0; L2 , we consider Dirichlet boundary conditions for the velocity field, and

Ref. [22]

the fluid is heated from the bottom wall (hh ) and cooled at the upper wall (hc ). The two side walls were thermally insulated (see Fig. A.8 (a)). Fig. A.8 (b) shows a discrete domain Xh ¼ x1;i ; x2;j j1 6 i 6 N 1 ; 1 6 j 6 N 2 for X with boundary @ Xh ¼ x1;1 ; x2;j ; x1;N1 ; x2;j ; x1;i ; x2;1 ; x1;i ; x2;N2 j1 6 i 6 N 1 ; 1 6 j 6 N 2 . Here, N 1 and N 2 are the grid points in the x1 - and x2 -directions in X, respectively. For convenience, we introduce the notations ip ¼ i þ 1 for 0 6 i 6 N 1 1; im ¼ i 1 for 1 6 i 6 N 1 ; jp ¼ j þ 1 for 0 6 j 6 N 2 1, and jm ¼ j 1 for 1 6 j 6 N 2 . As shown in Fig. A.8 (b), we define the intervals h1;i ; h2;j ; g 1;i , and g 2;j . Let n1 n1 un1;i;j ; un2;i;j ; hi;j 2 , and pi;j 2 be approximations of u1 x1;i ; x2;j þ 0:5h2;j ; nDtÞ; u2 x1;i þ 0:5h1;i ; x2;j ; nDt ; h x1;i þ 0:5h1;i ; x2;j þ 0:5h2;j ; n 12 Dt , respectively, for and p x1;i þ 0:5h1;i ; x2;j þ 0:5h2;j ; n 12 Dt , 0 6 i 6 N 1 and 0 6 j 6 N 2 in a staggered grid. Here, the discrete boundary conditions at time level n þ 1 are defined as,

unþ1 1;i;0 ¼ u1 x1;i ; 0; ðn þ 1ÞDt ; unþ1 2;i;1 ¼ u2 x1;i þ 0:5h1;i ; 0; ðn þ 1ÞDt ; unþ1 1;i;N 2 ¼ u1 x1;i ; L2 ; ðn þ 1ÞDt ; unþ1 2;i;N 2 ¼ u2 x1;i þ 0:5h1;i ; L2 ; ðn þ 1ÞDt ; unþ1 1;1;j ¼ u1 0; x2;j þ 0:5h2;j ; ðn þ 1ÞDt ; unþ1 2;0;j ¼ u2 0; x2;j ; ðn þ 1ÞDt ; unþ1 1;N 1 ;j ¼ u L1 ; x2;j þ 0:5h2;j ; ðn þ 1ÞDt ; unþ1 2;N 1 ;j ¼ u2 L1 ; x2;j ; ðn þ 1ÞDt ; nþ12 hi;0 ¼ h x1;i þ 0:5h1;i ; 0; n þ 12 Dt ; nþ1 hi;N22 ¼ h x1;i þ 0:5h1;i ; L2 ; n þ 12 Dt ; nþ1

@h0;j 2 @x1 nþ1 @hN 2;j 1

@x1

¼ ¼

@hð0;x2;j þ0:5h2;j ;ðnþ12ÞDtÞ @x1

;

@hðL1 ;x2;j þ0:5h2;j ;ðnþ12ÞDt Þ @x1

:

Fig. A.8. ðaÞ Schematic diagram and boundary conditions for Rayleigh–Bénard convection and ðbÞ discrete domain.

11

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677 nþ1

n

ebci;j ¼ 14 un2;i;jm g 1 hi;0 2 þ 12 2;jm

1 h2;j g 2;jm

nþ1

n

ebci;j ¼ 14 un2;i;jp g 1 hi;N22 þ 12 2;jp

nþ1

n

ebci;j ¼ 14 un1;i;j n

@h0;j 2 @x1

ebci;j ¼ 14 un1;ip;j

nþ1

hi;0 2 for j ¼ 1,

1 h2;j g 2;jp

nþ1

hi;N22 for j ¼ N 2 1,

nþ1

12

nþ1 @hN 2;j 1

@x1

1 h1;i

þ 12

@h0;j 2 @x1 1 h1;i

for i ¼ 1, and

nþ1 @hN 2;j 1

@x1

for i ¼ N 1 1.

Appendix B. Theoretical analysis of the proposed method In this appendix, we discuss the temporal accuracy and stability of the proposed method theoretically. The following hypotheses are required for the theoretical analysis.

Fig. A.9. Location for vertical velocity u2 . nþ12

We present the explicit form of mbc

T nþ1 nþ1 :¼ x1 bc 2 ; x2 bc 2 in

nþ1 nþ1 the momentum equation, where x1 bci;j 2 and x2 bci;j 2 at the same locations with un1;i;j and un2;i;j , respectively.

are defined nþ1

and linearized buoyancy terms at the ðn þ 1Þth time level are as follows:

rﬃﬃﬃﬃﬃﬃ 1 unþ1 1 n @unþ1 1 nþ1 @un2 @pnþ2 1 Pr @ 2 unþ1 2 2 2 u2 þ þ u2 þ 2 2 2 2 Ra @x2 Dt C @x2 @x2 C @x2 C C 1 n unþ1 1 nþ12 1 @u 2;i;j nþ1 2 h ¼ un unþ1 unþ1 þ 2;i;jm þ u2;i;j 2 4g 2;j 2;i;j 2;i;jp 2 Dt @x2 C C ! rﬃﬃﬃﬃﬃﬃ nþ1 nþ1 nþ1 nþ1 nþ1 nþ1 pi;j 2 pi;jm2 1 Pr u2;i;jp u2;i;j u2;i;j u2;i;jm þ 2g 2;j Ra g 2;j h2;j h2;jm nþ1

2 2 1 h2;jm hi;j þ h2;j hi;jm 2 2g 2;j

ðA:1Þ

Because the velocities at the boundary are known, we only need to nþ12

consider x2 bci;j for 1 6 i 6 N1 ; 2 6 j 6 N 2 , which corresponds to the interior points. Moreover, because we used the second-order central difference scheme in space, the boundary conditions do not directly affect the unknown variables unþ1 2;i;j for 3 6 j 6 N 2 1. Thus, nþ12

x2 bci;j

vanishes for 3 6 j 6 N2 1.

Therefore, from Eq. (A.1), we have qﬃﬃﬃﬃ unþ1 nþ1 1 Pr 2;i;1 x2 bci;j 2 ¼ 4g1 un2;i;j unþ1 for j ¼ 2 and 2;i;1 þ 2g 2;j Ra h2;jm 2;j qﬃﬃﬃﬃ unþ1 1 nþ 1 Pr 2;i;N2 for j ¼ N 2 1. x2 bci;j 2 ¼ 4g12;j un2;i;j unþ1 2;i;N2 þ 2g 2;j Ra h2;j Similarly, based on the corresponding spatial discretizations, we obtain qﬃﬃﬃﬃ unþ1 nþ1 Pr 1;1;j

x1 bci;j 2 ¼ 4g1 un1;i;j unþ1 þ 2g1 for i ¼ 2 and 1;1;j Ra h1;im 1;i 1;i qﬃﬃﬃﬃ unþ1 nþ12 1;N nþ1 Pr 1 ;j for i ¼ N 1 1. x1 bci;j ¼ 4g1 un1;i;j u1;N1 ;j þ 2g1 Ra h 1;i

nþ1 cbci;j

¼

nþ1 u2;i;j

h2;j

h2;j

/2 þ

unþ1 1;i;j hxi

/3

1;i

unþ1 1;ip;j hxi

/4 ,where

1

Hypothesis 3. We assume that C is a Lipschitz operator with conTherefore, given uniform bounded stant CC .

kuh1 k; kuh2 k; kv h1 k; kv h2 k; khh1 k, and khh2 k, there exists a constant C C > 0 such that

h h C u ; v C uh ; v h 6 C C uh þ uh v h v h 1 1 2 2 1 2 1 2 þ C C v h1 þ v h2 uh1 uh2 and

C uh1 ; hh1 C uh2 ; hh2 6 C C uh1 þ uh2 hh1 hh2 þ C C hh1 þ hh2 uh1 uh2 ; where uh1 ; uh2 ; v h1 ; v h2 2 H1 ðXh Þd and hh1 ; hh2 2 H1 ðXh Þ. Hypothesis 4. We assume that L :¼ DG is a Lipschitz operator with constant C L . In other words, given uniform bounded kuh1 k; kuh2 k; khh1 k, and khh2 k, there exists a constant C L > 0 such that

h Lu Luh 6 C L uh uh 1 2 1 2 6 C L hh1 hh2 ;

and

h Lh1 Lhh2

where uh1 ; uh2 2 H1 ðXh Þd and hh1 ; hh2 2 H1 ðXh Þ. B.1. Temporal accuracy To derive the temporal accuracy of MVDPM-STD, we first estimate the global error of MVCPM-STD solutions approximated to the numerical solutions for Eqs. (4)–(6) (i.e., v n ; nn2 , and qn2 at time nDt). We assume that all the solutions satisfy the uniform boundedness in Hypothesis 2. The differences between the MVCPM-STD solutions and the solutions for Eqs. (4)–(6) are defined as 1

1

1

1

1

n ¼ un v n ; and gn2 ¼ pn2 qn2 : 1

1

1

1; j ¼ 1

wn2 ¼ hn2 nn2 ;

0; otherwise 1; j ¼ N2 1

0; otherwise 1; i ¼ 1

According to the solution procedure of the MVCPM-STD shown in Eqs. (16)–(19), MVCPM-STD can be expressed as

0; otherwise 1; i ¼ N1 1

/1 ¼ /2 ¼ /3 ¼ /4 ¼

/1

1;i

nþ1 u2;i;jp

1

assume that kun k < C; khn2 k < C, and kpn2 k < C for nDt < T f . Here, C is a given positive constant.

Here, based

on the cell shown in Fig. A.9, the detailed expression of x2 bc 2 relating to unsteady, linearized convection, pressure, diffusion,

nþ1

Hypothesis 2. If we consider the process up to time T f , we can

0; otherwise

:

1 1 1 1 1 hnþ2 hn2 1 ¼ C un ; hnþ2 C un ; hn2 2 2 Dt 1 1 1 n þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Lhnþ2 þ Lhn2 þ ebc 2 RaPr

ðB:1Þ

12

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

1 unþ1 un 1 ¼ C un ; u;nþ1 C u;nþ1 ; un 2 rﬃﬃﬃﬃﬃﬃ 2 Dt 1 Pr ;nþ1 1 þ þ Lun Gpnþ2 Lu 2 Ra 1

nþ12

þ F hnþ2 þ mbc

nþ1 where I 2 ¼ DtkPk12 C v nþ1 ; v nþ1 12 Cðv n ; v n Þ þ 12 C v n ; v nþ1 nþ1 nþ1 12 C v nþ1 ; v n þ mbc 2 mbcCN2 k. Taking the summation of inequalities (B.7) and (B.10) enables us to obtain

ðB:2Þ

1 þ C 1 Dt nþ12 n12 w þ nþ1 6 w þ kn k 1 C 1 Dt

along with discrete continuity Eq. (9). Furthermore, a fully discrete form of the governing equations based on the Crank–Nicolson scheme is 1 1 1 1 1 nnþ2 nn2 1 ¼ C v n ; nnþ2 C v n ; nn2 2 2 Dt 1 1 n nþ12 þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Ln þ Lnn2 þ ebcCN ; 2 RaPr

Dv nþ1 ¼ cbcCN ; nþ1

v

v Dt

nþ1

n

ðB:3Þ ðB:4Þ

Subtracting Eq. (B.3) from Eq. (B.1) gives

1 1 1 1 Dt n nþ12 C v ;n wnþ2 ¼ wn2 þ þ C v n ;nn2 C un ; hnþ2 2 n 1 1 Dt n n n12 C u ;h þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Lwnþ2 þ Lwn2 þ Dt ebc ebcCN ; 2 RaPr ðB:6Þ n n where the last term Dt ebc ebcCN vanishes owing to Eq. (12). Applying the triangle inequality, the mean value theorem, and Hypotheses 3 and 4 to Eq. (B.6), we infer that

ðB:7Þ

Dt n nþ1 C v ;v þ C v nþ1 ; v n C un ; u;nþ1 2 rﬃﬃﬃﬃﬃﬃ ;nþ1 n Dt Pr ;nþ1 Lu ;u þ þ Lun Lv nþ1 Lv n C u 2 Ra 1 1 1 nþ1 nþ1 DtGgnþ2 þ Dt F hnþ2 F nnþ2 þ Dt mbc 2 mbcCN2 ;

nþ1 ¼ n þ

ðB:8Þ and subtracting Eq. (B.4) from Eq. (9) and applying Eq. (13) gives

ðB:9Þ

Referring to the proof of Theorem 2 in [12], we derive

rﬃﬃﬃﬃﬃﬃ ! nþ1 6 kn k þ Dt C C ku;nþ1 k þ kv nþ1 k þ 1 Pr C L kn k 2 Ra ! rﬃﬃﬃﬃﬃﬃ 1 1 Pr þ Dt C C ðkun k þ kv n kÞ þ C L nþ1 þ DtC F wnþ2 2 Ra rﬃﬃﬃﬃﬃﬃ ! nþ1 1 Pr n n C L u;nþ1 unþ1 þ I 2 ; þ Dt C C ðku k þ kv kÞ þ 2 Ra

ðB:11Þ

Lemma 1. Given h2 ¼ n2 and u0 ¼ v 0 , from property (B.11), we can conclude that

n12 w þ kn k 6

1

1

j2 n j C 1 Dt X u u;j þ I r nj 0 1 C 1 Dt j¼0 C 1 Dt

!

1

for all nDt < T f , where r0 ¼ ð1 þ C 1 Dt Þ=ð1 C 1 DtÞ and I 2 ¼ 0. Proof. For the detailed proof of the present lemma, readers can refer to that of Lemma 1 in [8]. Subsequently, we obtain the conclusion without further proof. h nþ1 I 2 From Lemma 1, we found that unþ1 u;nþ1 and C1 Dt should be 1 further estimated to obtain a complete estimate of wnþ2 þ nþ1 . We now estimate unþ1 u;nþ1 by applying Eq. (19):

nþ1 1 1 u u;nþ1 ¼ dunþ1 du;nþ1 ¼ Dt Gdpnþ2 6 Dt kGkdpnþ2 :

ðB:12Þ In addition, multiplying D by Eq. (B.2) along with Eq. (9) gives

Subtracting Eq. (B.5) from Eq. (B.2) produces

Dnþ1 ¼ 0:

1 ! 2 ;nþ1 Inþ C 1 Dt u ; unþ1 þ 1 C 1 Dt C 1 Dt

where C 1 is a positive constant that does not depend on either Dt or h. Here, we assume that C 1 Dt < 12. Consequently, we can provide the following lemma. 1

1 1 ¼ C v nþ1 ; v nþ1 Cðv n ; v n Þ 2 rﬃﬃﬃﬃﬃﬃ 2 1 1 Pr nþ1 1 nþ1 Lv þ Lv n Gqnþ2 þ F nnþ2 þ mbcCN2 : þ 2 Ra ðB:5Þ

1 CL nþ12 n12 Dt C C ðkun k þ kv n kÞ þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ wnþ2 w 6 w þ 2 RaPr 1 Dt C L C C ðkun k þ kv n kÞ þ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ wn2 þ 2 RaPr Dt nþ12 nþ12 n12 n12 þ C C n þ h þ n þ h kn k: 2

þ

nþ1 n n1 cbc 2cbc þ cbc 1 1 dpnþ2 ¼ ðDG Þ1 þ ðDG Þ1 D C un1 ; u;n Dt 2 þC u;n ; un1 C un ; u;nþ1 C u;nþ1 ; un rﬃﬃﬃﬃﬃﬃ 1 Pr ðDG Þ1 D Lu;nþ1 þ Lun Lu;n Lun1 þ 2 Ra 1 nþ1 n1 þðDG Þ1 D mbc 2 mbc 2 þ ðDGÞ1 D F hnþ2 1 ðB:13Þ F hn2 ¼ I 1 þ I2 þ I3 þ I 4 þ I5 : nþ1

I1 þ I2 þ I3 þ I4 6 C I1 þ C I2 þ C I3 þ C I4 Dt þ C I2 kun u;n k;

ðB:14Þ

where C I1 ; C I2 ; C I3 , and C I4 are positive constants. By the uniform boundedness of kðDGÞ1 k and kDk, the mean value theorem, and Hypothesis 1, we estimate I5 as

1 1 1 I5 6 kðDG Þ kkDkF hnþ2 F hn2 1 1 6 kðDG Þ1 kkDkC F hnþ2 hn2 6 C I5 Dt; 1

ðB:10Þ

nþ1

Under Hypothesis 3 and the assumption that mbc 2 and cbc have bounded discrete temporal second-order derivatives, referring to the proof of Theorem 4 in [12] for the estimate of I1 þ I2 þ I3 þ I4 gives us

ðB:15Þ 1

where C I5 is a positive constant and khnþ2 hn2 k ¼ OðDt Þ is applied. Given the sum of all the estimates along with inequality (B.12), there exists a positive constant C 2 such that

nþ1 u u;nþ1 6 C 2 DtðDt þ kun u;n kÞ:

ðB:16Þ

13

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

Hereafter, we use C 2 > 0 as a generic constant that may differ from line to line or even within the same line. Subsequently, analogous to Lemma 3 in [12], we obtain the following lemma from inequality (B.16).

We require the following lemma stated without proof, and we refer readers to Corollary 8 in [12] for further details. Lemma 3. Given a discrete domain, Xh , let uh and v h 2 H1 ðXh Þd , and hh 2 H1 ðXh Þ. If Duh ¼ 0 in Xh , and uh njCh ¼ 0, there exist constants

Lemma 2. Given u0 ¼ u;0 for all ðn þ 1ÞDt < T f , the following estimate holds

C uC ; C hC 2 R such that the discrete operator C satisfies

nþ1 u u;nþ1 6 C 2 Dt 2 ð1 þ C 2 Dt Þn :

nþ1

Considering the boundedness of

I 2 , C 1 Dt

for all ðn þ 1ÞDt < T f , we

derive nþ12

I 6 C I Dt 2 C 1 Dt

ðB:17Þ

nþ1 with a positive constant C I , where we apply I 2 ¼ O Dt 3 on account of the temporal second-order accuracy in Eq. (11). 1 Subsequently, the estimate of wnþ2 þ nþ1 can be finalized. nþ1

Theorem 1. Assuming that Hypotheses 1–4 hold, and mbc 2 and n ebc have bounded discrete second-order derivatives in time, starting from the same initial conditions, the MVCPM-STD solutions converge to the solutions for Eqs. (4)–(6) with second-order temporal accuracy:

nþ12 w ¼ O Dt2 ; nþ1 ¼ O Dt2 ;

1 and gnþ2 ¼ O Dt 2

vh ; C

uh ; v h

Proof. For a detailed proof of the present lemma, readers can refer to the proof of Theorem 1 in [8]. Therefore, we present the conclusion without further proof. h The above-mentioned discussion and arguments show that the MVCPM-STD solutions converge to the solutions for Eqs. (4)–(6) with second-order temporal accuracy. Thus, we obtain Theorem 2 as follows without further proof.

Remark 1. Theorems 1 and 2, along with the triangle inequality, ensure that the MVDPM-STD solutions are second-order accurate with respect to the solutions for Eqs. (4)–(6) in time. Hence, owing to the fact that the Crank–Nicolson discretization (4)–(6) gives a second-order temporal approximation to the continuous Eqs. (1)–(3), we conclude that both MVCPM-STD and MVDPM-STD solutions are O Dt2 of the exact solutions in terms 2

of a discrete l -norm. B.2. Stability analysis In the current appendix, we analyze the stability of the proposed MVDPM-STD in terms of energy evolution. Periodic or homogeneous boundary conditions are assumed for velocity and nþ1

n

n

nþ1

nþ12

mbcCN2 ; mbc

D

2 E 2 hh ; C uh ; hh ¼ C hC h hh ;

First, we consider the energy evolution of the linear system given in Eqs. (16)–(19), which ensure decoupling among temperature, velocity, and pressure. Theorem 3. Under Hypotheses 1–4, for the linear system (16)–(19), there exist constants s; B 2 Rþ such that

D E 1 1 1 umþ1 þ um umþ1 þ um 1 ; þ RaC 2F hmþ2 ; hmþ2 2 2 2 2 rﬃﬃﬃﬃﬃﬃ m Dt X Pr G unþ1 þ 2un þ un1 2 ðd 1Þ þ 16d n¼1 Ra ! rﬃﬃﬃﬃﬃﬃ 2 Ra nþ12 n12 G h þ h þð4d 1ÞC 2F 6B Pr

nþ1

; cbcCN ;

cbc ; ebcCN , and ebc vanish [26]. Divergence-free velocity fields are assumed to be orthogonal to the gradient field, and L and G are commutative [27].

and

Proof. We define kinetic energy as nþ12

Ek

¼

1 unþ1 þ un unþ1 þ un ; 2 2 2

at t ¼ n þ 12 Dt. The difference in kinetic energy for time levels n 12 and n þ 12 is denoted by nþ12

Ek

n12

Ek

1 unþ1 þ un un þ un1 unþ1 þ un un þ un1 þ ; 2 2 2 2 2 D E Dt nþ1 1 1 u þ 2un þ un1 ; Gpn2 þ Gpnþ2 ¼ 8 rﬃﬃﬃﬃﬃﬃ Dt Pr nþ1 u þ 2un þ un1 ; Lu;nþ1 þ Lun þ 16 Ra Dt nþ1 u þ 2un þ un1 ; C un ; u;nþ1 þLu;n þ Lun1 16 þC u;nþ1 ; un þ C un1 ; u;n þ C u;n ; un1 1 1 E Dt D nþ1 u þ 2un þ un1 ; F hnþ2 þ F hn2 þ 8

¼

:¼ J1 þ J 2 þ J 3 þ J 4 : Owing to the orthogonal property of the pressure gradient and velocity, it is clear that J1 vanishes. From Eq. (19), we can rewrite J2 as

rﬃﬃﬃﬃﬃﬃ Pr nþ1 u þ 2un þ un1 ; Lunþ1 þ Lun þ Lun þ Lun1 Ra rﬃﬃﬃﬃﬃﬃ E Dt Pr D nþ1 1 1 þ u þ 2un þ un1 ; L Gdpnþ2 þ L Gdpn2 16 Ra :¼ J 21 þ J 22 ;

J2 ¼

Subsequently,

and

where n is the outward unit normal vector on Ch .

Theorem 2. Under the same hypotheses as those of Theorem 1, because the application of the second-order approximate factorization technique, the MVDPM-STD solutions converge to the MVCPM-STD solutions with second-order temporal accuracy under the same initial conditions.

fields.

2 2 ¼ C uC h v h

when Dt 6 s for all ðm þ 1ÞDt < T f . Here, ðd 1Þ 2 Rþ ð4d 1Þ 2 Rþ for d ¼ 2 or 3.

for all ðn þ 1ÞDt < T f .

temperature

Dt 16

where J 22 vanishes owing to the commutative property of discrete operators L and G and the orthogonal property of the pressure gradient and velocity. By applying Hypothesis 4, we can estimate J 21 as

14

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

J 21

rﬃﬃﬃﬃﬃﬃ Dt Pr nþ1 G u ¼ þ 2un þ un1 ; G unþ1 þ 2un þ un1 16 Ra rﬃﬃﬃﬃﬃﬃ Dt Pr G unþ1 þ 2un þ un1 2 ; ¼ 16 Ra

From Hypothesis 3 and Lemma 3, J 5 can be written as

E 1 1 1 Dt D nþ12 h þ hn2 ; C un ; hnþ2 þ hn2 4 n h;nþ1 h;n1 o 1 1 2 Dt 2 6 max C C 2 ; C C 2 h hnþ2 þ hn2 4 2 nþ1 h Dt n1 1 2 1 2 2 6 C 5 hnþ2 þ hn2 ¼ C 5 h Dt Eh 2 þ Eh 2 ; 2

J5 ¼

where the periodic or homogeneous boundary conditions are considered for the discrete velocities. Using Eq. (19), we reformulate J3 as D n n1 unþ1 þun un þun1 u;nþ1 þu;n E unþ1 þ un þ u þu 2 ; ;C 2 2 2 D unþ1 þun un þun1 u;nþ1 þu;n un þun1 E Mt n n1 ; ;C 4 u þu þ 2 2 2 2 D n ;nþ1 ;nþ1 n n1 ;n E Mt unþ1 þ2un þun1 ; C u ; u ; u ;u þ C u þ C u þ C u;n ;un1 16 2 D E Mt unþ1 þ2un þun1 ; Cðun ; u;n Þ þ Cðu;n ; un Þ þ C un1 ; u;nþ1 þ C u;nþ1 ; un1 þ 16 2 D E D E n n1 unþ1 þun n n1 nþ1 n n1 n unþ1 þun ¼ Mt ; C u þu ; 2 ;C u 2þu ; u 2þu þ u þu 2 2 2 2 D E nþ1 n n n1 u;nþ1 þu;n ;nþ1 ;n n n1 un þun1 u 2þu ; C u þu ; C u 2þu ; u þu Mt 4 2 2 2 2 1 1 nþ1 n1 2 n n1 Gpnþ2 þGpn2 n n1 2 2 un1 þun unþ1 þun þ u þu Mt2 ;C u þu ; ; C Gp þGp ; 2 2 2 2 2 2 D E Mt unþ1 þ2un þun1 16 ; C un un1 ;u;nþ1 u;n þ C u;nþ1 u;n ; un un1 2

J 3 ¼ Mt 4

:¼ J 31 þ J 32 þ J 33 þ J 34 :

For the estimate of J31 þ J32 þ J33 , we refer to the proof of Theorem 9 in [12]:

n1 n1 nþ1 nþ1 2 J 31 þ J 32 þ J 33 6 C 31 h Dt Ek 2 þ Ek 2 þ C 32 Dt3 þ C 33 Dt3 3Ek 2 þ Ek 2 ;

h;nþ12

where C C

1

1

* !+ 1 1 Dt hnþ2 þ hn2 unþ1 þ 2un þ un1 ; 2F 8 2 1 1 Dt nþ1 6 ku þ 2un þ un1 kC F hnþ2 þ hn2 8 ! rﬃﬃﬃﬃﬃﬃ rﬃﬃﬃﬃﬃﬃ 2 Dt 1 Pr nþ1 1 2 Ra nþ12 n12 n n1 2 6 ku þ 2u þ u k þ C F h þ h 8 2 Ra 2 Pr rﬃﬃﬃﬃﬃﬃ Dt Pr nþ1 kG u þ 2un þ un1 k2 6 16d Ra rﬃﬃﬃﬃﬃﬃ 1 1 2 Dt 2 Ra CF þ G hnþ2 þ hn2 : 16d Pr

J4 ¼

Then, we consider the thermal energy at time level n þ 12 on the basis of the following definition: nþ1 Eh 2

1 D nþ12 nþ12 E : h ;h ¼ 2

The energy difference between time levels n 12 and n þ 12 is denoted by nþ12

Eh

n12

Eh

E 1 1 1 1 D nþ12 h þ hn2 ; hnþ2 hn2 2 E 1 1 1 Dt D nþ12 h þ hn2 ; C un ; hnþ2 þ C un ; hn2 ¼ 4 E 1 1 1 Dt 1 D nþ12 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ h þ hn2 ; Lhnþ2 þ Lhn2 :¼ J 5 þ J 6 : þ 4 RaPr

¼

h;n12

max jC C

06n
j. Fur-

ther, J6 is reformulated as

1 E 1 1 Dt 1 D nþ12 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ G h þ hn2 ; G hnþ2 þ hn2 4 RaPr 1 1 2 Dt 6 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ G hnþ2 þ hn2 : 4 RaPr

J6 ¼

Defining the combination of two energies nþ12

E

nþ12

¼ Ek

nþ12

þ RaC 2F Eh

at time level n þ 12 gives 1

rﬃﬃﬃﬃﬃﬃ Dt Pr G unþ1 þ 2un þ un1 2 ðd 1Þ 16d Ra rﬃﬃﬃﬃﬃﬃ 1 1 2 Dt Pr ð4d 1ÞC 2F G hnþ2 þ hn2 16d Ra 1 1 2 þ C T1 h Dt En2 þ Enþ2 þ C T2 Dt 3 1 1 þ C T3 Dt 3 3En2 þ Enþ2 ;

1

Enþ2 En2 6

J 34 6 C 34 Dt3 ;

DtGdpnþ2 un DtGdpn2 k ¼ OðDt Þ and kun un1 k ¼ OðDt Þ for estimating J34 . From Hypothesis 1, the mean value theorem, Young’s inequality, and the discrete Poincaré inequality, J4 can be estimated as

are the corresponding constants at time

levels n þ 12 and n 12, respectively, and C 51 ¼

and J34 is estimated to be

where C 31 ; C 32 ; C 33 , and C 34 denote generic positive constants that may differ for each time level; however, they are independent of both Dt and h. Here, we apply ku;nþ1 u;n k ¼ kunþ1 þ

h;n12

and C C

Dt ðd 16d T C 1 ; C T2 ,

where

1Þ and

Dt ð4d 16d

ðB:18Þ

1Þ are positive for d ¼ 2 and d ¼ 3.

C T3

and are generic positive constants. Here, Summing up inequality (B.18) for n from 1 to m, we obtain

rﬃﬃﬃﬃﬃﬃ m Dt X Pr G unþ1 þ 2un þ un1 2 E þ ð d 1Þ 16d n¼1 Ra ! rﬃﬃﬃﬃﬃﬃ 2 Ra nþ12 2 n12 þð4d 1ÞC F G h þ h Pr m m 1 X X 1 2 C T1 h þ 3C T3 Dt2 En2 6 E2 þ Dt 3 C T2 þ Dt mþ12

n¼1

n¼1

m þ1 1 X 1 2 þ C T1 h þ C T3 Dt2 Enþ2 6 b1 þ Dt Bn1 En2 ; 2

n¼2

2

1 P 2 T where b1 ¼ 1 þ C T1 h Dt þ 3C T3 Dt 3 E2 þ Dt 3 m n¼1 C 2 2 2 Bn1 ¼ 2 C T1 h þ 3C T3 Dt 2 . Assuming that Dt 6 s such 2

and that

Bn1 Dt < 1 for all 2 6 n 6 m þ 1, and applying the discrete 2

Gronwall’s inequality [28], we obtain

rﬃﬃﬃﬃﬃﬃ m Dt X Pr G unþ1 þ 2un þ un1 2 ð d 1Þ 16d n¼1 Ra ! rﬃﬃﬃﬃﬃﬃ 2 Ra 2 nþ12 n12 þð4d 1ÞC F G h þ h 6B Pr 1

Emþ2 þ

with B > 0. Thus, the proof of Theorem 3 is complete. On the basis of the energy evolution in Theorem 3, we can obtain the following theorem without proof. Theorem 4. Under Hypotheses 1–4, for MVDPM-STD, there exist constants sstd ; Bstd 2 Rþ such that

X. Pan et al. / International Journal of Heat and Mass Transfer 144 (2019) 118677

mþ1 D mþ1 mþ1 E m m 1 umþ1 1 std þ ustd ustd þ ustd ; þ RaC 2F hstd 2 ; hstd 2 2 2 2 2 rﬃﬃﬃﬃﬃﬃ m X Dt Pr n n1 2 ð d 1Þ G unþ1 þ std þ 2ustd þ ustd 16d n¼1 Ra ! rﬃﬃﬃﬃﬃﬃ 2 Ra nþ12 n12 2 þð4d 1ÞC F G hstd þ hstd 6 Bstd Pr

ðB:19Þ

Dt when Dt 6 sstd for all ðm þ 1ÞDt < T f . Here, 8d ð2d 1Þ 2 Rþ for d ¼ 2or3. Throughout this paper, the variables and operators accompanied by the subscript ’VD’ are associated with MVDPM-STD.

Theorem 4 can be proved similarly to Theorem 4 in [8]. Remark 2. Based on the assumption that Bstd Dt < 1, MVDPM-STD is stable in the sense of inequality (B.19). This assumption is reasonable because small values of h and Dt are generally required to perform accurate numerical simulations. Referring to the stable property of FDMPM-1P in Theorem 4 of [8], MVDPM-STD is more stable than FDMPM-1P for small Ra, and the stability of MVDPMSTD is comparable to that of FDMPM-1P as Ra increases. References [1] R.M. Kerr, Rayleigh number scaling in numerical convection, J. Fluid Mech. 310 (1996) 139–179. [2] H. Johnston, C.R. Doering, Comparison of turbulent thermal convection between conditions of constant temperature and constant flux, Phys. Rev. Lett. 102 (6) (2009) 064501. [3] J. Loewe, G. Lube, A projection-based variational multiscale method for largeeddy simulation with application to non-isothermal free convection problems, Math. Models Methods Appl. Sci. 22 (02) (2012) 1150011. [4] J. Deteix, A. Jendoubi, D. Yakoubi, A coupled prediction scheme for solving the Navier-Stokes and convection-diffusion equations, SIAM J. Numer. Anal. 52 (5) (2014) 2415–2439. [5] Y. Zhang, Y. Hou, J. Zhao, Error analysis of a fully discrete finite element variational multiscale method for the natural convection problem, Comput. Math. Appl. 68 (4) (2014) 543–567. [6] Y. Qian, T. Zhang, The second order projection method in time for the timedependent natural convection problem, Appl. Math. 61 (3) (2016) 299–315. [7] X. Pan, K. Kim, C. Lee, J.-I. Choi, A decoupled monolithic projection method for natural convection problems, J. Comput. Phys. 314 (2016) 160–166. [8] X. Pan, K. Kim, C. Lee, J.-I. Choi, Fully decoupled monolithic projection method for natural convection problems, J. Comput. Phys. 334 (2017) 582–606. [9] X.-J. Huang, Y.-R. Li, L. Zhang, Y.-P. Hu, Turbulent Rayleigh-Bénard convection in a cubical container filled with cold water near its maximum density, Int. J. Heat Mass Transf. 127 (2018) 21–31.

15

[10] H. Karatas, T. Derbentli, Natural convection in differentially heated rectangular cavities with time periodic boundary condition on one side, Int. J. Heat Mass Transf. 129 (2019) 224–237. [11] K. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Trans. Antennas Propag. 14 (3) (1966) 302–307. [12] X. Pan, C. Lee, K. Kim, J.-I. Choi, Analysis of velocity-components decoupled projection method for the incompressible Navier-Stokes equations, Comput. Math. Appl. 71 (8) (2016) 1722–1743. [13] O. Zerkak, T. Kozlowski, I. Gajev, Review of multi-physics temporal coupling methods for analysis of nuclear reactors, Ann. Nucl. Energy 84 (2015) 225– 233. [14] R.M. Beam, R. Warming, An implicit factored scheme for the compressible Navier-Stokes equations, AIAA J. 16 (4) (1978) 393–402. [15] K. Kim, S.-J. Baek, H.J. Sung, An implicit velocity decoupling procedure for the incompressible Navier-Stokes equations, Int. J. Numer. Meth. Fluids 38 (2) (2002) 125–138. [16] X. Pan, C. Lee, J.-I. Choi, Efficient monolithic projection method for timedependent conjugate heat transfer problems, J. Comput. Phys. 369 (2018) 191– 208. [17] J. Kim, P. Moin, Application of a fractional-step method to incompressible Navier-Stokes equations, J. Comput. Phys. 59 (2) (1985) 308–323. [18] G. Yang, J. Wu, Effects of natural convection, wall thermal conduction, and thermal radiation on heat transfer uniformity at a heated plate located at the bottom of a three-dimensional rectangular enclosure, Numer. Heat Tr. A-Appl. 69 (6) (2016) 589–606. [19] G. Oh, K.M. Noh, H. Park, J.-I. Choi, Extended synthetic eddy method to generate inflow data for turbulent thermal boundary layer, Int. J. Heat Mass Transf. 134 (2019) 1261–1267. [20] I. Otic´, G. Grötzbach, M. Wörner, Analysis and modelling of the temperature variance equation in turbulent natural convection for low-Prandtl-number fluids, J. Fluid Mech. 525 (2005) 237–261. [21] J. Niemela, L. Skrbek, K. Sreenivasan, R. Donnelly, Turbulent convection at very high Rayleigh numbers, Nature 404 (6780) (2000) 837. [22] A. Xu, L. Shi, H.-D. Xi, Lattice Boltzmann simulations of three-dimensional thermal convective flows at high Rayleigh number, Int. J. Heat Mass Transf. 140 (2019) 359–370. [23] M. Sussman, E.G. Puckett, A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows, J. Comput. Phys. 162 (2) (2000) 301–337. [24] J.A. Sethian, P. Smereka, Level set methods for fluid interfaces, Annu. Rev. Fluid Mech. 35 (1) (2003) 341–372. [25] A. Chertock, K. Fellner, A. Kurganov, A. Lorz, P. Markowich, Sinking, merging and stationary plumes in a coupled chemotaxis-fluid model: a high-resolution numerical approach, J. Fluid Mech. 694 (2012) 155–190. [26] J.B. Perot, An analysis of the fractional step method, J. Comput. Phys. 108 (1) (1993) 51–58. [27] E.P. Newren, A.L. Fogelson, R.D. Guy, R.M. Kirby, Unconditionally stable discretizations of the immersed boundary equations, J. Comput. Phys. 222 (2) (2007) 702–719. [28] J. Shen, On error estimates of the projection methods for the Navier-Stokes equations: second-order schemes, Math. Comp. 65 (215) (1996) 1039–1065.

Efficient monolithic projection method with staggered time discretization for natural convection problems

Efficient monolithic projection method with staggered time discretization for natural convection problems

Recommend Documents