diff --git a/Manuals/FDS_Technical_Reference_Guide/Appendices.tex b/Manuals/FDS_Technical_Reference_Guide/Appendices.tex index b9a19322077..2244c9f8d61 100644 --- a/Manuals/FDS_Technical_Reference_Guide/Appendices.tex +++ b/Manuals/FDS_Technical_Reference_Guide/Appendices.tex @@ -424,11 +424,11 @@ \chapter{Gas Phase Absorption Coefficients} By default in FDS, the pathlength, $L$, is 10~cm. It can also be specified by the user. If $T=T_{\rm rad}$ the intensity does not depend on $\kappa_{n,i,{\rm e}}$, and the value $\kappa_{n,i,{\rm e}}(T_{\rm rad})$ is therefore interpolated from the neighboring temperatures. -In cases with only one band ($N$=1), the smaller of the two absorption coefficients is used: +In cases of non-zero path length, the smaller of the two absorption coefficients is used: \be \kappa_{n,i}=\min \Big( \kappa_{n,i,{\rm P}}(Y_i,T),\kappa_{n,i,{\rm e}}(Y_i,T,L) \Big) \ee -If $N>1$ or $L=0$, $\kappa_{n,i}=\kappa_{n,i,{\rm P}}$. Note that the spectral data within RADCAL are used whenever the gas mixture contains water vapor, fuel or combustion products, regardless of the number of radiation bands $N$. +If $L=0$, $\kappa_{n,i}=\kappa_{n,i,{\rm P}}$. Note that the spectral data within RADCAL are used whenever the gas mixture contains water vapor, fuel or combustion products, regardless of the number of radiation bands $N$. \textbf{Note on wavenumber, wavelength, and frequency}: some confusion might arise when dealing with the various quantities describing the wave nature of radiation. These quantities are wavenumber $\om$, wavelength $\la$, and frequency denoted here $\nu$. Most users may be familiar with the frequency $\nu$, in units of \textit{hertz}, ${\rm Hz}$, representing the number of cycles per second. While this unit is preferred for radiation waves of low energy such as radio waves, wavelength and wavenumber are preferred for waves of higher energy. Wavenumber and wavelength are related to frequency through \cite{Penner:1959} \be diff --git a/Manuals/FDS_User_Guide/FDS_User_Guide.tex b/Manuals/FDS_User_Guide/FDS_User_Guide.tex index 383006ab3df..5b390ef7c52 100644 --- a/Manuals/FDS_User_Guide/FDS_User_Guide.tex +++ b/Manuals/FDS_User_Guide/FDS_User_Guide.tex @@ -306,9 +306,10 @@ \subsubsection{Devices and Output} \chapter{Getting Started} \label{info:gettingstarted} -FDS is a computer program that solves equations that describe the evolution of fire. It is a Fortran program that reads input parameters from a text file, computes a numerical solution to the governing equations, and writes user-specified output data to files. Smokeview is a companion program that reads FDS output files and produces animations on the computer screen. Smokeview has a simple menu-driven interface. FDS does not. However, there are various third-party programs that have been developed to generate the text file containing the input parameters needed by FDS. +FDS is a Fortran program that solves the equations governing fire or other thermally-driven, low Mach number flows. FDS reads input parameters from a text file, computes a numerical solution to the governing equations, and writes user-specified output data to files. Smokeview\footnote{A separate document~\cite{Smokeview_Users_Guide} describes how to use Smokeview.} is a companion program that reads FDS output files and produces animations on the computer screen. Smokeview has a simple menu-driven interface. FDS does not. However, there are various third-party programs that have been developed to generate the text file containing the input parameters needed by FDS. + +Some third-party applications, like \href{https://www.thunderheadeng.com/pyrosim/}{PyroSim from Thunderhead Engineering}, have graphical interfaces for specifying input parameters and for viewing results. FDS is the core solver of such a package, but Smokeview would be replaced by a program that is compatible with the geometry engine used to build the model. If you are using a package like PyroSim, FDS would be installed as part of the larger package and you can skip the next section which is intended for those who want to download and install FDS and Smokeview only. -This guide describes how to obtain FDS and Smokeview and how to use FDS. A separate document~\cite{Smokeview_Users_Guide} describes how to use Smokeview. \section{How to Acquire FDS and Smokeview} \label{info:acquire} @@ -322,7 +323,6 @@ \section{How to Acquire FDS and Smokeview} \noindent The typical FDS/Smokeview distribution consists of an installation package or compressed archive, which is available for MS Windows, macOS, and Linux. -Old versions can be kept by copying the old version's installation directly to another location so that it is not overwritten when installing a new version. \section{Computer Hardware Requirements} @@ -347,6 +347,8 @@ \section{Installation Testing} If you are running FDS under a quality assurance plan that requires installation testing, a test procedure is provided in Appendix B of the FDS Verification Guide~\cite{FDS_Verification_Guide}. This guide can be obtained from the FDS-SMV website. +If you are a first-time user or you simply want to ensure that FDS and Smokeview have been installed properly on your computer, run\footnote{On a Windows computer, FDS is typically installed in a directory that is read-only for most users. Copy the sample input file and move it to a directory for which you have write-privilege.} the case called \ct{Examples/Fire/simple_test.fds} following the instructions in the next chapter. This case requires only a few tens of minutes to complete and allows you to visualize the results of a very simple fire scenario. As you learn about new features in the chapters to follow, this case can serve as a useful testbed. The text file \ct{simple_test.fds} contains parameters organized by {\em namelist} groups that are listed alphabetically in Chapter~\ref{info:namelists}. All FDS input parameters are listed in these tables, and each has a hyperlink to a section in this guide that explains its use. + \chapter{Running FDS} @@ -10545,7 +10547,7 @@ \subsection{Thermocouples} \frac{D_{\rm TC}}{6} \, \rho_{\rm TC} \, c_{\rm TC} \, \frac{\d T_{\rm TC}}{\d t} = \epsilon_{\rm TC} \, (U/4 - \sigma T_{\rm TC}^4) + h(T_{\rm g} - T_{\rm TC}) \label{TC} \ee -where $\epsilon_{\rm TC}$ is the emissivity of the thermocouple, $U$ is the integrated radiative intensity, $T_{\rm g}$ is the true gas temperature, and $h$ is the heat transfer coefficient to a small sphere, $h=k \, \NU / D_{\rm TC}$. The bead \ct{DIAMETER}, \ct{EMISSIVITY}, \ct{DENSITY}, and \ct{SPECIFIC_HEAT} are given on the associated \ct{PROP} line. To over-ride the calculated value of the heat transfer coefficient, set \ct{HEAT_TRANSFER_COEFFICIENT} on the \ct{PROP} line (\si{W/(m.K)}). The default value for the bead diameter is 0.001~m. The default emissivity is 0.85. The default values for the bead density and specific heat are that of nickel; 8908~kg/m$^3$ and 0.44~kJ/kg/K, respectively. See the discussion on heat transfer to a water droplet in the Technical Reference Guide for details of the convective heat transfer to a small sphere. +where $\epsilon_{\rm TC}$ is the emissivity of the thermocouple, $U$ is the integrated radiative intensity, $T_{\rm g}$ is the true gas temperature, and $h$ is the heat transfer coefficient to a small sphere, $h=k \, \NU / D_{\rm TC}$. The bead \ct{DIAMETER}, \ct{EMISSIVITY}, \ct{DENSITY}, and \ct{SPECIFIC_HEAT} or \ct{SPECIFIC_HEAT_RAMP} are given on the associated \ct{PROP} line. To over-ride the calculated value of the heat transfer coefficient, set \ct{HEAT_TRANSFER_COEFFICIENT} on the \ct{PROP} line (\si{W/(m.K)}). The default value for the bead diameter is 0.001~m. The default emissivity is 0.85. The default values for the bead density and specific heat are that for a Type-K thermocouple; 8700~kg/m$^3$ and a linear ramp from 0.4515~kJ/kg/K at 20\,$^\circ$C to 0.6010~kJ/kg/K at 1200\,$^\circ$C, respectively. See the discussion on heat transfer to a water droplet in the Technical Reference Guide for details of the convective heat transfer to a small sphere. The above discussion is appropriate for a so-called ``bare bead'' thermocouple, but often thermocouples are shielded or sheathed in various ways to mitigate the effect of thermal radiation. In such cases, there is no obvious bead diameter and there may be multiple metals and air gaps in the construction. Usually, the manufacturer provides a time constant, which is defined as the time required for the sensor to respond to 63.2~\% of its total output signal when suddenly plunged into a warm air stream flowing at 20~m/s. The analysis and testing is typically done at relatively low temperature, in which case the radiation term in Eq.~(\ref{TC}) can be neglected and the time constant, $\tau$, can be defined in terms of effective thermal properties and an effective diameter: \be @@ -13188,7 +13190,7 @@ \section{\texorpdfstring{{\tt PROP}}{PROP} (Device Properties)} \ct{CALIBRATION_CONSTANT} & Real & Section~\ref{info:bidir_probe} & & 0.93 \\ \hline \ct{CHARACTERISTIC_VELOCITY} & Real & Section~\ref{info:pressure_coefficient} & m/s & 1. \\ \hline \ct{C_FACTOR} & Real & Section~\ref{info:sprinklers} & (m/s)$^{1/2}$ & 0. \\ \hline -\ct{DENSITY} & Real & Section~\ref{info:THERMOCOUPLE} & kg/m$^3$ & 8908. \\ \hline +\ct{DENSITY} & Real & Section~\ref{info:THERMOCOUPLE} & kg/m$^3$ & 8700. \\ \hline \ct{DIAMETER} & Real & Section~\ref{info:THERMOCOUPLE} & m & 0.001 \\ \hline \ct{EMISSIVITY} & Real & Section~\ref{info:THERMOCOUPLE} & & 0.85 \\ \hline \ct{FED_ACTIVITY} & Integer & Section~\ref{info:FED} & & 2 \\ \hline @@ -13232,7 +13234,8 @@ \section{\texorpdfstring{{\tt PROP}}{PROP} (Device Properties)} \ct{SMOKEVIEW_PARAMETERS(:)} & Char.~Array & Section~\ref{info:SMOKEVIEW_PARAMETERS} & & \\ \hline \ct{SPARK} & Logical & Section~\ref{info:ignition} & & \ct{F} \\ \hline \ct{SPEC_ID} & Character & Section~\ref{info:alternative_smoke} & & \\ \hline -\ct{SPECIFIC_HEAT} & Real & Section~\ref{info:THERMOCOUPLE} & $\si{kJ/(kg.K)}$ & 0.44 \\ \hline +\ct{SPECIFIC_HEAT} & Real & Section~\ref{info:THERMOCOUPLE} & $\si{kJ/(kg.K)}$ & \\ \hline +\ct{SPECIFIC_HEAT_RAMP} & Character & Section~\ref{info:THERMOCOUPLE} & & \\ \hline \ct{SPRAY_ANGLE(2,2)} & Real & Section~\ref{info:sprinklers} & degrees & 60.,75. \\ \hline \ct{SPRAY_PATTERN_BETA} & Real & Section~\ref{info:sprinklers} & degrees & 5. \\ \hline \ct{SPRAY_PATTERN_MU} & Real & Section~\ref{info:sprinklers} & degrees & 0. \\ \hline diff --git a/Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex b/Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex index e1d39d90dca..5c6eece5317 100644 --- a/Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex +++ b/Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex @@ -4423,6 +4423,15 @@ \subsection{Thermophoretic Settling and Deposition of Aerosols\\(\texorpdfstring This verification test consists of two test cases. The second case, \ct{aerosol\_thermophoretic\_deposition\_2}, reverses the temperature gradient. The case consists of a box 1~cm on a side with adiabatic, free-slip side walls and a 100 K temperature gradient over the height of the box. The box is filled with two gas species each having a molecular weight of 28.8~g/mol, a viscosity of $2\times10^{-5}$~\si{kg/(m.s)}, a thermal conductivity of 0.025~\si{W/(m.K)}, and specific heat of 1~\si{kJ/(kg.K)}, and zero diffusivity. One of the gas species is defined as an aerosol with a diameter of 1~$\mu$m, a solid phase density of 2000~\si{kg/m^3}, and a solid phase conductivity of 1~\si{W/(m.K)}. The initial mass fraction of the aerosol is $1\times10^{-5}$. The gas temperature is initialized to its steady-state temperature gradient. \ct{STRATIFICATION}, \ct{NOISE}, and all aerosol behaviors except for \ct{THERMOPHORETIC\_SETTLING} and \ct{THERMOPHORETIC\_DEPOSITION} are turned off. Thermophoretic settling rates are weakly dependent on the gas density. Since there is a temperature gradient, the settling rates are not uniform over the height of the box. Unlike the gravitational settling case, this means over long enough time periods the overall settling rate is not linear in time; however, for a short time period a near linear settling rate is expected and can be determined analytically. +Using the equation for the thermophoretic velocity, $u_{\rm th}$, the thermophoretic velocities for the gas cells adjacent to the hot and cold walls are respectively 2.25E-4~\si{m/s} and 2.15E-4~\si{m/s}. Adjacent to the cold wall, $u_{\rm th}$ is the rate at which the aerosol deposits onto the wall. Over the simulation time of 2~s, the soot will move 4.29E-4~m which is 42.9~\% of the cell size. That fraction of the aerosol in the cell will deposit which is 5.05E-9~\si{kg/m^2}. Adjacent to the hot wall, $u_{\rm th}$ is the rate at which the aerosol leaves the cell. Over the simulation time of 2~s, the soot will move 4.49E-4~m which is 44.9~\% of the cell size. One minus that fraction of the original aerosol density in the cell will remain which is 4.98E-6~\si{kg/m^3}. + +The thermophorertic velocity is given by +\be +u_{\rm th} = \frac{K_{\rm th} \nu}{T_{\rm g}} \; \frac{\d T}{\d x} +\ee +where $K_{\rm th}$ is a coefficient dependent on the particle size and solid conductivity, $\nu$ is the kinematic viscosity, $T_{\rm g}$ is the cell gas temperature, and $\frac{\d T}{\d x}$ is the temperature gradient. + + \begin{figure}[ht] \noindent \begin{tabular*}{\textwidth}{l@{\extracolsep{\fill}}r} diff --git a/Source/chem.f90 b/Source/chem.f90 index e8d36c9f618..424e8995d89 100644 --- a/Source/chem.f90 +++ b/Source/chem.f90 @@ -783,7 +783,8 @@ END FUNCTION DDTMP_TROE !> \brief cvode interface for ODE integrator. Call sundials cvode in serial mode. !> \param CC species concentration (kmol/m3) array !> \param ZZ_0 initial species mass fraction array (of the unbuned zone),needed for mixing+chem -!> \param TMP_IN is the temperature of the cell (unburned zone) +!> \param TMP_IN is the temperature of the mixed zone (For DNS it is cell temp) +!> \param TMP_UNMIX is the unmix zone temperature !> \param PR_IN is the pressure !> \param ZETA0 is the initial unmixed fraction !> \param TAU_MIX is Mixing timescale @@ -792,13 +793,13 @@ END FUNCTION DDTMP_TROE !> \param TEND is the end time in seconds !> \param RTOL is the relative error for all the species (REAL_EB) !> \param ATOL is the absolute error tolerance array for the species (REAL_EB) -!> \param TMP_OUT reactor calculated temperature at the end +!> \param TMP_OUT reactor calculated mixing zone temperature at the end !> \param CHEM_TIME Chemical time scale !> \param WRITE_SUBSTEPS Whether to write cvode substeps. Only write for first cfd step. !> \param CVODE_CALL_OPTION 1:CV_NORMAL, 2=CV_ONE_STEP !> \details This is the interface subroutine to the other modules. -SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,TEND, RTOL, ATOL, & +SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, TMP_UNMIX, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR,TEND, RTOL, ATOL, & TMP_OUT, CHEM_TIME, WRITE_SUBSTEPS, CVODE_CALL_OPTION) USE PHYSICAL_FUNCTIONS, ONLY : MOLAR_CONC_TO_MASS_FRAC, CALC_EQUIV_RATIO, GET_ENTHALPY, GET_MOLECULAR_WEIGHT USE COMP_FUNCTIONS, ONLY: GET_FILE_NUMBER @@ -813,7 +814,7 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR, USE FSUNDIALS_NVECTOR_MOD ! FORTRAN INTERFACE TO GENERIC N_VECTOR REAL(EB), INTENT(INOUT) :: CC(N_TRACKED_SPECIES) -REAL(EB), INTENT(IN) :: ZZ_0(N_TRACKED_SPECIES),TMP_IN,PR_IN,ZETA0,TAU_MIX,CELL_MASS,TCUR,TEND +REAL(EB), INTENT(IN) :: ZZ_0(N_TRACKED_SPECIES),TMP_IN,TMP_UNMIX,PR_IN,ZETA0,TAU_MIX,CELL_MASS,TCUR,TEND REAL(EB), INTENT(IN) :: ATOL(N_TRACKED_SPECIES) REAL(EB), INTENT(IN) :: RTOL REAL(EB), INTENT(OUT) :: TMP_OUT,CHEM_TIME @@ -841,7 +842,7 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR, TYPE(C_PTR) :: USERDATAPTR ! USER DATA CONTAINS MIXING INFORMATION REAL(EB) :: ZZ(N_TRACKED_SPECIES), EQUIV, H_IN -INTEGER :: CVODE_TASK, NS, NTRY, MAXTRY, SUBSTEP_COUNT +INTEGER :: CVODE_TASK, NS, NTRY, MAXTRY, SUBSTEP_COUNT, MAXTRY_FAC REAL(EB) :: H_G TYPE(USERDATA), TARGET :: USER_DATA LOGICAL :: ONLY_FIRST_STEP=.TRUE. ! Needed in CV_ONE_STEP @@ -948,7 +949,7 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR, USER_DATA%CELL_MASS = CELL_MASS ALLOCATE(USER_DATA%ZZ_0(N_TRACKED_SPECIES)) USER_DATA%ZZ_0 = ZZ_0 -CALL GET_ENTHALPY(ZZ_0,H_IN,TMP_IN) +CALL GET_ENTHALPY(ZZ_0,H_IN,TMP_UNMIX) USER_DATA%H_IN = H_IN USERDATAPTR = C_LOC(USER_DATA) IERR_C = FCVODESETUSERDATA(CVODE_MEM, USERDATAPTR) @@ -987,12 +988,14 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR, ENDIF END DO - ENDIF IF (IERR_C /= 0) THEN - MAXTRY = CVODE_MAX_TRY ! If all internal substeps are taken try two more times. This will allow larger CFD timestep. + ! Make MAXTRY at least CVODE_MAX_TRY or for larger timestep (>1E-3) scale it proportionaly. + MAXTRY_FAC = CEILING((TEND-TCUR)/1.0E-3_EB) + MAXTRY_FAC = MIN(MAX(MAXTRY_FAC,1),50) + MAXTRY = CVODE_MAX_TRY*MAXTRY_FAC IF (IERR_C == CV_TOO_MUCH_WORK) THEN !CV_TOO_MUCH_WORK == all internal substeps are taken NTRY = 0 DO WHILE (NTRY < MAXTRY) @@ -1008,8 +1011,8 @@ SUBROUTINE CVODE_SERIAL(CC,ZZ_0, TMP_IN, PR_IN, ZETA0, TAU_MIX, CELL_MASS, TCUR, IF (IERR_C .NE. CV_SUCCESS) THEN IF (IERR_C == CV_TOO_MUCH_WORK) THEN - WRITE(LU_ERR,'(A, 2E18.8, A)')" WARN: CVODE took all internal substeps. CUR_CFD_TIME, DT=", CUR_CFD_TIME, (TEND-TCUR), & - ". If the warning persists, reduce the timestep." + WRITE(LU_ERR,'(A, 2E18.8, I8, A)')" WARN: CVODE took all internal substeps. CUR_CFD_TIME, DT, MAXTRY=", CUR_CFD_TIME, & + (TEND-TCUR), MAXTRY, ". If the warning persists, reduce the timestep." ELSE WRITE(LU_ERR,'(A, I4, A, 2E18.8, A)')" WARN: CVODE didn't finish ODE solution with message code:", IERR_C, & " and CUR_CFD_TIME, DT=", CUR_CFD_TIME, (TEND-TCUR), ". If the warning persists, reduce the timestep." diff --git a/Source/cons.f90 b/Source/cons.f90 index cadac973bc0..e8bb81da934 100644 --- a/Source/cons.f90 +++ b/Source/cons.f90 @@ -291,7 +291,6 @@ MODULE GLOBAL_CONSTANTS LOGICAL :: WROTE_SMOKE3D=.FALSE. !< Flag indicating if a Smoke3D file was written during primary out LOGICAL :: WROTE_BNDF=.FALSE. !< Flag indicating if a BNDF file was written during primary out LOGICAL :: WROTE_PART=.FALSE. !< Flag indicating if a PART file was written during primary out -LOGICAL :: FLUX_LIMITER_SINGLE_COEF=.TRUE. !< Flag to base flux coefficients off of a single worst-case scalar gradient LOGICAL :: STORE_FIRE_ARRIVAL=.FALSE. !< Flag for tracking arrival of spreading fire front over a surface LOGICAL :: STORE_FIRE_RESIDENCE=.FALSE. !< Flag for tracking residence time of spreading fire front over a surface LOGICAL :: STORE_LS_SPREAD_RATE=.FALSE. !< Flag for outputting local level set spread rate magnitude @@ -942,14 +941,25 @@ END MODULE RADCONS MODULE CHEMCONS USE PRECISION_PARAMETERS +!> \brief Parameters associated with IGNITION_ZONES +TYPE IGNITION_ZONE_TYPE + REAL(EB) :: X1 !< Lower x bound of Ignition Zone + REAL(EB) :: X2 !< Upper x bound of Ignition Zone + REAL(EB) :: Y1 !< Lower y bound of Ignition Zone + REAL(EB) :: Y2 !< Upper y bound of Ignition Zone + REAL(EB) :: Z1 !< Lower z bound of Ignition Zone + REAL(EB) :: Z2 !< Upper z bound of Ignition Zone + INTEGER :: DEVC_INDEX=0 !< Index of device controlling the status of the zone + CHARACTER(LABEL_LENGTH) :: DEVC_ID='null' !< Name of device controlling the status of the zone +END TYPE IGNITION_ZONE_TYPE + INTEGER, ALLOCATABLE, DIMENSION(:) :: YP2ZZ REAL(EB) :: ODE_MIN_ATOL= -1._EB LOGICAL :: EQUIV_RATIO_CHECK = .FALSE. -REAL(EB) :: MIN_EQUIV_RATIO=0.0_EB +REAL(EB) :: MIN_EQUIV_RATIO=0.1_EB REAL(EB) :: MAX_EQUIV_RATIO=20.0_EB LOGICAL :: DO_CHEM_LOAD_BALANCE = .FALSE. INTEGER :: MAX_CVODE_SUBSTEPS=100000 -REAL(EB) :: MAX_CHEM_TIME=1.E-6_EB INTEGER :: CVODE_MAX_TRY=4 LOGICAL :: IS_EXPONENT_LT_1 = .FALSE. @@ -958,7 +968,18 @@ MODULE CHEMCONS REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: CVODE_SUBSTEP_DATA INTEGER :: TOTAL_SUBSTEPS_TAKEN -! FOR THICKENED FLAME MODEL -REAL(EB) :: FLAME_THICK_FACTOR=1.0_EB +! Adiabatic flame temperature calculation +CHARACTER(LABEL_LENGTH) :: FUEL_ID_FOR_AFT='null' +INTEGER :: I_FUEL,I_CO2,I_H2O,I_O2,I_N2 ! Store the index of the species in the ZZ array. +LOGICAL :: USE_MIXED_ZN_AFT_TMP = .FALSE. + +! Mixing +REAL(EB) :: ZETA_ARTIFICAL_MIN_LIMIT=0.99_EB +REAL(EB) :: ZETA_ARTIFICAL_MAX_LIMIT=0.9999_EB +REAL(EB) :: ZETA_FIRST_STEP_DIV=10._EB + +! IGNITION ZONES (mainly for premixed flame) +INTEGER :: N_IGNITION_ZONES = 0 +TYPE(IGNITION_ZONE_TYPE), DIMENSION(MAX_IGNITION_ZONES) :: IGNITION_ZONES !< Coordinates of ignition zones END MODULE CHEMCONS diff --git a/Source/data.f90 b/Source/data.f90 index b46fc0e5b66..c5813108a3b 100644 --- a/Source/data.f90 +++ b/Source/data.f90 @@ -677,6 +677,10 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES OUTPUT_QUANTITY(144)%SHORT_NAME = 'Y_ELEM' OUTPUT_QUANTITY(144)%ELEM_ID_REQUIRED = .TRUE. +OUTPUT_QUANTITY(145)%NAME='EQUIVALENCE RATIO' +OUTPUT_QUANTITY(145)%UNITS='' +OUTPUT_QUANTITY(145)%SHORT_NAME = 'equivRatio' + OUTPUT_QUANTITY(150)%NAME = 'SUM LUMPED VOLUME FRACTIONS' OUTPUT_QUANTITY(150)%UNITS = '' OUTPUT_QUANTITY(150)%SHORT_NAME = 'sum X' @@ -2148,7 +2152,7 @@ SUBROUTINE COLOR2RGB(RGB,COLOR) CASE ('NAVY');RGB = (/0,0,128/) CASE ('OLD LACE');RGB = (/253,245,230/) CASE ('OLIVE');RGB = (/128,128,0/) -CASE ('OLIVE DRAB');RGB = (/192,255,62/) +CASE ('OLIVE DRAB');RGB = (/107,142,35/) CASE ('OLIVE DRAB 1');RGB = (/179,238,58/) CASE ('OLIVE DRAB 2');RGB = (/154,205,50/) CASE ('OLIVE DRAB 3');RGB = (/105,139,34/) diff --git a/Source/devc.f90 b/Source/devc.f90 index e7de967e69c..58bc47dcdf6 100644 --- a/Source/devc.f90 +++ b/Source/devc.f90 @@ -10,17 +10,18 @@ MODULE DEVICE_VARIABLES !> \brief Derived type storing inputs for the PROP namelist group TYPE PROPERTY_TYPE - REAL(EB) :: DENSITY,DIAMETER,EMISSIVITY,HEAT_TRANSFER_COEFFICIENT,SPECIFIC_HEAT,RTI,TIME_CONSTANT, & + REAL(EB) :: DENSITY,DIAMETER,EMISSIVITY,HEAT_TRANSFER_COEFFICIENT,RTI,TIME_CONSTANT, & ACTIVATION_TEMPERATURE,ACTIVATION_OBSCURATION, CALIBRATION_CONSTANT,& ALPHA_E,ALPHA_C,BETA_E,BETA_C,CHARACTERISTIC_VELOCITY,PARTICLE_VELOCITY,MASS_FLOW_RATE,FLOW_RATE,FLOW_TAU, & GAUGE_EMISSIVITY,GAUGE_TEMPERATURE,INITIAL_TEMPERATURE,K_FACTOR,C_FACTOR,OPERATING_PRESSURE,OFFSET,& SPRAY_ANGLE(2,2),P0=0._EB,PX(3)=0._EB,PXX(3,3)=0._EB,VIEW_ANGLE,PROBE_DIAMETER INTEGER :: PDPA_M=0,PDPA_N=0,N_SMOKEVIEW_PARAMETERS=0,N_SMOKEVIEW_IDS=0,N_INSERT,I_VEL=0,PARTICLES_PER_SECOND LOGICAL :: PDPA_INTEGRATE=.TRUE.,PDPA_NORMALIZE=.TRUE.,HISTOGRAM_NORMALIZE=.TRUE.,HISTOGRAM=.FALSE., & - HISTOGRAM_CUMULATIVE=.FALSE.,SPARK=.FALSE.,TC=.TRUE. + HISTOGRAM_CUMULATIVE=.FALSE.,SPARK=.FALSE.,TC=.TRUE.,IGNITION_ZONE=.FALSE. REAL(EB) :: PDPA_START=0._EB,PDPA_END=1.E6_EB,PDPA_RADIUS=0.1_EB - REAL(EB), ALLOCATABLE, DIMENSION(:) :: TABLE_ROW, V_FACTOR - INTEGER :: PART_INDEX=-1,FLOW_RAMP_INDEX,SPRAY_PATTERN_INDEX,Z_INDEX=-999,Y_INDEX=-999,PRESSURE_RAMP_INDEX + REAL(EB), ALLOCATABLE, DIMENSION(:) :: TABLE_ROW, V_FACTOR,SPECIFIC_HEAT + INTEGER :: PART_INDEX=-1,FLOW_RAMP_INDEX,SPRAY_PATTERN_INDEX,Z_INDEX=-999,Y_INDEX=-999,PRESSURE_RAMP_INDEX,& + SPECIFIC_HEAT_RAMP_INDEX=-1 CHARACTER(LABEL_LENGTH) :: SMOKEVIEW_ID(SMOKEVIEW_OBJECTS_DIMENSION),PART_ID,ID,QUANTITY,TABLE_ID,SPEC_ID='null', & SMOKEVIEW_PARAMETERS(SMOKEVIEW_OBJECTS_DIMENSION)='null' REAL(EB), ALLOCATABLE, DIMENSION(:) :: SPRAY_LON_CDF,SPRAY_LON,SPRAY_LAT diff --git a/Source/divg.f90 b/Source/divg.f90 index 8133546dd7d..6c4d2dd524e 100644 --- a/Source/divg.f90 +++ b/Source/divg.f90 @@ -30,7 +30,6 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) USE COMPLEX_GEOMETRY, ONLY : CC_CGSC, CC_UNKZ, CC_SOLID, CC_CUTCFE USE CC_SCALARS, ONLY : ADD_CUTCELL_PSUM,ADD_LINKEDCELL_PSUM,SET_EXIMDIFFLX_3D,SET_EXIMRHOHSLIM_3D,& SET_EXIMRHOZZLIM_3D,CC_DIVERGENCE_PART_1,CC_VELOCITY_FLUX,CFACE_PREDICT_NORMAL_VELOCITY -USE CHEMCONS, ONLY : FLAME_THICK_FACTOR INTEGER, INTENT(IN) :: NM REAL(EB), INTENT(IN) :: T,DT @@ -150,10 +149,6 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) ENDIF ENDIF - IF (FLAME_THICK_FACTOR > 1._EB) THEN - RHO_D = RHO_D*FLAME_THICK_FACTOR - ENDIF - ! Manufactured solution IF (PERIODIC_TEST==7) RHO_D = DIFF_MMS @@ -591,7 +586,7 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) CONST_GAMMA_IF_1: IF (.NOT.CONSTANT_SPECIFIC_HEAT_RATIO) THEN - CALL ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) ! Compute u dot grad rho h_s + CALL ENTHALPY_ADVECTION_NEW(U_DOT_DEL_RHO_H_S) ! Compute u dot grad rho h_s DO K=1,KBAR DO J=1,JBAR @@ -640,11 +635,7 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) CONST_GAMMA_IF_2: IF (.NOT.CONSTANT_SPECIFIC_HEAT_RATIO) THEN - IF (FLUX_LIMITER_SINGLE_COEF) THEN - CALL SPECIES_ADVECTION_PART_1_NEW - ELSE - CALL SPECIES_ADVECTION_PART_1 ! Compute and store face values of (rho Z_n) - ENDIF + CALL SPECIES_ADVECTION_PART_1_NEW ! Compute and store face values of (rho Z_n) DO N=1,N_TRACKED_SPECIES @@ -794,14 +785,14 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM) END SUBROUTINE DIVERGENCE_PART_1 -SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) +SUBROUTINE ENTHALPY_ADVECTION_NEW(U_DOT_DEL_RHO_H_S) USE PHYSICAL_FUNCTIONS, ONLY: GET_SENSIBLE_ENTHALPY -USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE +USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_COEF,GET_SCALAR_FACE_VALUE_NEW REAL(EB), POINTER, DIMENSION(:,:,:) :: FX_H_S,FY_H_S,FZ_H_S,RHO_H_S_P,U_DOT_DEL_RHO_H_S REAL(EB) :: UN,UN_P,TMP_F_GAS,DU_P,DU_M,DV_P,DV_M,DW_P,DW_M,DU,H_S REAL(EB), ALLOCATABLE, DIMENSION(:) :: ZZ_GET -REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP +REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP,B_TEMP REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP INTEGER :: IC,I,J,K,IW TYPE(WALL_TYPE), POINTER :: WC @@ -814,6 +805,10 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) FZ_H_S=>WORK4 ; FZ_H_S = 0._EB U_DOT_DEL_RHO_H_S=>WORK6 ; U_DOT_DEL_RHO_H_S=0._EB +BFX = 1.E10_EB +BFY = 1.E10_EB +BFZ = 1.E10_EB + ! Compute and store rho*h_s !$OMP PARALLEL PRIVATE(ZZ_GET, H_S) @@ -834,9 +829,13 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) ! Compute scalar face values -CALL GET_SCALAR_FACE_VALUE(UU,RHO_H_S_P,FX_H_S,1,IBM1,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) -CALL GET_SCALAR_FACE_VALUE(VV,RHO_H_S_P,FY_H_S,1,IBAR,1,JBM1,1,KBAR,2,I_FLUX_LIMITER) -CALL GET_SCALAR_FACE_VALUE(WW,RHO_H_S_P,FZ_H_S,1,IBAR,1,JBAR,1,KBM1,3,I_FLUX_LIMITER) +CALL GET_SCALAR_FACE_COEF(UU,RHO_H_S_P,BFX,0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) +CALL GET_SCALAR_FACE_COEF(VV,RHO_H_S_P,BFY,1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) +CALL GET_SCALAR_FACE_COEF(WW,RHO_H_S_P,BFZ,1,IBAR,1,JBAR,0,KBAR,3,I_FLUX_LIMITER) + +CALL GET_SCALAR_FACE_VALUE_NEW(UU,RHO_H_S_P,FX_H_S,BFX,0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) +CALL GET_SCALAR_FACE_VALUE_NEW(VV,RHO_H_S_P,FY_H_S,BFY,1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) +CALL GET_SCALAR_FACE_VALUE_NEW(WW,RHO_H_S_P,FZ_H_S,BFZ,1,IBAR,1,JBAR,0,KBAR,3,I_FLUX_LIMITER) ALLOCATE(ZZ_GET(1:N_TRACKED_SPECIES)) @@ -876,7 +875,8 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) IF ((UU(BC%II+1,BC%JJ,BC%KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II+1,BC%JJ,BC%KK))%WALL_INDEX(+1)>0)) THEN Z_TEMP(0:2,1,1) = (/RHO_H_S_P(BC%II+1,BC%JJ,BC%KK),RHO_H_S_P(BC%II+1:BC%II+2,BC%JJ,BC%KK)/) U_TEMP(1,1,1) = UU(BC%II+1,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFX(BC%II+1,BC%JJ,BC%KK) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) FX_H_S(BC%II+1,BC%JJ,BC%KK) = F_TEMP(1,1,1) ENDIF CASE(-1) OFF_WALL_SELECT_1 @@ -886,35 +886,40 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) IF ((UU(BC%II-2,BC%JJ,BC%KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II-1,BC%JJ,BC%KK))%WALL_INDEX(-1)>0)) THEN Z_TEMP(1:3,1,1) = (/RHO_H_S_P(BC%II-2:BC%II-1,BC%JJ,BC%KK),RHO_H_S_P(BC%II-1,BC%JJ,BC%KK)/) U_TEMP(1,1,1) = UU(BC%II-2,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFX(BC%II-2,BC%JJ,BC%KK) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) FX_H_S(BC%II-2,BC%JJ,BC%KK) = F_TEMP(1,1,1) ENDIF CASE( 2) OFF_WALL_SELECT_1 IF ((VV(BC%II,BC%JJ+1,BC%KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ+1,BC%KK))%WALL_INDEX(+2)>0)) THEN Z_TEMP(1,0:2,1) = (/RHO_H_S_P(BC%II,BC%JJ+1,BC%KK),RHO_H_S_P(BC%II,BC%JJ+1:BC%JJ+2,BC%KK)/) U_TEMP(1,1,1) = VV(BC%II,BC%JJ+1,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFY(BC%II,BC%JJ+1,BC%KK) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) FY_H_S(BC%II,BC%JJ+1,BC%KK) = F_TEMP(1,1,1) ENDIF CASE(-2) OFF_WALL_SELECT_1 IF ((VV(BC%II,BC%JJ-2,BC%KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ-1,BC%KK))%WALL_INDEX(-2)>0)) THEN Z_TEMP(1,1:3,1) = (/RHO_H_S_P(BC%II,BC%JJ-2:BC%JJ-1,BC%KK),RHO_H_S_P(BC%II,BC%JJ-1,BC%KK)/) U_TEMP(1,1,1) = VV(BC%II,BC%JJ-2,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFY(BC%II,BC%JJ-2,BC%KK) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) FY_H_S(BC%II,BC%JJ-2,BC%KK) = F_TEMP(1,1,1) ENDIF CASE( 3) OFF_WALL_SELECT_1 IF ((WW(BC%II,BC%JJ,BC%KK+1)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ,BC%KK+1))%WALL_INDEX(+3)>0)) THEN Z_TEMP(1,1,0:2) = (/RHO_H_S_P(BC%II,BC%JJ,BC%KK+1),RHO_H_S_P(BC%II,BC%JJ,BC%KK+1:BC%KK+2)/) U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK+1) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK+1) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) FZ_H_S(BC%II,BC%JJ,BC%KK+1) = F_TEMP(1,1,1) ENDIF CASE(-3) OFF_WALL_SELECT_1 IF ((WW(BC%II,BC%JJ,BC%KK-2)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ,BC%KK-1))%WALL_INDEX(-3)>0)) THEN Z_TEMP(1,1,1:3) = (/RHO_H_S_P(BC%II,BC%JJ,BC%KK-2:BC%KK-1),RHO_H_S_P(BC%II,BC%JJ,BC%KK-1)/) U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-2) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) + B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK-2) + CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) FZ_H_S(BC%II,BC%JJ,BC%KK-2) = F_TEMP(1,1,1) ENDIF END SELECT OFF_WALL_SELECT_1 @@ -971,117 +976,7 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S) ENDDO !$OMP END PARALLEL DO -END SUBROUTINE ENTHALPY_ADVECTION - - -SUBROUTINE SPECIES_ADVECTION_PART_1 - -USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE -USE PHYSICAL_FUNCTIONS, ONLY: GET_MOLECULAR_WEIGHT -REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P -REAL(EB), POINTER, DIMENSION(:,:,:,:) :: FX_ZZ,FY_ZZ,FZ_ZZ -REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP -REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP -INTEGER :: I,J,K,IW,N -TYPE(WALL_TYPE), POINTER :: WC -TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC -TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 - -FX_ZZ=>SWORK1 -FY_ZZ=>SWORK2 -FZ_ZZ=>SWORK3 -RHO_Z_P=>WORK_PAD - -! Species face values - -SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS - - !$OMP PARALLEL DO - DO K=-1,KBP1+1 - DO J=-1,JBP1+1 - DO I=-1,IBP1+1 - RHO_Z_P(I,J,K) = RHOP(I,J,K)*ZZP(I,J,K,N) - ENDDO - ENDDO - ENDDO - !$OMP END PARALLEL DO - - ! Compute scalar face values - - CALL GET_SCALAR_FACE_VALUE(UU,RHO_Z_P,FX_ZZ(:,:,:,N),0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE(VV,RHO_Z_P,FY_ZZ(:,:,:,N),1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE(WW,RHO_Z_P,FZ_ZZ(:,:,:,N),1,IBAR,1,JBAR,0,KBAR,3,I_FLUX_LIMITER) - - !$OMP PARALLEL DO PRIVATE(IW,WC,BC,B1,U_TEMP,Z_TEMP,F_TEMP) - WALL_LOOP_2: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS - WC=>WALL(IW) - IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_LOOP_2 - BC=>BOUNDARY_COORD(WC%BC_INDEX) - B1=>BOUNDARY_PROP1(WC%B1_INDEX) - - ! Overwrite first off-wall advective flux if flow is away from the wall and if the face is not also a wall cell - - OFF_WALL_IF_2: IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY .AND. WC%BOUNDARY_TYPE/=OPEN_BOUNDARY) THEN - - OFF_WALL_SELECT_2: SELECT CASE(BC%IOR) - CASE( 1) OFF_WALL_SELECT_2 - ! ghost FX/UU(II+1) - ! /// II /// II+1 | II+2 | ... - ! ^ WALL_INDEX(II+1,+1) - IF ((UU(BC%II+1,BC%JJ,BC%KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II+1,BC%JJ,BC%KK))%WALL_INDEX(+1)>0)) THEN - Z_TEMP(0:2,1,1) = (/RHO_Z_P(BC%II+1,BC%JJ,BC%KK),RHO_Z_P(BC%II+1:BC%II+2,BC%JJ,BC%KK)/) - U_TEMP(1,1,1) = UU(BC%II+1,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - FX_ZZ(BC%II+1,BC%JJ,BC%KK,N) = F_TEMP(1,1,1) - ENDIF - CASE(-1) OFF_WALL_SELECT_2 - ! FX/UU(II-2) ghost - ! ... | II-2 | II-1 /// II /// - ! ^ WALL_INDEX(II-1,-1) - IF ((UU(BC%II-2,BC%JJ,BC%KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II-1,BC%JJ,BC%KK))%WALL_INDEX(-1)>0)) THEN - Z_TEMP(1:3,1,1) = (/RHO_Z_P(BC%II-2:BC%II-1,BC%JJ,BC%KK),RHO_Z_P(BC%II-1,BC%JJ,BC%KK)/) - U_TEMP(1,1,1) = UU(BC%II-2,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - FX_ZZ(BC%II-2,BC%JJ,BC%KK,N) = F_TEMP(1,1,1) - ENDIF - CASE( 2) OFF_WALL_SELECT_2 - IF ((VV(BC%II,BC%JJ+1,BC%KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ+1,BC%KK))%WALL_INDEX(+2)>0)) THEN - Z_TEMP(1,0:2,1) = (/RHO_Z_P(BC%II,BC%JJ+1,BC%KK),RHO_Z_P(BC%II,BC%JJ+1:BC%JJ+2,BC%KK)/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ+1,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - FY_ZZ(BC%II,BC%JJ+1,BC%KK,N) = F_TEMP(1,1,1) - ENDIF - CASE(-2) OFF_WALL_SELECT_2 - IF ((VV(BC%II,BC%JJ-2,BC%KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ-1,BC%KK))%WALL_INDEX(-2)>0)) THEN - Z_TEMP(1,1:3,1) = (/RHO_Z_P(BC%II,BC%JJ-2:BC%JJ-1,BC%KK),RHO_Z_P(BC%II,BC%JJ-1,BC%KK)/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ-2,BC%KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - FY_ZZ(BC%II,BC%JJ-2,BC%KK,N) = F_TEMP(1,1,1) - ENDIF - CASE( 3) OFF_WALL_SELECT_2 - IF ((WW(BC%II,BC%JJ,BC%KK+1)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ,BC%KK+1))%WALL_INDEX(+3)>0)) THEN - Z_TEMP(1,1,0:2) = (/RHO_Z_P(BC%II,BC%JJ,BC%KK+1),RHO_Z_P(BC%II,BC%JJ,BC%KK+1:BC%KK+2)/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK+1) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - FZ_ZZ(BC%II,BC%JJ,BC%KK+1,N) = F_TEMP(1,1,1) - ENDIF - CASE(-3) OFF_WALL_SELECT_2 - IF ((WW(BC%II,BC%JJ,BC%KK-2)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(BC%II,BC%JJ,BC%KK-1))%WALL_INDEX(-3)>0)) THEN - Z_TEMP(1,1,1:3) = (/RHO_Z_P(BC%II,BC%JJ,BC%KK-2:BC%KK-1),RHO_Z_P(BC%II,BC%JJ,BC%KK-1)/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-2) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - FZ_ZZ(BC%II,BC%JJ,BC%KK-2,N) = F_TEMP(1,1,1) - ENDIF - END SELECT OFF_WALL_SELECT_2 - - ENDIF OFF_WALL_IF_2 - - ENDDO WALL_LOOP_2 - !$OMP END PARALLEL DO - -ENDDO SPECIES_LOOP - -END SUBROUTINE SPECIES_ADVECTION_PART_1 +END SUBROUTINE ENTHALPY_ADVECTION_NEW SUBROUTINE SPECIES_ADVECTION_PART_1_NEW diff --git a/Source/dump.f90 b/Source/dump.f90 index 5eafda355dd..d11f90bb2cb 100644 --- a/Source/dump.f90 +++ b/Source/dump.f90 @@ -8482,7 +8482,7 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z USE PHYSICAL_FUNCTIONS, ONLY: GET_MASS_FRACTION,FED,FIC,GET_SPECIFIC_HEAT,RELATIVE_HUMIDITY, & GET_CONDUCTIVITY,GET_MOLECULAR_WEIGHT,GET_MASS_FRACTION_ALL,GET_ENTHALPY,GET_SENSIBLE_ENTHALPY, & GET_VISCOSITY,GET_POTENTIAL_TEMPERATURE,GET_SPECIFIC_GAS_CONSTANT,& - SURFACE_DENSITY + SURFACE_DENSITY, CALC_EQUIV_RATIO USE COMP_FUNCTIONS, ONLY : CURRENT_TIME,SYSTEM_MEM_USAGE USE RADCONS, ONLY: WL_LOW, WL_HIGH, RADTMP USE RAD, ONLY: BLACKBODY_FRACTION @@ -8500,7 +8500,7 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z EXPON,Y_SPECIES,MEC,Y_SPECIES2,Y_H2O,R_Y_H2O,R_DN,SGN,Y_ALL(N_SPECIES),H_S,D_Z_N(0:I_MAX_TEMP),& DISSIPATION_RATE,S11,S22,S33,S12,S13,S23,DUDX,DUDY,DUDZ,DVDX,DVDY,DVDZ,DWDX,DWDY,DWDZ,ONTHDIV,SS,ETA,DELTA,R_DX2,& UVW,UODX,VODY,WODZ,XHAT,ZHAT,BBF,GAMMA_LOC,VC,VOL,PHI,GAS_PHASE_OUTPUT_CC,& - GAS_PHASE_OUTPUT_CFA,CFACE_AREA,VELOCITY_COMPONENT(1:3),ATOTV(1:3),TMP_F,R_D,MW,RHO_AIR,PROBE_TMP + GAS_PHASE_OUTPUT_CFA,CFACE_AREA,VELOCITY_COMPONENT(1:3),ATOTV(1:3),TMP_F,R_D,MW,PROBE_TMP INTEGER :: N,I,J,K,NN,IL,III,JJJ,KKK,IP,JP,KP,FED_ACTIVITY,IP1,JP1,KP1,IM1,JM1,KM1,IIM1,JJM1,KKM1,NR,NS,RAM,& ICC,JCC,NCELL,AXIS,ICF,NFACE,JCF,JCC_LO,JCC_HI,PDPA_FORMULA,IC REAL(FB) :: RN @@ -9009,24 +9009,27 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z GAS_PHASE_OUTPUT_RES = FED(ZZ_GET,RSUM(II,JJ,KK),FED_ACTIVITY) CASE(110) ! THERMOCOUPLE - TMP_G = TMP(II,JJ,KK) - IF (PY%HEAT_TRANSFER_COEFFICIENT<0._EB) THEN - UU = U(II,JJ,KK) - VV = V(II,JJ,KK) - WW = W(II,JJ,KK) - VEL2 = UU**2+VV**2+WW**2 - ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES) - CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP(II,JJ,KK)) - CALL GET_CONDUCTIVITY(ZZ_GET,K_G,TMP(II,JJ,KK)) - RE_D = RHO(II,JJ,KK)*SQRT(VEL2)*PY%DIAMETER/MU_G - CALL FORCED_CONVECTION_MODEL(NUSSELT,RE_D,PR_ONTH,SURF_SPHERICAL) - H_TC = NUSSELT*K_G/PY%DIAMETER - ELSE - H_TC = PY%HEAT_TRANSFER_COEFFICIENT + IF (T > T_BEGIN) THEN + TMP_G = TMP(II,JJ,KK) + IF (PY%HEAT_TRANSFER_COEFFICIENT<0._EB) THEN + UU = U(II,JJ,KK) + VV = V(II,JJ,KK) + WW = W(II,JJ,KK) + VEL2 = UU**2+VV**2+WW**2 + ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES) + CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP(II,JJ,KK)) + CALL GET_CONDUCTIVITY(ZZ_GET,K_G,TMP(II,JJ,KK)) + RE_D = RHO(II,JJ,KK)*SQRT(VEL2)*PY%DIAMETER/MU_G + CALL FORCED_CONVECTION_MODEL(NUSSELT,RE_D,PR_ONTH,SURF_SPHERICAL) + H_TC = NUSSELT*K_G/PY%DIAMETER + ELSE + H_TC = PY%HEAT_TRANSFER_COEFFICIENT + ENDIF + FAC = 6._EB/(PY%DENSITY*PY%SPECIFIC_HEAT(NINT(DV%TMP_L))*PY%DIAMETER)*DT + DV%TMP_L = (DV%TMP_L + FAC*(H_TC*(TMP_G-0.5_EB*DV%TMP_L) + & + PY%EMISSIVITY*(0.25_EB*UII(II,JJ,KK)+SIGMA*DV%TMP_L**4))) / & + (1._EB + FAC*(0.5_EB*H_TC+2._EB*PY%EMISSIVITY*SIGMA*DV%TMP_L**3)) ENDIF - RHS = (6._EB/(PY%DENSITY*PY%SPECIFIC_HEAT*PY%DIAMETER))* & - ( H_TC*(TMP_G-DV%TMP_L) + PY%EMISSIVITY*(0.25_EB*UII(II,JJ,KK)-SIGMA*DV%TMP_L**4) ) - IF (T>T_BEGIN) DV%TMP_L = DV%TMP_L + DT*RHS GAS_PHASE_OUTPUT_RES = DV%TMP_L - TMPM CASE(111:113) ! ENTHALPY FLUX @@ -9059,17 +9062,17 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z VV = V(II,JJ,KK) WW = W(II,JJ,KK) VEL2 = UU**2+VV**2+WW**2 + VEL = SQRT(VEL2) DP = 0.5_EB*VEL2*RHO(II,JJ,KK) COSTHETA = (UU*ORIENTATION_VECTOR(1,DV%ORIENTATION_INDEX)+VV*ORIENTATION_VECTOR(2,DV%ORIENTATION_INDEX)+& - WW*ORIENTATION_VECTOR(3,DV%ORIENTATION_INDEX))/SQRT(VEL2) - FAC = -2.308_EB*ABS(COSTHETA)**3 + 2.533_EB*ABS(COSTHETA)**2 + 0.7847_EB*ABS(COSTHETA) - VEL = FAC*SQRT(VEL2) + WW*ORIENTATION_VECTOR(3,DV%ORIENTATION_INDEX))/VEL ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES) CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP(II,JJ,KK)) RE_D = MIN(3800._EB,MAX(40._EB,RHO(II,JJ,KK)*VEL*PY%PROBE_DIAMETER/MU_G)) + FAC = MAX(0._EB,-2.308_EB*ABS(COSTHETA)**3 + 2.533_EB*ABS(COSTHETA)**2 + 0.7847_EB*ABS(COSTHETA) - 0.0097_EB) + VEL = FAC*VEL FAC = 1.533_EB-0.001366_EB*RE_D+0.000001688_EB*RE_D**2-0.0000000009706_EB*RE_D**3+& 0.0000000000002555_EB*RE_D**4-2.484E-17_EB*RE_D**5 - RHO_AIR = 350.9736_EB/PROBE_TMP !350 is 0.0288 101325/ 8.314472 GAS_PHASE_OUTPUT_RES = SIGN(1._EB,COSTHETA)*VEL*PY%CALIBRATION_CONSTANT*FAC CASE(130) ! EXTINCTION @@ -9140,6 +9143,10 @@ REAL(EB) RECURSIVE FUNCTION GAS_PHASE_OUTPUT(T,DT,NM,II,JJ,KK,IND,IND2,Y_INDEX,Z GAS_PHASE_OUTPUT_RES = GAS_PHASE_OUTPUT_RES + & ZZ_GET(NS)*SPECIES_MIXTURE(NS)%ATOMS(ELEM_INDX)*ELEMENT(ELEM_INDX)%MASS/SPECIES_MIXTURE(NS)%MW ENDDO + CASE(145) ! EQUIVALENCE RATIO + ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES) + GAS_PHASE_OUTPUT_RES = 0._EB + CALL CALC_EQUIV_RATIO(ZZ_GET(1:N_TRACKED_SPECIES), GAS_PHASE_OUTPUT_RES) CASE(150) ! SUM LUMPED VOLUME FRACTIONS ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(II,JJ,KK,1:N_TRACKED_SPECIES) CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW) @@ -13138,10 +13145,10 @@ SUBROUTINE DUMP_CVODE_SUBSTEPS() DO COLI = 1, NCOLS IF (COLI == NCOLS) THEN ! Writing the last column without a trailing comma - WRITE(LU_CVODE_SUBSTEPS, '(F18.5)') CVODE_SUBSTEP_DATA(ROWI, COLI) + WRITE(LU_CVODE_SUBSTEPS, '(ES12.5)') CVODE_SUBSTEP_DATA(ROWI, COLI) ELSE ! Writing columns with commas - WRITE(LU_CVODE_SUBSTEPS, '(F18.5, A)', ADVANCE="NO") CVODE_SUBSTEP_DATA(ROWI, COLI), ',' + WRITE(LU_CVODE_SUBSTEPS, '(ES12.5, A)', ADVANCE="NO") CVODE_SUBSTEP_DATA(ROWI, COLI), ',' ENDIF ENDDO ENDDO diff --git a/Source/fire.f90 b/Source/fire.f90 index 146f5917281..ef6d70462a6 100644 --- a/Source/fire.f90 +++ b/Source/fire.f90 @@ -53,7 +53,7 @@ SUBROUTINE COMBUSTION_LOAD_BALANCED(T,DT) ! Set CVODES options IF (.NOT. COMBUSTION_INIT) THEN COMBUSTION_INIT = .TRUE. - NVAR_TO_SEND = N_TRACKED_SPECIES +6 ! NS, TEMP, RHO, PRES, MU, DELTA, VOL + NVAR_TO_SEND = N_TRACKED_SPECIES +7 ! NS, TEMP, RHO, PRES, MU, DELTA, VOL, IGN_ZN NVAR_TO_RECEIVE = N_TRACKED_SPECIES +4 ! NS, Q_OUT, MIX_TIME_OUT, CHI_R_OUT, CHEM_SUBIT_TMP_OUT ENDIF @@ -107,6 +107,7 @@ SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED(T,DT) TYPE (CC_CUTCELL_TYPE), POINTER :: CC LOGICAL :: IS_CHEM_ACTIVE REAL(EB) :: MIX_TIME_OUT, Q_OUT, CHI_R_OUT, MYTEMP, MYRHO, MYMU, DELTA, VOL, TNOW2 +INTEGER :: IGN_ZN Q_EXISTS = .FALSE. @@ -125,12 +126,13 @@ SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED(T,DT) DO I=1,IBAR ZZ_GET = ZZ(I,J,K,1:N_TRACKED_SPECIES) PRES = PBAR(K,PRESSURE_ZONE(I,J,K)) + RHO(I,J,K)*(H(I,J,K)-KRES(I,J,K)) - CALL CHECK_CHEMICALLY_ACTIVE_STATE(ZZ_GET, TMP(I,J,K), I, J, K, .FALSE., IS_CHEM_ACTIVE) + CALL CHECK_CHEMICALLY_ACTIVE_STATE(ZZ_GET, TMP(I,J,K), I, J, K, .FALSE., IS_CHEM_ACTIVE, IGN_ZN) IF (STOP_STATUS/=NO_STOP) RETURN IF (.NOT. IS_CHEM_ACTIVE) CYCLE NCHEM_ACTIVE_CELLS_AND_CC = NCHEM_ACTIVE_CELLS_AND_CC + 1 NCHEM_ACTIVE_CELLS = NCHEM_ACTIVE_CELLS +1 CHEM_ACTIVE_CELLS(NCHEM_ACTIVE_CELLS,1:3)=[I, J ,K] + IF(IGN_ZN > 0) CHEM_ACTIVE_CELLS(NCHEM_ACTIVE_CELLS,4) = IGN_ZN ENDDO ENDDO ENDDO @@ -158,12 +160,13 @@ SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED(T,DT) IF ( ABS(CUT_CELL(ICC)%VOLUME(JCC)/VCELL) < 1.E-12_EB ) CYCLE ZZ_GET = CUT_CELL(ICC)%ZZ(1:N_TRACKED_SPECIES,JCC) PRES = PBAR(K,PRESSURE_ZONE(I,J,K)) + RHO(I,J,K)*(H(I,J,K)-KRES(I,J,K)) - CALL CHECK_CHEMICALLY_ACTIVE_STATE(ZZ_GET, CUT_CELL(ICC)%TMP(JCC), I, J, K, .TRUE., IS_CHEM_ACTIVE) + CALL CHECK_CHEMICALLY_ACTIVE_STATE(ZZ_GET, CUT_CELL(ICC)%TMP(JCC), I, J, K, .TRUE., IS_CHEM_ACTIVE, IGN_ZN) IF (STOP_STATUS/=NO_STOP) RETURN IF (.NOT. IS_CHEM_ACTIVE) CYCLE NCHEM_ACTIVE_CELLS_AND_CC = NCHEM_ACTIVE_CELLS_AND_CC + 1 NCHEM_ACTIVE_CC = NCHEM_ACTIVE_CC +1 CHEM_ACTIVE_CC(NCHEM_ACTIVE_CC,1:2)=[ICC, JCC] + IF(IGN_ZN > 0) CHEM_ACTIVE_CC(NCHEM_ACTIVE_CC,3) = IGN_ZN ENDDO ENDDO ENDIF @@ -171,6 +174,8 @@ SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED(T,DT) MESHES(NM)%NCHEM_ACTIVE_CC=NCHEM_ACTIVE_CC ENDDO + + !------ !STEP2: Solve chemistry by distributing the chemistry load accross all MPI processes for parallel runs. ! For 1 MPI process (serial) distribution is not needed. @@ -204,12 +209,12 @@ SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED(T,DT) MYMU = REALS_TO_RECV_ARRAY(N_TRACKED_SPECIES+4,INDX) DELTA = REALS_TO_RECV_ARRAY(N_TRACKED_SPECIES+5,INDX) VOL = REALS_TO_RECV_ARRAY(N_TRACKED_SPECIES+6,INDX) - + IGN_ZN = INT(REALS_TO_RECV_ARRAY(N_TRACKED_SPECIES+7,INDX)) !*************************************************************************************** ! Call combustion integration routine for Cartesian cell (I,J,K) CALL COMBUSTION_MODEL( T,DT,ZZ_GET,Q_OUT,MIX_TIME_OUT,CHI_R_OUT,& CHEM_SUBIT_TMP_OUT,REAC_SOURCE_TERM_TMP,Q_REAC_TMP,& - MYTEMP,MYRHO,PRES,MYMU,DELTA,VOL,ZETA_0) + MYTEMP,MYRHO,PRES,MYMU,DELTA,VOL,ZETA_0,IGN_ZN) !*************************************************************************************** RESULTS_TO_SEND_ARRAY(1:N_TRACKED_SPECIES,INDX) = ZZ_GET @@ -249,6 +254,7 @@ SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED(T,DT) DO NC=1,NCHEM_ACTIVE_CELLS I=CHEM_ACTIVE_CELLS(NC,1); J=CHEM_ACTIVE_CELLS(NC,2); K= CHEM_ACTIVE_CELLS(NC,3) + IGN_ZN=CHEM_ACTIVE_CELLS(NC,4) ZZ_GET = ZZ(I,J,K,1:N_TRACKED_SPECIES) PRES = PBAR(K,PRESSURE_ZONE(I,J,K)) + RHO(I,J,K)*(H(I,J,K)-KRES(I,J,K)) DZZ = ZZ_GET ! store old ZZ for divergence term @@ -257,7 +263,7 @@ SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED(T,DT) CALL COMBUSTION_MODEL( T,DT,ZZ_GET,Q(I,J,K),MIX_TIME(I,J,K),CHI_R(I,J,K),& CHEM_SUBIT_TMP,REAC_SOURCE_TERM_TMP,Q_REAC_TMP,& TMP(I,J,K),RHO(I,J,K),PRES, MU(I,J,K),& - LES_FILTER_WIDTH(I,J,K),DX(I)*DY(J)*DZ(K),ZETA_0,IIC=I,JJC=J,KKC=K ) + LES_FILTER_WIDTH(I,J,K),DX(I)*DY(J)*DZ(K),ZETA_0,IGN_ZN,IIC=I,JJC=J,KKC=K ) !*************************************************************************************** !IF (STOP_STATUS/=NO_STOP) RETURN IF (OUTPUT_CHEM_IT) CHEM_SUBIT(I,J,K) = CHEM_SUBIT_TMP @@ -269,6 +275,7 @@ SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED(T,DT) !$OMP DO SCHEDULE(DYNAMIC) DO NC=1,NCHEM_ACTIVE_CC ICC= CHEM_ACTIVE_CC(NC,1); JCC= CHEM_ACTIVE_CC(NC,2) + IGN_ZN=CHEM_ACTIVE_CC(NC,3) CC => CUT_CELL(ICC); I = CC%IJK(IAXIS); J = CC%IJK(JAXIS); K = CC%IJK(KAXIS) CC%CHI_R(JCC) = 0._EB ZZ_GET = CC%ZZ(1:N_TRACKED_SPECIES,JCC) @@ -281,7 +288,7 @@ SUBROUTINE COMBUSTION_GENERAL_LOAD_BALANCED(T,DT) CC%CHI_R(JCC),& CHEM_SUBIT_TMP,REAC_SOURCE_TERM_TMP,Q_REAC_TMP,& CC%TMP(JCC),CC%RHO(JCC),PRES,MU(I,J,K),& - LES_FILTER_WIDTH(I,J,K),CC%VOLUME(JCC),ZETA_0,IIC=I,JJC=J,KKC=K) + LES_FILTER_WIDTH(I,J,K),CC%VOLUME(JCC),ZETA_0,IGN_ZN,IIC=I,JJC=J,KKC=K) !*************************************************************************************** CALL SET_SPECIES_SOURCE_TERM_CUTCELL(DT, ICC, JCC, ZZ_GET, DZZ, REAC_SOURCE_TERM_TMP, Q_REAC_TMP) END DO @@ -349,27 +356,32 @@ SUBROUTINE CHECK_REACTION (ZZ_GET, DO_REACTION) END SUBROUTINE CHECK_REACTION -SUBROUTINE CHECK_CHEMICALLY_ACTIVE_STATE (ZZ_GET, TMP, I, J , K, IS_CUT_CELL, CHEM_ACTIVE) +SUBROUTINE CHECK_CHEMICALLY_ACTIVE_STATE (ZZ_GET, TMP, I, J , K, IS_CUT_CELL, CHEM_ACTIVE, IGN_ZN) +USE DEVICE_VARIABLES, ONLY: DEVICE USE PHYSICAL_FUNCTIONS, ONLY: IS_REALIZABLE, CALC_EQUIV_RATIO USE COMPLEX_GEOMETRY, ONLY : CC_CGSC, CC_GASPHASE -USE CHEMCONS, ONLY: EQUIV_RATIO_CHECK, MIN_EQUIV_RATIO, MAX_EQUIV_RATIO +USE CHEMCONS, ONLY: EQUIV_RATIO_CHECK,MIN_EQUIV_RATIO,MAX_EQUIV_RATIO,N_IGNITION_ZONES,IGNITION_ZONES,USE_MIXED_ZN_AFT_TMP + REAL(EB), INTENT(IN) :: ZZ_GET(1:N_TRACKED_SPECIES) REAL(EB), INTENT(IN) :: TMP INTEGER, INTENT(IN), OPTIONAL :: I, J, K LOGICAL, INTENT(IN) :: IS_CUT_CELL LOGICAL, INTENT(INOUT) :: CHEM_ACTIVE -LOGICAL :: DO_REACTION, REALIZABLE +INTEGER, INTENT(INOUT) :: IGN_ZN +LOGICAL :: DO_REACTION, REALIZABLE, LES_EQUIV_CHECK REAL(EB) :: EQUIV +INTEGER :: IZ CHEM_ACTIVE = .TRUE. +IGN_ZN = -1 +LES_EQUIV_CHECK = .FALSE. IF (CELL(CELL_INDEX(I,J,K))%SOLID) THEN CHEM_ACTIVE = .FALSE. RETURN ENDIF - IF (CC_IBM .AND. .NOT. IS_CUT_CELL) THEN ! Check if cell is regular gas phase cell, if not return. IF (CCVAR(I,J,K,CC_CGSC) /= CC_GASPHASE) THEN CHEM_ACTIVE = .FALSE. @@ -377,9 +389,42 @@ SUBROUTINE CHECK_CHEMICALLY_ACTIVE_STATE (ZZ_GET, TMP, I, J , K, IS_CUT_CELL, CH ENDIF ENDIF -IF (.NOT.ALL(REACTION%FAST_CHEMISTRY) .AND. TMP < FINITE_RATE_MIN_TEMP) THEN - CHEM_ACTIVE = .FALSE. - RETURN +IF (.NOT.ALL(REACTION%FAST_CHEMISTRY)) THEN + ! Check for ignition zones for detailed chemistry + IF (COMBUSTION_ODE_SOLVER==CVODE_SOLVER) THEN + IF (N_IGNITION_ZONES > 0) THEN + DO IZ=1,N_IGNITION_ZONES + IF (IGNITION_ZONES(IZ)%DEVC_INDEX>0) THEN + IF (.NOT.DEVICE(IGNITION_ZONES(IZ)%DEVC_INDEX)%CURRENT_STATE) THEN + CYCLE + ENDIF + ENDIF + IF (XC(I)>=IGNITION_ZONES(IZ)%X1 .AND. XC(I)<=IGNITION_ZONES(IZ)%X2 .AND. & + YC(J)>=IGNITION_ZONES(IZ)%Y1 .AND. YC(J)<=IGNITION_ZONES(IZ)%Y2 .AND. & + ZC(K)>=IGNITION_ZONES(IZ)%Z1 .AND. ZC(K)<=IGNITION_ZONES(IZ)%Z2) THEN + CALL CALC_EQUIV_RATIO(ZZ_GET(1:N_TRACKED_SPECIES), EQUIV) + ! No need to solve for no fuel (phi < MIN ) case and only fuel (phi is rich > MAX) case. + IF (EQUIV .GT. MIN_EQUIV_RATIO .AND. EQUIV .LT. MAX_EQUIV_RATIO) THEN + IGN_ZN = IZ + EXIT + ENDIF + ENDIF + ENDDO + ENDIF + ENDIF + + IF (IGN_ZN .LE. 0 .AND. SIM_MODE .NE. DNS_MODE .AND. USE_MIXED_ZN_AFT_TMP) THEN ! LES and other modes. Check equiv ratio + CALL CALC_EQUIV_RATIO(ZZ_GET(1:N_TRACKED_SPECIES), EQUIV) + ! No need to solve for no fuel (phi < MIN ) case and only fuel (phi is rich > MAX) case. + IF (EQUIV .GT. MIN_EQUIV_RATIO .AND. EQUIV .LT. MAX_EQUIV_RATIO) THEN + LES_EQUIV_CHECK = .TRUE. ! Chemically active + ENDIF + ENDIF + + IF (IGN_ZN .LE. 0 .AND. .NOT. LES_EQUIV_CHECK .AND. TMP < FINITE_RATE_MIN_TEMP) THEN + CHEM_ACTIVE = .FALSE. + RETURN + ENDIF ENDIF IF (CHECK_REALIZABILITY) THEN @@ -393,7 +438,6 @@ SUBROUTINE CHECK_CHEMICALLY_ACTIVE_STATE (ZZ_GET, TMP, I, J , K, IS_CUT_CELL, CH RETURN ENDIF ENDIF - CALL CHECK_REACTION (ZZ_GET, DO_REACTION) IF (.NOT.DO_REACTION) THEN CHEM_ACTIVE = .FALSE. @@ -408,7 +452,6 @@ SUBROUTINE CHECK_CHEMICALLY_ACTIVE_STATE (ZZ_GET, TMP, I, J , K, IS_CUT_CELL, CH RETURN ENDIF ENDIF - END SUBROUTINE CHECK_CHEMICALLY_ACTIVE_STATE @@ -569,6 +612,7 @@ SUBROUTINE DISTRIBUTE_CELLS_ACCROSS_MPI_PROCESSES (NRECEIVE_CELLS) REALS_TO_SEND_ARRAY(N_TRACKED_SPECIES+4,INDX) = MU(I,J,K) REALS_TO_SEND_ARRAY(N_TRACKED_SPECIES+5,INDX) = LES_FILTER_WIDTH(I,J,K) REALS_TO_SEND_ARRAY(N_TRACKED_SPECIES+6,INDX) = DX(I)*DY(J)*DZ(K) + REALS_TO_SEND_ARRAY(N_TRACKED_SPECIES+7,INDX) = REAL(CHEM_ACTIVE_CELLS(NC,4),EB) ! Is ignition zone RANK_TO=RANK_TO+1; IF(RANK_TO==N_MPI_PROCESSES+1) RANK_TO=1 ENDDO IF (CC_IBM) THEN @@ -590,6 +634,7 @@ SUBROUTINE DISTRIBUTE_CELLS_ACCROSS_MPI_PROCESSES (NRECEIVE_CELLS) REALS_TO_SEND_ARRAY(N_TRACKED_SPECIES+4,INDX) = MU(I,J,K) REALS_TO_SEND_ARRAY(N_TRACKED_SPECIES+5,INDX) = LES_FILTER_WIDTH(I,J,K) REALS_TO_SEND_ARRAY(N_TRACKED_SPECIES+6,INDX) = CC%VOLUME(JCC) + REALS_TO_SEND_ARRAY(N_TRACKED_SPECIES+7,INDX) = REAL(CHEM_ACTIVE_CC(NC,3),EB) RANK_TO=RANK_TO+1; IF(RANK_TO==N_MPI_PROCESSES+1) RANK_TO=1 ENDDO ENDIF @@ -718,11 +763,12 @@ END SUBROUTINE GATHER_CELLS_FROM_MPI_PROCESSES SUBROUTINE COMBUSTION_MODEL(T,DT,ZZ_GET,Q_OUT,MIX_TIME_OUT,CHI_R_OUT,CHEM_SUBIT_OUT,REAC_SOURCE_TERM_OUT,Q_REAC_OUT,& - TMP_IN,RHO_IN,PRES_IN,MU_IN,DELTA,CELL_VOLUME,ZETA_0_IN,IIC,JJC,KKC) + TMP_IN,RHO_IN,PRES_IN,MU_IN,DELTA,CELL_VOLUME,ZETA_0_IN,IGN_ZN,IIC,JJC,KKC) USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP USE PHYSICAL_FUNCTIONS, ONLY: GET_REALIZABLE_MF USE COMP_FUNCTIONS, ONLY: SHUTDOWN -USE CHEMCONS, ONLY: ODE_MIN_ATOL, FLAME_THICK_FACTOR +USE CHEMCONS, ONLY: ODE_MIN_ATOL +INTEGER, INTENT(IN) :: IGN_ZN INTEGER, INTENT(IN), OPTIONAL :: IIC,JJC,KKC REAL(EB), INTENT(IN) :: T,DT,RHO_IN,PRES_IN,MU_IN,DELTA,CELL_VOLUME,ZETA_0_IN REAL(EB), INTENT(OUT) :: Q_OUT,MIX_TIME_OUT,CHI_R_OUT,REAC_SOURCE_TERM_OUT(N_TRACKED_SPECIES),Q_REAC_OUT(N_REACTIONS) @@ -881,9 +927,9 @@ SUBROUTINE COMBUSTION_MODEL(T,DT,ZZ_GET,Q_OUT,MIX_TIME_OUT,CHI_R_OUT,CHEM_SUBIT_ DO NS =1,N_TRACKED_SPECIES ATOL(NS) = DBLE(SPECIES_MIXTURE(NS)%ODE_ABS_ERROR) ENDDO - IF (1==2) WRITE(LU_ERR,*) PRES_IN ! To avoid unused variable error. + IF (1==2) WRITE(LU_ERR,*) PRES_IN, IGN_ZN ! To avoid unused variable error. #ifdef WITH_SUNDIALS - CALL CVODE(ZZ_MIXED,TMP_IN,PRES_IN,ZETA_0, ZETA,TAU_MIX,CELL_MASS,T,DT_SUB,GLOBAL_ODE_REL_ERROR, ATOL) + CALL CVODE(ZZ_MIXED,TMP_IN,PRES_IN,ZETA_0, ZETA,TAU_MIX,CELL_MASS,IGN_ZN,T,DT_SUB,GLOBAL_ODE_REL_ERROR, ATOL) #endif ZETA_0 = ZETA Q_REAC_SUB = 0._EB @@ -926,12 +972,6 @@ SUBROUTINE COMBUSTION_MODEL(T,DT,ZZ_GET,Q_OUT,MIX_TIME_OUT,CHI_R_OUT,CHEM_SUBIT_ ENDDO INTEGRATION_LOOP - -! Consider thickened flame model consideration here - IF (COMBUSTION_ODE_SOLVER == CVODE_SOLVER .AND. FLAME_THICK_FACTOR > 1._EB) THEN - ZZ_GET = ZZ_0 + (ZZ_GET-ZZ_0)/(FLAME_THICK_FACTOR) - ENDIF - ! Compute heat release rate Q_OUT = -RHO_IN*SUM(SPECIES_MIXTURE%H_F*(ZZ_GET-ZZ_0))/DT ! FDS Tech Guide (5.47) @@ -998,27 +1038,45 @@ END SUBROUTINE COMBUSTION_MODEL !> \param ATOL is the absolute error tolerance array for the species (REAL_EB) #ifdef WITH_SUNDIALS -SUBROUTINE CVODE(ZZ, TMP_IN, PRES_IN, ZETA_IN, ZETA_OUT, TAU_MIX, CELL_MASS, T_CFD, DT, GLOBAL_ODE_REL_ERROR, ATOL) +SUBROUTINE CVODE(ZZ, TMP_IN, PRES_IN, ZETA_IN, ZETA_OUT, TAU_MIX, CELL_MASS, IGN_ZN, T_CFD, DT, GLOBAL_ODE_REL_ERROR, ATOL) USE PHYSICAL_FUNCTIONS, ONLY : GET_MOLECULAR_WEIGHT -USE CHEMCONS, ONLY: WRITE_CVODE_SUBSTEPS, MAX_CHEM_TIME - +USE CHEMCONS, ONLY: WRITE_CVODE_SUBSTEPS, ZETA_ARTIFICAL_MAX_LIMIT,ZETA_ARTIFICAL_MIN_LIMIT, ZETA_FIRST_STEP_DIV,& + USE_MIXED_ZN_AFT_TMP REAL(EB), INTENT(INOUT) :: ZZ(N_TRACKED_SPECIES) REAL(EB), INTENT(IN) :: ATOL(N_TRACKED_SPECIES) REAL(EB), INTENT(IN) :: TMP_IN,PRES_IN,ZETA_IN,TAU_MIX,CELL_MASS,T_CFD,DT,GLOBAL_ODE_REL_ERROR REAL(EB), INTENT(OUT) ::ZETA_OUT +INTEGER, INTENT(IN) :: IGN_ZN -REAL(EB) :: CC(N_TRACKED_SPECIES), CC_CHEM_TIME(N_TRACKED_SPECIES) -REAL(EB) :: MW, RHO_IN, RHO_OUT, T1, T2, ZETA0, TMP_IN_MOD, TMP_OUT, CHEM_TIME, DT_MOD +REAL(EB) :: CC(N_TRACKED_SPECIES), CC_CHEM_TIME(N_TRACKED_SPECIES), ZZ_IN(N_TRACKED_SPECIES) +REAL(EB) :: MW, RHO_IN, RHO_OUT, T1, T2, ZETA0, ZETA_IN_MOD, TMP_IN_MOD, TMP_OUT, & + CHEM_TIME, DT_MOD, ZETA_ARTF, ZETA_FINAL, ZETA_MAX_LIMIT, ZETA_MIN_LIMIT,AFT INTEGER :: NS, CVODE_CALL_OPTION -LOGICAL :: WRITE_SUBSTEPS +LOGICAL :: WRITE_SUBSTEPS, CALL_CHEM_AGAIN CVODE_CALL_OPTION = 1 ! CV_NORMAL WRITE_SUBSTEPS = .FALSE. DT_MOD = DT +ZETA_IN_MOD = ZETA_IN TMP_IN_MOD = TMP_IN +CALL_CHEM_AGAIN = .FALSE. + +IF(IGN_ZN > 0) THEN + CALL CALC_ADIABATIC_FLAME_TEMPERATURE(ZZ,TMP_IN,AFT) + TMP_IN_MOD = MAX(TMP_IN,AFT) + ZETA_IN_MOD = 0.0_EB +ELSE + IF (SIM_MODE .NE. DNS_MODE .AND. USE_MIXED_ZN_AFT_TMP) THEN + CALL CALC_ADIABATIC_FLAME_TEMPERATURE(ZZ,TMP_IN,AFT) + TMP_IN_MOD = MAX(TMP_IN,AFT) ! TO DO: Ideal would be a equilibrium temperature based on ZZ(:) and TMP_IN + ENDIF +ENDIF +ZETA_OUT = ZETA_IN_MOD*EXP(-DT/TAU_MIX) +IF(ZETA_OUT > ZETA_ARTIFICAL_MAX_LIMIT) RETURN ! The mixing can be ignored due to large mixing time. + CC = 0._EB -! Get the initial mass and concentration of mixed zone +! Calculate RHO based on actual temperrature, such that RHO and Concentration comes out to be same as in the cell. CALL GET_MOLECULAR_WEIGHT(ZZ,MW) RHO_IN = PRES_IN*MW/R0/TMP_IN ! [PR]= Pa, [MW] = g/mol, [R0]= J/K/kmol, [TMP]=K, [RHO]= kg/m3 DO NS =1,N_TRACKED_SPECIES @@ -1027,38 +1085,51 @@ SUBROUTINE CVODE(ZZ, TMP_IN, PRES_IN, ZETA_IN, ZETA_OUT, TAU_MIX, CELL_MASS, T_ WHERE(CC<0._EB) CC=0._EB ! Get the initial mass and concentration of mixed zone -IF(ZETA_IN > ONE_M_EPS) THEN +IF(ZETA_IN_MOD > ONE_M_EPS) THEN !With ZETA_IN =1 the CVODE ODE become too stiff to solve. Hence, performing negligible !artifical mixing based on a chemical time scale. - !1. First find a reasonable chemical time scale. - !2. Do cemistry on the chemical time scale - !3. Set the final zeta0 based on the chemical time + !1. Do cemistry and find a chemical time scale with a cvode call (one substep). + !2. Set the final zeta0 based on the chemical time T1 = 0._EB T2 = DT ZETA0 = 0._EB ! Assume completely mixed. CVODE_CALL_OPTION = 2 ! CV_ONE_STEP CUR_CFD_TIME = T_CFD ! Set current cfd time in cvode, for logging purpose. CC_CHEM_TIME(1:N_TRACKED_SPECIES)=CC(1:N_TRACKED_SPECIES) - CALL CVODE_SERIAL(CC_CHEM_TIME,ZZ, TMP_IN, PRES_IN, ZETA0, TAU_MIX, CELL_MASS, T1,T2, & + CALL CVODE_SERIAL(CC_CHEM_TIME,ZZ, TMP_IN_MOD, TMP_IN, PRES_IN, ZETA0, TAU_MIX, CELL_MASS, T1,T2, & GLOBAL_ODE_REL_ERROR, ATOL, TMP_OUT, CHEM_TIME, WRITE_SUBSTEPS, CVODE_CALL_OPTION) ! Find the chem time scale - IF (CHEM_TIME > MAX_CHEM_TIME) THEN ! Limit artifical mixing. + ! Check 1) too much artificial mixing (low zeta after first step) or + ! 2) too little mixing (high zeta after first step) to avoid steep problem in cvode. + ZETA_ARTF = ZETA_IN_MOD*EXP(-CHEM_TIME/TAU_MIX) ! Zeta after CHEM_TIME + ZETA_FINAL = ZETA_IN_MOD*EXP(-DT/ZETA_FIRST_STEP_DIV/TAU_MIX) ! Allowed final zeta after artifical mixing substep + ZETA_MIN_LIMIT = MAX(ZETA_FINAL,ZETA_ARTIFICAL_MIN_LIMIT) + ZETA_MAX_LIMIT = MAX(ZETA_FINAL,ZETA_ARTIFICAL_MAX_LIMIT) + IF (ZETA_ARTF < ZETA_MIN_LIMIT) THEN ! Limit too much artifical mixing to ZETA_ARTIFICAL_MAX_LIMIT. + T2 = TAU_MIX*(LOG(ZETA_IN_MOD) - LOG(ZETA_MIN_LIMIT)) + CALL_CHEM_AGAIN = .TRUE. + ELSEIF (ZETA_ARTF > ZETA_MAX_LIMIT) THEN ! Check too liitle mixing, then make little more mixing such + ! that cvode is not stiff in subsequent calls. + T2 = TAU_MIX*(LOG(ZETA_IN_MOD) - LOG(ZETA_MAX_LIMIT)) + CALL_CHEM_AGAIN = .TRUE. + ENDIF + + IF (CALL_CHEM_AGAIN) THEN T1 = 0._EB - T2 = MAX_CHEM_TIME + !T2 = T2 !Set above ZETA0 = 0._EB ! Assume completely mixed. CVODE_CALL_OPTION = 1 ! CV_NORMAL CUR_CFD_TIME = T_CFD ! Set current cfd time in cvode, for logging purpose. - CALL CVODE_SERIAL(CC,ZZ, TMP_IN, PRES_IN, ZETA0, TAU_MIX, CELL_MASS, T1,T2, & + CALL CVODE_SERIAL(CC,ZZ, TMP_IN_MOD, TMP_IN, PRES_IN, ZETA0, TAU_MIX, CELL_MASS, T1,T2, & GLOBAL_ODE_REL_ERROR, ATOL, TMP_OUT, CHEM_TIME, WRITE_SUBSTEPS, CVODE_CALL_OPTION) - CHEM_TIME = MAX_CHEM_TIME + CHEM_TIME = T2 ELSE CC(1:N_TRACKED_SPECIES)=CC_CHEM_TIME(1:N_TRACKED_SPECIES) ENDIF TMP_IN_MOD = TMP_OUT - ZETA0 = ZETA_IN*EXP(-CHEM_TIME/TAU_MIX) - !WRITE(LU_ERR,*)"T_CFD,TMP_IN, TMP_OUT=",T_CFD,TMP_IN, TMP_OUT + ZETA0 = ZETA_IN_MOD*EXP(-CHEM_TIME/TAU_MIX) DT_MOD = DT_MOD - CHEM_TIME ELSE - ZETA0 = ZETA_IN + ZETA0 = ZETA_IN_MOD ENDIF CVODE_CALL_OPTION = 1 @@ -1071,18 +1142,97 @@ SUBROUTINE CVODE(ZZ, TMP_IN, PRES_IN, ZETA_IN, ZETA_OUT, TAU_MIX, CELL_MASS, T_ T1 = 0._EB T2 = DT_MOD CUR_CFD_TIME = T_CFD ! Set current cfd time in cvode, for logging purpose. -CALL CVODE_SERIAL(CC,ZZ, TMP_IN_MOD, PRES_IN, ZETA0, TAU_MIX, CELL_MASS, T1,T2, GLOBAL_ODE_REL_ERROR, ATOL, & +ZZ_IN = ZZ +CALL CVODE_SERIAL(CC,ZZ, TMP_IN_MOD, TMP_IN, PRES_IN, ZETA0, TAU_MIX, CELL_MASS, T1,T2, GLOBAL_ODE_REL_ERROR, ATOL, & TMP_OUT, CHEM_TIME, WRITE_SUBSTEPS, CVODE_CALL_OPTION) + + ! Convert back to mass fraction (Check for negative concentration) WHERE(CC<0._EB) CC=0._EB ZZ(1:N_TRACKED_SPECIES) = CC(1:N_TRACKED_SPECIES)*SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%MW RHO_OUT = SUM(ZZ) -ZZ = ZZ/RHO_OUT - -ZETA_OUT = ZETA_IN*EXP(-DT/TAU_MIX) +ZZ = ZZ/(RHO_OUT+TWO_EPSILON_EB) END SUBROUTINE CVODE + +!> \brief Constant pressure adiabatic flame temperature calculation +!> \param ZZ species mass fraction array +!> \param TMP_IN is the temperature +!> \param AFT is the adiabatic flame temperature (output) +SUBROUTINE CALC_ADIABATIC_FLAME_TEMPERATURE(ZZ,TMP_IN,AFT) +USE PHYSICAL_FUNCTIONS, ONLY: CALC_EQUIV_RATIO,GET_ENTHALPY,GET_SPECIFIC_HEAT,GET_TEMPERATURE +USE CHEMCONS, ONLY: I_FUEL +REAL(EB), INTENT(INOUT) :: ZZ(N_TRACKED_SPECIES) +REAL(EB), INTENT(IN) :: TMP_IN +REAL(EB), INTENT(OUT) :: AFT +REAL(EB) :: ZZ_REAC(N_TRACKED_SPECIES),ZZ_PROD(N_TRACKED_SPECIES) +REAL(EB) :: HS_IN + +IF (I_FUEL <= 0) THEN ! No FUEL_ID_FOR_AFT specified. + AFT = TMP_IN + RETURN +ENDIF +CALL CALC_AFT_REAC_AND_PROD(ZZ,ZZ_REAC,ZZ_PROD) +CALL GET_ENTHALPY(ZZ,HS_IN,TMP_IN) +AFT = TMP_IN +CALL GET_TEMPERATURE(AFT,HS_IN,ZZ_PROD) + +END SUBROUTINE CALC_ADIABATIC_FLAME_TEMPERATURE + +! Calculate Reactants and products +SUBROUTINE CALC_AFT_REAC_AND_PROD(ZZ,ZZ_REAC,ZZ_PROD) +USE PHYSICAL_FUNCTIONS, ONLY: CALC_EQUIV_RATIO +USE CHEMCONS, ONLY: I_FUEL,I_CO2,I_H2O,I_O2,I_N2 +REAL(EB), INTENT(IN) :: ZZ(N_TRACKED_SPECIES) +REAL(EB), INTENT(OUT) :: ZZ_REAC(N_TRACKED_SPECIES),ZZ_PROD(N_TRACKED_SPECIES) +REAL(EB) :: EQUIV, X,Y,Z,A,B,C,D,E, SUM_ZZ + + +CALL CALC_EQUIV_RATIO(ZZ(1:N_TRACKED_SPECIES), EQUIV) + +! Based on CxHyOz + a(O2+3.76N2) = bCO2 + cH2O + dCxHyOz + eO2 + 3.76aN2 +X=SPECIES_MIXTURE(I_FUEL)%ATOMS(6) !C +Y=SPECIES_MIXTURE(I_FUEL)%ATOMS(1) !H +Z=SPECIES_MIXTURE(I_FUEL)%ATOMS(8) !O +A=0.5_EB*((2.0_EB*X+0.5_EB*Y)/EQUIV -Z) +IF(ABS(EQUIV - 1.0_EB) <1E-3_EB) THEN ! Stoich + B = X + C = 0.5_EB*Y + D = 0._EB ! No fuel + E = 0._EB ! No O2 +ELSEIF (EQUIV > 1) THEN ! Rich + D = (2.0_EB*X+0.5_EB*Y-2.0_EB*A-Z)/(2.0_EB*X+0.5_EB*Y-Z) + B = (1._EB-D)*X + C = 0.5_EB*(1._EB-D)*Y + E = 0._EB ! No O2 +ELSE !EQUIV < 1, Lean + B = X + C = 0.5_EB*Y + D = 0._EB ! No fuel + E = A + 0.5*Z - B - 0.5_EB*C +ENDIF + +! Setup reactants +ZZ_REAC = 0.0_EB +ZZ_REAC(I_FUEL)=1._EB*SPECIES_MIXTURE(I_FUEL)%MW +ZZ_REAC(I_O2)=A*SPECIES_MIXTURE(I_O2)%MW +ZZ_REAC(I_N2)=3.76*A*SPECIES_MIXTURE(I_N2)%MW +SUM_ZZ=SUM(ZZ_REAC) +ZZ_REAC = ZZ_REAC/SUM_ZZ + +! Setup products +ZZ_PROD = 0.0_EB +ZZ_PROD(I_CO2)=B*SPECIES_MIXTURE(I_CO2)%MW +ZZ_PROD(I_H2O)=C*SPECIES_MIXTURE(I_H2O)%MW +ZZ_PROD(I_FUEL)=D*SPECIES_MIXTURE(I_FUEL)%MW +ZZ_PROD(I_O2)=E*SPECIES_MIXTURE(I_O2)%MW +ZZ_PROD(I_N2)=3.76_EB*A*SPECIES_MIXTURE(I_N2)%MW +SUM_ZZ=SUM(ZZ_PROD) +ZZ_PROD = ZZ_PROD/SUM_ZZ + + +END SUBROUTINE CALC_AFT_REAC_AND_PROD #endif SUBROUTINE CHECK_AUTO_IGNITION(EXTINCT,TMP_IN,AIT,IIC,JJC,KKC,REAC_INDEX) diff --git a/Source/func.f90 b/Source/func.f90 index 199cb5f6825..5459232ac2e 100644 --- a/Source/func.f90 +++ b/Source/func.f90 @@ -2187,8 +2187,10 @@ SUBROUTINE CALC_EQUIV_RATIO (ZZ, EQUIV) DENOM = DENOM + ZZ(NS)*SPECIES_MIXTURE(NS)%OXA ENDDO -IF (DENOM < TWO_EPSILON_EB) THEN +IF (NUMER <= 0.0_EB) THEN EQUIV= 0._EB; +ELSEIF (DENOM <= 0.0_EB) THEN + EQUIV= 1000._EB; ! Set a high value ELSE EQUIV = NUMER/ DENOM; ENDIF @@ -4130,6 +4132,9 @@ SUBROUTINE PACK_BOUNDARY_COORD(NM,IC,RC,OS,BC_INDEX,UNPACK_IT,COUNT_ONLY) IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%II,UNPACK_IT) IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%JJ,UNPACK_IT) IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%KK,UNPACK_IT) +IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%II2,UNPACK_IT) +IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%JJ2,UNPACK_IT) +IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%KK2,UNPACK_IT) IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%IIG,UNPACK_IT) IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%JJG,UNPACK_IT) IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),BC%KKG,UNPACK_IT) diff --git a/Source/geom.f90 b/Source/geom.f90 index 939d8727b5f..7742e62c009 100644 --- a/Source/geom.f90 +++ b/Source/geom.f90 @@ -1762,7 +1762,7 @@ SUBROUTINE SET_CUTCELLS_3D DO ICC=1,M%N_CUTCELL_MESH SUM_CC = SUM_CC + CC%NCELL ENDDO - ALLOCATE(M%CHEM_ACTIVE_CC(SUM_CC,2)) + ALLOCATE(M%CHEM_ACTIVE_CC(SUM_CC,3)) M%CHEM_ACTIVE_CC=-1 ENDDO MAIN_MESH_LOOP_4 diff --git a/Source/init.f90 b/Source/init.f90 index a9e56f9289d..aa4ddfb6ead 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -572,7 +572,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM) ALLOCATE(M%CHEM_SUBIT(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','CHEM_SUBIT',IZERO) ; M%CHEM_SUBIT=0._EB ENDIF -ALLOCATE(M%CHEM_ACTIVE_CELLS(IBP1*JBP1*KBP1,3),STAT=IZERO) ; CALL ChkMemErr('INIT','CHEM_ACTIVE_CELLS',IZERO) +ALLOCATE(M%CHEM_ACTIVE_CELLS(IBP1*JBP1*KBP1,4),STAT=IZERO) ; CALL ChkMemErr('INIT','CHEM_ACTIVE_CELLS',IZERO) M%CHEM_ACTIVE_CELLS=-1 ALLOCATE(M%Q(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','Q',IZERO) @@ -595,18 +595,6 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM) ALLOCATE(M%RSUM(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','RSUM',IZERO) ; M%RSUM=RSUM0 -! Allocate scalar face values - -ALLOCATE( M%FX(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FX',IZERO) ; M%FX=0._EB -ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FY',IZERO) ; M%FY=0._EB -ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FZ',IZERO) ; M%FZ=0._EB - -IF (FLUX_LIMITER_SINGLE_COEF) THEN - ALLOCATE( M%BFX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFX',IZERO) ; M%BFX=0._EB - ALLOCATE( M%BFY(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFY',IZERO) ; M%BFY=0._EB - ALLOCATE( M%BFZ(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFZ',IZERO) ; M%BFZ=0._EB -ENDIF - ! Allocate storage for scalar total fluxes IF (STORE_SPECIES_FLUX) THEN @@ -831,6 +819,23 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM) M%ZZS(:,:,:,N) = INITIAL_UNMIXED_FRACTION ENDDO +! Allocate scalar face values + +ALLOCATE( M%FX(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FX',IZERO) ; M%FX=0._EB +ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FY',IZERO) ; M%FY=0._EB +ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO); CALL ChkMemErr('INIT','FZ',IZERO) ; M%FZ=0._EB +DO N=1,N_TRACKED_SPECIES + M%FX(:,:,:,N)=RHOA*SPECIES_MIXTURE(N)%ZZ0 + M%FY(:,:,:,N)=RHOA*SPECIES_MIXTURE(N)%ZZ0 + M%FZ(:,:,:,N)=RHOA*SPECIES_MIXTURE(N)%ZZ0 +ENDDO + +! Allocate flux limiter coefficients + +ALLOCATE( M%BFX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFX',IZERO) ; M%BFX=0._EB +ALLOCATE( M%BFY(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFY',IZERO) ; M%BFY=0._EB +ALLOCATE( M%BFZ(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO); CALL ChkMemErr('INIT','BFZ',IZERO) ; M%BFZ=0._EB + ! Allocate and Initialize Mesh-Dependent Radiation Arrays M%QR = 0._EB diff --git a/Source/main.f90 b/Source/main.f90 index a8baada7bd1..49a5f7bb4ae 100644 --- a/Source/main.f90 +++ b/Source/main.f90 @@ -725,11 +725,7 @@ PROGRAM FDS COMPUTE_FINITE_DIFFERENCES_1: DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX CALL INSERT_ALL_PARTICLES(T,NM) IF (.NOT.SOLID_PHASE_ONLY .AND. .NOT.FREEZE_VELOCITY) CALL COMPUTE_VISCOSITY(NM,APPLY_TO_ESTIMATED_VARIABLES=.FALSE.) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - CALL MASS_FINITE_DIFFERENCES_NEW(NM) - ELSE - CALL MASS_FINITE_DIFFERENCES(NM) - ENDIF + CALL MASS_FINITE_DIFFERENCES_NEW(NM) ENDDO COMPUTE_FINITE_DIFFERENCES_1 ! Estimate quantities at next time step, and decrease/increase time step if necessary based on CFL condition @@ -910,11 +906,7 @@ PROGRAM FDS COMPUTE_FINITE_DIFFERENCES_2: DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX IF (.NOT.SOLID_PHASE_ONLY .AND. .NOT.FREEZE_VELOCITY) CALL COMPUTE_VISCOSITY(NM,APPLY_TO_ESTIMATED_VARIABLES=.TRUE.) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - CALL MASS_FINITE_DIFFERENCES_NEW(NM) - ELSE - CALL MASS_FINITE_DIFFERENCES(NM) - ENDIF + CALL MASS_FINITE_DIFFERENCES_NEW(NM) CALL DENSITY(T,DT,NM) IF (LEVEL_SET_MODE>0) CALL LEVEL_SET_FIRESPREAD(T,DT,NM) ENDDO COMPUTE_FINITE_DIFFERENCES_2 diff --git a/Source/mass.f90 b/Source/mass.f90 index dfb7dae2e4c..f92260f7728 100644 --- a/Source/mass.f90 +++ b/Source/mass.f90 @@ -9,7 +9,7 @@ MODULE MASS IMPLICIT NONE (TYPE,EXTERNAL) PRIVATE -PUBLIC MASS_FINITE_DIFFERENCES,MASS_FINITE_DIFFERENCES_NEW,DENSITY +PUBLIC MASS_FINITE_DIFFERENCES_NEW,DENSITY CONTAINS @@ -17,172 +17,6 @@ MODULE MASS !> \brief Compute spatial differences for mass transport equations !> \param NM Mesh index -SUBROUTINE MASS_FINITE_DIFFERENCES(NM) - -USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE -USE PHYSICAL_FUNCTIONS, ONLY: GET_MOLECULAR_WEIGHT - -INTEGER, INTENT(IN) :: NM -REAL(EB) :: TNOW -REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP -REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP -REAL(EB), PARAMETER :: DUMMY=0._EB -INTEGER :: I,J,K,N,IOR,IW,IIG,JJG,KKG,II,JJ,KK,IC -REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,RHOP -REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZP -REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P -TYPE(WALL_TYPE), POINTER :: WC -TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC -TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1 - -IF (SOLID_PHASE_ONLY) RETURN - -TNOW=CURRENT_TIME() -CALL POINT_TO_MESH(NM) - -IF (PREDICTOR) THEN - UU => U - VV => V - WW => W - RHOP => RHO - ZZP => ZZ -ELSE - UU => US - VV => VS - WW => WS - RHOP => RHOS - ZZP => ZZS -ENDIF - -! Reset counter for CLIP_RHOMIN, CLIP_RHOMAX -! Done here so DT_RESTRICT_COUNT will persist until WRITE_DIAGNOSTICS is called - -IF (PREDICTOR) DT_RESTRICT_COUNT = 0 - -! Species face values - -SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS - - RHO_Z_P=>WORK_PAD - - !$OMP PARALLEL DO - DO K=-1,KBP1+1 - DO J=-1,JBP1+1 - DO I=-1,IBP1+1 - RHO_Z_P(I,J,K) = RHOP(I,J,K)*ZZP(I,J,K,N) - ENDDO - ENDDO - ENDDO - !$OMP END PARALLEL DO - - ! Compute scalar face values - - CALL GET_SCALAR_FACE_VALUE(UU,RHO_Z_P,FX(:,:,:,N),0,IBAR,1,JBAR,1,KBAR,1,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE(VV,RHO_Z_P,FY(:,:,:,N),1,IBAR,0,JBAR,1,KBAR,2,I_FLUX_LIMITER) - CALL GET_SCALAR_FACE_VALUE(WW,RHO_Z_P,FZ(:,:,:,N),1,IBAR,1,JBAR,0,KBAR,3,I_FLUX_LIMITER) - - !$OMP PARALLEL DO PRIVATE(IW,WC,BC,B1,II,JJ,KK,IIG,JJG,KKG,IOR,IC,U_TEMP,Z_TEMP,F_TEMP) - WALL_LOOP_2: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS - WC=>WALL(IW) - IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_LOOP_2 - BC=>BOUNDARY_COORD(WC%BC_INDEX) - B1=>BOUNDARY_PROP1(WC%B1_INDEX) - - II = BC%II - JJ = BC%JJ - KK = BC%KK - IIG = BC%IIG - JJG = BC%JJG - KKG = BC%KKG - IOR = BC%IOR - IC = CELL_INDEX(II,JJ,KK) - - IF (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. .NOT.CELL(IC)%SOLID .AND. .NOT.CELL(IC)%EXTERIOR) THEN - SELECT CASE(IOR) - CASE( 1); FX(IIG-1,JJG,KKG,N) = 0._EB - CASE(-1); FX(IIG,JJG,KKG,N) = 0._EB - CASE( 2); FY(IIG,JJG-1,KKG,N) = 0._EB - CASE(-2); FY(IIG,JJG,KKG,N) = 0._EB - CASE( 3); FZ(IIG,JJG,KKG-1,N) = 0._EB - CASE(-3); FZ(IIG,JJG,KKG,N) = 0._EB - END SELECT - ELSE - SELECT CASE(IOR) - CASE( 1); FX(IIG-1,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-1); FX(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE( 2); FY(IIG,JJG-1,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-2); FY(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE( 3); FZ(IIG,JJG,KKG-1,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-3); FZ(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - END SELECT - ENDIF - - ! Overwrite first off-wall advective flux if flow is away from the wall and if the face is not also a wall cell - - OFF_WALL_IF_2: IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY .AND. WC%BOUNDARY_TYPE/=OPEN_BOUNDARY) THEN - - OFF_WALL_SELECT_2: SELECT CASE(IOR) - CASE( 1) OFF_WALL_SELECT_2 - ! ghost FX/UU(II+1) - ! /// II /// II+1 | II+2 | ... - ! ^ WALL_INDEX(II+1,+1) - IF ((UU(II+1,JJ,KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II+1,JJ,KK))%WALL_INDEX(+1)>0)) THEN - Z_TEMP(0:3,1,1) = (/RHO_Z_P(II+1,JJ,KK),RHO_Z_P(II+1:II+2,JJ,KK),DUMMY/) - U_TEMP(1,1,1) = UU(II+1,JJ,KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - FX(II+1,JJ,KK,N) = F_TEMP(1,1,1) - ENDIF - CASE(-1) OFF_WALL_SELECT_2 - ! FX/UU(II-2) ghost - ! ... | II-2 | II-1 /// II /// - ! ^ WALL_INDEX(II-1,-1) - IF ((UU(II-2,JJ,KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II-1,JJ,KK))%WALL_INDEX(-1)>0)) THEN - Z_TEMP(0:3,1,1) = (/DUMMY,RHO_Z_P(II-2:II-1,JJ,KK),RHO_Z_P(II-1,JJ,KK)/) - U_TEMP(1,1,1) = UU(II-2,JJ,KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - FX(II-2,JJ,KK,N) = F_TEMP(1,1,1) - ENDIF - CASE( 2) OFF_WALL_SELECT_2 - IF ((VV(II,JJ+1,KK)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ+1,KK))%WALL_INDEX(+2)>0)) THEN - Z_TEMP(1,0:3,1) = (/RHO_Z_P(II,JJ+1,KK),RHO_Z_P(II,JJ+1:JJ+2,KK),DUMMY/) - U_TEMP(1,1,1) = VV(II,JJ+1,KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - FY(II,JJ+1,KK,N) = F_TEMP(1,1,1) - ENDIF - CASE(-2) OFF_WALL_SELECT_2 - IF ((VV(II,JJ-2,KK)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ-1,KK))%WALL_INDEX(-2)>0)) THEN - Z_TEMP(1,0:3,1) = (/DUMMY,RHO_Z_P(II,JJ-2:JJ-1,KK),RHO_Z_P(II,JJ-1,KK)/) - U_TEMP(1,1,1) = VV(II,JJ-2,KK) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - FY(II,JJ-2,KK,N) = F_TEMP(1,1,1) - ENDIF - CASE( 3) OFF_WALL_SELECT_2 - IF ((WW(II,JJ,KK+1)>0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ,KK+1))%WALL_INDEX(+3)>0)) THEN - Z_TEMP(1,1,0:3) = (/RHO_Z_P(II,JJ,KK+1),RHO_Z_P(II,JJ,KK+1:KK+2),DUMMY/) - U_TEMP(1,1,1) = WW(II,JJ,KK+1) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - FZ(II,JJ,KK+1,N) = F_TEMP(1,1,1) - ENDIF - CASE(-3) OFF_WALL_SELECT_2 - IF ((WW(II,JJ,KK-2)<0._EB) .AND. .NOT.(CELL(CELL_INDEX(II,JJ,KK-1))%WALL_INDEX(-3)>0)) THEN - Z_TEMP(1,1,0:3) = (/DUMMY,RHO_Z_P(II,JJ,KK-2:KK-1),RHO_Z_P(II,JJ,KK-1)/) - U_TEMP(1,1,1) = WW(II,JJ,KK-2) - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - FZ(II,JJ,KK-2,N) = F_TEMP(1,1,1) - ENDIF - END SELECT OFF_WALL_SELECT_2 - - ENDIF OFF_WALL_IF_2 - - ENDDO WALL_LOOP_2 - !$OMP END PARALLEL DO - -ENDDO SPECIES_LOOP - -T_USED(3)=T_USED(3)+CURRENT_TIME()-TNOW -END SUBROUTINE MASS_FINITE_DIFFERENCES - - SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE_NEW,GET_SCALAR_FACE_COEF @@ -294,6 +128,7 @@ SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) IC = CELL_INDEX(II,JJ,KK) IF (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. .NOT.CELL(IC)%SOLID .AND. .NOT.CELL(IC)%EXTERIOR) THEN + ! thin obstruction SELECT CASE(IOR) CASE( 1); FX(IIG-1,JJG,KKG,N) = 0._EB CASE(-1); FX(IIG,JJG,KKG,N) = 0._EB @@ -303,14 +138,16 @@ SUBROUTINE MASS_FINITE_DIFFERENCES_NEW(NM) CASE(-3); FZ(IIG,JJG,KKG,N) = 0._EB END SELECT ELSE - SELECT CASE(IOR) - CASE( 1); FX(IIG-1,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-1); FX(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE( 2); FY(IIG,JJG-1,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-2); FY(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - CASE( 3); FZ(IIG,JJG,KKG-1,N) = B1%RHO_F*B1%ZZ_F(N) - CASE(-3); FZ(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) - END SELECT + IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) THEN + SELECT CASE(IOR) + CASE( 1); FX(IIG-1,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) + CASE(-1); FX(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) + CASE( 2); FY(IIG,JJG-1,KKG,N) = B1%RHO_F*B1%ZZ_F(N) + CASE(-2); FY(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) + CASE( 3); FZ(IIG,JJG,KKG-1,N) = B1%RHO_F*B1%ZZ_F(N) + CASE(-3); FZ(IIG,JJG,KKG,N) = B1%RHO_F*B1%ZZ_F(N) + END SELECT + ENDIF ENDIF ! Overwrite first off-wall advective flux if flow is away from the wall and if the face is not also a wall cell diff --git a/Source/mesh.f90 b/Source/mesh.f90 index ab04b28ffe5..2daabf458cf 100644 --- a/Source/mesh.f90 +++ b/Source/mesh.f90 @@ -53,8 +53,8 @@ MODULE MESH_VARIABLES REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: D_SOURCE!< Source terms in the expression for the divergence REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: CSD2 !< \f$ C_s \Delta^2 \f$ in Smagorinsky turbulence expression REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: CHEM_SUBIT !< Number of chemistry sub-iterations - INTEGER, ALLOCATABLE, DIMENSION(:,:) :: CHEM_ACTIVE_CELLS !< I , J ,K info of chemically active cells. - INTEGER, ALLOCATABLE, DIMENSION(:,:) :: CHEM_ACTIVE_CC !< ICC, JCC of chemically active cells. + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: CHEM_ACTIVE_CELLS !< I , J ,K, igntion_zone info of chemically active cells. + INTEGER, ALLOCATABLE, DIMENSION(:,:) :: CHEM_ACTIVE_CC !< ICC, JCC, igntion_zone of chemically active cells. REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: MIX_TIME !< Mixing-controlled combustion reaction time (s) REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: STRAIN_RATE !< Strain rate \f$ |S|_{ijk} \f$ (1/s) REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: D_Z_MAX !< \f$ \max D_\alpha \f$ diff --git a/Source/pois.f90 b/Source/pois.f90 index e51f035cb6b..0ae9d7a710d 100644 --- a/Source/pois.f90 +++ b/Source/pois.f90 @@ -4511,7 +4511,7 @@ SUBROUTINE SSWAP(N,SX,INCX,SY,INCY) INTEGER :: N INTEGER :: INCX INTEGER :: INCY -REAL(EB) SX(1),SY(1),STEMP1,STEMP2,STEMP3 +REAL(EB) SX(*),SY(*),STEMP1,STEMP2,STEMP3 INTEGER :: I, IX, IY, M, MP1, NS !***FIRST EXECUTABLE STATEMENT SSWAP diff --git a/Source/prec.f90 b/Source/prec.f90 index 3f2ab1cb9ca..18cadc2bd51 100644 --- a/Source/prec.f90 +++ b/Source/prec.f90 @@ -24,6 +24,7 @@ MODULE PRECISION_PARAMETERS INTEGER, PARAMETER :: MAX_INPUT_ID=40 !< Maximum number of CTRL INPUT_IDs INTEGER, PARAMETER :: N_ZONE_POINTS=100 !< Maximum number of declared ZONE points (deprecated) INTEGER, PARAMETER :: MAX_AIT_EXCLUSION_ZONES=10 !< Maximum number of AUTO_IGNITION_TEMPERATURE exclusion zones +INTEGER, PARAMETER :: MAX_IGNITION_ZONES=10 !< Maximum number of Ignition zones INTEGER, PARAMETER :: SMOKEVIEW_OBJECTS_DIMENSION=20 !< Number of parameters that can be passed to Smokeview to describe objects INTEGER, PARAMETER :: LABEL_LENGTH=60 !< Maximum length of most labels INTEGER, PARAMETER :: MESSAGE_LENGTH=200 !< Maximum length of error and warning labels diff --git a/Source/radi.f90 b/Source/radi.f90 index eb14f65d8da..a342e288794 100644 --- a/Source/radi.f90 +++ b/Source/radi.f90 @@ -3291,10 +3291,10 @@ SUBROUTINE INIT_RADIATION PARTIAL_PRESSURES_ATM(14,1) = 1._EB END SELECT CALL SUB_RADCAL(AMEAN,AP0,RADIANCE,TRANSMISSIVITY) - IF (NSB==1 .AND. PATH_LENGTH > 0.0_EB) THEN - RADCAL_SPECIES2KAPPA(NS,J,K,1) = MIN(AMEAN,AP0) - ELSE - RADCAL_SPECIES2KAPPA(NS,J,K,IBND) = AMEAN/BBF + IF (PATH_LENGTH > 0.0_EB) THEN + RADCAL_SPECIES2KAPPA(NS,J,K,1) = MIN(AMEAN,AP0)/BBF + ELSE ! zero path length + RADCAL_SPECIES2KAPPA(NS,J,K,1) = AP0/BBF ENDIF END DO RADCAL_SPECIES_LOOP ENDDO Y_LOOP_Z diff --git a/Source/read.f90 b/Source/read.f90 index f792194b30a..4873f60ef98 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -1743,7 +1743,7 @@ SUBROUTINE READ_MISC CFL_MAX,CFL_MIN,CFL_VELOCITY_NORM,CHECK_FO,CHECK_HT,CHECK_VN, & CNF_CUTOFF,CONSTANT_SPECIFIC_HEAT_RATIO,& C_SMAGORINSKY,C_VREMAN,C_WALE,DEPOSITION,EXTERNAL_FILENAME,& - FIXED_LES_FILTER_WIDTH,FLUX_LIMITER,FLUX_LIMITER_SINGLE_COEF,& + FIXED_LES_FILTER_WIDTH,FLUX_LIMITER,& FREEZE_VELOCITY,FYI,GAMMA,GRAVITATIONAL_DEPOSITION,& GRAVITATIONAL_SETTLING,GVEC,H_F_REFERENCE_TEMPERATURE,& HUMIDITY,HVAC_LOCAL_PRESSURE,HVAC_MASS_TRANSPORT_CELL_L,HVAC_PRES_RELAX,HVAC_QFAN,IBLANK_SMV,I_MAX_TEMP,& @@ -3548,8 +3548,21 @@ SUBROUTINE READ_SPEC ENDDO DO N1 = 1, N_TRACKED_SPECIES - SPECIES_MIXTURE(N1)%OXR = 2*(SPECIES_MIXTURE(N1)%ATOMS(1) * 0.25_EB + SPECIES_MIXTURE(N1)%ATOMS(6)) & - * ELEMENT(8)%MASS / SPECIES_MIXTURE(N1)%MW + SPECIES_MIXTURE(N1)%OXR = (SPECIES_MIXTURE(N1)%ATOMS(1) * 0.5_EB + & ! H (H2O) + SPECIES_MIXTURE(N1)%ATOMS(6) * 2.0_EB + & ! C (CO2) + SPECIES_MIXTURE(N1)%ATOMS(16) * 2.0_EB + & ! S (SO2) + SPECIES_MIXTURE(N1)%ATOMS(15) * 2.5_EB + & ! P (P4O10) + SPECIES_MIXTURE(N1)%ATOMS(26) * 1.5_EB + & ! Fe (Fe2O3) + SPECIES_MIXTURE(N1)%ATOMS(13) * 1.5_EB + & ! Al (Al2O3) + SPECIES_MIXTURE(N1)%ATOMS(14) * 2.0_EB + & ! Si (SiO2) + SPECIES_MIXTURE(N1)%ATOMS(22) * 2.0_EB + & ! Ti (TiO2) + SPECIES_MIXTURE(N1)%ATOMS(29) * 1.0_EB + & ! Cu (CuO) + SPECIES_MIXTURE(N1)%ATOMS(30) * 1.0_EB + & ! Zn (ZnO) + SPECIES_MIXTURE(N1)%ATOMS(11) * 0.5_EB + & ! Na (Na2O) + SPECIES_MIXTURE(N1)%ATOMS(19) * 0.5_EB + & ! K (K2O) + SPECIES_MIXTURE(N1)%ATOMS(20) * 1.0_EB ) & ! Ca (CaO) + * ELEMENT(8)%MASS / SPECIES_MIXTURE(N1)%MW + SPECIES_MIXTURE(N1)%OXA = SPECIES_MIXTURE(N1)%ATOMS(8) * ELEMENT(8)%MASS / SPECIES_MIXTURE(N1)%MW ENDDO @@ -4607,14 +4620,15 @@ END SUBROUTINE PROC_SPEC_2 !> \brief Read the COMBustion namelist line SUBROUTINE READ_COMB -USE CHEMCONS, ONLY: ODE_MIN_ATOL, EQUIV_RATIO_CHECK, MIN_EQUIV_RATIO, MAX_EQUIV_RATIO, DO_CHEM_LOAD_BALANCE,FLAME_THICK_FACTOR +USE CHEMCONS, ONLY: ODE_MIN_ATOL, EQUIV_RATIO_CHECK, MIN_EQUIV_RATIO, MAX_EQUIV_RATIO, DO_CHEM_LOAD_BALANCE, & + FUEL_ID_FOR_AFT,USE_MIXED_ZN_AFT_TMP REAL(EB) :: ODE_REL_ERROR CHARACTER(LABEL_LENGTH) :: RAMP_ZETA_0='null' NAMELIST /COMB/ CHECK_REALIZABILITY,COMPUTE_ADIABATIC_FLAME_TEMPERATURE, DO_CHEM_LOAD_BALANCE, EQUIV_RATIO_CHECK, & - EXTINCTION_MODEL,FINITE_RATE_MIN_TEMP, FIXED_MIX_TIME,FLAME_THICK_FACTOR,FREE_BURN_TEMPERATURE,& - INITIAL_UNMIXED_FRACTION, MAX_CHEMISTRY_SUBSTEPS, MAX_EQUIV_RATIO, MIN_EQUIV_RATIO, & + EXTINCTION_MODEL,FINITE_RATE_MIN_TEMP, FIXED_MIX_TIME,FREE_BURN_TEMPERATURE, & + FUEL_ID_FOR_AFT,INITIAL_UNMIXED_FRACTION, MAX_CHEMISTRY_SUBSTEPS, MAX_EQUIV_RATIO, MIN_EQUIV_RATIO, & N_FIXED_CHEMISTRY_SUBSTEPS, ODE_MIN_ATOL,ODE_REL_ERROR,ODE_SOLVER,SUPPRESSION,TAU_CHEM, & - TAU_FLAME,RAMP_ZETA_0,ZZ_MIN_GLOBAL + TAU_FLAME,RAMP_ZETA_0,USE_MIXED_ZN_AFT_TMP,ZZ_MIN_GLOBAL ODE_SOLVER = 'null' ODE_REL_ERROR = -1._EB @@ -4644,6 +4658,13 @@ SUBROUTINE READ_COMB 26 IF (IOS>0) THEN ; CALL SHUTDOWN('ERROR(101): Problem with COMB line.') ; RETURN ; ENDIF ENDDO COMB_LOOP2 25 REWIND(LU_INPUT) ; INPUT_FILE_LINE_NUMBER = 0 + IF(USE_MIXED_ZN_AFT_TMP) THEN + IF (TRIM(FUEL_ID_FOR_AFT) == 'null') THEN + WRITE(MESSAGE,'(A)') 'ERROR(*): FUEL_ID_FOR_AFT must be specified when USE_MIXED_ZN_AFT_TMP is true (default).' & + // NEW_LINE('A') // 'Either set USE_MIXED_ZN_AFT_TMP=FALSE, or provide FUEL_ID_FOR_AFT in the COMB line.' + CALL SHUTDOWN(MESSAGE) ; RETURN + ENDIF + ENDIF ENDIF ! Extinction Model @@ -4668,6 +4689,8 @@ SUBROUTINE READ_COMB ENDIF ENDIF + + FINITE_RATE_MIN_TEMP = FINITE_RATE_MIN_TEMP + TMPM ! Convert C to K for EXTINCTION 1 temperature cut-off @@ -5183,7 +5206,7 @@ END SUBROUTINE READ_REAC SUBROUTINE PROC_REAC_1 -USE CHEMCONS, ONLY: YP2ZZ, IS_EXPONENT_LT_1 +USE CHEMCONS, ONLY: YP2ZZ,IS_EXPONENT_LT_1,FUEL_ID_FOR_AFT,I_FUEL,I_CO2,I_H2O,I_O2,I_N2 USE PROPERTY_DATA, ONLY : PARSE_EQUATION, SHUTDOWN_ATOM REAL(EB) :: MASS_PRODUCT,MASS_REACTANT,REACTION_BALANCE(118),NU_Y(N_SPECIES) INTEGER :: NS,NS2, NR @@ -5621,6 +5644,30 @@ SUBROUTINE PROC_REAC_1 DZZ_CLIP = ZZ_MIN_GLOBAL * GLOBAL_ODE_REL_ERROR ENDIF +! Compute the index of key species for adiabatic flame temperature calculation. +IF(TRIM(ODE_SOLVER)=='CVODE') THEN + I_FUEL = -1 + I_O2 = -1 + I_N2 = -1 + I_CO2 = -1 + I_H2O = -1 + + DO NS=1,N_TRACKED_SPECIES + IF(TRIM(FUEL_ID_FOR_AFT)==TRIM(SPECIES_MIXTURE(NS)%ID) .OR. TRIM(FUEL_ID_FOR_AFT)==TRIM(SPECIES_MIXTURE(NS)%ALT_ID)) THEN + I_FUEL = NS + ELSEIF(TRIM(SPECIES_MIXTURE(NS)%ID) == "O2" .OR. TRIM(SPECIES_MIXTURE(NS)%ID) == "OXYGEN") THEN + I_O2 = NS + ELSEIF(TRIM(SPECIES_MIXTURE(NS)%ID) == "N2" .OR. TRIM(SPECIES_MIXTURE(NS)%ID) == "NITROGEN") THEN + I_N2 = NS + ELSEIF(TRIM(SPECIES_MIXTURE(NS)%ID) == "CO2" .OR. TRIM(SPECIES_MIXTURE(NS)%ID) == "CARBON DIOXIDE") THEN + I_CO2 = NS + ELSEIF(TRIM(SPECIES_MIXTURE(NS)%ID) == "H2O" .OR. TRIM(SPECIES_MIXTURE(NS)%ID) == "WATER VAPOR") THEN + I_H2O = NS + ENDIF + ENDDO + +ENDIF + END SUBROUTINE PROC_REAC_1 @@ -6600,24 +6647,25 @@ SUBROUTINE READ_PROP C_FACTOR,CHARACTERISTIC_VELOCITY,ORIFICE_DIAMETER,EMISSIVITY, & PARTICLE_VELOCITY,FLOW_RATE,FLOW_TAU,GAUGE_EMISSIVITY,GAUGE_TEMPERATURE,INITIAL_TEMPERATURE,K_FACTOR,& LENGTH,SPRAY_ANGLE(2,2),OFFSET,OPERATING_PRESSURE,RTI,PDPA_START,PDPA_END,PDPA_RADIUS,MASS_FLOW_RATE,& - SPRAY_PATTERN_MU,SPRAY_PATTERN_BETA,HISTOGRAM_LIMITS(2),P0,PX(3),PXX(3,3),TIME_CONSTANT,VIEW_ANGLE,PROBE_DIAMETER + SPRAY_PATTERN_MU,SPRAY_PATTERN_BETA,HISTOGRAM_LIMITS(2),P0,PX(3),PXX(3,3),TIME_CONSTANT,VIEW_ANGLE, & + PROBE_DIAMETER INTEGER ::I,N,NN,PDPA_M,PDPA_N,PARTICLES_PER_SECOND,VELOCITY_COMPONENT,HISTOGRAM_NBINS,FED_ACTIVITY -LOGICAL :: PDPA_INTEGRATE,PDPA_NORMALIZE,HISTOGRAM_NORMALIZE,HISTOGRAM,HISTOGRAM_CUMULATIVE,SPARK,TC +LOGICAL :: PDPA_INTEGRATE,PDPA_NORMALIZE,HISTOGRAM_NORMALIZE,HISTOGRAM,HISTOGRAM_CUMULATIVE,SPARK,TC,IGNITION_ZONE CHARACTER(LABEL_LENGTH) :: SMOKEVIEW_ID(SMOKEVIEW_OBJECTS_DIMENSION),QUANTITY='null',PART_ID='null',FLOW_RAMP='null', & SPRAY_PATTERN_TABLE='null',SPEC_ID='null',& PRESSURE_RAMP='null',SMOKEVIEW_PARAMETERS(SMOKEVIEW_OBJECTS_DIMENSION), & - SPRAY_PATTERN_SHAPE='GAUSSIAN' + SPRAY_PATTERN_SHAPE='GAUSSIAN',SPECIFIC_HEAT_RAMP TYPE (PROPERTY_TYPE), POINTER :: PY NAMELIST /PROP/ ACTIVATION_OBSCURATION,ACTIVATION_TEMPERATURE,ALPHA_C,ALPHA_E,BETA_C,BETA_E,CALIBRATION_CONSTANT,& CHARACTERISTIC_VELOCITY,C_FACTOR,DENSITY,DIAMETER,EMISSIVITY,FED_ACTIVITY,FLOW_RAMP,FLOW_RATE,FLOW_TAU, & GAUGE_EMISSIVITY,GAUGE_TEMPERATURE,HEAT_TRANSFER_COEFFICIENT,HISTOGRAM,HISTOGRAM_CUMULATIVE, & - HISTOGRAM_LIMITS,HISTOGRAM_NBINS,HISTOGRAM_NORMALIZE,ID, & + HISTOGRAM_LIMITS,HISTOGRAM_NBINS,HISTOGRAM_NORMALIZE,ID, IGNITION_ZONE, & INITIAL_TEMPERATURE,K_FACTOR,LENGTH,MASS_FLOW_RATE,OFFSET,OPERATING_PRESSURE,ORIFICE_DIAMETER,P0,& PARTICLES_PER_SECOND,PARTICLE_VELOCITY,PART_ID,PDPA_END,& PDPA_INTEGRATE,PDPA_M,PDPA_N,PDPA_NORMALIZE,PDPA_RADIUS,& PDPA_START,PRESSURE_RAMP,PROBE_DIAMETER,PX,PXX,QUANTITY,RTI,SMOKEVIEW_ID,SMOKEVIEW_PARAMETERS,SPARK,& - SPEC_ID,SPECIFIC_HEAT,SPRAY_ANGLE,& + SPEC_ID,SPECIFIC_HEAT,SPECIFIC_HEAT_RAMP,SPRAY_ANGLE,& SPRAY_PATTERN_BETA,SPRAY_PATTERN_MU,SPRAY_PATTERN_SHAPE,SPRAY_PATTERN_TABLE,TC,TIME_CONSTANT,VELOCITY_COMPONENT,& VIEW_ANGLE @@ -6666,7 +6714,17 @@ SUBROUTINE READ_PROP PY%DIAMETER = DIAMETER PY%EMISSIVITY = EMISSIVITY PY%HEAT_TRANSFER_COEFFICIENT= HEAT_TRANSFER_COEFFICIENT - PY%SPECIFIC_HEAT = SPECIFIC_HEAT*1000._EB/TIME_SHRINK_FACTOR + ALLOCATE(PY%SPECIFIC_HEAT(0:I_MAX_TEMP)) + IF(SPECIFIC_HEAT > 0._EB) THEN + PY%SPECIFIC_HEAT = SPECIFIC_HEAT*1000._EB + ELSE + ! Type-K CP(20 C)=0.4515 kJ/kg, CP(1200 C)=0.6010 kJ/kg + DO J = 0,I_MAX_TEMP + PY%SPECIFIC_HEAT(J) = 451.5_EB+0.001_EB*(REAL(MAX(293,MIN(1493,J)),EB)-293._EB)*(601.0_EB-451.5_EB) + ENDDO + ENDIF + PY%SPECIFIC_HEAT = PY%SPECIFIC_HEAT/TIME_SHRINK_FACTOR + IF (SPECIFIC_HEAT_RAMP /= 'null') CALL GET_RAMP_INDEX(SPECIFIC_HEAT_RAMP,'TEMPERATURE',PY%SPECIFIC_HEAT_RAMP_INDEX) PY%C_FACTOR = C_FACTOR PY%CHARACTERISTIC_VELOCITY = CHARACTERISTIC_VELOCITY PY%GAUGE_EMISSIVITY = GAUGE_EMISSIVITY @@ -6710,6 +6768,7 @@ SUBROUTINE READ_PROP IF (PY%SMOKEVIEW_PARAMETERS(I)/='null') PY%N_SMOKEVIEW_PARAMETERS = PY%N_SMOKEVIEW_PARAMETERS + 1 ENDDO PY%SPARK = SPARK + PY%IGNITION_ZONE = IGNITION_ZONE PY%SPEC_ID = SPEC_ID IF (PART_ID/='null' .AND. SPRAY_PATTERN_TABLE /= 'null') THEN CALL GET_TABLE_INDEX(SPRAY_PATTERN_TABLE,SPRAY_PATTERN,PY%SPRAY_PATTERN_INDEX) @@ -6907,11 +6966,12 @@ SUBROUTINE SET_PROP_DEFAULTS BETA_C = -1.0_EB BETA_E = -1.0_EB CALIBRATION_CONSTANT = 0.93_EB -DENSITY = 8908._EB ! kg/m3 (Nickel) +DENSITY = 8700._EB ! kg/m3 (Type-K) DIAMETER = 0.001 ! m EMISSIVITY = 0.85_EB HEAT_TRANSFER_COEFFICIENT= -1._EB ! W/m2/K -SPECIFIC_HEAT = 0.44_EB ! kJ/kg/K (Nickel) +SPECIFIC_HEAT = -1._EB ! kJ/kg/K +SPECIFIC_HEAT_RAMP = 'null' C_FACTOR = 0.0_EB CHARACTERISTIC_VELOCITY = 1.0_EB ! m/s PARTICLE_VELOCITY = 0._EB ! m/s @@ -6954,6 +7014,7 @@ SUBROUTINE SET_PROP_DEFAULTS SMOKEVIEW_ID = 'null' SMOKEVIEW_PARAMETERS = 'null' SPARK = .FALSE. +IGNITION_ZONE = .FALSE. SPEC_ID = 'null' SPRAY_ANGLE(1,1) = 60._EB ! degrees SPRAY_ANGLE(2,1) = 75._EB ! degrees @@ -6977,6 +7038,7 @@ END SUBROUTINE READ_PROP SUBROUTINE PROC_PROP USE DEVICE_VARIABLES +USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP REAL(EB) :: TOTAL_FLOWRATE, SUBTOTAL_FLOWRATE INTEGER :: N,NN,N_V_FACTORS,ILPC LOGICAL :: TABLE_NORMED(1:N_TABLE) @@ -7075,7 +7137,21 @@ SUBROUTINE PROC_PROP PY%V_FACTOR = PY%PARTICLE_VELOCITY/SQRT(PY%OPERATING_PRESSURE) ENDIF ENDIF + + ! Check units of specific heat + IF (PY%SPECIFIC_HEAT_RAMP_INDEX > 0) THEN + IF (.NOT.RAMPS(PY%SPECIFIC_HEAT_RAMP_INDEX)%DEP_VAR_UNITS_CONVERTED) THEN + RAMPS(PY%SPECIFIC_HEAT_RAMP_INDEX)%INTERPOLATED_DATA(:) = & + RAMPS(PY%SPECIFIC_HEAT_RAMP_INDEX)%INTERPOLATED_DATA(:)*1000._EB/TIME_SHRINK_FACTOR + RAMPS(PY%SPECIFIC_HEAT_RAMP_INDEX)%DEP_VAR_UNITS_CONVERTED = .TRUE. + ENDIF + IF (RAMPS(PY%SPECIFIC_HEAT_RAMP_INDEX)%DEPENDENT_DATA(1) > 10._EB .AND. MY_RANK==0) & + WRITE(LU_ERR,'(A,A)') 'WARNING: SPECIFIC_HEAT units are kJ/kg/K check PROP ',TRIM(PY%ID) + DO I=0,I_MAX_TEMP + PY%SPECIFIC_HEAT(I)=EVALUATE_RAMP(REAL(I,EB),PY%SPECIFIC_HEAT_RAMP_INDEX) + ENDDO + ENDIF ENDDO PROP_LOOP END SUBROUTINE PROC_PROP @@ -10401,6 +10477,8 @@ SUBROUTINE READ_RAMP CTRL_ID = 'null' X = -1.E6_EB Z = -1.E6_EB + T = 0._EB + F = 1._EB CALL CHECKREAD('RAMP',LU_INPUT,IOS) ; IF (STOP_STATUS==SETUP_STOP) RETURN IF (IOS==1) EXIT SEARCH_LOOP2 READ(LU_INPUT,RAMP) @@ -13587,7 +13665,7 @@ END SUBROUTINE READ_ZONE SUBROUTINE READ_DEVC USE DEVICE_VARIABLES, ONLY: DEVICE_TYPE,SUBDEVICE_TYPE,DEVICE,N_DEVC,N_DEVC_TIME,N_DEVC_LINE,MAX_DEVC_LINE_POINTS,& - DEVC_PIPE_OPERATING + DEVC_PIPE_OPERATING, PROPERTY USE GEOMETRY_FUNCTIONS, ONLY: TRANSFORM_COORDINATES,SEARCH_OTHER_MESHES INTEGER :: N,NN,NM,MESH_NUMBER=0,N_DEVC_READ,IOR,TRIP_DIRECTION,VELO_INDEX,POINTS,I_POINT,PIPE_INDEX,ORIENTATION_INDEX, & N_INTERVALS,MOVE_INDEX,IIG,JJG,KKG,NOM,ERROR_CODE,LP_TAG @@ -13612,6 +13690,7 @@ SUBROUTINE READ_DEVC TIME_AVERAGED,TIME_HISTORY,TIME_PERIOD,TRIP_DIRECTION,UNITS,VELO_INDEX,XB,XBP,XYZ,X_ID,Y_ID,Z_ID,XYZ_UNITS INTEGER, ALLOCATABLE, DIMENSION(:,:) :: MESH_DEVICE_ARRAY INTEGER, ALLOCATABLE, DIMENSION(:) :: MESH_DEVICE +INTEGER :: PROP_INDEX ! Read the input file and count the number of DEVC lines @@ -13619,6 +13698,7 @@ SUBROUTINE READ_DEVC N_DEVC_READ = 0 N_DEVC_TIME = 0 N_DEVC_LINE = 0 +PROP_INDEX = -1 REWIND(LU_INPUT) ; INPUT_FILE_LINE_NUMBER = 0 COUNT_DEVC_LOOP: DO @@ -13726,6 +13806,11 @@ SUBROUTINE READ_DEVC .OR. QUANTITY=='LAYER HEIGHT' .OR. QUANTITY=='UPPER TEMPERATURE' .OR. QUANTITY=='LOWER TEMPERATURE') XB_XYZ= .FALSE. IF (CTRL_ID/='null' .OR. DEVC_ID/='null' .OR. DUCT_ID/='null' .OR. NODE_ID(1)/='null' .OR. INIT_ID/='null') XB_XYZ= .FALSE. IF (SPATIAL_STATISTIC/='null' .AND. SPATIAL_STATISTIC/='INTERPOLATION') XB_XYZ= .FALSE. + ! Check for ignition zone + IF (PROP_ID/='null') THEN + CALL GET_PROPERTY_INDEX(PROP_INDEX,'DEVC',PROP_ID) + IF (PROPERTY(PROP_INDEX)%IGNITION_ZONE) XB_XYZ= .FALSE. + ENDIF IF (POINTS==1 .AND. SPATIAL_STATISTIC=='null' .AND. ALL(XB>-1E6_EB) .AND. ALL(XYZ<=-1E6_EB) .AND. XB_XYZ) THEN XYZ(1) = 0.5_EB*(XB(1)+XB(2)) XYZ(2) = 0.5_EB*(XB(3)+XB(4)) @@ -14215,8 +14300,7 @@ SUBROUTINE READ_DEVC DV%QUANTITY(1)=='GAUGE HEAT FLUX GAS' .OR. & DV%QUANTITY(1)=='RADIANCE' .OR. & DV%QUANTITY(1)=='ADIABATIC SURFACE TEMPERATURE GAS' .OR. & - DV%QUANTITY(1)=='RADIOMETER GAS' .OR. & - DV%QUANTITY(1)=='BI-DIRECTIONAL PROBE') THEN + DV%QUANTITY(1)=='RADIOMETER GAS') THEN IF (DV%ORIENTATION_INDEX==0) THEN WRITE(MESSAGE,'(3A)') 'ERROR(887): DEVC ',TRIM(ID),' must have an ORIENTATION.' CALL SHUTDOWN(MESSAGE) ; RETURN @@ -14229,6 +14313,12 @@ SUBROUTINE READ_DEVC INIT_RESERVED(N_INIT_RESERVED)%DEVC_INDEX = N_DEVC INIT_RESERVED(N_INIT_RESERVED)%N_PARTICLES = POINTS ENDIF + ENDIF + + ! Special case for QUANTITY='BI-DIRECTIONAL PROBE' + IF (DV%QUANTITY(1)=='BI-DIRECTIONAL PROBE' .AND. DV%ORIENTATION_INDEX==0) THEN + WRITE(MESSAGE,'(3A)') 'ERROR(887): DEVC ',TRIM(ID),' must have an ORIENTATION.' + CALL SHUTDOWN(MESSAGE) ; RETURN ENDIF ! Miscellaneous actions taken based on specific device attributes @@ -14915,6 +15005,7 @@ END SUBROUTINE PROC_OBST SUBROUTINE PROC_DEVC USE COMP_FUNCTIONS, ONLY: CHANGE_UNITS +USE CHEMCONS, ONLY: N_IGNITION_ZONES,IGNITION_ZONES,FUEL_ID_FOR_AFT USE GEOMETRY_FUNCTIONS, ONLY: SEARCH_OTHER_MESHES USE PHYSICAL_FUNCTIONS, ONLY: GET_VISCOSITY,GET_CONDUCTIVITY USE CONTROL_VARIABLES @@ -14986,6 +15077,40 @@ SUBROUTINE PROC_DEVC ENDDO ENDIF + IF (PROPERTY(DV%PROP_INDEX)%IGNITION_ZONE) THEN + N_IGNITION_ZONES = N_IGNITION_ZONES + 1 + IF (N_IGNITION_ZONES > MAX_IGNITION_ZONES) THEN + WRITE(MESSAGE,'(3A,I3,A)') 'ERROR: DEVC ',TRIM(DV%ID),': Total ', MAX_IGNITION_ZONES,' ignition zones are allowed.' + CALL SHUTDOWN(MESSAGE) ; RETURN + ENDIF + IF (DV%X1 > -1E6_EB) THEN !XB + IGNITION_ZONES(N_IGNITION_ZONES)%X1 = DV%X1 + IGNITION_ZONES(N_IGNITION_ZONES)%X2 = DV%X2 + IGNITION_ZONES(N_IGNITION_ZONES)%Y1 = DV%Y1 + IGNITION_ZONES(N_IGNITION_ZONES)%Y2 = DV%Y2 + IGNITION_ZONES(N_IGNITION_ZONES)%Z1 = DV%Z1 + IGNITION_ZONES(N_IGNITION_ZONES)%Z2 = DV%Z2 + ELSEIF (DV%X > -1E6_EB .AND. DV%Y > -1E6_EB .AND. DV%Z > -1E6_EB .AND. PROCESS(DV%MESH)==MY_RANK) THEN !XYZ + M => MESHES(DV%MESH) + CALL GET_IJK(DV%X,DV%Y,DV%Z,DV%MESH,XI,YJ,ZK,II,JJ,KK) + IGNITION_ZONES(N_IGNITION_ZONES)%X1 = M%X(II-1) + IGNITION_ZONES(N_IGNITION_ZONES)%X2 = M%X(II) + IGNITION_ZONES(N_IGNITION_ZONES)%Y1 = M%Y(JJ-1) + IGNITION_ZONES(N_IGNITION_ZONES)%Y2 = M%Y(JJ) + IGNITION_ZONES(N_IGNITION_ZONES)%Z1 = M%Z(KK-1) + IGNITION_ZONES(N_IGNITION_ZONES)%Z2 = M%Z(KK) + ENDIF + IGNITION_ZONES(N_IGNITION_ZONES)%DEVC_INDEX = N + IGNITION_ZONES(N_IGNITION_ZONES)%DEVC_ID = DV%ID + ENDIF + ENDIF + + IF(TRIM(ODE_SOLVER)=='CVODE' .AND. N_IGNITION_ZONES >= 1) THEN + IF (TRIM(FUEL_ID_FOR_AFT) == 'null') THEN + WRITE(MESSAGE,'(A)') 'ERROR(*): FUEL_ID_FOR_AFT must be specified when IGNITION_ZONE(s) are speciefied.' & + // NEW_LINE('A') // 'Provide FUEL_ID_FOR_AFT in the COMB line.' + CALL SHUTDOWN(MESSAGE) ; RETURN + ENDIF ENDIF ! Check if the output QUANTITY exists and is appropriate @@ -15180,7 +15305,7 @@ SUBROUTINE PROC_DEVC PY => PROPERTY(DV%PROP_INDEX) IF (PY%TIME_CONSTANT>0._EB) THEN ! Convert a specified TIME_CONSTANT into an equivalent bead DIAMETER IF (PY%HEAT_TRANSFER_COEFFICIENT>0._EB) THEN ! Calculate effective diameter directly - PY%DIAMETER = 6._EB*PY%HEAT_TRANSFER_COEFFICIENT*PY%TIME_CONSTANT/(PY%DENSITY*PY%SPECIFIC_HEAT) + PY%DIAMETER = 6._EB*PY%HEAT_TRANSFER_COEFFICIENT*PY%TIME_CONSTANT/(PY%DENSITY*PY%SPECIFIC_HEAT(293)) ELSE ! Calculate effective diameter implicitly ZZ_G = 0._EB ZZ_G(1) = 1._EB @@ -15190,15 +15315,16 @@ SUBROUTINE PROC_DEVC TOL = 1._EB RE = RHOA*UU*PY%DIAMETER/MU_G ! Make initial estimate of Re and Nu based on default bead diameter NU = 2._EB + 0.6_EB*SQRT(RE)*PR_AIR**ONTH - PY%DIAMETER = SQRT(6._EB*K_G*NU*PY%TIME_CONSTANT/(PY%DENSITY*PY%SPECIFIC_HEAT)) + PY%DIAMETER = SQRT(6._EB*K_G*NU*PY%TIME_CONSTANT/(PY%DENSITY*PY%SPECIFIC_HEAT(293))) DO WHILE(TOL>1.E-5_EB) RE = RHOA*UU*PY%DIAMETER/MU_G NU = 2._EB + 0.6_EB*SQRT(RE)*PR_AIR**ONTH - F = 6._EB*K_G*NU*PY%TIME_CONSTANT - PY%DENSITY*PY%SPECIFIC_HEAT*PY%DIAMETER**2 + F = 6._EB*K_G*NU*PY%TIME_CONSTANT - PY%DENSITY*PY%SPECIFIC_HEAT(293)*PY%DIAMETER**2 DFDD = 1.8_EB*K_G*PY%TIME_CONSTANT*PR_AIR**ONTH*RE**(-0.5)*RHOA*UU/MU_G - & - 2._EB*PY%DENSITY*PY%SPECIFIC_HEAT*PY%DIAMETER + 2._EB*PY%DENSITY*PY%SPECIFIC_HEAT(293)*PY%DIAMETER PY%DIAMETER = PY%DIAMETER - F/DFDD TOL = ABS(F/DFDD) + WRITE(*,*)'C',TRIM(PY%ID),PY%DIAMETER ENDDO ENDIF ENDIF diff --git a/Source/type.f90 b/Source/type.f90 index 35c4e727473..802ae18fe86 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -187,9 +187,12 @@ MODULE TYPES TYPE BOUNDARY_COORD_TYPE - INTEGER :: II,II2 !< Ghost cells x index - INTEGER :: JJ,JJ2 !< Ghost cells y index - INTEGER :: KK,KK2 !< Ghost cells z index + INTEGER :: II !< First ghost cell inside wall in the x direction + INTEGER :: JJ !< First ghost cell inside wall in the y direction + INTEGER :: KK !< First ghost cell inside wall in the z direction + INTEGER :: II2 !< Second ghost cell, x direction + INTEGER :: JJ2 !< Second ghost cell, y direction + INTEGER :: KK2 !< Second ghost cell, z direction INTEGER :: IIG !< Gas cell x index INTEGER :: JJG !< Gas cell y index INTEGER :: KKG !< Gas cell z index diff --git a/Source/wall.f90 b/Source/wall.f90 index 0b51fb2d7c3..57a0b5c9aaa 100644 --- a/Source/wall.f90 +++ b/Source/wall.f90 @@ -270,13 +270,16 @@ END SUBROUTINE WALL_BC SUBROUTINE ASSIGN_GHOST_VALUE(WALL_INDEX,BC,B1) USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT +USE COMPLEX_GEOMETRY, ONLY : CC_CGSC, CC_SOLID INTEGER, INTENT(IN) :: WALL_INDEX REAL(EB) :: ARO,ZZ_GET(1:N_TRACKED_SPECIES),RHO_OTHER,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS),& - RSUM_TMP,RHO_OTHER_2,RHO_ZZ_OTHER_2(1:N_TOTAL_SCALARS) -INTEGER :: IIO,JJO,KKO,II2,JJ2,KK2 + RSUM_TMP,RHO_OTHER_2,RHO_ZZ_OTHER_2(1:N_TOTAL_SCALARS),PBAR_P_2 +INTEGER :: IIO,JJO,KKO,II2,JJ2,KK2,ICG,ICO +LOGICAL :: CC_SOLID_FLAG,SECOND_ORDER_INTERPOLATED_BOUNDARY REAL(EB), POINTER, DIMENSION(:,:,:,:) :: OM_ZZP REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_RHOP TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC +TYPE(WALL_TYPE), POINTER :: WC TYPE(OMESH_TYPE), POINTER :: OM TYPE(MESH_TYPE), POINTER :: MM TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC @@ -331,13 +334,44 @@ SUBROUTINE ASSIGN_GHOST_VALUE(WALL_INDEX,BC,B1) CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM(BC%II,BC%JJ,BC%KK)) TMP(BC%II,BC%JJ,BC%KK) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM(BC%II,BC%JJ,BC%KK)*RHOP(BC%II,BC%JJ,BC%KK)) -! Assign RHOP and ZZP to the second ghost cell +! Assign RHOP and ZZP to the second ghost cell (if it is not a solid) -RHOP(BC%II2,BC%JJ2,BC%KK2) = RHO_OTHER_2 -ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TOTAL_SCALARS) = MAX(0._EB,MIN(1._EB,RHO_ZZ_OTHER_2(1:N_TOTAL_SCALARS)/RHO_OTHER_2)) -ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TRACKED_SPECIES) -CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_TMP) -TMP(BC%II2,BC%JJ2,BC%KK2) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM_TMP*RHOP(BC%II2,BC%JJ2,BC%KK2)) +SECOND_ORDER_INTERPOLATED_BOUNDARY = .FALSE. +WC => WALL(WALL_INDEX) +INTERPOLATED_IF: IF (WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY) THEN + ! Determine if there are 4 equally sized cells spanning the interpolated boundary + CC_SOLID_FLAG = .FALSE. + IF (ABS(EWC%AREA_RATIO-1._EB)<0.01_EB) THEN + IIO = EWC%IIO_MIN + JJO = EWC%JJO_MIN + KKO = EWC%KKO_MIN + SELECT CASE(BC%IOR) + CASE( 1) ; ICG = CELL_INDEX(BC%IIG+1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO-1,JJO,KKO) + CASE(-1) ; ICG = CELL_INDEX(BC%IIG-1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO+1,JJO,KKO) + CASE( 2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG+1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO-1,KKO) + CASE(-2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG-1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO+1,KKO) + CASE( 3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG+1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO-1) + CASE(-3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG-1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO+1) + END SELECT + IF (CC_IBM) THEN ! Test if one of surrounding cells is CC_SOLID. + IF(CCVAR(BC%IIG,BC%JJG,BC%KKG,CC_CGSC)==CC_SOLID .OR. CCVAR(BC%II,BC%JJ,BC%KK,CC_CGSC)==CC_SOLID) CC_SOLID_FLAG = .TRUE. + ENDIF + IF (.NOT.CELL(ICG)%SOLID.AND..NOT.MM%CELL(ICO)%SOLID.AND..NOT.CC_SOLID_FLAG) SECOND_ORDER_INTERPOLATED_BOUNDARY = .TRUE. + ENDIF +ENDIF INTERPOLATED_IF + +IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN + RHOP(BC%II2,BC%JJ2,BC%KK2) = RHO_OTHER_2 + ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TOTAL_SCALARS) = MAX(0._EB,MIN(1._EB,RHO_ZZ_OTHER_2(1:N_TOTAL_SCALARS)/RHO_OTHER_2)) + ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TRACKED_SPECIES) + CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_TMP) + PBAR_P_2 = 2._EB*PBAR_P(BC%KK,B1%PRESSURE_ZONE) - PBAR_P(BC%KKG,B1%PRESSURE_ZONE) + TMP(BC%II2,BC%JJ2,BC%KK2) = PBAR_P_2/(RSUM_TMP*RHOP(BC%II2,BC%JJ2,BC%KK2)) +ELSE + RHOP(BC%II2,BC%JJ2,BC%KK2) = RHOP(BC%II,BC%JJ,BC%KK) + ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TOTAL_SCALARS) = ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS) + TMP(BC%II2,BC%JJ2,BC%KK2) = TMP(BC%II,BC%JJ,BC%KK) +ENDIF END SUBROUTINE ASSIGN_GHOST_VALUE @@ -447,26 +481,22 @@ END SUBROUTINE NEAR_SURFACE_GAS_VARIABLES SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX) USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,INTERPOLATE1D_UNIFORM,GET_SCALAR_FACE_VALUE,GET_SCALAR_FACE_VALUE_NEW,GET_SCALAR_FACE_COEF -USE COMPLEX_GEOMETRY, ONLY : CC_CGSC, CC_SOLID USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT,GET_VISCOSITY,GET_MOLECULAR_WEIGHT USE DEVICE_VARIABLES, ONLY : PROPERTY,PROPERTY_TYPE -USE CHEMCONS, ONLY : FLAME_THICK_FACTOR REAL(EB), INTENT(IN) :: T INTEGER, INTENT(IN) :: NM INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX -REAL(EB) :: ARO,QNET,RAMP_FACTOR,RHO_G_2,RSUM_F,PBAR_F,TSI,UN, & +REAL(EB) :: ARO,QNET,RAMP_FACTOR,RSUM_F,PBAR_F,TSI,UN, & RHO_ZZ_F(1:N_TOTAL_SCALARS),ZZ_GET(1:N_TRACKED_SPECIES),& - RHO_OTHER,RHO_OTHER_2,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS),RHO_ZZ_OTHER_2,RHO_ZZ_G,RHO_ZZ_G_2, & - DDO,PBAR_G,PBAR_OTHER,DENOM,D_Z_N(0:I_MAX_TEMP),D_Z_G,D_Z_OTHER,TMP_OTHER, & - MU_DNS_G,MU_DNS_OTHER,MU_OTHER,RHO_D,RHO_D_TURB,RHO_D_DZDN,RHO_D_DZDN_OTHER,RSUM_OTHER,& + RHO_OTHER, PBAR_OTHER,D_Z_N(0:I_MAX_TEMP),D_Z_G,D_Z_OTHER,TMP_OTHER,DDO,DENOM,PBAR_G, & + MU_DNS_G,MU_DNS_OTHER,MU_OTHER,RHO_D,RHO_D_TURB,RHO_D_DZDN,RHO_D_DZDN_OTHER,RSUM_OTHER, & BBB,CCC,PPP,QQQ,RRR,UUU,YYY,WWW,HTC_OLD,RSC_LOC,DTMP,RSUM_G,MU_G -LOGICAL :: SECOND_ORDER_INTERPOLATED_BOUNDARY,ATMOSPHERIC_INTERPOLATION,CC_SOLID_FLAG -INTEGER :: IIO,JJO,KKO,N,ADCOUNT,ICG,ICO +LOGICAL :: ATMOSPHERIC_INTERPOLATION +INTEGER :: IIO,JJO,KKO,N,ADCOUNT +REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_RHOP REAL(EB), POINTER, DIMENSION(:,:,:,:) :: OM_ZZP REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_MUP -REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP,B_TEMP -REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC TYPE(OMESH_TYPE), POINTER :: OM TYPE(MESH_TYPE), POINTER :: MM @@ -478,7 +508,6 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC -REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_RHOP TYPE(PROPERTY_TYPE), POINTER :: PY IF (PRESENT(WALL_INDEX)) THEN @@ -668,6 +697,7 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I EWC => EXTERNAL_WALL(WALL_INDEX) OM => OMESH(EWC%NOM) + MM => MESHES(EWC%NOM) IF (PREDICTOR) THEN OM_RHOP => OM%RHOS OM_ZZP => OM%ZZS @@ -676,36 +706,9 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I OM_ZZP => OM%ZZ ENDIF - MM => MESHES(EWC%NOM) - RHO_G_2 = B1%RHO_G ! first-order - RHO_OTHER = RHOP(BC%II,BC%JJ,BC%KK) - RHO_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK) - RHO_ZZ_OTHER(1:N_TOTAL_SCALARS) = ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS)*RHO_OTHER - - ! Determine if there are 4 equally sized cells spanning the interpolated boundary - - SECOND_ORDER_INTERPOLATED_BOUNDARY = .FALSE. - CC_SOLID_FLAG = .FALSE. - IF (ABS(EWC%AREA_RATIO-1._EB)<0.01_EB) THEN - IIO = EWC%IIO_MIN - JJO = EWC%JJO_MIN - KKO = EWC%KKO_MIN - SELECT CASE(BC%IOR) - CASE( 1) ; ICG = CELL_INDEX(BC%IIG+1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO-1,JJO,KKO) - CASE(-1) ; ICG = CELL_INDEX(BC%IIG-1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO+1,JJO,KKO) - CASE( 2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG+1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO-1,KKO) - CASE(-2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG-1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO+1,KKO) - CASE( 3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG+1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO-1) - CASE(-3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG-1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO+1) - END SELECT - IF (CC_IBM) THEN ! Test if one of surrounding cells is CC_SOLID. - IF(CCVAR(BC%IIG,BC%JJG,BC%KKG,CC_CGSC)==CC_SOLID .OR. CCVAR(BC%II,BC%JJ,BC%KK,CC_CGSC)==CC_SOLID) CC_SOLID_FLAG = .TRUE. - ENDIF - IF (.NOT.CELL(ICG)%SOLID.AND..NOT.MM%CELL(ICO)%SOLID.AND..NOT.CC_SOLID_FLAG) SECOND_ORDER_INTERPOLATED_BOUNDARY = .TRUE. - ENDIF - - ! Density + ! interp or extrap RHO_OTHER for jump in vertical grid resolution, linear in temperature to match heat flux in divg + RHO_OTHER = RHOP(BC%II,BC%JJ,BC%KK) ATMOSPHERIC_INTERPOLATION = .FALSE. DDO = 1._EB KKO = EWC%KKO_MIN @@ -717,7 +720,6 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ATMOSPHERIC_INTERPOLATION = .TRUE. IF (ATMOSPHERIC_INTERPOLATION) THEN - ! interp or extrap RHO_OTHER for jump in vertical grid resolution, linear in temperature to match heat flux in divg PBAR_G = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) PBAR_OTHER = EVALUATE_RAMP(MM%ZC(EWC%KKO_MIN),I_RAMP_P0_Z) ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,MIN(1._EB,ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS))) @@ -727,199 +729,9 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I TMP(BC%II,BC%JJ,BC%KK) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM(BC%II,BC%JJ,BC%KK)*RHOP(BC%II,BC%JJ,BC%KK)) ENDIF - SELECT CASE(BC%IOR) - CASE( 1) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG+1,BC%JJG,BC%KKG) - RHO_OTHER_2 = RHOP(BC%II-1,BC%JJ,BC%KK) - ENDIF - Z_TEMP(0:3,1,1) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) - U_TEMP(1,1,1) = UU(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFX(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - CASE(-1) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG-1,BC%JJG,BC%KKG) - RHO_OTHER_2 = RHOP(BC%II+1,BC%JJ,BC%KK) - ENDIF - Z_TEMP(0:3,1,1) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) - U_TEMP(1,1,1) = UU(BC%II-1,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFX(BC%II-1,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - CASE( 2) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG,BC%JJG+1,BC%KKG) - RHO_OTHER_2 = RHOP(BC%II,BC%JJ-1,BC%KK) - ENDIF - Z_TEMP(1,0:3,1) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFY(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - CASE(-2) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG,BC%JJG-1,BC%KKG) - RHO_OTHER_2 = RHOP(BC%II,BC%JJ+1,BC%KK) - ENDIF - Z_TEMP(1,0:3,1) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ-1,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFY(BC%II,BC%JJ-1,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - CASE( 3) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG+1) - RHO_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK-1) - ENDIF - Z_TEMP(1,1,0:3) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - CASE(-3) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG-1) - RHO_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK+1) - ENDIF - Z_TEMP(1,1,0:3) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-1) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK-1) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ENDIF - B1%RHO_F = F_TEMP(1,1,1) - END SELECT - - ! Species and temperature + ! Density, Species and temperature SINGLE_SPEC_IF: IF (N_TOTAL_SCALARS > 1) THEN - SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS - - RHO_ZZ_G = B1%RHO_G*B1%ZZ_G(N) - RHO_ZZ_G_2 = RHO_ZZ_G ! first-order (default) - RHO_ZZ_OTHER_2 = RHO_ZZ_OTHER(N) - - SELECT CASE(BC%IOR) - CASE( 1) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG+1,BC%JJG,BC%KKG)*ZZP(BC%IIG+1,BC%JJG,BC%KKG,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II-1,BC%JJ,BC%KK)*ZZP(BC%II-1,BC%JJ,BC%KK,N) - ENDIF - Z_TEMP(0:3,1,1) = (/RHO_ZZ_OTHER_2,RHO_ZZ_OTHER(N),RHO_ZZ_G,RHO_ZZ_G_2/) - U_TEMP(1,1,1) = UU(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFX(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) - CASE(-1) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG-1,BC%JJG,BC%KKG)*ZZP(BC%IIG-1,BC%JJG,BC%KKG,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II+1,BC%JJ,BC%KK)*ZZP(BC%II+1,BC%JJ,BC%KK,N) - ENDIF - Z_TEMP(0:3,1,1) = (/RHO_ZZ_G_2,RHO_ZZ_G,RHO_ZZ_OTHER(N),RHO_ZZ_OTHER_2/) - U_TEMP(1,1,1) = UU(BC%II-1,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFX(BC%II-1,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) - CASE( 2) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG+1,BC%KKG)*ZZP(BC%IIG,BC%JJG+1,BC%KKG,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II,BC%JJ-1,BC%KK)*ZZP(BC%II,BC%JJ-1,BC%KK,N) - ENDIF - Z_TEMP(1,0:3,1) = (/RHO_ZZ_OTHER_2,RHO_ZZ_OTHER(N),RHO_ZZ_G,RHO_ZZ_G_2/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFY(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) - CASE(-2) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG-1,BC%KKG)*ZZP(BC%IIG,BC%JJG-1,BC%KKG,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II,BC%JJ+1,BC%KK)*ZZP(BC%II,BC%JJ+1,BC%KK,N) - ENDIF - Z_TEMP(1,0:3,1) = (/RHO_ZZ_G_2,RHO_ZZ_G,RHO_ZZ_OTHER(N),RHO_ZZ_OTHER_2/) - U_TEMP(1,1,1) = VV(BC%II,BC%JJ-1,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFY(BC%II,BC%JJ-1,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) - CASE( 3) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG+1)*ZZP(BC%IIG,BC%JJG,BC%KKG+1,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK-1)*ZZP(BC%II,BC%JJ,BC%KK-1,N) - ENDIF - Z_TEMP(1,1,0:3) = (/RHO_ZZ_OTHER_2,RHO_ZZ_OTHER(N),RHO_ZZ_G,RHO_ZZ_G_2/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = (PBAR_P(BC%KK,B1%PRESSURE_ZONE)*DZ(BC%KKG) + PBAR_P(BC%KKG,B1%PRESSURE_ZONE)*DZ(BC%KK)) / & - (DZ(BC%KK)+DZ(BC%KKG)) - CASE(-3) - IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN - RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG-1)*ZZP(BC%IIG,BC%JJG,BC%KKG-1,N) - RHO_ZZ_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK+1)*ZZP(BC%II,BC%JJ,BC%KK+1,N) - ENDIF - Z_TEMP(1,1,0:3) = (/RHO_ZZ_G_2,RHO_ZZ_G,RHO_ZZ_OTHER(N),RHO_ZZ_OTHER_2/) - U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-1) - IF (FLUX_LIMITER_SINGLE_COEF) THEN - B_TEMP(1,1,1) = BFZ(BC%II,BC%JJ,BC%KK-1) - CALL GET_SCALAR_FACE_VALUE_NEW(U_TEMP,Z_TEMP,F_TEMP,B_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ELSE - CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER) - ENDIF - RHO_ZZ_F(N) = F_TEMP(1,1,1) - PBAR_F = (PBAR_P(BC%KK,B1%PRESSURE_ZONE)*DZ(BC%KKG) + PBAR_P(BC%KKG,B1%PRESSURE_ZONE)*DZ(BC%KK)) / & - (DZ(BC%KK)+DZ(BC%KKG)) - END SELECT - ENDDO SPECIES_LOOP - - ! face value of temperature IF (ATMOSPHERIC_INTERPOLATION) THEN B1%TMP_F = (TMP(BC%II,BC%JJ,BC%KK)*DZ(BC%KKG) + TMP(BC%IIG,BC%JJG,BC%KKG)*DZ(BC%KK)) / (DZ(BC%KK)+DZ(BC%KKG)) @@ -927,11 +739,29 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ZZP(BC%IIG,BC%JJG,BC%KKG,1:N_TOTAL_SCALARS)*DZ(BC%KK)) / (DZ(BC%KK)+DZ(BC%KKG)) ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F) + SELECT CASE(ABS(BC%IOR)) + CASE DEFAULT + PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) + CASE(3) + ! for very large cells with stratification the pressure variation in z can be significant + PBAR_F = (PBAR_P(BC%KK,B1%PRESSURE_ZONE)*DZ(BC%KKG) + PBAR_P(BC%KKG,B1%PRESSURE_ZONE)*DZ(BC%KK)) / & + (DZ(BC%KK)+DZ(BC%KKG)) + END SELECT B1%RHO_F = PBAR_F/(RSUM_F*B1%TMP_F) ELSE + SELECT CASE(BC%IOR) + CASE( 1); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FX(BC%II ,BC%JJ,BC%KK,1:N_TOTAL_SCALARS) + CASE(-1); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FX(BC%II-1,BC%JJ,BC%KK,1:N_TOTAL_SCALARS) + CASE( 2); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FY(BC%II,BC%JJ ,BC%KK,1:N_TOTAL_SCALARS) + CASE(-2); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FY(BC%II,BC%JJ-1,BC%KK,1:N_TOTAL_SCALARS) + CASE( 3); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FZ(BC%II,BC%JJ,BC%KK ,1:N_TOTAL_SCALARS) + CASE(-3); RHO_ZZ_F(1:N_TOTAL_SCALARS) = FZ(BC%II,BC%JJ,BC%KK-1,1:N_TOTAL_SCALARS) + END SELECT + B1%RHO_F = SUM(RHO_ZZ_F(1:N_TRACKED_SPECIES)) B1%ZZ_F(1:N_TOTAL_SCALARS) = MAX(0._EB,MIN(1._EB,RHO_ZZ_F(1:N_TOTAL_SCALARS)/B1%RHO_F)) ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES) CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F) + PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE) B1%TMP_F = PBAR_F/(RSUM_F*B1%RHO_F) ENDIF @@ -971,9 +801,6 @@ SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_I ENDIF RHO_D = 0.5_EB*( RHO_OTHER*D_Z_OTHER + B1%RHO_G*D_Z_G ) + RHO_D_TURB END SELECT MODE_SELECT - IF (FLAME_THICK_FACTOR > 1._EB) THEN - RHO_D = RHO_D*FLAME_THICK_FACTOR - ENDIF SELECT CASE(BC%IOR) CASE( 1) ARO = MM%DY(JJO)*MM%DZ(KKO)/(DY(BC%JJ)*DZ(BC%KK)) @@ -1043,7 +870,6 @@ END SUBROUTINE SURFACE_HEAT_TRANSFER !> \details Calculate the diffusion coefficient, rho*D (kg/m/s) SUBROUTINE CALCULATE_RHO_D_F(B1,BC,WALL_INDEX,CFACE_INDEX) -USE CHEMCONS, ONLY : FLAME_THICK_FACTOR INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX REAL(EB) :: MU_G @@ -1077,9 +903,6 @@ SUBROUTINE CALCULATE_RHO_D_F(B1,BC,WALL_INDEX,CFACE_INDEX) B1%RHO_D_F(N) = B1%RHO_F*D_Z(ITMP,N) ENDDO END SELECT -IF (FLAME_THICK_FACTOR > 1._EB) THEN - B1%RHO_D_F = B1%RHO_D_F*FLAME_THICK_FACTOR -ENDIF END SUBROUTINE CALCULATE_RHO_D_F @@ -1500,7 +1323,7 @@ SUBROUTINE CALCULATE_ZZ_F(T,DT,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX) ENDIF IF (PREDICTOR) B1%U_NORMAL_S = -UN - IF (CORRECTOR) B1%U_NORMAL = -UN + IF (CORRECTOR) B1%U_NORMAL = -UN END SELECT METHOD_OF_MASS_TRANSFER @@ -1509,6 +1332,7 @@ SUBROUTINE CALCULATE_ZZ_F(T,DT,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX) IF (PRESENT(WALL_INDEX)) THEN IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS .AND. .NOT.CELL(IC)%SOLID .AND. .NOT.CELL(ICG)%SOLID) & ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES) = 2._EB*B1%ZZ_F(1:N_TRACKED_SPECIES) - B1%ZZ_G(1:N_TRACKED_SPECIES) + ZZP(BC%II2,BC%JJ2,BC%KK2,1:N_TRACKED_SPECIES) = 2._EB*B1%ZZ_F(1:N_TRACKED_SPECIES) - B1%ZZ_G(1:N_TRACKED_SPECIES) ENDIF END SUBROUTINE CALCULATE_ZZ_F @@ -1674,7 +1498,10 @@ SUBROUTINE CALCULATE_RHO_F(BC,B1,WALL_INDEX,CFACE_INDEX) ENDIF IF (PRESENT(WALL_INDEX)) THEN - IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS .AND. BOUNDARY_TYPE/=OPEN_BOUNDARY) RHOP(BC%II,BC%JJ,BC%KK) = 2._EB*B1%RHO_F - B1%RHO_G + IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS .AND. BOUNDARY_TYPE/=OPEN_BOUNDARY) THEN + RHOP(BC%II ,BC%JJ, BC%KK ) = 2._EB*B1%RHO_F - B1%RHO_G + RHOP(BC%II2,BC%JJ2,BC%KK2) = 2._EB*B1%RHO_F - B1%RHO_G + ENDIF ENDIF END SUBROUTINE CALCULATE_RHO_F @@ -3344,7 +3171,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F ! Compute the Spalding B number - SUM_Y_SV = SUM(Y_SV) + SUM_Y_SV = SUM(Y_SV(1:N_MATS)) SUM_Y_SV_SMIX = 0._EB MATERIAL_LOOP_1: DO N=1,N_MATS IF (.NOT.LIQUID(N)) CYCLE MATERIAL_LOOP_1 @@ -3352,7 +3179,7 @@ SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F ENDDO MATERIAL_LOOP_1 IF (SUM_Y_SV