Steady-State Convection in a Saturated Porous Medium Heated from Below

Note

To view this project in FLAC3D, use the menu command Help ‣ Examples…. The project’s main data files are shown at the end of this example.

The problem of steady-state convection in a porous layer heated from below (the Horton-Rogers-Lapwood Convection Problem) is a classical convection application. We consider a saturated layer of finite thickness and lateral extent. The temperature is fixed at the top and the base of the layer. The boundaries are impermeable to fluid flow, and the lateral boundaries are adiabatic. The problem solution depends on the value of the Rayleigh number, \(Ra\), which is defined for this problem as

(1)\[Ra = {k \over k^T} (\rho_0 c_w) \cdot (\rho_0 g H) \beta_f (T_1 - T_0)\]

where \(H\) is the height of the layer, and \(T_1\), \(T_0\) are the temperatures at bottom and top (see, for example, Lapwood 1948). For Rayleigh numbers below the minimum critical value of \(4 \pi^2\) (or about 39.48), the problem has a trivial steady-state solution, which corresponds to zero specific discharge (the conduction solution). For Rayleigh numbers above \(4 \pi^2\), non-trivial solutions may also exist in the form of convection cells. A spectrum of different steady convection modes can be realized physically; we consider various geometrical aspect ratios and Rayleigh numbers, and we seek some of these solutions in this section.

Saturated Porous Cube

The first example concerns two-dimensional convection in a cube. With the properties and dimensions adopted in the model, the Rayleigh number for this problem is approximately 42. This quantity is larger than the critical value of \(4 \pi^2\). Therefore, non-trivial solutions can be expected. We will illustrate one of these solutions: a two-dimensional roll, which we compare to an analytic solution derived from a linear stability analysis. First, we start with an illustration of the computational methodology. The data file for this example is “HRL-Convection-Cube-LowRa.dat”.

We take, as initial conditions, the conduction solution (see Equation (1) and (3) in natural convection). In FLAC3D, the transient equations are solved for the fluid, and a value for the bulk modulus of the fluid must be provided. In the present algorithm, this value affects the magnitude of stable timestep adopted for the simulation. We select a value of 2 × 10⁵ Pa, which means that, for the property values and zone size adopted in the explicit simulation, the fluid timestep is one order of magnitude smaller than the thermal timestep.

For an explicit mode of calculation, we fix the thermal timestep at 5.5 × 10¹⁰ seconds and the fluid timestep at 5.5 × 10⁹ seconds. We then design a “supersolve” command with the FISH language, which takes one thermal step and 10 fluid-flow steps for each super-step in the calculation. In this way, the thermal and fluid times are synchronized after execution of each super-step.

The explicit mode of calculation is very inefficient for this problem. Alternatively, to reduce the computation time, we use the implicit mode of calculation with a timestep ten times larger for both the thermal and fluid simulations.

The settings for the implicit calculation are

zone fluid implicit on
zone thermal implicit on
model fluid   timestep fix 5.5e10
model thermal timestep fix 5.5e11

The implicit mode commands are listed below.

The FLAC3D grid is shown in Figure 1, together with the location of monitoring points. The initial pore pressure contours and temperature contours corresponding to the conduction solution for a fluid with temperature-dependent density are plotted in Figure 2 and Figure 3.

Figure 1: FLAC3D grid and location of monitoring points.

With this type of problem, there is more than one possible solution (one conduction and one or more advection solutions, depending on the Rayleigh number). The conduction solution is usually obtained, and a small perturbation is often needed to make the system bifurcate to one of the possible convection solutions (e.g., a random distribution of a property, or tilting the gravity slightly). In this example, we were able to obtain a two-dimensional roll by tilting the angle of gravity by 0.06° in the \(x\)-direction. The perturbation is applied for a finite number of steps and then removed. Note that if the same example is run for a long enough time without any perturbation, there seems to be enough numerical noise to cause the solution to bifurcate, but this time to a three-dimensional convection cell.

Figure 2: Initial temperature contours—conduction solution.

Figure 3: Initial pore pressure contours—conduction solution.

Starting with the conduction solution, we take a number of calculation super-steps. For a Rayleigh number close to critical (as in this example), the number of super-steps that needs to be taken before noticeable changes are detected can be quite large. After 100 super-steps, hardly any change is seen in the temperature in the model; however, a convection cell starts to form. This may be seen in Figure 4, where flow vectors are plotted on a background of temperature contours.

The steady-state solution has not been reached. As the simulation proceeds, the temperature begins to evolve in the model, and the convection cells develop further into a two-dimensional roll. Figure 5 and Figure 6 show the temperature contours and evolution of temperature with time at the control points after 700 super-steps.

Figure 4: Temperature contours and flow vectors after 100 super-steps (Ra = 42).

Figure 5: Temperature contours after 700 super-steps.

Figure 6: Evolution of temperature at 5 control points after 700 super-steps.

The simulation is continued for another 7000 super-steps. From the evolution of temperature (see Figure 7), it appears that the system has reached a steady state. The temperature contours at the end of the simulation are plotted in Figure 8. A contour plot of pore pressure at the end of the simulation is sketched in Figure 9.

Figure 7: Evolution of temperature at 5 control points after 7700 super-steps.

Figure 8: Temperature contours on a plane perpendicular to the x-axis after 7700 super-steps.

Figure 9: Pore pressure contours and flow vectors after 7700 super-steps.

The temperature contours shown in Figure 8 (2D roll in a cube at a Rayleigh number close to critical) may be compared to the two-dimensional analytical temperature contour values for mode 1, derived from linear perturbation theory (see, for example, Zhao et al. 1997) and plotted in Figure 10. The 2D analytical solution has the form

(2)\[{{T - T_0} \over {T_1 - T_0}} = C_1 {\rm cos} \bigl(\pi {x^{\ast} \over H} \bigr) {\rm sin} \bigl(\pi {y^{\ast} \over H} \bigr) + \bigl(1 - {z^{\ast} \over H} \bigr)\]

where \(x^{\ast} = x + {H \over 2}\), \(y^{\ast} = y + {H \over 2}\), \(z^{\ast} = z + {H \over 2}\) and \(C_1\) is a constant, defined here by requiring that numerical and analytical temperature values be equal at \(x^{\ast}\) = 0, \(z^{\ast} / H\) = 0.5.

Figure 10: Temperature contours, analytical steady-state solution, Rayleigh = \(4\pi^2\).

Saturated Porous Cube with a High Rayleigh Number

For this simulation at a high Rayleigh number, we use data similar to “HRL-Convection-Cube-LowRa.dat”, except for a higher value of the fluid thermal expansion property. A small perturbation, made by tilting the gravity vector slightly in the \(x\)-direction, is done to trigger the convection cell formation.

The resulting Rayleigh number for this model is approximately 381. The multiple three-dimensional cells obtained at steady state are shown in Figure 11. The data file is listed in “HRL-Convection-Cube-HighRa.dat”.

Note that the explicit calculation mode is used at this high Rayleigh number. The implicit calculation is not as beneficial in this case.

Figure 11: Steady-state temperature contours and flow vectors for a 1 × 1 × 1 box, Ra = 381.

Long Saturated Porous Box

The previous simulation is repeated, this time for a box with an aspect ratio of 8 × 1 × 1. A perturbation is done on the gravity vector, in the direction of the long dimension of the box. The cell pattern obtained at steady state is shown in Figure 12 and Figure 13. The pattern corresponds to a series of 2D rolls. The data file is listed in “HRL-Convection-LongBox.dat”.

Figure 12: Steady-state temperature contours and flow vectors for 8 × 1 × 1 box, Ra = 42.

Figure 13: Close-up view of temperature contours and flow vectors for 8 × 1 × 1 box, Ra = 42.

Grid Sensitivity Analysis

To test the influence of mesh density on the steady-state results, the convection simulation is repeated for a porous box with aspect ratio 4 × 1 × 1, at a Rayleigh number of 42, using three different regular grid resolutions: coarse (24 × 6 × 6 zones); medium (48 × 12 × 12 zones); and fine (96 × 24 × 24 zones). We used the same perturbation in gravity (long box direction) for all three cases. Two-dimensional rolls are obtained for all three grid resolutions. The temperature contours and flow vectors are shown in Figure 14 to Figure 16 for coarse, medium, and fine resolution, respectively. The data file for the simulations is provided in “HRL-Convection-LongBox-Coarse.dat” to “HRL-Convection-LongBox-Fine.dat”.

Figure 14: Steady-state temperature contours and flow vectors—coarse grid.

Figure 15: Steady-state temperature contours and flow vectors—medium grid.

Figure 16: Steady-state temperature contours and flow vectors—fine grid.

References

Lapwood, E. R. “Convection of a Fluid in a Porous Medium,” Proc. Cambridge Philos. Soc., 44, 508-521 (1948).

Zhao, C., H. B. Mühlhaus and B. E. Hobbs. “Finite Element Analysis of Steady-State Natural Convection Problems in Fluid-Saturated Porous Media Heated from below,” Int. J. Num. & Analy. Method. Geomech., 21, 863-881 (1997).

Data Files

HortonRogersLapwoodConvection.dat

program call 'HRL-Convection-Cube-LowRa'
program call 'HRL-Convection-Cube-HighRa'
program call 'HRL-Convection-LongBox'
program call 'HRL-Convection-LongBox-Coarse'
program call 'HRL-Convection-LongBox-Medium'
program call 'HRL-Convection-LongBox-Fine'
program call 'check'

HRL-Convection-Cube-LowRa.dat

; 3D convection in a porous layer
model new
model title '3D convection in a porous layer - Ra = 42'
fish automatic-create off
model configure fluid thermal
program call 'HRL-Convection-FISH' suppress
[setup]
; --- geometry ---
zone create brick size 12 12 12 point 0 (-5e3,-5e3,-5e3) ...
                                point 1 ( 5e3,-5e3,-5e3) ...
                                point 2 (-5e3, 5e3,-5e3) ...
                                point 3 (-5e3,-5e3, 5e3)
zone face skin
; --- thermal, fluid and mechanical models ---
zone thermal cmodel assign advection-conduction
zone fluid cmodel assign isotropic
zone cmodel assign elastic
; --- properties ---
;    (mechanical)
zone property density 2000 bulk=[_bulk] shear=[_shear]
;    (thermal)
zone thermal property conductivity=[_condu] specific-heat=[_mspec]
zone thermal property conductivity-fluid=[_fcond] specific-heat-fluid=[_fspec]
zone thermal property expansion-fluid=[_fthex] expansion=[_thexp]
zone thermal property temperature-reference=20.0
;    (fluid)
zone fluid property porosity=[_poros] permeability=[_mobil]
zone fluid property biot=1 undrained-thermal-coefficient=0
zone fluid biot off
zone gridpoint initialize fluid-modulus=2e5 
zone gridpoint initialize fluid-tension=-1e10
zone initialize fluid-density=[_fdens]
; --- initial conditions ---
zone gridpoint initialize pore-pressure 0.5e8 gradient 0 0 -1e4
zone initialize stress-xx -1.05e8 gradient 0 0 2.1e4
zone initialize stress-yy -1.05e8 gradient 0 0 2.1e4
zone initialize stress-zz -1.05e8 gradient 0 0 2.1e4
; --- boundary conditions ---
zone face apply velocity (0,0,0)  range group 'Bottom'
zone face apply velocity-x 0      range group 'West' or 'East'
zone face apply velocity-y 0      range group 'North' or 'South'
zone face apply temperature [_T1] range group 'Bottom'
zone face apply temperature [_T0] range group 'Top'
; --- settings ---
model gravity 0.01 0 [-_grav]
; --- for explicit calculations, use this:
model fluid   timestep fix 5.5e9
model thermal timestep fix 5.5e10
; --- for implicit calculations, use this:
;zone fluid implicit on
;zone thermal implicit on
;model fluid   timestep fix 5.5e10
;model thermal timestep fix 5.5e11
model fluid active off thermal on mechanical on
; --- histories ---
history interval 5000
model history fluid time-total
zone history pore-pressure position -4e3 -4e3 -4e3
zone history pore-pressure position -2e3 -2e3 -2e3
zone history pore-pressure position 0 0 0
zone history pore-pressure position 2e3 2e3 2e3
zone history pore-pressure position 4e3 4e3 4e3
model history thermal time-total
zone history temperature   position -4e3 -4e3 -4e3
zone history temperature   position -2e3 -2e3 -2e3
zone history temperature   position 0 0 0
zone history temperature   position 2e3 2e3 2e3
zone history temperature   position 4e3 4e3 4e3
; --- test ---
model mechanical active off
zone gridpoint fix velocity
; --- conduction solution ---
zone gridpoint initialize temperature 120 gradient 0 0 -2e-2
[ini_pp]
model save 'cube-lowRa-cd'
; --- convection ---
; --- Note: reset to vertical in this run ---
[supersolve(250)]
model save 'cube-lowRa-hf'
model gravity 0 0 [-_grav]       
zone thermal property specific-heat-fluid=[_fspec]
[supersolve(250)]
model save 'cube-lowRa-he'
[supersolve(250)]
[supersolve(250)]
model save 'cube-lowRa-h4'   ; 1000 supersteps, implicit calculation:stop here
[supersolve(2000)]             
;model save 'cube-lowRa-h6'  ; 3000 supersteps
[supersolve(2000)]             
;model save 'cube-lowRa-h8'  ; 5000 supersteps
[supersolve(2000)]             
model save 'cube-lowRa-h10'  ; 7000 supersteps
[supersolve(2000)]             
;model save 'cube-lowRa-h12' ; 9000 supersteps
[supersolve(2000)]             
;model save 'cube-lowRa-h14' ; 11000 supersteps
[supersolve(6000)]             
;model save 'cube-lowRa-h20' ; 17000 supersteps
[supersolve(10000)]            
;model save 'cube-lowRa-h30' ; 27000 supersteps
[supersolve(50000)]            
model save 'cube-lowRa-h40'  ; 77000 supersteps
program return

HRL-Convection-Cube-HighRa.dat

model new
model title '3D convection in a porous layer - Ra = 381'
model configure fluid thermal
program call 'HRL-Convection-FISH' suppress
[setup]
[global _fthex=1.4e-3]  ; increase fluid thermal expansion coefficient
; --- geometry ---
zone create brick size 18 18 18 point 0 (-5e3,-5e3,-5e3) ...
                                point 1 ( 5e3,-5e3,-5e3) ...
                                point 2 (-5e3, 5e3,-5e3) ...
                                point 3 (-5e3,-5e3, 5e3)
zone face skin
; --- thermal, fluid and mechanical models ---
zone thermal cmodel assign advection-conduction
zone fluid cmodel assign isotropic
zone cmodel assign elastic
; --- properties ---
;    (mechanical)
zone property density 2000 bulk=[_bulk] shear=[_shear]
;    (thermal)
zone thermal property conductivity=[_condu] specific-heat=[_mspec]
zone thermal property conductivity-fluid=[_fcond] specific-heat-fluid=[_fspec]
zone thermal property expansion-fluid=[_fthex] expansion=[_thexp]
zone thermal property temperature-reference=20.0
;    (fluid)
zone fluid property porosity=[_poros] permeability=[_mobil]
zone fluid property biot=1 undrained-thermal-coefficient=0
zone fluid biot off
zone gridpoint initialize fluid-modulus=2e5 
zone gridpoint initialize fluid-tension=-1e10
zone initialize fluid-density=[_fdens]
; --- initial conditions ---
zone gridpoint initialize pore-pressure 0.5e8 gradient 0 0 -1e4
zone initialize stress-xx -1.05e8 gradient 0 0 2.1e4
zone initialize stress-yy -1.05e8 gradient 0 0 2.1e4
zone initialize stress-zz -1.05e8 gradient 0 0 2.1e4
; --- boundary conditions ---
zone face apply velocity (0,0,0) range group 'Bottom'
zone face apply velocity-x 0     range group 'West' or 'East'
zone face apply velocity-y 0     range group 'North' or 'South'
zone face apply temperature [_T1] range group 'Bottom'
zone face apply temperature [_T0] range group 'Top'
; --- settings ---
model gravity 0.01 0 [-_grav]
; --- for explicit calculations, use this:
model fluid   timestep fix 2.4e9
model thermal timestep fix 2.4e10
model fluid active off thermal on mechanical on
; --- test ---
model mechanical active off
zone gridpoint fix velocity
; --- conduction solution ---
zone gridpoint initialize temperature 120 gradient 0 0 -2e-2
[ini_pp]
; --- convection ---
; perturbation
model gravity 0.01 0. [-_grav]
[supersolve(500)]
model save 'cube-highRa-hf'
; --- reset gravity to vertical
model gravity 0 0 [-_grav] 
[supersolve(500)]
model save 'cube-highRa-he'
[supersolve(500)]
[supersolve(500)]
;model save 'cube-highRa-h4'
[supersolve(2000)]
;model save 'cube-highRa-h6'
[supersolve(2000)]
;model save 'cube-highRa-h8'
[supersolve(2000)]
model save 'cube-highRa-h10'
program return

HRL-Convection-LongBox.dat

model new
model title '3D convection in a long porous box - Ra = 42 ? medium grid'
fish automatic-create off
model configure fluid thermal
program call 'HRL-Convection-FISH' suppress
[setup]
; --- geometry ---
zone create brick size 96 12 12 point 0 -4e4 -5e3 -5e3 ...
                                point 1  4e4 -5e3 -5e3 ...
                                point 2 -4e4  5e3 -5e3 ...
                                point 3 -4e4 -5e3  5e3
zone face skin 
; --- thermal, fluid and mechanical models ---
zone thermal cmodel assign advection-conduction
zone fluid cmodel assign isotropic
zone cmodel assign elastic
; --- properties ---
;    (mechanical)
zone property density 2000 bulk=[_bulk] shear=[_shear]
;    (thermal)
zone thermal property conductivity=[_condu] specific-heat=[_mspec]
zone thermal property conductivity-fluid=[_fcond] specific-heat-fluid=[_fspec]
zone thermal property expansion-fluid=[_fthex] expansion=[_thexp]
zone thermal property temperature-reference=20.0
;    (fluid)
zone fluid property porosity=[_poros] permeability=[_mobil]
zone fluid property biot=1 undrained-thermal-coefficient=0
zone fluid biot off
zone gridpoint initialize fluid-modulus=2e5 
zone gridpoint initialize fluid-tension=-1e10
zone initialize fluid-density=[_fdens]
; --- initial conditions ---
zone gridpoint initialize pore-pressure 0.5e8 gradient 0 0 -1e4
zone initialize stress-xx -1.05e8 gradient 0 0 2.1e4
zone initialize stress-yy -1.05e8 gradient 0 0 2.1e4
zone initialize stress-zz -1.05e8 gradient 0 0 2.1e4
; --- boundary conditions ---
zone face apply velocity (0,0,0) range group 'Bottom'
zone face apply velocity-x 0     range group 'West' or 'East'
zone face apply velocity-y 0     range group 'North' or 'South'
zone face apply temperature [_T1] range group 'Bottom'
zone face apply temperature [_T0] range group 'Top'
; --- settings ---
model gravity 0 0 [-_grav]           
;
zone fluid implicit on
zone thermal implicit on
model fluid timestep   fix 5.5e10 
model thermal timestep fix 5.5e11
; --- histories ---
model history fluid time-total
history interval 20
zone history pore-pressure position -4e3 -4e3 -4e3
zone history pore-pressure position -2e3 -2e3 -2e3
zone history pore-pressure position 0 0 0
zone history pore-pressure position 2e3 2e3 2e3
zone history pore-pressure position 4e3 4e3 4e3
model history thermal time-total
zone history temperature position -4e3 -4e3 -4e3
zone history temperature position -2e3 -2e3 -2e3
zone history temperature position 0 0 0
zone history temperature position 2e3 2e3 2e3
zone history temperature position 4e3 4e3 4e3
; --- test ---
model mechanical active off
zone gridpoint fix velocity
; --- conduction solution ---
zone gridpoint initialize temperature 120 gradient 0 0 -2e-2
[ini_pp]
model save 'longBox0-cd'
; --- convection ---
; perturbation
model gravity 0.01 0 [-_grav]
;
[supersolve(25)]
model save 'longBox0-hf'
; reset gravity to vertical
model gravity 0 0 [-_grav]  
[supersolve(25)]
model save 'longBox0-he'
[supersolve(25)]
[supersolve(25)]
;model save 'longBox0-h4'
[supersolve(200)]
;model save 'longBox0-h6'
[supersolve(200)]
;model save 'longBox0-h8'
[supersolve(200)]
;model save 'longBox0-h10'
[supersolve(200)]
;model save 'longBox0-h12'
[supersolve(200)]
;model save 'longBox0-h14'
[supersolve(600)]
;model save 'longBox0-h20'
[supersolve(1000)]
;model save 'longBox0-h30'
[supersolve(1000)]
model save 'longBox0-h40'
program return

HRL-Convection-LongBox-Coarse.dat

model new
model title '3D convection in a porous 4x1x1 box - Ra = 42 - coarse'
fish automatic-create off
model configure fluid thermal
program call 'HRL-Convection-FISH' suppress
[setup]
; --- geometry ---
zone create brick size 24 6 6 point 0 -2e4 -5e3 -5e3 point 1  2e4 -5e3 -5e3 ...
                              point 2 -2e4  5e3 -5e3 point 3 -2e4 -5e3  5e3
zone face skin                              
; --- thermal, fluid and mechanical models ---
zone thermal cmodel assign advection-conduction
zone fluid cmodel assign isotropic
zone cmodel assign elastic
; --- properties ---
;    (mechanical)
zone property bulk=[_bulk] shear=[_shear]
zone initialize density 2000 
;    (thermal)
zone thermal property conductivity=[_condu] specific-heat=[_mspec]
zone thermal property conductivity-fluid=[_fcond] specific-heat-fluid=[_fspec]
zone thermal property expansion-fluid=[_fthex] expansion=[_thexp]
zone thermal property temperature-reference=20.0
;    (fluid)
zone fluid property porosity=[_poros] permeability=[_mobil]
zone fluid property biot=1 undrained-thermal-coefficient=0
zone fluid biot off
zone gridpoint initialize fluid-modulus=2e5 
zone gridpoint initialize fluid-tension=-1e10
zone initialize fluid-density=[_fdens]
; --- initial conditions ---
zone gridpoint initialize pore-pressure 0.5e8 gradient 0 0 -1e4
zone initialize stress-xx -1.05e8 gradient 0 0 2.1e4
zone initialize stress-yy -1.05e8 gradient 0 0 2.1e4
zone initialize stress-zz -1.05e8 gradient 0 0 2.1e4
; --- boundary conditions ---
zone face apply temperature [_T1] range group 'Bottom'
zone face apply temperature [_T0] range group 'Top'
; --- settings ---
model gravity 0 0 [-_grav]          
zone fluid implicit on
zone thermal implicit on
model fluid timestep   fix 2.2e11 
model thermal timestep fix 2.2e12
;
model fluid active off thermal on mechanical on
; --- test ---
model mechanical active off
zone gridpoint fix velocity
; --- conduction solution ---
zone gridpoint initialize temperature 120 gradient 0 0 -2e-2
[ini_pp]
; --- convection ---
; perturbation
model gravity 0.01 0 [-_grav] 
[supersolve(10)]
model save 'longBox-coarse-hf'
; reset gravity to vertical
model gravity 0 0 [-_grav]
[supersolve(10)]
model save 'longBox-coarse-he'
[supersolve(10)]
[supersolve(10)]
;model save 'longBox-coarse-h4'
[supersolve(200)]
;model save 'longBox-coarse-h6'
[supersolve(200)]
;model save 'longBox-coarse-h8'
[supersolve(200)]
;model save 'longBox-coarse-h10'
[supersolve(200)]
;model save 'longBox-coarse-h12'
[supersolve(200)]
;model save 'longBox-coarse-h14'
[supersolve(600)]
;model save 'longBox-coarse-h20'
[supersolve(1000)]
model save 'longBox-coarse-h30'
program return

HRL-Convection-LongBox-Medium.dat

model new
model title '3D convection in a 4x1x1 porous box - Ra = 42 - medium'
fish automatic-create off
model configure fluid thermal
program call 'HRL-Convection-FISH' suppress
[setup]
; --- geometry ---
zone create brick size 48 12 12 point 0 -2e4 -5e3 -5e3 ...
                                point 1 2e4 -5e3 -5e3  ...
                                point 2 -2e4 5e3 -5e3  ...
                                point 3 -2e4 -5e3 5e3
zone face skin
; --- thermal, fluid and mechanical models ---
zone thermal cmodel assign advection-conduction
zone fluid cmodel assign isotropic
zone cmodel assign elastic
; --- properties ---
;    (mechanical)
zone property bulk=[_bulk] shear=[_shear]
zone initialize density 2000 
;    (thermal)
zone thermal property conductivity=[_condu] specific-heat=[_mspec]
zone thermal property conductivity-fluid=[_fcond] specific-heat-fluid=[_fspec]
zone thermal property expansion-fluid=[_fthex] expansion=[_thexp]
zone thermal property temperature-reference=20.0
;    (fluid)
zone fluid property porosity=[_poros] permeability=[_mobil]
zone fluid property biot=1 undrained-thermal-coefficient=0
zone fluid biot off
zone gridpoint initialize fluid-modulus=2e5 
zone gridpoint initialize fluid-tension=-1e10
zone initialize fluid-density=[_fdens]
; --- initial conditions ---
zone gridpoint initialize pore-pressure 0.5e8 gradient 0 0 -1e4
zone initialize stress-xx -1.05e8 gradient 0 0 2.1e4
zone initialize stress-yy -1.05e8 gradient 0 0 2.1e4
zone initialize stress-zz -1.05e8 gradient 0 0 2.1e4
; --- boundary conditions ---
zone face apply temperature [_T1] range group 'Bottom'
zone face apply temperature [_T0] range group 'Top'
; --- settings ---
model gravity 0 0 [-_grav]          
zone fluid implicit on
zone thermal implicit on
model fluid timestep   fix 5.5e10 
model thermal timestep fix 5.5e11
;
model fluid active off thermal on mechanical on
; --- test ---
model mechanical active off
zone gridpoint fix velocity
; --- conduction solution ---
zone gridpoint initialize temperature 120 gradient 0 0 -2e-2
[ini_pp]
; --- convection ---
; perturbation
model gravity 0.01 0 [-_grav]
[supersolve(25)]
model save 'longBox-medium-hf'
; reset gravity to vertical
model gravity 0 0 [-_grav]
model save 'longBox-medium-he'
[supersolve(25)]
[supersolve(25)]
;model save 'longBox-medium-h4'
[supersolve(200)]
;model save 'longBox-medium-h6'
[supersolve(200)]
;model save 'longBox-medium-h8'
[supersolve(200)]
;model save 'longBox-medium-h10'
[supersolve(200)]
;model save 'longBox-medium-h12'
[supersolve(200)]
;model save 'longBox-medium-h14'
[supersolve(600)]
;model save 'longBox-medium-h20'
[supersolve(1000)]
model save 'longBox-medium-h30'
program return

HRL-Convection-LongBox-Fine.dat

model new
model title '3D convection in a porous 4x1x1 box - Ra = 42 - fine'
fish automatic-create off
model configure fluid thermal
program call 'HRL-Convection-FISH' suppress
[setup]
; --- geometry ---
zone create brick size 96 24 24 point 0 -2e4 -5e3 -5e3 ...
                                point 1 2e4 -5e3 -5e3  ...
                                point 2 -2e4 5e3 -5e3  ...
                                point 3 -2e4 -5e3 5e3
zone face skin
; --- thermal, fluid and mechanical models ---
zone thermal cmodel assign advection-conduction
zone fluid cmodel assign isotropic
zone cmodel assign elastic
; --- properties ---
;    (mechanical)
zone property bulk=[_bulk] shear=[_shear]
zone initialize density 2000
;    (thermal)
zone thermal property conductivity=[_condu] specific-heat=[_mspec]
zone thermal property conductivity-fluid=[_fcond] specific-heat-fluid=[_fspec]
zone thermal property expansion-fluid=[_fthex] expansion=[_thexp]
zone thermal property temperature-reference=20.0
;    (fluid)
zone fluid property porosity=[_poros] permeability=[_mobil]
zone fluid property biot=1 undrained-thermal-coefficient=0
zone fluid biot off
zone gridpoint initialize fluid-modulus=2e5
zone gridpoint initialize fluid-tension=-1e10
zone initialize fluid-density=[_fdens]
; --- initial conditions ---
zone gridpoint initialize pore-pressure 0.5e8 gradient 0 0 -1e4
zone initialize stress-xx -1.05e8 gradient 0 0 2.1e4
zone initialize stress-yy -1.05e8 gradient 0 0 2.1e4
zone initialize stress-zz -1.05e8 gradient 0 0 2.1e4
; --- boundary conditions ---
zone face apply temperature [_T1] range group 'Bottom'
zone face apply temperature [_T0] range group 'Top'
; --- settings ---
model gravity 0 0 [-_grav]
zone fluid implicit on
zone thermal implicit on
model fluid timestep   fix 1.3e10
model thermal timestep fix 1.3e11
;
model fluid active off thermal on mechanical on
; --- test ---
model mechanical active off
zone gridpoint fix velocity
; --- conduction solution ---
zone gridpoint initialize temperature 120 gradient 0 0 -2e-2
[ini_pp]
; --- convection ---
; perturbation
model gravity 0.01 0 [-_grav]
[supersolve(50)]
model save 'longBox-fine-hf'
; reset gravity to vertical
model gravity 0 0 [-_grav]
[supersolve(50)]
model save 'longBox-fine-he'
[supersolve(50)]
[supersolve(50)]
;model save 'longBox-fine-h4'
[supersolve(200)]
;model save 'longBox-fine-h6'
[supersolve(200)]
;model save 'longBox-fine-h8'
[supersolve(200)]
;model save 'longBox-fine-h10'
[supersolve(200)]
;model save 'longBox-fine-h12'
[supersolve(200)]
;model save 'longBox-fine-h14'
[supersolve(600)]
;model save 'longBox-fine-h20'
[supersolve(400)]
;model save 'longBox-fine-h30'
[supersolve(1000)]
;model save 'longBox-fine-h40'
[supersolve(1000)]
;model save 'longBox-fine-h50'
[supersolve(1000)]
;model save 'longBox-fine-h60'
[supersolve(1000)]
model save 'longBox-fine-h70'
program return