3.10. Physical Models
Last updated
Was this helpful?
Last updated
Was this helpful?
Function: The momentum conservation equations solved by GASFLOW-MPI include body forces. The most common body force is that due to gravity. Because the user can orient the computational mesh arbitrarily with respect to the gravity vector, the code by default sets the body force acceleration term to zero in all three directions.
Input Array: the user has to define the values of the following variables in NAMELIST group xput:
gx Acceleration due to gravity in the i- (x- or r-) direction (cm/s^2). Default = 0.
gy Acceleration due to gravity in the j- (y- or θ-) direction (cm/s^2). Default = 0.
gz Acceleration due to gravity in the k- (z-) direction (cm/s^2). Default = 0.
Besides specifying the gravity term, the code offers the option of starting a calculation with a pressure gradient in the fluid that is in equilibrium with its own body weight. This option is specified through the ihystat variable in NAMELIST group xput:
ihystat Option flag for imposing an initial hydrostatic pressure gradient in the fluid cells according to the acceleration components gx, gy, and gz: 1 means ON; 0 means OFF (default).
When these options are selected, the pressure that the user has provided with the gasdef input statements will not be used. Instead, the analytical hydrostatic solution for the pressure field with be set in all computational cells and gasdef statements at the beginning of the simulation according to
The solution of this equation requires a reference point to which the integration may be associated, an initial pressure at an exact location. We provide this information in the XPUT namelist variables reference pressure (pamb0) and locations of reference point (xamb0, yamb0, zamb0). The default values of these variables are as follows:
pamb0 = 1.013253e+06,
xamb0 = 0.0,
yamb0 = 0.0,
zamb0 = 0.0,
If an initial hydrostatic pressure gradient is required, it is important the user investigates the rami-fications of using these initial variables. The reason is that the initial hydrostatic pressure is determined by an analytical evaluation of the above equation and not the numerical approximation to the entire momentum equations (Navier-Stokes equations) that GASFLOW-MPI solves, and there is likely to be a difference between the analytical and numerical methods.
Definition: Transport due to gradient diffusion is modeled by Fick's law-type fluxes in the conservation equations of mass, energy, and momentum.
Input Array: The user must explicitly specify them using the following input variables in NAMELIST group xput.
idiffme Option flag for mass and energy diffusion: 1 means ON, 0 means OFF (default).
idiffmom Option flag for momentum diffusion: 1 means ON, 0 means OFF (default).
Note that these option flags apply to the diffusion terms, which include both molecular and turbulent transport. However, by default, the code does not calculate turbulent diffusion.
For momentum diffusion, if idiffmom is set to 1 and no turbulence model is activated (tmodel = 'none'), there is an additional option to calculate the molecular diffusion coefficient. The user can specify whether the input dynamic viscosity, μ (in units of poise or g/cm·s), or the input kinematic viscosity, ν = μ/ρ (cm2/s), is to be used for calculating the diffusional momentum fluxes or viscous stresses. This is done via the following variables in NAMELIST xput:
muoption Option specifying whether the input variable nu and cmug is to be used for viscous stress calculation:
1 means nu is used (default)
2 means cmug is used.
nu Kinematic viscosity, ν = μ/ρ (cm2/s). Default = 0.15.
cmug Dynamic viscosity, μ (g/cm·s). Default = 1.8e-04.
For the diffusion of mass species and energy, the code uses the input variables schmidt and prandtl, respectively, if the user has selected the itopt = 0 option, to determine the appropriate diffusion coefficient.
The Navier-Stokes equation are nonlinear, and when the flow speed exceeds certain criteria, they become unstable in the sense that the solution exhibits oscillatory or fluctuating behavior.
Therefore, the equations solved in GASFLOW-MPI are actually Reynolds averaged, with the dependent variables being mean quantities after the turbulent fluctuations are removed. However, the fluctuations cannot be simply removed. The equations retain terms that are averages of products like u'v' and u'φ', where u', v' are the fluctuating velocity components and φ' is the fluctuating scalar quantity such as internal energy or a mass species concentration.
In GASFLOW-MPI, all turbulent fluxes are modeled like a molecular diffusion process. Turbulent momentum fluxes (stresses) are modeled as
where δui/δxj is the velocity gradient tensor. This is the so-called Boussinesq’s hypothesis, which states that turbulent transport can be modeled as a diffusion process with an eddy viscosity, μt. In GASFLOW, we further assume that the turbulent diffusion coefficients for mass and energy are both the same as that for momentum and are given by μt/ρ, where ρ is the fluid density.
Input Array. In GASFLOW-MPI, two turbulent models are available.
The choice can be specified via the input character variable tmodel in NAMELIST group xput:
tmodel Symbol designating the type of turbulence model to be used:
'none', no turbulence model, i.e., only molecular diffusion;
'alg', algebraic turbulence model;
'ke', κ-ε turbulence model;
‘ko’, κ-ω turbulence model;
‘sstko’, SST κ-ω turbulence model;
‘des’, detached eddy simulation;
‘les’, large eddy simulation.
Note that the turbulence model invoked via tmodel is only in effect when diffusion calculations have also been specified, i.e., when idiffmom and/or idiffme has been set to 1.
Wall Functions. When rigid no-slip conditions are specified and turbulence is activated, the user may wish to use wall functions rather than resolving the boundary layers. The options are as follows:
iwallfunc = 0 ; default value
; no wall functions are active
iwallfunc = 1 ; no-slip conditions must to active
; assumes smooth walls
iwallfunc = 2 ; no-slip conditions must to active
; assumes rough walls and krough must be specified
The initial turbulent conditions in the ambient can be defined by using two of the following variables,
tkeamb0 Initial turbulent kinetic energy in the ambient. Default = 10 (cm2/s2).
epsamb0 Initial turbulent dissipation rate in the ambient. Default = 10000 (cm2/s3).
sclamb0 Initial turbulent scale in the ambient.
The algebraic turbulence model used in GASFLOW assumes that the turbulent viscosity can be written as
where κ is the turbulent kinetic energy per unit mass, ρ is the fluid density, 𝑙 is a turbulent length scale, and cμ is a constant coefficient. In this model, it is assumed that the turbulent kinetic energy is a specified fraction of the fluid mean kinetic energy, i.e.,
The parameters to be input for the algebraic turbulence model are cμ, 𝑙, and 𝑓, which are represented by the following variables in NAMELIST group xput:
cmu Constant cμ for the algebraic turbulence model. Default = 0.09.
clength Length scale, 𝑙, for the algebraic turbulence model (cm). Default = 30.48.
fractke Fraction of the mean kinetic energy, 𝑓, used in the algebraic turbulence model to determine the turbulent kinetic energy. Default = 0.1
In the k-ε model, it is assumed that the turbulent viscosity can be written as:
where k is the turbulent kinetic energy (per mass unit), ρ is the fluid density, ε is the rate of dissipation of turbulent kinetic energy, cμ is a constant coefficient. The turbulent kinetic energy is determined from solution of its transport equation:
The terms on the right hand side are, in order, diffusional transport, production by viscous stresses, production by buoyancy, viscous dissipation, and generation from other sources. Similarly, the rate of dissipation, ε, is determined from solution of its transport equation:
The following variables in NAMELIST group xput are used to set the κ-ε model parameters:
cmuke cμ used in the κ-ε turbulence model. Default = 0.09.
c1ke c1 used in the κ-ε turbulence model. Default = 1.44.
c2ke c2 used in the κ-ε turbulence model. Default = 1.92.
sigmak σκ used in the κ-ε turbulence model. Default = 1.
sigmae σε used in the κ-ε turbulence model. Default = 1.3.
The turbdef array described below is used for defining initial and boundary conditions for the κ-ε model and is in the NAMELIST group xput.
turbdef(1,*) Beginning i mesh index (cell face number).
turbdef(2,*) Ending i mesh index (cell face number).
turbdef(3,*) Beginning j mesh index (cell face number).
turbdef(4,*) Ending j mesh index (cell face number).
turbdef(5,*) Beginning k mesh index (cell face number).
turbdef(6,*) Ending k mesh index (cell face number).
turbdef(7,*) Block number (must be 1 for GASFLOW-MPI).
turbdef(8,*) Integer pointer to location in tkeval array for value of turbulent kinetic energy.
turbdef(9,*) Integer pointer to location in epsval array for value of turbulence dissipation rate.
turbdef(10,*) Integer pointer to location in sclval array for value of turbulence length scale.
turbdef(11,*) Time (s) at which “turbulence definition” begins.
turbdef(12,*) Time (s) at which “turbulence definition” ends.
Example.
The user must define a valid pointer for the array tkeval in turbdef(8,*), and a valid pointer for the array epsval in turbdef(9,*). The input in turbdef(10,*) will be ignored. The following input illustrates the use of turbdef to define an initial condition and a boundary condition for solving the κ and ε transport equations:
turbdef(1:12,1) = 1, 50, 1, 32, 1, 45, 1, 2, 1, 0, 0.0, 0.0,
turbdef(1:12,2) = 0, 1, 2, 5, 4, 6, 1, 1, 2, 0, 0.0, 1.e9,
tkeval = 1.0e3, 0.0,
epsval = 0.0, 800.0,
The first turbdef defines an initial condition for κ and ε. The values of both κ and ε are initially 0. The second turbdef specifies a boundary condition at the –i boundary, between j-indices of 2 and 5, and between k-indices of 4 and 6.
The DES turbulent model, a hybrid of LES and RANS, has proven to be a successful tool in various industrial applications. GASFLOW-MPI utilizes a k-ε based DES model to capture turbulent behavior, incorporating equations for turbulent kinetic energy and dissipation rate. Unlike the traditional k-ε model, the DES model dynamically selects between the standard k-ε and LES models based on the turbulent state. In regions with larger mesh sizes, the DES model behaves similarly to the standard k-ε model.
The turbulent kinetic energy is determined from solution of its transport equation:
The terms on the right hand side are, in order, diffusional transport, production by viscous stresses, production by buoyancy, viscous dissipation, and generation from other sources.
Similarly, the rate of dissipation, ε, is determined from solution of its transport equation:
The following variables in NAMELIST group xput are used to set the des model parameters:
cmuke cμ used in the κ-ε turbulence model. Default = 0.09.
c1ke c1 used in the κ-ε turbulence model. Default = 1.44.
c2ke c2 used in the κ-ε turbulence model. Default = 1.92.
sigmak σκ used in the κ-ε turbulence model. Default = 1.
sigmae σε used in the κ-ε turbulence model. Default = 1.3.
The Smagorinsky model is employed in GASFLOW-MPI to calculate the SGS turbulent viscosity due to its simplicity and practicality. In the Smagorinsky model, turbulent viscosity is proposed by assuming that the turbulent energy transfers from the large resolved scales to the small subgrid scales and the turbulent energy dissipation by the small subgrid scales are in equilibrium. And the turbulent viscosity could be expressed as
C_les is a constant in Smagorinsky model and Δ is the LES filter width. Different from RANS model, the turbulent viscosity is a function of the mesh size in LES model. The filter width Δ is computed from the geometry size of the computational cell,
The following variables in NAMELIST group xput are used to set the des model parameters:
cles C_les used in the les turbulence model. Default = 0.1.
The following variables in NAMELIST group xput are used to set the SST κ-ω model parameters
sigmak3 σ_k3 used in the SST κ-ω turbulence model.
sigmao3 σ_ω3 used in the SST κ-ω turbulence model.
alphao3 α_3 used in the SST κ-ω turbulence model.
betao3 β_3 used in the SST κ-ω turbulence model.
The coefficient of the standard k-ω model are multiplied by blending F_1 and the coefficients of the transformed κ-ε equations are multiplied by a function 1-F1, the closure coefficients for SST_ κ-ω turbulent model are calculated by
sigmas σ_k1 used in the κ-ω turbulence model. Default = 0.5.
sigmao σ_ω1 used in the κ-ω turbulence model. Default = 0.5.
alphao α_1 used in the κ-ω turbulence model. Default = 5/9.
betao β_1 used in the κ-ω turbulence model. Default = 0.075.
sigmak σ_k2 used in the transformed κ-ε turbulence model. Default = 1.
sigmaeo σ_ω2 used in the transformed κ-ε turbulence model. Default = 0.856.
c1ke α_2 used in the transformed κ-ε turbulence model. Default = 0.44.
c2ke β_2 used in the transformed κ-ε turbulence model. Default = 0.0828.
The turbdef array used for defining initial and boundary conditions for the SST κ-ω model is identical as the turbdef array of κ-ε model. The user must define two valid pointers for the arrays tkeval in turdef(8,*), epsval in turdef(9,*) and sclval in turdef(10,*). Two of these three arrays must be defined. For the SST κ-ω model, the value of omega, omgval, can be calculated by
omgval = epsval/(cmuke*tkeval), when tkeval and epsval are defined;
omgval = (cmuke*sclval^2)/epsval, when epsval and sclval are defined;
omgval = tkeval^0.5/sclval, when tkeval and sclval are defined.
GASFLOW-MPI models the combustion of hydrogen in air with the following single, global reaction:
where the coefficients a, b, c, d, e, f, and g are generally referred to as the stoichiometric coefficients. The reaction rate is determined from a modified Arrhenius law, where
and a = 2, b = 1, and e = 2.
The progress variable is defined as
The forward and backward reaction rates are, respectively,
When the mixture is fuel rich, then the molar density change for the oxidizer is calculated first, and the molar densities are determined as
The reaction rate parameters and heat of reaction are fixed in the code. Therefore, the user only has to specify the following variable, in NAMELIST group xput, in order to turn on the hydrogen combustion model in the calculation:
iburn = +1, Hydrogen combustion model on (forward reaction only)
= 0, Hydrogen combustion model off (default)
= -1, Hydrogen combustion model on (both forward and reverse reactions)
Note: ignitior model should be only used together with the combustion model in the section 3.10.4.1.
At user-specified locations GASFLOW-MPI first checks the gas composition to determine if the mixture is combustible, based on the Shapiro Diagram for hydrogen, oxygen, and steam mixtures.
The locations where the ignitor model is applied are defined with the input array ignitdef, which is in the xput NAMELIST group. The parameter ignitdef is described below:
ignitdef(1,*) i cell index.
ignitdef(2,*) j cell index.
ignitdef(3,*) k cell index.
ignitdef(4,*) Block number (must be 1 for GASFLOW-MPI).
ignitdef(5,*) Ignitor type:
= 0, continuous ignitor;
>0, spark ignitor.
When the ignitdef(5,*) is greater than zero, then the ignitor types are defined in the input array spxider and the input variable spxigdt, which are also in the xput NAMELIST group. Input array spxidef and input variable spxigdt are defined below:
spxigdef(1,igtyp) Time when ignitor is first turned on.
spxigdef(2,igtyp) Time when ignitor is turned off.
spxigdef(3,igtyp) Time interval between sparks.
spxigdt Time of the spark duration. Default = 0.001 s.
In the spxigdef definitions, igtyp is the ignitor type number defined by ignitdef(5,*). The maximum number of ignitor definitions (i.e., ignitdef(*,n)) and the maximum number of ignitor types are both currently set with a parameter statement to 300.
If the user wishes to allow auto ignition to occur, they must activate the following model by setting ignitaut > 0 in the NAMELIST xput. The auto ignition temperature, autotemp, can also be set to a value other than it's default value of autotemp = 800.0 (Kelvin).
We have incorporated the Kumar flammability limits for hydrogen-oxygen-diluent mixtures. The gas concentration criteria proposed for this model is shown in the Table below
ignitaut = 0,
GASFLOW simulates hydrogen combustion in the limit of the Kumar criterion using the Arrhenius correlation with the physical temperature of each computational cell. The reference temperature T**_**ref applied in the Arrhenius correlation is set to the physical temperature in each computational cell.
ignitaut = 1
Hydrogen burn is simulated only for gas temperatures above auto ignition with the parameter ignitaut=1 in cells that fulfill the Kumar criterion. The temperature T_ref (m) for cell m in the Arrhenius correlation is then set different from the cell temperature. In a loop over all cells the code first determines if the cell temperature is below the input value autotemp for the auto ignition temperature. In this case the temperature Tref(m) for cell m in the Arrhenius correlation is set to 1 K only. The cell temperature Tcell (m) defines T_ref (m) only when it is above autotemp.
ignitaut = 2
It prevents auto ignition to occur PAR cells defined by rcombdef. The reference temperature for PAR cells are always set to 1 when the users sets the ignitaut=2.
Cells with ignitors
In active ignitor cells defined by ignitdef , spxigdef, and spxigdt the code resets the value of Tref to 2000 K. A value of Tref (m) >1 then denotes all cells which fulfill necessary conditions for ignition. But hydrogen burn is only simulated in these cells if their hydrogen and oxygen concentrations are within the lean and rich flammability limits as defined by the Kumar criterion.
A general "combustion progress variable" transport equation has been solved model the premixed turbulent combustion in GASFLOW-MPI.
In order to turn on these hydrogen combustion models, the user has to first specify iburn = 4 in NAMELIST group xput.
Since the key to this modeling approach is the source term, the user then has to specify isourcexi in NAMELIST group xput to tell the code which model will be used to obtain the source term.
isourcexi
= 0, source term off;
= 1, Arrhenius rate model;
= 2, model based on progress variable gradient (default);
= 3, under development;
= 4, under development;
= 5, Eddy dissipation model.
Please note currently the turbulent combustion models (isourcexi = 2, 3, 4 and 5) are coupled with the κ-ε, DES and LES turbulence model. Therefore,these turbulence models must be turned on to calculate the source term in the reaction progress variable transport equation. The initial conditions of the turbulent properties could strongly influence the combustion process and the time step of the calculation. Therefore, the user should estimate the initial conditions of experiments or simulations and give reasonable initial values to the turbulence model. There are two ways to set up the initial turbulent conditions.
The first way is to use two of the variables: tkeamb0, epsamb0 and sclamb0. For example, to set up the initial turbulent kinetic energy and dissipation rate in the whole computational domain, we could define
tkeamb0 = 1.0e2,
epsamb0 = 8.0e3,
Another more flexible way is to use turbdef which allows the user to define various conditions in different sub-domains. For example,
turbdef(1:12,1) = 2, ‘im1’, 2, ‘jm1’, 2, ‘km1’, 1, 1, 1, 0, 0.0, 0.0,
turbdef(1:12,2) = 2, 3, 2, 3, 2, 3, 1, 2, 2, 0, 0.0, 0.0,
tkeval = 1.0e2, 2.5e3,
epsval = 8.0e3, 2.0e4,
When isourcexi equals 2, the users has to define iturbflame in NAMELIST group xput to select a correlation for calculation of turbulent flame speed.
iturbflame
= 0, use laminar flame speed;
= 2, use Kawanabe correlation;
= 3, use Peters correlation;
= 4, use Zimont correlation (default);
= 5, use Zimont-Mesheriakov correlation;
= 6, use Schmidt correlation.
To obtain the reference initial pressure, density and thermal diffusivity to calculate the Damkoehler number for the models above (when isourcexi = 2), the user has to define a reference point to represent the unburnt mixture using the variables below.
iref_unburnt i cell index for reference point of unburnt mixture;
jref_unburnt j cell index for reference point of unburnt mixture;
kref_unburnt k cell index for reference point of unburnt mixture;
iblkref_unburnt Block number for reference point of unburnt mixture.
Note: this ignitor model should only be used together with the combustion models in Chapter 3.10.4.3.
The cells where the xi_ignitor model is applied are defined with the input array xi_ignitdef, which is in the xput NAMELIST group. Basically, the xi_ignitor model describes how long it takes for the reaction progress variable, xi, in the specified ignition cell to change from 0 to 1.
The parameter, xi_ignitdef, is described below:
xi_ignitdef(1,*) beginning i mesh index (cell face number).
xi_ignitdef(2,*) ending i mesh index (cell face number).
xi_ignitdef(3,*) beginning j mesh index (cell face number).
xi_ignitdef(4,*) ending j mesh index (cell face number).
xi_ignitdef(5,*) beginning k mesh index (cell face number).
xi_ignitdef(6,*) ending k mesh index (cell face number).
xi_ignitdef(7,*) block number (must be 1 for GASFLOW-MPI)
xi_ignitdef(8,*) time when ignitors are first turned on (s).
xi_ignitdef(9,*) time of the ignition duration (s). It means the duration for xi to increase from zero to one. A value close to the current time step or less is recommended if the exact ignition duration is unknown.
The asterisk (*) should be replaced by an integer that identifies the particular xi_ignitdef definition. You need to check xi_ignitdef(8,*) and xi_ignitdef(9,*) if the ignition process is too slow or no ignition occurs.
In many cases, users may have to restart a combustion simulation based on a distribution simulation (iburn = 0). When restarting from iburn > 0, memory has to be allocated for the combustion simulation from the beginning of the distribution simulation by defining a parameter in $xput
iburn_malloc = 1,
Heat transfer between the fluid and any solid structures are, by default, not calculated. Therefore, in a problem that, for instance, involves hydrogen combustion (iburn > 0) but the heat-transfer model is not explicitly requested, a temperature rise will still be calculated, which is only valid if the process is adiabatic or if the problem time scale is fast compared with the heat-transfer time scale.
To activate heat-transfer and steam-condensation calculations, using models based on the argument of analogies between momentum, heat, and mass transfer, the user must specify the following variable in NAMELIST group rheat:
ihtflag Option flag to turn on heat-transfer and steam-condensation calculations:
1 means ON; 0 means OFF (default).
When ihtflag > 0 other input values control heat transfer options through the rheat NAMELIST input stream for activating specific heat transfer relationships.
When the default model is used (ihtcoef = 4), the multiplication factor cons1 is preset.
cons1 Multiplication factor for the heat-transfer coefficient (applies to convection and condensation/evaporation). Default = 3.6.
There are other heat and mass transfer parameters that the user can control:
crelax Relaxation factor for depleting liquid H2O. Default = 0.01 (1/s)
icond Flag for condensation on structures: 1 means ON; 0 means OFF. Default = 1
ienh Film-model enhancement from Bird Stewart Lightfoot: 1 means ON; 0 means OFF; so far only for film condensation. Default = 1.
ifcca Flag for Chilton-Colburn analogy: 1 means ON; 0 means OFF. Default = 1.
These are the wall function options:
Within the general heat transfer model, when ihtflag > 0, the user may activate a radiation heat transfer model. The radiation model is activated by the rheat NAMELIST variable irad > 0. There are a number of options:
There are several other rheat NAMELIST variables connected with the radiation heat transfer model that the user has some control:
emismat(*) Material emissivity. Each material in the heat transfer material data base has a default emismat(*) = 1.
clamda radiation mean free path (cm).
If clamda > 0, then the user is specifying a constant radiation mean free path;
otherwise, the mean free path is computed using the method of Lechner, “Spectral and Total Emissivity of Water Vapor and Carbon Dioxide”, Combustion and Flame, Vol. 19, pp. 33-48, 1972.
itmaxr maximum number of iterations allowed for the solution of the radiation heat transfer equation (default = 500).
epsrad convergence criteria for the solution of the radiation heat transfer equation (default = 10^-6).
The particulate materials are characterized by class and size. The material class corresponds to a material density. Each material class may have multiple particle size.
The input parameter pmass(lc,ls) is the total mass of each particle class and size that is to be input over the total problem time. The class material density is set in prhoin(lc) and the particle diameters are specified for each class and size by pdiamin(lc,ls). From this input data the number of real particles, npreal(lc,ls), for each class and size and the total number of real particles, npsum, are computed. Next, the volumes, volp(lc), within the computational mesh into which each particle class is to be set are determined from the input parameters specifying the maximum and minimum x-, y-, and z-coordinates of the particle rectangular volumes including xpe(lc), xpw(lc), ypn(lc), yps(lc), zpt(lc), and zpb(lc). In some cases, e.g. to define a source volume for spray water droplets, it would be more convenient to enable a definition of a cylindrical or ring shaped volume. The user can turn on the option iparscyl(lc) = 1 together with the input parameters, rpsmin(lc), rpsmax(lc), thetapsmin(lc) and thetapsmax(lc). rpsmin(lc) and rpsmax(lc) are the minimum and maximum radius of the ring. The center of the ring is located at ((xpe(lc)+xpw(lc))/2,(ypn(lc)+yps(lc))/2). The z-coordinate of the cylinder is between zpt(lc), and zpb(lc). thetapsmin(lc) and thetapsmax(lc) are the minimum and maximum angle of the ring.
When the particle model is activated, the GASFLOW random number generator is also activated with a seed based on the day and time of execution. The input parameter init_random allows the initialization of the seed.
The user selects the total number of simulation particles for each class and size to be input and transported by the code. This input parameter is npinpt(lc,ls) (number particles input total for each class and size). From the total number of simulation particles input, npinpt(lc,ls), the number of real particles for each class and size, npreal(lc,ls), and the total number of real particles, npsum, the code computes the number of simulation particles for each class and size npsim(lc,ls). Using the information available, the real mass represented by each simulation particle, rmpp(lc,ls), is computed.
The input parameter twpinp(lc) designates the initial time when simulation particles of each class are to be input into the system. The total number of particles for each class may all be set in the mesh initially or may be injected as a function of time. This choice is governed by the input parameters tinjt(lc) (total injection time) and pinpdt(lc) (particle input Δt). To inject the particles as a function of time, set the initial time when particles are to be injected, twpinp(lc), the total injection time, tinjt(lc), and the time interval between particle injections, pinpdt(lc). From this input data the code computes the number of injection times and the number of simulation particles for each class and size, npinj(lc,ls), to be injected at the appropriate intervals. To inject a single volume of particles, input the time of injection, twpinp(lc), and values for tinjt(lc) and pinpdt(lc) such that the ratio tinjt(lc)/pinpdt(lc) is less than 1.0. The code will set the number of injection times, ninjts(lc), to 1.0. The value of pinpdt(lc) also must be greater than the problem time.
When restarting a problem from a previous particle run, set the parameter twpinp(lc) to a value greater than the final problem time to inhibit the injection of additional particles. If injection of particles is wanted at the time determined in the previous run, do not set a value for twpinp(lc) in the restart input parameter file. The correct value will be read from the restart tape file.
Input Array: The variables itpcl, itpsz, and msp are included in the xput NAMELIST group and are used to set array sizes in the particle model. Therefore, whenever solatype is greater than zero, these three inputs must be supplied in the xput NAMELIST group. A calculation can be started with solatype equal to zero and then restarted with solatype greater than zero. The restart must include the three inputs itpcl, itpsz, and msp. Once these three inputs have been used, all subsequent restarts must include the same values for these input variables. Initial particle velocities could be specified to each particle class with the part NAMELIST variable partvel(:,lc).
The following aerosol input parameters are in the xput NAMELIST group:
itpcl Total number of particle classes.
itpsz(lc) Number of initial particle sizes for each particle class.
msp Total number of simulation particles allowed in this calculation.
solatype
= 0.0 implies a hydrodynamic calculation without the particle model.
= 1.0 implies a hydrodynamic calculation with the particle model.
= 2.0 implies a particle-model-only calculation with gas velocities fixed in time.
init_random Particle random number generator option:
init_random=0 (default) implies the initialization of the seed.
init_random=n implies n is the seed which allows the user to initialize the random number generator with the same seed regardless of the time of execution. This option is helpful when reproducing the same results for each separate calculation is necessary.
ipblkin(lc) Mesh block in which particle class, index lc, is initially located (must be 1 for GASFLOW-MPI).
ipclin(lc) Particle class number input.
npinpt(lc,ls) Total number of simulation particles to be input for each class and size.
partvel(:,lc) Flag for specifying initial particle velocity (cm/s).
= 0; no initial particle velocity specified for particle class lc (default).
partvel(1,lc) = 1; particle velocity specified for particle class lc and signaling GASFLOW to assign particle class initial velocities partvel(2:4,lc). partvel(2,lc) specifies particle class lc u-component initial velocity (default = 0.0). partvel(3,lc) particle class lc v-component initial velocity (default = 0.0). partvel(4,lc) particle class lc w-component initial velocity (default = 0.0).
partvel(1,lc) = 2; particle velocity specified for particle class lc and signaling GASFLOW to assign particle class initial velocities the same as the local gas velocity. partvel(2:4,lc) will not be used.
partvel(1,lc) = 3; particle velocity specified for particle class lc and signaling GASFLOW to assign particle class initial velocities using the parameters defined by partvel(2:4,lc). This option provides an way to inject a particle with an random angle relative to the axial direction of the particle injection.
partvel(2,lc) specifies particle class lc injection axial direction: 1: X direction; 2: Y direction; 3: Z direction (default = 0.0).
partvel(3,lc) specifies particle class lc injection velocity magnitude and direction (default = 0.0).
partvel(4,lc) specifies particle class lc maximum random angle relative to the axial direction (default = 0.0).
pinpdt(lc) Time interval between particle inputs for each particle class (s).
pdiamin(lc,ls) Particle diameter input for each class and size (cm).
pmass(lc,ls) Total real particle mass of each class and size (g).
prhoin(lc) Material density input for each particle class (g/cm^3).
pcsubpin(lc) Material specific heat capacity input for each particle class (erg/(g-K)).
ptcondin(lc) Material thermal conductivity input for each particle class (erg/(s-cm-K)).
ptempin(lc) Temperature input for each particle class (K).
tinjt(lc) Total injection time for each particle class (s).
twpinp(lc) Time when to input particles of each class (s).
iparscyl(lc) The volume shape of the particle injection source in the Cartesian coordinate:
0: rectangular;
1: cylindrical.
Used in Cartesian coordinate only (default = 0).
rpsmin(lc) Minimum radius of initial particle ring for each particle class (cm). The ring center is located at ((xpe(lc)+xpw(lc))/2,(ypn(lc)+yps(lc))/2) . Used only when iparscyl(lc) = 1.
rpsmax(lc) Maximum radius of initial particle ring for each particle class (cm). ). The ring center is located at ((xpe(lc)+xpw(lc))/2,(ypn(lc)+yps(lc))/2) . Used only when iparscyl(lc) = 1.
thetapsmin(lc) Minimum angle (between 0 and 360) of initial particle block for each particle class. Used only when iparscyl(lc) = 1.
thetapsmax(lc) Maximum angle (between 0 and 360) of initial particle block for each particle class. Used only when iparscyl(lc) = 1.
xpe(lc) Maximum x-coordinate of initial particle block for each particle class (cm).
xpw(lc) Minimum x-coordinate of initial particle block for each particle class (cm).
ypn(lc) Maximum y-coordinate of initial particle block for each particle class (cm).
yps(lc) Minimum y-coordinate of initial particle block for each particle class (cm).
zpb(lc) Minimum z-coordinate of initial particle block for each particle class (cm).
zpt(lc) Maximum z-coordinate of initial particle block for each particle class (cm).
Particle material density, prhoin(lc), particle size, pdiamin(lc,ls), and total mass of particulate material, pmass(lc,ls).
Particle input volume and location: iparscyl(lc), rpsmin(lc), rpsmax(lc) (Note: used only when iparscyl(lc) = 1), thetapsmin(lc), thetapsmax(lc) (Note: used only when iparscyl(lc) = 1), xpw(lc), xpe(lc), yps(lc), ypn(lc), zpb(lc), zpt(lc).
Particle injection time and rate: twpinp(lc) (Initial time when to input particles), tinjt(lc) (Total injection time), pinpdt(lc) (Time increment between particle injections).
Total number of simulation particles to be input over the entire injection time: npinpt(lc,ls).
Determine the total number of simulation particles to be simulated, msp, the number of classes, itpcl, and the number of different sizes for each class, itpsz(lc).
Define the location of the particles to be input for each class and size. The mesh block, ipblkin(lc), in which the particles of class lc are placed, has a default value of 1. Specify the rectangular volume in the computational mesh into which the simulation particles are to be placed by inputting the maximum and minimum x-, y-, and z-coordinates of the particle volume for each class; these are, respectively, xpe(lc), xpw(lc), ypn(lc), yps(lc), zpt(lc), and zpb(lc). init_random allows to select the way of initializing the seed for particle random number generator.
Define the number of aerosols to be modeled. Set the number of different particle materials, i.e., particle classes, itpcl. Specify the number of particle sizes to be modeled in each of these classes, itpsz(lc). The default value for each of these is 1. The particle class number, ipclin(lc), is currently set by default and does not need to be changed.
Characterize the real particle material to be modeled. Specify the total material mass of each class and size, pmass(lc,ls). Set the density of the material of each class, prhoin(lc), and specify the particle diameter for each class and size, pdiamin(lc,ls). (These same particle densities and diameters are used for calculating the simulation particle dynamics in the code particle transport model.)
Define the simulation particle injection number and mode. Set the total number of simulation particles, which includes all the particles of each class and size to be input over all time, npinpt(lc,ls). Select the initial time and duration for which particles of each class will be injected by setting the time when to input particles, twpinp(lc), and the total injection time for each class tinjt(lc). Set the time interval between particle injections during this period, pinpdt(lc). To inject a single volume of particles, i.e, a one-time injection, input the time of injection, twpinp(lc). Set the value of pinpdt(lc) to a value greater than the final problem time. Select a value for tinjt(lc) such that the ratio tinjt(lc)/pinpdt(lc) is less than 1.0.
Determine the input volume dimension in the inflow velocity direction to be compatible with Uin . dt (cm).
Compute the particle number concentrations, cpinj(lc,ls), for each particle class and size (no. of particles/cm^3).
Compute particle input volume, volp(lc), from xpe(lc), etc., for each class (cm^3).
Compute the number of particles injected each injection time for each class and size, npinj(lc,ls) = cpinj(lc,ls) ·volp(lc).
Select a maximum injection time, tinjt(lc) (s), and the particle injection time interval, pinpdt(lc) (s).
Compute the number of injection times anticipated in the total problem time, ninjts(lc) = tinjt(lc) / pinpdt(lc).
Compute the number of particles for each class and size, npsim(lc,ls)= npinj(lc,ls)·ninjts(lc).
Sum the total number of particles in each class and size. This gives the input parameter npinpt(lc,ls), which is the total number of particles to be injected.
Compute the mass of each simulation particle class and size injected. The product of this and the number of particles for each class and size is the input parameter pmass(lc,ls), pmass(lc,ls)= [(pi/6) • **pdiamin(lc,ls)^**3 • prhoin(lc) ] • npsim(lc,ls) (g).
If imarker = 1, the particles input are used as marker particles and moved at the local fluid velocity. Most of the particle physics is bypassed. The default value of imarker is 0. In this case, the complete particle dynamics algorithm is executed.
The particle transport model includes the inertial, drag, and gravitational forces. In cylindrical coordinates, the centrifugal force also is modeled. The drag coefficient in the Stokes flow region is a function of the flow Reynolds number and, consequently, of the relative particle-fluid velocity. This requires an iterative procedure (the Newton-Raphson method is used) to simultaneously solve for the drag coefficient and the particle velocity. niterp is input to set the maximum number of iterations. The default value for niterp is set to 30. The slip correction factor, scfacin(mxpcl,mxpsz), reduces the drag force for small particles and should be used for particles with diameters less than about 1.0 μm.The default value for scfacin(mxpcl,mxpsz) is 1.0.
The turbulent diffusion of the particles is characterized by the particle turbulent diffusion coefficient, tdcp. The numerical model for turbulent diffusion computes the diffusion velocity components for each particle. In this model, the fluid velocity components at the locations of each particle are the sum of the mean velocity components and the random turbulent velocity component.
The user can specify the model used to move each particle by the parts NAMELIST variable fluid_velocity. The user has the option to interpolate two or four fluid velocities for each coordinate direction by specifying fluid_velocity.
The coefficient of kinematic viscosity, visf, is set by default to the value of nu input in NAMELIST group xput and typically does not need to be changed.
In the evaluation of the particle Reynolds Number, the user now has the option to input constant values for fluid density and kinetic viscosity by inputting rhogas (default = 0.0012) in the $parts namelist, and nu (default = 0.1535) in the xput namelist.
In the two-way momentum coupling, the continuous fluid phase is impacted by the discrete particle phase (and vice versa). To activate the two-way coupling option in GASFLOW, solatype must be equal to 1 which allows the hydrodynamic calculation with the particle model. The parts NAMELIST variable fpcpmom should be set to 2 to activate two-way momentum coupling. The parts NAMELIST variable fpcpme should be set to 2 to activate two-way mass and energy coupling. Note: when fpcpme is set to 2, fpcpmom must be set to 2.
To enable the phase change (condensation and evaporation) of the particle, the user has to set iparph2o(lc) = 1. Currently, phase change of the water droplets are supported. Ordinary differential equations of mass and energy for each water droplet are then solved iteratively. cflpmass needs to be set to ensure the numerical stability when the diameter of the water droplet approaches zero during evaporation. In this case, the time step will be limited due to the water droplet evaporation:
pdiammin is the minimum diameter of water droplet given by the user. When the diameter of water droplet is less than pdiammin, the water droplet is considered to be completely vaporized and the particle information will be erased. Note: pdiammin is used for water droplet evaporation only.
The following input parameters for the particle transport model are in the parts NAMELIST group.
cflpmass Time step control number due to the mass change of water droplet. Default = 0.25.
fluid_velocity Flag for methodology of computing local particle velocity in fluid cell.
= 0, uses local particle cell fluid velocities for computing the fluid velocity at the particle location.
= 1, uses all fluid velocities for computing the fluid velocity at the particle location.
fpcpmom Flag for two-way momentum coupling.
=1, one-way momentum coupling
=2, two-way momentum coupling
fpcpme Flag for two-way mass and energy couplings.
= 1, one-way mass and energy coupling ,
= 2, two-way mass and energy coupling (used for water droplets only in the current version).
fpcpturb Flag for two-way turbulence couplings.
= 1, one-way turbulence coupling ,
= 2, two-way turbulence coupling (reqiures fpcpmom = 2 and tdcp /= 0).
iparph2o(lc) Droplets with heat and mass transfer for the class lc.
= 0, solid particle or water droplet with constant diameter without phase change. Default = 0.
= 1, water droplets with heat and mass transfer. Computing the evaporation of the water droplets and steam condensation on the surface of the water droplets (reqiures fpcpme = 2 and ipdep(lc) = 2).
= 2, computing the diameter change of the aerosol droplets by Köhler equation without affecting the gas phase. It is valid for droplets with diameter less than 20 μm.
iparwfilm Water droplets as the mass and energy sources to the water film when the water droplets adhere on the solid surfaces. Default = 1.
= 0, disable.
= 1, enable (reqiures fpcpme = 2, ipdep(lc) = 2 and iparpclass(lc) = 1).
imarker Marker particle flag.
= 0, the complete particle dynamics algorithm is executed. Default = 0.
= 1, particles are moved at the local fluid velocity (most of the particle physics is bypassed).
inpvol Input parameter to choose actual or fictitious volume assigned to each particle.
= 0, actual volume is computed;
= 1, δx ·δy ·δz of cell in which particle is located.
local_rhogas
= 0 (default), the constant values for rhogas is used.
> 0, actual spatial and time dependent local fluid density is used.
local_visf
= 0 (default), the constant values for visf is used.
> 0, actual spatial and time dependent local fluid kinetic viscosity is used.
niterp Number of iterations for the Newton iteration cycle for determining drag-induced particle velocity. Default = 30.
pdiammin The minimum diameter during water droplet evaporation. Default = 3e-4 cm.
scfacin(lc, lc) Cunningham slip correction factor for each particle class and size. Default = 1.0.
tdcp Turbulent diffusion coefficient for particles (cm^2/s).
> 0, a constant particle turbulent diffusion coefficient.
= –1.0, then μ/(ρ·S·cp) is used for the turbulent diffusion coefficient for particles, where S·cp is the particle Schmidt Number and μ is the total viscosity of the gas mixture including the turbulent viscosity.
= –2.0, turbulent kinetic energy is used to determine the turbulent diffusion coefficient. In this case, if tmodel = 'ke' or tmodel = 'alg', the turbulent viscosity is used to determine the turbulent diffusion coefficient.
stokes(lc,ls)
= 0, solve the equations of particle motion.
> 0, simulate the Stokes flows.
visf Coefficient of kinematic viscosity of fluid.
wmax Maximum argument of error function used in inverse error function table. Default = 2.0.
Set up particle input as described in Section 3.10.6.1.
Select the value of imarker. The default value of 0 utilizes the full particle dynamics algorithm.
If the particle diameters are 1.0 μm or less, set the slip correction factor, scfacin(mxpcl,mxpsz). The default value is 1.0.
If particle turbulent diffusion is to be modeled, set the coefficient tdcp to an estimated value. The default value is 0.0 cm2/s.
Choose a value for niterp, typically between 10 and 40, depending on the particle size and accuracy required (default = 30).
Use the default values for inpvol and wmax. The values of rhogas and visf can be set as constant or dynamic by specifying local_rhogas and local_visf.
The default value of the particle deposition flag, ipdep(lc), is 0. If no particle adhesion upon impact for particle class lc is to be modeled, no further input is necessary. At the other extreme, if all particles adhere upon impact, an ipdep(lc) value of 2 is input. A deposition model based on the particle material properties (and some stochastic considerations) is invoked if ipdep(lc) is set equal to 1.
Example: if the first class of particle material (class number 1) is a liquid mist, then ipdep(1) = 2 is appropriate because the impacting material will most likely not bounce. Then, if the second particle class is a harder material and the deposition model is to be used to choose whether a particle bounces or adheres on impaction, ipdep(2) = 1 is input.
The deposition model uses the Dahneke theory of rebound to compute a critical rebound velocity, at which there is a 50% probability of bounce. The particle energy is a function of the Hamaker constant, hca, the particle diameter, pdiamin(lc,ls), the separation distance of the particle and substrate, sdz, and particle bulk mechanical properties, bmck(lc), which are a function of the material Poisson ratio, poisrt(lc), and Young’s modulus, yngmod(lc).
The coefficient of restitution, core, has been experimentally determined for many particle-surface combinations. The values typically range from 0.90 to 0.99 for a hard particle impinging onto a hard surface. The input parameter core has a default value of 0.96. Experimental data shows the trend for the coefficient of restitution to reach this maximum value at a given incident velocity and then almost immediately the ratio of rebound velocity to incident velocity begins to decrease as the incident velocity increases. This phenomenon is programmed into the code and does not need any input data.
Deposition is in a real sense a stochastic process that follows the general trend of the theoretical and empirical models developed from and compared with experimental data. Because of this, it is a reasonable assumption that some small, unknown percentage of the particles that impact a surface will adhere. To account for this, the parameter depperc(lc) is input. The default value is set to 5.0%. If the substrate is dirty or there are other reasons to believe a smaller or larger percentage of the particles will adhere on impact, this value can be changed.
Input Array: The following input parameters are in the parts NAMELIST group:
core Coefficient of restitution of particle material.
depperc(lc) Percentage of particles that immediately deposit upon impact.
hca Hamaker constant (ergs).
ipdep(lc) Particle deposition flag input for each class:
if ipdep(lc)= -1, no adhesion with elastic rebound;
= 0, no adhesion with a computed rebound velocity;
= 1, adhesion determined by deposition model;
= 2, all impacting particles adhere.
ndxpd Number of longitudinal sections in x-direction in which particle deposition (only on bottom of mesh) is monitored: if - 1, no plot.
pdiamin(lc,ls) Particle diameter input for each class and size (cm).
poisrt(lc) Poisson ratio of particle class material.
prhoin(lc) Material density input for each particle class (g/cm^3).
sdz Separation distance of particle and substrate (cm).
yngmod(lc) Young’s modulus of particle material (dynes/cm^2).
Set up particle input as described in the latter sections
Select a value of ipdep(lc) for each particle class being modeled. There is no deposition of particles when the default value of 0 is set.
Use the default values for the Hamaker constant, hca, the particle-substrate separation distance, sdz, and the coefficient of restitution, core. These input parameters may be changed if there is a reason to do so.
Input appropriate material property values, Poisson ratio, poisrt(lc), and Young’s modulus, yngmod(lc), for all particle classes for which the deposition model is used, that is, when ipdep(lc) = 1. The default values are those for steel.
Use the default value (5%) for the percentage of particles that deposit on impact, depperc(lc); or, considering the condition of the substrate and any other pertinent information, set it to a better guess.
If a plot is wanted of the total particle mass deposited on the “floor” of the mesh (i.e., the bottom face of the k = 2 layer of cells) as a function of distance in the x-direction, set the input parameter ndxpd to the number of sections into which the distance from the mesh minimum x-coordinate to the maximum x-coordinate is to be divided and in which the deposition is to be monitored. If this information is not wanted, do not set this parameter, and the default value is such that this computation and plot will be skipped.
If time-history plots are wanted for the computational mesh cell faces onto which particles have deposited, input the necessary information according to instructions in Section 9.1.5, Graphics and Tabular Particle Data Output.
The default value for the particle entrainment flag, intrn, is 0. If particle entrainment is not to be modeled, no further input is required. If particle entrainment is to be modeled, set intrn to 1, and the entrainment model is used to suspend deposited particles when the fluid velocity equals or exceeds the computed particle threshold suspension velocity.
The adhesive force is the van der Waals inter-surface molecular force and is a function of the Hamaker constant, hca, the particle diameter, pdiamin(lc,ls), the separation distance of the particle and the substrate, sdz, and particle bulk mechanical properties, bmck(lc), which is a function of the material Poisson ratio, poisrt(lc), and Young’s modulus, yngmod(lc). The drag and lift forces depend on the gas stream velocity. The friction force is proportional to the coefficient of sliding friction and is set to a value of 0.45, which is the value suggested from experimental studies.
If the input parameter dbl is set to a value less than 0, the code will compute a boundary layer thickness as a function of a Reynolds number and the distance along the wall from the point at which the turbulent fluid initially contacts the wall. It is assumed that this is the point at which the fluid enters the computational mesh. However, if this model is not compatible with the problem setup, an estimate of a constant boundary layer thickness may be set by the parameter dbl.
The input parameters cdrf and sdzrf are multipliers for the computed drag coefficient and the usual surface-substrate separation distance, respectively, that allow a qualitative representation of the roughness of the particle and substrate surfaces. The default values are 10.0 and 100.0, respectively. These values were determined by comparison of computations with one set of experimental data and require further evaluation.
Input Array. The following input parameters are in the parts NAMELIST group:
intrn Input parameter to choose particle entrainment option:
= 0, no particle pickup (default and no more input is required.);
= 1, pickup option is on. dbl Boundary layer thickness. Default = 10 cm.
cdrf Aerodynamic drag coefficient roughness factor. Default =10.
sdzrf Particle-substrate separation distance roughness factor. Default = 100.
Setup Procedure
Set up particle input and deposition input as described in the latter sections.
Select a value of intrn. There is no entrainment of particles for the default value of 0. Input a value of 1 to turn on the entrainment model.
Set dbl to –1.0 for the code to compute the boundary layer thickness, or set an estimate of a constant boundary layer thickness. The default value is 10.0 cm.
Use the default values of cdrf and sdzrf. These may be changed if better estimates are known.
The default value of the particle cloud model flag, icloud, is 0. If the particle cloud model is not to be used, no further input is necessary. If the cloud model is to be used (icloud = 1), then a value for the particle cloud diffusion coefficient, pcdc, is chosen. The default value is 0. An estimated value for this input parameter is the value of the turbulent diffusion coefficient, tdcp.
The cloud model is designed to be used with particle monitors at selected locations in the computational domain. The total number of monitors is specified by the input parameter ntmntr, which has a default value of 0. The selected number of monitors must not be greater than the maximum number of monitors, mntrmx, which is set to the default value 100. The monitor locations are specified by the input parameters xm, ym, and zm for each monitor.
Input Array. The input parameters are in the $parts NAMELIST group:
icloud Particle cloud model flag.
= 0, cloud model off (default);
= 1, cloud model is on.
pcdc Particle cloud diffusion coefficient. Default = 0.
If pcdc is greater than zero, it is used for the particle cloud diffusion coefficient.
If pcdc = –1.0, then the particle cloud diffusion coefficient is equal to μ /(ρ*Sc), where μ is the total viscosity of the gas mixture including the turbulent viscosity, and Sc is the Schmidt number.
xm x-coordinate of each particle cloud monitor.
ym y-coordinate of each particle cloud monitor.
zm z-coordinate of each particle cloud monitor.
Setup Procedure
Choose icloud=1 to turn on the particle cloud model.
Set pcdc to an estimated value.
Choose the total number of particle cloud monitors, ntmntr, and their locations, xm, ym, and zm.
Specify in the graphics input parameter pthp 'pmntr' for time-history plots of the real particle cloud mass at each monitor.
Setting the NAMELIST group xput input parameter tddt to the appropriate time interval between tape dumps will generate a GASFLOW-MPI restart file (see next Section) file from which to restart a run. A common procedure in running the aerosol model is to set solatype, also in NAMELIST group xput, to a value of 0.0, for a fluid solution only, and generating a steady-state fluid flow field. Then, solatype = 2.0 is set for particle transport only, and particles are injected into this steady flow field. Parameters that may need changing when restarting from a tape dump include: twfin, maxcyc, solatype, nrsdump, pthpt0 and twpinp.
GASFLOW-MPI allows the user to specify a sump, which is basically a water film that has a constant thickness and a time-dependent temperature. The relevant phenomenon is shown in Figure below.
We see that the sump mass is determined by five sources:
Droplet depletion or rainout of the liquid component of the fluid field directly into the sump;
Liquid films draining from slabs and walls directly into the sump;
Phase change (condensation and/or evaporation) at the sump surface;
Sump to Sump mass exchange, and
Direct mass addition.
Input Array. To activate the sump model, the input parameter ihtflag in NAMELIST RHEAT must be none zero and positive, and the liquid (h2ol) and vapor (h2o) components of water must be specified with the mat input variable in NAMELIST XPUT. The sump model or the old constant film thickness model is activated by inputting information in the NAMELIST RHEAT input stream with the cfilmdef input variable array. This input variable array, cfilmdef(*,n), is defined for the nth surface as:
cfilmdef(1,n) beginning i mesh index (cell face number).
cfilmdef(2,n) ending i mesh index (cell face number).
cfilmdef(3,n) beginning j mesh index (cell face number).
cfilmdef(4,n) ending j mesh index (cell face number).
cfilmdef(5,n) k mesh index.
cfilmdef(6,n) k mesh index.
cfilmdef(7,n) block number (must be 1 for GASFLOW-MPI)
cfilmdef(8,n)
> 0, constant film thickness (must be less than 100 cm)
< 0, sump number (must be a negative integer)
The asterisk (*) should be replaced by an integer that identifies the particular cfilmdef definitions. GASFLOW supports 500 definitions for cfilmdef.
The temperature of the sump, sumptemp, as a function of time, sumptime, can be controlled, and input in either degrees Celsius or Kelvin. This can be accomplished via the input variable nsumppts and the array variables sumptemp and sumptime in NAMELIST group rheat. The definitions are as follows:
nsumppts Number of sump temperatures in the sumptemp and sumptime arrays. Absolute value of nsumppts must be less than 100.
< 0, sumptemp values are in °C.
>0, sumptemp values are in K.
sumptemp Sump temperatures < 100 values.
sumptime Time for the sump temperatures < 100 values.
When specifying a sump, the user must
Define a horizontal wall (walls) or obstacles (mobs) that coincides exactly with the sump location.
Define a constant film thickness for the sump with the cfilmdef input variable that exactly coincides with the defined wall or obstacle from 1.
Make certain that 'h2o' is listed in the problem composition, mat, and additionally 'h2ol' if the two-phase homogeneous equilibrium model (HEM) is to be used.
Input a sump temperature as a function of time through the paired single-dimensioned arrays sumptemp and sumptime.
Example
An example is as follows:
We assume the same geometry as the example from Figure below.
The sump is located at cell edge k = 2 and covers the entire bottom of the computation domain, i.e., 1 < i < 11 and 1 < j < 11, and a wall is constructed to indicate this sump.
The HEM is activated.
A constant film thickness is established to exactly coincide with the wall from 2.
Eight paired temperature and time values controlling the sump temperature in degrees Celsius or Kelvin are defined.
The input stream is shown here:
$xput
…
mobs = 3, 10, 1, 11, 5, 8, 1, 1, ; solid obstacle walls
2, 10, 1, 11, 3, 3, 1, 2, ; horizontal wall
2, 2, 1, 11, 6, 10, 1, 2, ; vertical wall
1, 11, 1, 11, 2, 2, 1, 2, ; sump surface
holes = 5, 7, 4, 7, 6, 8, 1, 0, 0, 0, 0, 0, 1, ; top hole
8, 9, 5, 6, 5, 8, 1, 0, 0, 0, 0, 1, 1, ; thru hole
8, 9, 5, 6, 2, 4, 1, 1, 1, 1, 1, 1, 1, ; wall hole
mat = 'h2', 'n2', 'o2', 'h2o', 'h2ol', ; components -> HEM
cfilmdef = 1, 11, 1, 11, 2, 2, 1, 5.0, ; sump film thickness
…
$end
$rheat
...
nsumppts = -8, ; sump paired temperature and time values in Celsius
sumptemp = 25.0, 37.0, 50.0, 45.0, ; sump temperature
60.0, 62.0, 63.0, 63.0, ; sump temperature
sumptime = 0.0, 2800.0, 6200.0, 9480.0, ; sump time
12600.0, 14400.0, 18000.0, 30000.0, ; sump time
...
$end
cfilmdef(8,n) = -1, we must provide additional sump parameters or characteristics. This is accomplished with the sumpchar(*,n) input variable array in the NAMELIST RHEAT input stream, where n is the sump number (1 for this version). So this input variable array, sumpchar(*,n), is defined for the nth sump as:
sumpchar(1,n) initial sump temperature
>0 for degrees Kelvin and
< 0 for degrees Celsius.
sumpchar(2,n) initial sump depth (cm).
sumpchar(3,n) maximum sump depth (cm).
Additional important input in the NAMELIST RHEAT input stream is:
filmth
>0, maximum film thickness (default is 0.1 cm). When the film exceeds this value, the excess mass and energy is transferred to the sump.
< 0, a liquid film is initialized on all slab, wall, and sink surfaces with a thickness of |filmth|.
Several tables of sump sources can also be input in the NAMELIST RHEAT input stream. These include the direct source of mass (sumpmas) and energy (sumpengy) and the energy source due to radio nuclide decay energy (sumprn). These tables can be specified as functions of time with input variable array (sumptime). Each of these input sump source variable arrays may accept up to 99 values where the scalar variable nsumppts defines the actual number of table entries. For completeness we define these variables as:
sumptime number of input values for each of the sump source tables in seconds
>0, sumpengy is in ergs/s
< 0, sumpengy is in ergs/(g·s)
sumpengy energy source associated with sumpmas, which is the direct energy source in ergs/s or ergs/(g·s), depending upon the sign of sumptime.
sumpmas mass source associated with sumpengy, which is the direct mass source in g/s.
sumprn energy source associated with the decay of fission products, radio nuclides, in ergs/s.
Example
$rheat
…
nsumppts = 4,
sumptime = 0.0, 750.0, 1250.0, 3600.0,
sumpmas = 0.0, 1000.0, 2000.0, 0.0,
sumpengy = 0.0, 1.2e+10, 2.0e+10, 0.0,
sumprn = 0.0, 6.0e+12, 8.0e+12, 1.0e+13,
...
$end
The location in the i, j plane for each recombiner is given in the rcombdef array described below. The rcombdef array is in the rheat NAMELIST group.
rcombdef(1,*) Beginning i mesh index (cell face number).
rcombdef(2,*) Ending i mesh index (cell face number).
rcombdef(3,*) Beginning j mesh index (cell face number).
rcombdef(4,*) Ending j mesh index (cell face number).
rcombdef(5,*) Beginning k mesh index (cell face number).
rcombdef(6,*) Ending k mesh index (cell face number).
rcombdef(7,*) Block number (must be 1 for GASFLOW-MPI).
rcombdef(8,*) Actual flow area into the recombiner irrespective of the computational mesh for recombiner type 1-4 and 10-18. NOTE: new features in GASFLOW-MPI. For recombiner type 5-9 and 19, rcombdef(8,*) is the recombiner capability defined by the user. The default value is 1.0 meaning the recombiner is fully used. 0.5 means the recombination capability is only 50% of the recombiner.
rcombdef(9,*) Recombiner type.
= 1, implies NIS recombiner model with chemistry based on compositions at the location identified by the offset level, and forced volumetric flow through the recombiner based on the correlation for the NIS recombiner;
= 2, implies Siemens recombiner model with chemistry based on compositions at the location identified by the offset level, and forced volumetric flow through the recombiner based on the correlation for the Siemens recombiner;
= 3, implies Siemens recombiner model with chemistry based on compositions in the catalytic reaction volume according to the correlations for Siemens Type FR 90/1 Type 3 allows early downflow from containment convection and flow reversal to upflow from catalytic recombination
= 4, implies GRS recombiner foil model with chemistry based on compositions in the fluid cells with common faces to the specified foil
= 5, implies Siemens correlation for box type FR90/1–750, this box was used also in GX tests
= 6, implies Siemens correlation for box type FR90/1–150
= 7, implies Siemens correlation for box type FR90/1-320
= 8, implies Siemens correlation for box type FR90/1-960
= 9, implies Siemens correlation for box type FR90/1-1500
= 10-15 , implies old NIS PARs from Biblis A LOCA simulation for reruns of archive case
= 16, PAR model for Reinecke net recombiner no longer used
= 17, temperature correction to Areva PAR correlation no longer used
= 18, New NIS PAR validated with test THAI HR 14
= 19, Areva PAR FR380 from THAI HR Tests
rcombdef(10,*) Offset. It can be only set to -1 due to the limitation of layers of ghost cells in the parallel code.
rcombdef(11,*) Threshold value for hydrogen concentration. Hydrogen concentration must be above this value before recombiner will operate. User can input this choice. Default is 0.03.
If rcombdef(11,*) is zero and rcombdef(9,*) is one, then rcombdef(11,*) will be set to 0.008.
If rcombdef(11,*) is zero and rcombdef(9,*) is two, then rcombdef(11,*) will be set to 0.031.
If rcombdef(11,*) is zero and recombdef(9,*) is >=5 and <=9 then recombdef(11,*) will be set to 0.02.
If rcombdef(11,*) is zero, and rcombdef(9,*) is >=10 and <=15 rcombdef(11,*) will be set to 0.008.
If rcombdef(11,*) is zero and rcombdef(9,*) is 18 rcombdef(11,*) will be set to 0.015 as in THAI test HR14.
If rcombdef(11,*) is zero and rcombdef(9,*) is 19 rcombdef(11,*) will be set to 0.02 as in THAI HR3.
If rcombdef(11,*) is less than zero, then rcombdef(11,*) is the time in seconds when the recombinder will begin to operate.
rcombdef(12,*) The user can input any time constant above zero. Default value is set to 100s if the user inputs zero. If the time constant is input as zero for recombiner type 1, and 10 to 15 then the time constant is set to 1800.0 s. The time constant is not used for recombiner type 4.
The asterisk (*) should be replaced by an integer that identifies the particular recombiner definition. GASFLOW-MPI supports 500 definitions for recombiners.
The time history of the total recombination rate of all PARs is stored as variable sumh2loss in the plothist.nc file. The time history of the recombined hydrogen mass itemized for each PAR is stored as variable recmass(ncomb) in the plothist.nc file with ncomb being the total number of PARs. ncomb is automatically determined from the input rcombdef.
When specifying a recombiner, the user must
Define at least hydrogen, oxygen, and water vapor in the problem composition.
Define the recombiner box or foil.
Define the catalytic reaction volume, the type of recombiner, offset for the NIS and Siemens forced volumetric flow rate models, and other recombiner characteristics defined above.
Example
We assume the same geometry as the example from Figure above, but without the internal obstacles and walls other than the recombiner geometries.
Introduce a NIS recombiner box type 1 (see 1 in Figure below).
Introduce a Siemens recombiner box type 2 (see 2 in Figure below).
Introduce a Siemens Type 3 FR90/1 recombiner box (see 3 in Figure below).
Introduce a catalytic zone type 4 GRS foil on each side of a defined wall.
The input stream is shown here:
$xput
…
walls = 9, 10, 6, 6, 6, 9, 1, 2, ;NIS north wall (type 1)
10, 10, 5, 6, 6, 9, 1, 2, ;NIS east wall (type 1)
9, 9, 5, 6, 6, 9, 1, 2, ;NIS west wall (type 1)
9, 10, 5, 5, 6, 9, 1, 2, ;NIS south wall (type 1)
3, 4, 6, 6, 2, 8, 1, 2, ;Siemens north wall(type 2)
4, 4, 5, 6, 2, 8, 1, 2, ;Siemens east wall (type 2)
3, 3, 5, 6, 2, 7, 1, 2, ;Siemens west wall (type 2)
3, 4, 5, 5, 2, 8, 1, 2, ;Siemens south wall (type 2)
3, 4, 5, 6, 8, 8, 1, 2, ;Siemens top wall(type 2)
8, 9, 6, 6, 2, 5, 1, 2, ;Siemens north wall(type 3)
9, 9, 5, 6, 2, 5, 1, 2, ;Siemens east wall (type 3)
8, 8, 5, 6, 2, 5, 1, 2, ;Siemens west wall (type 3)
8, 9, 5, 5, 2, 5, 1, 2, ;Siemens south wall (type 3)
6, 6, 4, 8, 7,10, 1, 2, ;GRS foil model (type 4)
mat = 'h2','n2','o2','h2o','h2ol',
…
$end
$rheat
…
rcombdef = 9, 10, 5, 6, 7, 8, 1, 1.0e+04, 1, -2, 0, 0, ; NIS (type 1)
3, 4, 5, 6, 4, 5, 1, 2.0e+02, 2, -3, 0, 0, ; Siemens (type 2)
8, 9, 5, 6, 3, 4, 1, 9.6e+03, 3, 0, 0, 0, ; Siemens (type 3)
5, 7, 4, 8, 7, 10, 1, 5.0e+04, 4, 0, 0, 0, ;GRS foil (type 4)
matcomb = 2, ; par material for foil type 4
walldef = 2, 0.1, 0.0, 0.0, 0.0, ; GRS foil material 2
…
$end
For all Siemens type recombiners (rcombdef(9,*) = 5-9), new input variables are now available for users to define the coefficients, k1 and k2, in the equation (2-134) in the Theory Manual. The input variables can be defined by the users in $heat in ingf file
rcomb5_k1 k1 used in recombiner type 5. Default = 6.10e-6.
rcomb5_k2 K2 used in recombiner type 5. Default = 7.4.
rcomb6_k1 k1 used in recombiner type 6. Default = 0.48e-6.
rcomb6_k2 k2 used in recombiner type 6. Default = 0.58.
rcomb7_k1 k1 used in recombiner type 7. Default = 1.0e-6.
rcomb7_k2 k2 used in recombiner type 7. Default = 1.2.
rcomb8_k1 k1 used in recombiner type 8. Default = 3.1e-6.
rcomb8_k2 k2 used in recombiner type 8. Default = 3.7.
rcomb9_k1 k1 used in recombiner type 9. Default = 13.7e-6.
rcomb9_k2 k2 used in recombiner type 9. Default = 16.7.
GASFLOW-MPI has a Xenon gas component, which when activated allows decay energy to be added to the gas mixture in the containment atmosphere as a function of time and space according to the convective behavior of the Xenon gas component.
Input Array. In the NAMELIST input block XPUT, the following variables have been added:
xecoef Xenon decay coefficient [ergs/(g·m·s)] (Default = 0.0)
xet0 Xenon decay time shift [s] (Default = 0.0)
xetau0 Xenon decay time constant [s] (Default = 1.0)
xepowert Power (-1<=xepowert<=1) for the (t-xet0) term in Equ. below. Note that the absolute value is actually used in the evaluation of the equation. The reason for this is discussed below. (Default = 0.0).
>0, for specific decay energy rate input, ergs/(g·s).
<0, for total decay energy rate input, ergs/s. For this condition the input function becomes:
It is convenient to allow tabular input to define the decay heating rate, either as a specific energy rate or as a total energy rate. We have added this capability by using xetau0 as an input flag in the following fashion:
xetau0
> 0, for function input.
< 0, for tabular input.
When xetau0 < 0, the tabular input is accomplished through the xecoef input variable.
Example A
Tabular total decay energy rate, ergs/s, input:
The NAMELIST XPUT input stream is given here:
xetau0 = -5.68828e+03, ; Tabular
xepowert = -0.267, ; (ergs/s) input
xecoef ; VALUE TIME
= 0000.000e+10, 0.0000e+03,
0170.650e+10, 0.2010e+03,
1685.150e+10, 1.5210e+03,
2062.000e+10, 2.0010e+03,
1239.500e+10, 3.8010e+03,
0533.490e+10, 1.8200e+04,
Example B
Tabular specific decay energy rate, ergs/(g·s), input:
Consider the time dependent Xenon mass in the system as shown in the Figure below
Dividing the tabular input of Example A by this Xenon mass data yields the specific decay energy rate results, and presented as the NAMELIST XPUT input stream given here:
xetau0 = -5.68828e+03, ; Tabular
xepowert = +0.267, ; (ergs/s) input
xecoef ; VALUE TIME
= 0000.000e+04, 0.0000e+03,
2171.000e+04, 0.2010e+03,
2832.000e+04, 1.5210e+03,
2635.000e+04, 2.0010e+03,
1584.000e+04, 3.8010e+03,
0682.000e+04, 1.8200e+04,
Example C
Tabular total decay energy rate, ergs/s, input:
The NAMELIST XPUT input stream is give here:
xetau0 = +4.64684e+03, ; Defined function
xepowert = -0.431, ; (ergs/s) input
xecoef = 120.214e+10, ; Coefficient
xet0 = 0.0e+00, ; Time shift
To activate the generalized fan options, one must place fans within the 3-dimensional mesh by using the fandef statement to define the n-th fan within the NAMELIST XPUT:
fandef( 1,n) Beginning i mesh index (cell face number).
fandef( 2,n) Ending i mesh index (cell face number).
fandef( 3,n) Beginning j mesh index (cell face number).
fandef( 4,n) Ending j mesh index (cell face number).
fandef( 5,n) Beginning k mesh index (cell face number).
fandef( 6,n) Ending k mesh index (cell face number).
fandef( 7,n) Block number (must be 1 for GASFLOW-MPI).
fandef( 8,n) ITFAN table number for fan head versus volumetric flow (see fantb below). Maximum number allowed is 20.
fandef( 9,n) ITSPD table number for fan speed versus time (see fspdtb below) Maximum number allowed is 20.
fandef( 10,n) SPDR, rated fan blower speed, revolutions/s.
fandef( 11,n) QR, rated fan volumetric flow rate, cm3/s.
>0, fan is directed in the positive coordinate direction.
< 0, fan is directed in the negative coordinate direction.
fandef( 12,n) HR, rated fan head
>0, fan table (fantb) is in dynes/cm2.
< 0, fan table (fantb) is in cm H2O.
fandef(13,n) SPD, fan speed at time = 0.0, revolutions/s.
There is a current limitation of 600 fandef statements.
The fan characteristics table, which is usually provided by the fan manufacture, is a table of paired values in the NAMELIST XPUT for the array fantb(1:2,maxp,maxtb), where maxp is limited to 30 paired values and maxtb is limited to 20 tables. The input is defined:
fantb(1,i,j) Time of table pair i for table number j, s.
fantb(2,i,j) Volumetric flow of table pair i for table number j, cm^3/s.
The fan performance table, which is problem dependent, is provided to GASFLOW with a table of paired values in the NAMELIST XPUT for the array fspdtb(1:2,maxp,maxtb), where again maxp is limited to 30 paired values and maxtb is limited to 20 tables. The input is defined:
fspdtb(1,i,j) Time of table pair i for table number j, s.
fspdtb(2,i,j) Fan speed of table pair i for table j, revolutions/s.
In the following example we specify two fans to fully operate for the first 40 seconds with then the first of the fans to fail and not be active beyond 40 seconds.
$xput
…
fandef(1:13,1) = 25, 25, 43, 46, 5, 8, 1, 1, 2, 2603.0, 2.469e+05, 4000.0, 2603.0,
fandef(1:13,2) = 25, 25, 43, 46, 10, 13, 1, 1, 1, 2603.0, 2.469e+05, 4000.0, 2603.0,
fantb(1:2, 1,1) = 0.000e+00, 5.8300e+03,
fantb(1:2, 2,1) = 3.740e+04, 5.5000e+03,
fantb(1:2, 3,1) = 7.516e+04, 5.3000e+03,
fantb(1:2, 4,1) = 1.117e+05, 5.1000e+03,
fantb(1:2, 5,1) = 1.473e+05, 4.9000e+03,
fantb(1:2, 6,1) = 2.012e+05, 4.5000e+03,
fantb(1:2, 7,1) = 2.469e+05, 4.0000e+03,
fantb(1:2, 8,1) = 2.818e+05, 3.5000e+03,
fantb(1:2, 9,1) = 3.114e+05, 2.9990e+03,
fantb(1:2, 10,1) = 3.657e+05, 1.9990e+03,
fantb(1:2, 11,1) = 3.910e+05, 1.4990e+03,
fantb(1:2, 12,1) = 4.131e+05, 0.9990e+03,
fantb(1:2, 13,1) = 4.349e+05, 0.4989e+03,
fantb(1:2, 14,1) = 4.551e+05, 0.000,
fspdtb(1:2,1,1) = 0.0e+0, 2603.0,
fspdtb(1:2,2,1) = 4.0e+1, 2603.0,
fspdtb(1:2,3,1) = 1.0e+5, 2603.0,
fspdtb(1:2,1,2) = 0.0e+0, 2603.0,
fspdtb(1:2,2,2) = 4.0e+1, 2603.0,
fspdtb(1:2,3,2) = 4.1e+1, 0000.0,
fspdtb(1:2,4,2) = 4.1e+9, 0000.0,
...
$end
We have generalized and expanded the existing energy source input, esdef, to a maximum of 500 NAMELIST XPUT statements and changed the actual input energy units of power, i.e., ergs/s where 1 Watt = 10^7ergs/s. This NAMELIST XPUT statements are now defined as follows:
esdef( 1,*) Beginning i mesh index (cell face number).
esdef( 2,*) Ending i mesh index (cell face number).
esdef( 3,*) Beginning j mesh index (cell face number).
esdef( 4,*) Ending j mesh index (cell face number).
esdef( 5,*) Beginning k mesh index (cell face number).
esdef( 6,*) Ending k mesh index (cell face number).
esdef( 7,*) Block number (must be 1 for GASFLOW-MPI).
esdef( 8,*) Total power in the computational volume defined by the first 7 elements of this array, ergs/s.
esdef( 9,*) Time (s) at which "esdef" begins.
esdef(10,*) Time (s) at which "esdef" ends.
In the way of an example, we add statements to the NAMELIST XPUT of the above example (example from Section 3.10.8), to demonstrate how one can easily model fans in, for example CPU Boxes with heat sources:
fandef(1:13, 3) = 19, 19, 09, 10, 04, 06, 1, 2, 3, 3300.0, +15.00e+03, 0000.0, 3300.0, ;Box edit fans
fandef(1:13, 4) = 09, 09, 09, 10, 07, 08, 1, 2, 3, 3300.0, +15.00e+03, 0000.0, 3300.0, ; Box internal fans
fandef(1:13, 5) = 09, 09, 09, 10, 09, 10, 1, 2, 3, 3300.0, +15.00e+03, 0000.0, 3300.0, ; Box internal fans
fantb(1:2, 1,2) = 0.000e+00, 6.0000e+02,
fantb(1:2, 2,2) = 2.500e+03, 5.2000e+02,
fantb(1:2, 3,2) = 5.000e+03, 3.7500e+02,
fantb(1:2, 4,2) = 7.500e+03, 2.8000e+02,
fantb(1:2, 5,2) = 10.000e+03, 2.4000e+02,
fantb(1:2, 6,2) = 12.500e+03, 1.8000e+02,
fantb(1:2, 7,2) = 15.000e+03, 0.0000e+00,
fspdtb(1:2,1,3) = 000.0e+0, 3300.0,
fspdtb(1:2,2,3) = 004.0e+1, 3300.0,
fspdtb(1:2,3,3) = 001.0e+5, 3300.0,
and the power deposited in the gas from each CPU (120 W each), the CPU power supply (10 W), and the motherboard power (40 W):
esdef(1:10,1) = 10, 11, 09, 10, 07, 08, 1, 1.20e+9, 5.0, 1.0e+10, ; right side CPU
esdef(1:10,2) = 10, 11, 09, 10, 09, 10, 1, 1.20e+9, 5.0, 1.0e+10, ; left side CPU
esdef(1:10,3) = 13, 17, 09, 10, 04, 06, 1, 1.00e+8, 5.0, 1.0e+10, ; CPU power supply
esdef(1:10,4) = 13, 18, 09, 10, 07, 12, 1, 3.00e+8, 5.0, 1.0e+10, ; Motherboard
Warning: The GASFLOW-MPI Spray model is experimental and under development as time permits; and therefore, should be used with extreme caution.
GASFLOW-II solves the classical two-phase Homogenous Equilibrium Model (HEM) with the assump-tion that the liquid droplets are dispersed in a gaseous medium. This means that there is thermal and mechanical equilibrium (equal temperatures and velocities, respectively) between the phases and that the volume fraction of the gaseous components is close to unity.
The first version of the GASFLOW-II spray model assumes that there is still mechanical equilibrium between the phases (equal velocities) but non-equilibrium temperatures. This requires the solution of a specific internal energy equation for the liquid droplet phase where the gas temperature is now computed by subtracting the liquid droplet energy from the total energy equation and inverting this relationship for the gas temperature and also inverting the specific internal energy function to obtain the liquid temperature. The pressure field is determined from the gaseous components only. Convective heat transfer and phase change mass transfer is modeled between liquid and vapor components to obtain the appropriate coupling phenomena.
Droplet depletion is also modeled to provide the effect of liquid sinks from the fluid field seen in droplet sedimentation or rainout and with a mechanistic impaction model.
ispray = 1, ; Spray model activated
The liquid component, ‘h2ol’, must be last in the mat input list: mat = 'air', 'h2o', 'h2ol',
There must be some liquid "seed" in all gasdef lists
Use only one cell to define a spray source in spraydef definitions.
The cell defined in spraydef for spray source shall not be surrounded by walls.
Droplet Advection Algorithm Options:
impsprayd = 0 ; Explicit droplet advection algorithm (default);
impsprayd = 1 ; Implicit droplet advection algorithm.
When using a gasdef input variable to define droplet temperatures and diameters, the user can use the spraygdef input array with the standard gasdef input array to specify droplet temperature and average diameter:
gasdef(1:18,1) = 1,'im1', 1,'jm1', 1, 'km1', 1, 1.015e+06, 300.00, 2, 0.0, 0.0, 'air', 0.999999, 'h2o', 0.0, 'h2ol', 0.000001,
spraygdef(1:3,1) = 300.0, 0.00001, 0.0,
The above Initial Condition example sets both the gas and liquid phases at 300 K with the droplet seed having an initial diameter equaling 0.1 microns.
gasdef(1:18,2) = 6, 7, 6, 7, 0, 1, 1, 1.015e+06, 600.0, 2, 0.0, 50.0, 'air', 0.0, 'h2o', 0.999, 'h2ol', 0.001,
spraygdef(1:3,2) = 300.0, 0.1, 0.0,
In this example, a typical boundary condition is active from 0.0 < time < 50.0 s, the gas is at 600 K while the liquid droplets are at 300 K with a diameter of 1 mm.
The temperature of the liquid droplets can be coupled with the sortam file when spraygdef(1,*) < 0, spraygdef(3,*) > 0 and sortami > 0.
gasdef(1:18,2) = 6, 7, 6, 7, 0, 1, 1, 1.015e+06, 600.0, 2, 0.0, 50.0, 'air', 0.0, 'h2o', 0.999, 'h2ol', 0.001,
spraygdef(1:3,2) = -5, 0.1, 1,
In this example, a typical boundary condition is active from 0.0 < time < 50.0 s. The gas temperature is 600 K. The transient temperature of the liquid droplets is read from the 5th column of the sortam file. The diameter of the liquid droplets is 1 mm.
When spraygdef(1,*) < 0 and spraygdef(3,*) <= 0, the temperature of the liquid droplets equals the surrounding gas temperature at the same cell.
There is a maximum of 500 spraygdef's possible, i.e., spraygdef(1:3,500) is the maximum allowed. The default values of spraygdef(1:3,500) are -9.0e+30.
The user can define a liquid spray source any place in the computational mesh (for example, to locate a spray head or nozzle) by using the spraydef input arrays:
spraydef(1:12,1) = 3, 4, 6, 7, 9, 10, 1, 2.0e+04, 300.0, 1.0e-02, 50.0, 99999.0,
spraydef(1:12,2) = 9, 10, 6, 7, 9, 10, 1, 2.0e+04, 300.0, 1.0e-02, 50.0, 99999.0,
spraydef(1:12,3) = 6, 7, 3, 4, 9, 10, 1, 2.0e+04, 300.0, 1.0e-02, 50.0, 99999.0,
spraydef(1:12,4) = 6, 7, 9, 10, 9, 10, 1, 2.0e+04, 300.0, 1.0e-02, 50.0, 99999.0,
The action of the sprays, coupling the above spraydef input with the traditional mbc and mvalue input arrays, is shown here with the following example:
mbc = 3, 4, 6, 7, 9, 9, 1, 1, 50.0, 99999.0,
9, 10, 6, 7, 9, 9, 1, 2, 50.0, 99999.0,
6, 7, 3, 4, 9, 9, 1, 3, 50.0, 99999.0,
6, 7, 9, 10, 9, 9, 1, 4, 50.0, 99999.0,
mvalue = -2.000e+04, -2.000e+04,
-2.000e+04, -2.000e+04,
Note: negative mvalue is required to obtain a negative velocity when the water droplets flow downwards. For instance, if gz = -980 cm/s2, in order to obtain the correct velocity at the spray nozzle a negative mvalue is needed. The magnitude of mvalue must be consistent with spraydef(8,*). spraydef(8,*) is needed to calculate the mass balance, and it is a positive number. mvalue is used to calculate the injection velocity at the spray nozzle, and it could be positive or negative.
There is a maximum of 100 spraydef's possible, i.e., spraydef(1:12,100) is the maximum allowed in the current version.
This model can be activated with the input variable in the $RHEAT NAMELIST group iliq.
= 1, use the correlation ;
