Skip to content

ziolai/software

BranchesTags

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

140 Commits
figures
figures
 
 
milling
milling
 
 
my-simple-geometry-1D
my-simple-geometry-1D
 
 
my-simple-geometry-2D
my-simple-geometry-2D
 
 
my-simple-geometry-3D
my-simple-geometry-3D
 
 
myLaplacianFoam
myLaplacianFoam
 
 
pitzDaily_Reynolds_stress_R_manual
pitzDaily_Reynolds_stress_R_manual
 
 
reactive-flow-course
reactive-flow-course
 
 
thermal-oxidizer
thermal-oxidizer
 
 
Sandia_FlameD.ipynb
Sandia_FlameD.ipynb
 
 
chemical_reactor_networks.ipynb
chemical_reactor_networks.ipynb
 
 
counterflow-diffusion-flame.ipynb
counterflow-diffusion-flame.ipynb
 
 
easyMillLaminarSimpleFoam_residual_norm.png
easyMillLaminarSimpleFoam_residual_norm.png
 
 
fgm-table-generation.ipynb
fgm-table-generation.ipynb
 
 
helmholtz-unit-interval-fem.ipynb
helmholtz-unit-interval-fem.ipynb
 
 
helmholtz-unit-interval-gmres-fdm.ipynb
helmholtz-unit-interval-gmres-fdm.ipynb
 
 
helmholtz-unit-square-fdm.ipynb
helmholtz-unit-square-fdm.ipynb
 
 
intro_ode.ipynb
intro_ode.ipynb
 
 
p1-radiation-model-interval-fdm.ipynb
p1-radiation-model-interval-fdm.ipynb
 
 
poisson-unit-bar.ipynb
poisson-unit-bar.ipynb
 
 
read-foamLog.ipynb
read-foamLog.ipynb
 
 
readme.ipynb
readme.ipynb
 
 
readme.md
readme.md
 
 
simplefoam_resconv_relaxation.png
simplefoam_resconv_relaxation.png
 
 
swi-groningen-2023.ipynb
swi-groningen-2023.ipynb
 
 

Repository files navigation

OpenFoam for Turbulent Reactive Flow

Section 1:/ Motivation for these Notes

  • develop backbone for lectures and further research on turbulent reactive flows with emphasis on (non, partially)-premixed (laminar,turbulent) combustion using OpenFOAM;
  • develop model for turbulence, chemistry and radiation in the freeboard coupled with heat conduction in the lining through a conjugate heat transfer model;
  • develop understand of OpenFOAM code separated into source code and application ;

Section 2:/ Introduction

  • OpenFOAM: cell-centered finite volume method on unstructured 3D grids;
  • object oriented in C++ using templates: definition in header file, implementation in corresponding C-file (could give example here);
  • parallel by mesh decomposition and parallel vector, matrix and solver architecture;
  • examples: mesh = object + operations (check on consistency, non-orthogonality), vector, matrix, ODE solver;
  • class (primitive mesh: mesh and operations on the mesh: which ones?) and superclass (polymesh: definition of patches and operations on patches: which ones?);
  • question: what kind of volumetric and surface patches can OpenFoam represent?
  • tensor calculus: templated classes to support tensors with various types of entries (integer, real or complex);
  • fvc and fvm: calculus (explicit) and method (implicit); using namespaces; attempt to keep operators diagonally dominant (beware of the initial values used to initialize non-linear iterations); all operators in either implicit and explicit versions;
  • Rhie-Chow interpolation to suppress spurious oscillations in the pressure field (outer corrections);
  • OpenFOAM distribution divided into source and application;

Books

Papers

Videos

Section 3: Not in these Notes - Disclaimer

Section 4: OpenFoam Code

  • fvCFH.H: main include file;
  • object registry: you can think of an objectRegistry as being a database for holding objects (objects of type regIOobject to be specific). Since the objectRegistry itself is derived from a regIOobject, an objectRegistry can thus contain other object registries, and can potentially have a parent registry as well. Time is a special case of an objectRegistry that is also used to define the time coordination for multiple registries. It probably helps to think in terms of a file-system. The objectRegistry is analogous to a directory, the regIOobject is analogous to a file/directory. See also what-object-registry Gerhard Holzinger 2021 and openfoam wiki
  • examples of lookupObject include the examples below; const volScalarField& nut = mesh_.lookupObject("nut"); const volScalarField& alphat = mesh_.lookupObject("alphat"); const volScalarField& T = mesh_.lookupObject("T");
  • autoPtr and tmp: smart pointers: resource acquisition is initialization : if I wrap my newly allocated piece of data in a proper class with a destructor, it will get cleaned up automatically, as soon as the destructor of this class is called: understanding-autoptr tomislav_maric 2012 and tmp OpenFoam wiki;

Section 5: Install, compile, link, run and debug OpenFOAM

UnnamedMoose Tutorials

UnnamedMoose Tutorials

OpenFoam and Combustion Simulation Webinar

OpenFoam and Combustion Simulation Webinar

Parametric Studies

  • parametric studies using parametric studies PyFoam;

Parallel Runs

  • parallel processing using decomposePar (by a simple decomposition or by calling a mesh decomposition package) and reconstructPar;

  • what posssibly goes wrong in very large models: conda-forge/mpi4py-feedstock#15 see https://www.cfd-online.com/Forums/openfoam-solving/232470-cleaning-out-timesteps-processor-directories.html I would not delete processor* directories if you do not change the mesh.

  • If you want to replace the processor*/0 directory from the new 0.orig directory: -- rm -rf processor*/<time> && restore0Dir -processor -- You need to source /bin/tools/RunFunctions to use restore0Dir.

  • If you want to decompose fields only, and not the existing mesh: -- decomposePar -fields

  • If you want to perform decompositions in parallel: -- mpirun -np X redistributePar -decompose -parallel

Section 6: Mesh Generation and Handling

- lduMesh: provides LDU addressing: https://github.com/OpenFOAM/OpenFOAM-dev/blob/3db8158b7b90c204f1c33f3bb458cbf73a5585ff/src/OpenFOAM/meshes/lduMesh/lduMesh.H

- fvMesh: mesh data needed to do the Finite Volume discretisation: https://github.com/OpenFOAM/OpenFOAM-dev/blob/3db8158b7b90c204f1c33f3bb458cbf73a5585ff/src/finiteVolume/fvMesh/fvMesh.H - the fvPatch class has a method delta() which returns the cell-centre to the face-centre vector. One can get the magnitude of those vectors to get the distance to the nearest wall and compute an estimate for y+; - conservative interpolation between meshes; - multiple region and processors interconnected by interfaces; - mesh conversion from other sources; - format used by Metis to decompose the mesh;

Section 6.1: Mesh Generation using blockMesh

Section 6.2: Mesh Generation using classyBlocks

Section 6.3: Mesh Generation using cfMesh

run cartesianMesh with dictionary provided in system/meshDict for only the gas domain

Section 6.4: Mesh Manipulation

Section 7: Post Processing

Settings in system/controlDict

writeControl in case that variable time step is used; Talice uses writeControl adjustable: need to look into manual what this means;

Function Object

function objects for post-processing: https://www.openfoam.com/documentation/guides/latest/doc/openfoam-guide-post-processing.html

computing averages of fields: https://www.openfoam.com/documentation/guides/latest/doc/guide-fos-field-field-average.html

probe plots: https://cfd.direct/openfoam/user-guide/v8-graphs-monitoring/#x33-2680006.3 and https://www.youtube.com/watch?v=w_CPMRykzJM&feature=youtu.be&fbclid=IwAR0ZdYZGNKVC2lbCQkXdPs1W7YHj2Iu6bQQopElewuQiIeGhTf107mhiRQw

yPlus: simpleFoam -postProcess -latestTime -func yPlus

Mach number: using -postProcess -func "MachNo(U)” or force a function object to be executed via settings in system/controlDict and system/MachNo. Implementation is in https://github.com/OpenFOAM/OpenFOAM-6/blob/master/src/functionObjects/field/MachNo/MachNo.H writes a volScalarField that represents the Mach number;

wall heat flux: as before, writes a wallHeatFlux that can be visualized in Paraview: make volumetric plots of the wallheat flux. This will lower the value of the heat flux due to interpolation with zero values on the boundary. Later you can modify the plot to visualize on patches only.

Q-criterium https://www.cfd-online.com/Forums/openfoam-solving/232098-order-execution-post-processing-functions.html

Paraview and Python for Paraview

Paraview: Hakan Nillson ParaFoam tutorial:

pyFoam: slides Gschaider 11th OpenFoam Workshop on Portugal

pickled format in Python (used in connection with pyFoam): https://pythontips.com/2013/08/02/what-is-pickle-in-python/

basic information on Python for Paraview: https://www.paraview.org/Wiki/ParaView_and_Python more detailed information on Python for Paraview: https://www.paraview.org/Wiki/ParaView/Python_Scripting To run Python for Paraview outside of the ParaView GUI, please use pvpython (pvpython.exe on Windows). This executable is part of ParaView installation.

Residual Norm Plotting

foamLog and foamMonitor: see tutorial https://www.youtube.com/watch?v=ZWEaMNPZpYI&feature=share&fbclid=IwAR3k9wGz23r6BKWEadXexVGQWDQD51tiVV3rlf-IAJi1bbae_blwhAxYjzc

Plotting the residual norm of various fields vs. iteration number
write residual norms of the various fields to file instead of to screen. Do so using so-called unix pipe command, e.g., simpleFoam >& log.simpleFoam or use the runApplication utility instead (typically inside of Allrun);
use foamLog utility to generate the the logs-directory in the case directory; foamLog runs a bunch of scripting commands to read data from the log-file and reorders the information in separate files; use gnuplot to load and plot the data in the logs-directory: in separate terminal run ./myplot. See sample myplot file below

Sample myplot file to plot residual norm vs. iteration number #!/bin/bash

Requires running foamLog first!!

gnuplot -persist >/dev/null 2>&2 << EOF set logscale y plot "logs/UxFinalRes_0" with lines,
"logs/UyFinalRes_0" with lines,
"logs/UzFinalRes_0" with lines,
"logs/pFinalRes_0" with lines,
"logs/kFinalRes_0" with lines,
"logs/epsilonFinalRes_0" with lines

Residual Field Plotting

Changes to be made in the system directory

1.1/ In system-folder/controlDict-file

… writeInterval 10; …

functions { #includeFunc residualFieldDict }

In system-folder/residualFieldDict-file

type solverInfo; // specific for version-1906 libs ("libutilityFunctionObjects.so");

writeControl runTime; // timeStep; writeInterval 10; // make sure this to be equal as the writeInterval specified in the controlDict file

writeResidualFields yes;

fields (".*"); //

Remarks

How to make sure that writeInterval in system-folder/controlDict-file and system-folder/residualFieldDict-file

References

in OpenFoam.com: using solverInfo: https://develop.openfoam.com/Development/openfoam/blob/master/src/functionObjects/utilities/solverInfo/solverInfo.H residual.cfg: https://develop.openfoam.com/Development/openfoam/blob/master/etc/caseDicts/postProcessing/numerical/residuals.cfg

Producing plots as *.png file while the solver is running (no graphical terminal required) process in two stages in Stage1: write data to file as solver is running by extending the functions dictionary in the controlDict file; in Stage2: produce *.png files from the data (without required graphical terminal) the trick is to produce *.png files from data written to file without requiring a terminal with graphical capabilities; above can be extended for probes and surface plots; for more information and examples, see https://cfd.direct/openfoam/user-guide/v6-graphs-monitoring/

Stage1: Write data to file at values-of-time-instances-for-which-you-would-like-to-see-the-graphs in system-folder/controldict-file, include the line functions { .. // other functions #includeFunc singleGraph_1 #includeFunc singleGraph_2 #includeFunc singleGraph_3

} Here singleGraph_1, singleGraph_2 and singleGraph_3 (and possibly other) are files in system-folder/. Each of these files has the outline given below

/--------------------------------- C++ -----------------------------------
========= | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | \ / O peration | Version: v1806 |
| \ / A nd | Web: www.OpenFOAM.com | \/ M anipulation |

Description Writes graph data for specified fields along a line, specified by start and end points.

*---------------------------------------------------------------------------*/

start (0 2.7001 0); // start coordinate in xyz-space end (0 40 0); // end coordinate in xyzzy-space fields (U k epsilon T p); // fields to plot region gas; // region to sample

// Sampling and I/O settings #includeEtc "caseDicts/postProcessing/graphs/sampleDict.cfg"

// Override settings here, e.g. // setConfig { type midPoint; }

// Must be last entry #includeEtc "caseDicts/postProcessing/graphs/graph.cfg"

// ************************************************************************* //

Running the solver with this setting will produce files of the form case-folder/postProcessing/singleGraph_1/region/time-stamp/line.xy

Stage 2: Use Gnuplot in batch mode to create *.png file from the data at values-of-time-instances-for-which-you-would-like-to-see-the-graphs
run gnuplot in batch mode using the gnuplot_export script
the gnuplot_export script takes data from postprocessing subfolder and writes the figure as *.png file into the gnuplot_expo folder. Make sure to clean up the postprocessing folder before starting the simulations again;
run “gnuplot -c gnuplot_export “. The option -c allows to run gnuplot in script mode. No interactive session is required with graphical output is required.

Below is a sample gnuplot_export script

reset set xrange [:] set yrange [:] set term png size 1305,900 set output 'gnuplot_export/singleGraph_1/Ux.png' set style data linespoints set lmargin 15 set bmargin 5 set xrange [0:37.3] set tics font "Times,20" set key font "Times,20" set title "Ux" font "Times,20" set ylabel 'Ux [m/s]' font "Times,20" offset -3,0 set xlabel 'y [m]' font "Times,20" offset 0,-0.5 set xtics ("2.7" 0, "5" 2.3, "10" 7.3, "15" 12.3, "20" 17.3, "25" 22.3, "30" 27.3, "35" 32.3, "40" 37.3) offset 0,-0.5 plot "postProcessing/singleGraph_1/gas/18000/line_U.xy" u 1:2 title "18000", "postProcessing/singleGraph_1/gas/17995/line_U.xy" u 1:2 title "17995", "postProcessing/singleGraph_1/gas/17000/line_U.xy" u 1:2 title "17000"

Ideally we would like the time stamp as run time option for these scripts.

Averaging of there field quantities average the fields over the flame/jet has been developed (from iteration 2k or so onwards)

Plotting psi in OpenFoam in controlDict place { type writeRegisteredObject; functionObjectLibs ("libIOFunctionObjects.so"); objectNames ("bananas"); outputControl outputTime; outputInterval 1; } Here, "bananas" will give you the list of all registered object, replace by the one you want. See https://www.cfd-online.com/Forums/openfoam-post-processing/100307-how-output-other-properties.html

POD decomposition of the Flow Results https://github.com/IllinoisRocstar/AccelerateCFD_CE

Dynamic Mode Decomposition Function Object https://www.openfoam.com/news/main-news/openfoam-v20-06/post-processing

What triggers function objects to be executed https://www.cfd-online.com/Forums/openfoam-programming-development/144913-what-trigger-function-objects.html

Section 8: Linear Solvers

Describe here:

  1. vectors (or fields): field, scalarField (good example of object oriented programming), vectorField, surfaceField, volField; operations on vectors and matrices with a parallel layout;

  2. operations on vectors (schemes)

  3. matrices: lduMatrix and fvScalarMatrix; parallel assembly and parallel solve;

  4. linear solvers; solving an equation typically involves going through there following four steps:

a. define the equation with fvScalarMatrix as an output; b. relax the equation (as oppposed to relaxing the field) by calling the relax member function of the class; The objectivite here is to put as much weight on the diagonal of the matrix as possible; c. apply the constraints by calling the constraint (or correct ?) member function (Dirichlet and Neumann boundary conditions); d. solve the equations by calling the solve member function;

definition of relax and correct as member functions of fvScalarMatrix are given in https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

linear system residual computation: https://www.cfd-online.com/Forums/openfoam-solving/95540-residuals-calculation.html

nice documentation of foamacademy.com: https://www.foamacademy.com/wp-content/uploads/2018/03/OF2018matrix.pdf

Relaxation

AMGX

  1. CDF-Online dialogue: https://www.cfd-online.com/Forums/openfoam-programming-development/163162-has-any-one-tried-amgx-openfoam.html

PETSc4Foam

  1. video by Wolf Dynamics: https://www.youtube.com/watch?v=_TWYCwpRTlY

Section 9: ODE Solvers

Section 10: Basic Solvers

  • laplacianFOAM: stationary and transient diffusion; possibly variable diffusion coefficient; solves for T (volScalarField) with a diffusion coefficient DT (dimensionScalar); solve placed inside a SIMPLE loop; accuracy of the SIMPLE loop governed by the residual; the SIMPLE loop steers non-orthogonal corrections and checks the convergence;

  • scalartransportFOAM: stationary and in stationary convection-diffusion; velocity field U read from file and thus given by a previous solve; as laplacianFOAM except for added convective term;

  • potentialFOAM: potential flow solved for pressure potential Phi (volScalarField); reconstruct velocity U from the flux phi of the computed solution; better correspondence with analytical solution in case that orthogonal meshes are employed;

Temporal Discretization

  • more details on the discretization of the time-derivative term (see Moukaled-Darwish);

  • Euler: first order backward Euler;

  • backward: three-time-steps backward difference method;

  • CranckNicolson: modified Crank-Nicolson scheme where an off-center parameter with a value between 0 and 1 has been introduced to control the weight toward either standard Crank-Nicolson scheme or the implicit Euler scheme;

  • bound on the time-step (maxDeltaT) and on the Courant number (maxCo). Both are set in the PIMPLE dictionary in the system/fvSolutions file. Courant number is governing the time step even though an implicit time-stepping method is choosen;

  • local time-stepping method; explained on CFD-Online Local Time Stepping; Quote: compute "timesteps" for each cell based off a max courant number, smooth this timestep field (in space) to improve stability, and then damp this timestep field (in time, as in limit changes) to further improve stability: end quote. Smoothing is controlled by the parameters rDeltaTSmoothingCoeff and rDeltaTDampingCoeff set in the PIMPLE dictionary in the system/fvSolutions file; SLTS phi rho 1.0 in ddt (fvSchemes, details missing);

  • time integration for wave equation in OpenFoam: paper On the use of Euler and Crank-Nicolson time-stepping schemes for seakeeping simulations in OpenFOAM

  • second order time discretization: paper Consistent second-order time-accurate non-iterative PISO-algorithm

Spatial Discretization

Coded Sources and fvOptions

in the file fvOptions options { energySource { type scalarSemiImplicitSource;

    timeStart       0.2;
    duration        0.3;
    selectionMode   cellSet;
    cellSet         ignition;
    volumeMode      specific;

    injectionRateSuSp
    {
        h          (2e7 0); // kg/m/s^3
    }
}

Section 11: Boundary Conditions

general information on boundary conditions: https://www.openfoam.com/documentation/guides/latest/doc/openfoam-guide-boundary-conditions.html and https://www.cfd-online.com/Forums/openfoam-solving/88762-updating-boundary-conditions-each-iteration.html

difference between imposed and calculated boundary conditions. Imposed boundary conditions remain fixed throughout computations. Calculated boundary conditions maybe be updated throughout computations. Example is the computations of nut on the boundary depending on k and eps values. Imposed fixed (as opposed to calculated) conditions on density proved to be essential to reach convergence in both the rotary kiln and anode baking case;

valuable discussiion updateCoeff() and evaluate(): https://www.cfd-online.com/Forums/openfoam-programming-development/123913-updatecoeffs-evaluate.html

use of “value” in some boundary conditions: https://www.cfd-online.com/Forums/openfoam-pre-processing/234124-clarification-k-omega-settings.html

Inlet Boundary Conditions

k: turbulentIntensityKineticEnergyInlet: boundary condition for turbulent kinetic energy on an inlet patch: for details, see https://www.openfoam.com/documentation/guides/latest/doc/guide-bcs-inlet-turbulent-k-turbulent-intensity-kinetic-energy.html

inlet condition for k requires a turbulent intensity: typically set to 0.05 or 5 percent;

epsilon: turbulentMixingLengthDissipationRateInlet: for details, see https://www.openfoam.com/documentation/guides/latest/doc/guide-bcs-inlet-turbulent-epsilon-turbulent-mixing-length-dissipation-rate.html

inlet condition for epsilon requires a turbulent length scale; typically set to 0.038 times the hydraulic parameter; hydraulic parameter is defined as 4 * area / perimeter: see https://en.wikipedia.org/wiki/Hydraulic_diameter

Outlet Boundary Conditions

for both k and epsilon: inletOutletConditions: prevents reverse flow on a boundary patch. See https://www.cfd-online.com/Forums/openfoam-solving/60337-questions-about-inletoutlet-outletinlet-boundary-conditions.html

Wall Boundary Conditions

for k: kqRWallFunction with value same as initial value specified by internalField

for epsilon: epsilonWallFunction interesting discussion of implementation of wall functions in OpenFoam: Item #6 in https://www.cfd-online.com/Forums/openfoam-solving/229743-les-simulation-wall-modeling.html

Wave Transmissive Conditions

wave transmissive boundary conditions: from the users guide https://www.openfoam.com/documentation/guides/latest/api/waveTransmissiveFvPatchField_8H_source.html or https://openfoamwiki.net/index.php/HowTo_Using_the_WaveTransmissive_Boundary_condition (confusing to read)

Section 12: Transport Models

Example of a small class, private and public data members, private and public member functions.

Section 13: Thermodynamical Models

internal energy e, kinetic energy K and enthalpy h thermo dynamics library (JANAF, polynomial tables, possible time consuming to evaluate) based on rho or psi; equation for energy (internal energy or the enthalpy) EEqn; switch between e or h; see e.g. https://www.cfd-online.com/Forums/openfoam-solving/241020-hconst-econst.html definition of p = p_rgh + rho*gh; static pressure, total pressure and isentropic pressure: https://www.openfoam.com/documentation/guides/latest/doc/guide-fos-field-pressure.html rhoSimpleFoam does not allow to the source terms of chemistry and radiation into account; computing temperature from internal sensible enthalpy: https://www.cfd-online.com/Forums/openfoam-programming-development/216577-temperature-calculation-sensible-internal-energy.html compressibility: psi = \partial \rho / \partial p (change of density due to pressure variations): barotropic relationship thermodynamics class: https://wiki.openfoam.com/Thermodynamic_class_by_Isabelle_Choquet

Thermophysical properties

constant-folder/thermophysicalTransport-file: example of cavity for rhoPimpleFoam: https://github.com/OpenFOAM/OpenFOAM-8/blob/master/tutorials/compressible/rhoPimpleFoam/RAS/cavity/constant/thermophysicalTransport example: set non-unit Lewis number:
LES { model nonUnityLewisEddyDiffusivity; Prt 0.85; Sct 0.7; }

Thermodynamics in function of rho (density) or in terms of psi (compressibility)

RhoThermo: thermo.rho() returns a stored rho field that is calculated from pressure and temperature fields according to the selected thermophysical model in thermo.correct(); rhothermo: basic thermodynamic properties based on density; for fluids with low compressibility like water where the density change is mostly due to temperature variation; huge pressure jumps will not be represented correctly; solvers using rhothermo (like the buoyant solvers) use a simplified pressure equation where psi is accounted for explicitly;

Code for pEqn (as used for instance in heattransfer/buoyantPimpleFoam)

$$ fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) + fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)

fvOptions(psi, p_rgh, rho.name()) $$

PsiThermo: basic thermodynamic properties based on compressibility; for gas; combustion solvers and transsonic solvers use the psithermo model. They have another pressure equation where the change in time of pressure is treated implicit by dit(psi, p); larger pressure jumps are now possible. Such solvers are used when the flow is mainly driven by pressure changes and temperature change is only an effect of large compression and expansion;

Code for pEqn (as used for instance in combustion/reactingFoam)

$$ fvm::ddt(psi, p) + fvc::div(phiHbyA) - fvm::laplacian(rhorAUf, p)

fvOptions(psi, p, rho.name()) $$

note that alternatives above are essentially different in the ddt-terms. The code for pEqn in heat transfer/bouyantSimpleFoam contains fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)

rhoSimpleFoam now instantiates the lower-level fluidThermo which instantiates either a psiThermo or rhoThermo according to the 'type' specification in thermophysicalProperties. For details, see https://develop.openfoam.com/Development/OpenFOAM-plus/commit/655fc7874808927d14916307a2230a8965bdb860

References:

Implementation of the Thermodynamics: The case of hePsiThermo

definition of the class: https://github.com/OpenFOAM/OpenFOAM-7/blob/109ba3c8d53a8fa1e97854b168a16970aaeea951/src/thermophysicalModels/basic/psiThermo/hePsiThermo.H: inherits from heThermo public member data this->psi_ private member function calculate: called from public member function correct; loops over cells public member function correct(): called from EEqn.H in solver;

Section 14: Turbulence Models

wiki on turbulence modeling: https://en.wikipedia.org/wiki/Turbulence_modeling

cfd-online on the realizable k-epsilon model: https://www.cfd-online.com/Wiki/Realisable_k-epsilon_model

overview of turbulence models in OpenFoam (very nice): https://www.cfd-online.com/Forums/openfoam-solving/218879-latest-developed-turbulence-models.html

overview on turbulence modeling by Wolf Dynamics: slides

boundary condition settings: https://www.openfoam.com/documentation/guides/latest/doc/guide-turbulence-ras-k-epsilon.html

effect of compressibility on the turbulence modeling: dilation dissipation;

References: Fluent Manual Section 14

Wall functions for k (turbulent kinetic energy), eps (turbulent dissipation), nut (turbulent kinematic viscostity) and T (temperature)

wall function for k in OpenFOAM: two options are available. (1/2 )A low Reynolds number in kLowReWallFunction and (2/2) an alternative kqRWallFunction known to work for both low and high Reynolds numbers. Low or high Re is here to be intended actually in reference to the grid size, thus referring to the y+ of the near wall cells. Low Re means that the near wall cell is fine enough to give low y+ values (resolve all flow features at the wall, so that the linear interpolation is enough. High Re (or NON Low Re) means that the near wall cell is not fine enough to do the above (and the WF are needed). Low Reynolds is based on zero-gradient assumption used in master thesis Ali Kadar and Dimas Fikri. This is correct as both in the logarithmic region (y+ above 30 or so) and very close to the wall (k goes to 0 at wall with a larger than 1 power)): https://github.com/OpenFOAM/OpenFOAM-dev/tree/8ca408bf6c83999c71db180f5b2f306472f20b5c/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kLowReWallFunction The k wall function implements a zero-gradient boundary condition. The kqRWallFunction requires setting a value of k on the wall. In the Pitz-Daily test case this value is set equal to 0.375 m^2/s^2; this is based on an inlet/outle conditions; in case of outlet at the inlet, and fixed inlet value is set;

wall function for epsilon in OpenFOAM: two options are available: see https://github.com/OpenFOAM/OpenFOAM-dev/tree/8ca408bf6c83999c71db180f5b2f306472f20b5c/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction In the Pitz-Daily test case the epsilonWallFunction is choosen and the value equal to 14.855 m^2/s^3 is set as value on the boundary;

wall function for nut (eight options, tabulated requires large coding effort leaving 7, two rough wall options not applicable leaving 5, Spalding is an all y+ wall function, while the remaining ones are piecewise valid. Based on the Pitz-Dally tutorial, we will use nutkwallfunction): https://github.com/OpenFOAM/OpenFOAM-dev/tree/8ca408bf6c83999c71db180f5b2f306472f20b5c/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions The basic difference between nutUWallFunction and nutkWallfunction is the computation of y+. In nutkWallFunction it is computed from the turbulent kinetic energy and the constant c_mu, whereas in nutUWallFunction y+ is simply read in.

wall function for T (temperature) (two options of which only one uses the Jayatilleke model. We will use that one): https://github.com/OpenFOAM/OpenFOAM-dev/tree/8ca408bf6c83999c71db180f5b2f306472f20b5c/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions

general info on wall functions in OpenFOAM can be found here: https://github.com/OpenFOAM/OpenFOAM-dev/tree/8ca408bf6c83999c71db180f5b2f306472f20b5c/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions

Section 15: Incompressible Flow Solvers - Stationary Case

Section 1.15: Mathematical Model

Stationary incompressible laminar or turbulent flow. See Navier-Stokes equastions.

Section 2.15: Spatial Discretization, Non-Linear and Linear Solve

SIMPLE iteration (repeated substitution for linearization (requires relaxation) and Schur complement for pressure-velocity solve):

  1. wiki SIMPLE algorithm
  2. Chapter xx of Malaseekera and Versteeg

Section 3.15: Implementation in OpenFoam

  1. Chapter 15 of Darwish and Mangani
  2. OpenFoamwiki SimpleFoam: incompressible solver for pressure (pEqn), velocity (UEqn) and turbulent quantities (kEqn and epsEqn) inside of runtime loop; one pressure andf one velociy solve per SIMPLE iteration simpleFoam

pressure-velocity coupling in simpleFoam: Rhie-Chow interpolation (stable formulation that avoids checkerboarding): user has no immediate access to it: PhD thesis of Jassak: H-operator and report by Fabian Peng Karrholm at Chalmers: off-diagonal blocks in the pressure-velocity coupling; A-operator: diagonal blocks in the pressure-velocity coupling; https://www.cfd-online.com/Wiki/Velocity-pressure_coupling

pressure in simpleFoam (pEqn):

velocity in simpleFoam (UEqn):

turbulent quantities model (kEqn) and (epsEqn):

Section 16: Incompressible Flow Solvers - Transient Case

PIMPLE iteration (SIMPLE and PISO combined): Chapter 15 of Moukalled-Darwish-Mangani: confusion between PRIME and SIMPLE algorithms

pimpleFoam: transient solver; both laminar and turbulent; outer runtime loop; one velocity solve in case that momentum predictor is set to true (default); number of correctors pressure solves;

pisoFoam: equal pimpleFoam with one (single) outer corrector; pimpleFoam can be run essentially as pisoFoam by setting nOuterCorrectors to 1. This alleviates any coding problem. pimpleFoam supports adjustable timestep. https://www.cfd-online.com/Forums/openfoam-solving/130436-using-adjustabletimestep-pisofoam.html

Does increasing the number of outer iterations aid in the converge by allowing a larger step (or in any other way)?

For the first part of the question - pimple or piso - well, since you are running LES, you'll be ensuring Co<1, probably Co~0.3. In which case you don't get any benefit from the additional pressure correctors that allow the pimple solver to stretch the Co number to 10 and beyond. So, you basically want to be running piso. - Having said that - you can run pimpleFoam in "piso mode" by just limiting nOuterCorrectors to 1 (ie one pimple iteration per time step) and nCorrectors to 2 (ie two inner loop pressure correctors). That will then give you access to the other features of the pimpleFoam solver. See https://www.cfd-online.com/Forums/openfoam-solving/233906-large-eddy-simulation-pisofoam-pimplefoam.html

momentumPredictor predictor setting inside of SIMPLE and PIMPLE

what is the momentum correction? The momentum predictor is a first approximation of the velocity field obtained from the solution of a momentum equation assembled with the previous time-step pressure field. This velocity doesn't fulfill the restriction given by the continuity equation which is achieved later by the PISO method. Sometimes you can avoid the solution of the momPred saving time, but it can lead to a bad convergence in the PISO loop setting equal yes (computationally expensive, more stable in convergence): in general momentumPredictor should be set to no, especially for "challenging" test cases (i.e if there is problems with convergence or sharp volume fraction gradient). In this case the momentum equation is only constructed and used to obtain fluxes and solve the pressure equation;
setting equal no (computationally cheap, less stable in convergence): in special cases: in non-challenging steady-state calculations solving the momentum equation might speed up convergence;

nCorrectors vs. nOuterCorrectors

  • explained on the openfoamwiki on PimpleFom;
  • does overrelaxation play any role in case that nOuterCorrectors >= 2 is used?
  • if nOuterCorrectors = 1, then PIMPLE = PISO, and maxCo <= 1 should be used. If nOuterCorrectors >= 2, then maxCo >=1 can be used. Balance the time step and the number of pressure solves

Pitz - Dally pimpleFoam Test Case

Section 17: Radiative Heat Transfer

General Considerations

RTE equation, emission, absorption, scattering, grey model, optical thickness, introduction by Modest-Haworth References: radiation models own OpenFoam: https://www.openfoam.com/documentation/guides/latest/api/group__grpRadiationModels.html and

References:

P1 Model Description

in general: first harmonic in a spherical decomposition of the radiative intensity; works well for combustion applications where optical thickness, tau is large, i.e. tau = a*L > 3 (a = absorption coefficient, L = distance between objects)

domain equation: elliptic equation for the total emissivity G (need to write down the definition of G and the physical units of G; need to write down diffusion equation for G the equation assuming no scattering thus \sigma_s=0);

boundary condition:

  • for walls: Marshak boundary condition (need to write down the equation assuming unit refractive index thus n=1) (need to write down limit cases in which the boundary emissivity is equal to zero and one, respectively);
  • for inlet and outlet: black body absorption: wall emissivity equals one
  • for symmetry plane: symmetry boundary condition: homogeneous Neumann boundary conditions (zero flux)for G;

radiative flux term computation: computation of radiation::Sh() using Ru() and Rp(). Detailed expolanation is here Discussion between Tobias Holtzman and Chris85.

References:

P1 Model in OpenFoam

case set-up in OpenFoam:

  • constant folder: see folder-constant/radiationProperties for domain absorption coefficient and solverFreq and folder-constant/file-boundaryRadiationProperties for boundary emissivity;
  • system folder: see fvSchemes and fvSchemes for settings of G
  • zero folder: boundary conditions for G on inlet, outlet, walls and symmetry plane;

Default-File: the wall emissivity is set as follows
boundaryradiationProperies;
".*" { type greyDiffusiveRadiation; value uniform0; }

boundaryRadiationProperties: ".*" { type lookup; emissivity 1; absorptivity 0; }

Line 179 of the file https://develop.openfoam.com/Development/openfoam/blob/OpenFOAM-v2012/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C shows that the greyDiffusiveRadiation takes emissivity values from boundaryRadiationProperties; implementation:

post processing: radiative flux q_r = \Gamma \nabla G;

References:

fvDOM Model Description

in the fvDOM model one defines the number of rays to discretize a quarter of a sphere (discretization of the solid angle); nPhi specifies the number of rays over 90° (pi/2); nTheta defines the number of rays over 180° (pi) (quarter of a sphere); when you want an even discretization over the quarter sphere you should define nTheta=2nPhi. For a good discretization over the hole sphere, keep the symmetry in mind; the complete number of rays results as nRays=4nThetanPhi; nTheta = 2nPhi (depending on the dimension)

transport equation for individual rays, following a ray as it emits, absorbs and scatters in the medium (need to write down the equation for a ray assuming no scattering thus \sigma_s=0);

boundary condition for inlet, outlet, walls and symmetry condition (requires further study);

  • for walls;
  • for inlet and outlet: black body absorption: wall emissivity equals one
  • for symmetry plane;

post processing for the radiative flux by a weighted sum of the individual rays: what kind of quadrature is used?

References:

fvDOM Model in OpenFoam

case set-up in OpenFoam:

  • constant folder: see folder-constant/radiationProperties for domain absorption coefficient and solverFreq and folder-constant/file-boundaryRadiationProperties for boundary emissivity;
  • system folder: see fvSchemes and fvSchemes for settings of IDefault;
  • zero folder: boundary conditions for Default on inlet, outlet, walls and symmetry plane (Qr is unsed in the BC in the viewFactor model?);

Question on usage in OpenFoam: does fvDOM with symmetry boundary conditions work: Bug report;

References:

Radiation Submodels for absorption coefficient

greyMeanAbsorptionEmission radiation submodel :

Averaging over various grey models in the mixture. For each of these gray components, either P1 or fvDOM needs to be applied;

In this model the absorption coefficient is a volScalarField named as aCont. See L174 of this header file and L203 of this header file. What does this model implement? How is thne averaging implemented?

wideBandAbsorptionEmission radiation submodel : wideBandAbsorptionEmission.H provides the example of "bandLimits (1.0e-6 2.63e-6);" if 1.0e-6 is in metre, the band Limits would be 1 to 2.63 micrometres. I am not quite sure if it is a good estimative for the absorption bands of CO2 and H2O, but I think it is OK.

Turbulence-Radiation Interaction

What if hypothesis of grey medium no longer suffices? Extension of PDF model (for thermo-dynamical properties of the gas-mixture) to radiative heat transfer. See papers by Dirk.

References

Section 18: Compressible Flow Solvers

code in pEqn.H in compressible/RhoPimpleFOAM fvc::ddt(rho) + psi*correction(fvm::ddt(p)

  • dvc::div(phiHbyA) == fvOptions(psi, p, rho.name())

Initializing compressible solver from the incompressible solver

  1. Change the p file. You have to change the dimension, look at one p file in tutorial. You have also to add the atmosferic pressure.
  2. Add the T file
  3. Add transonic yes (if it is transonic) in your FvOption - Simple
  4. Add a small value for the relaxation factor of rho.
  5. You have to change the numerics in fvSchemes. Try and read the error in openfoam. see also: convert incompressible to compressible: https://www.cfd-online.com/Forums/openfoam-solving/99346-quick-way-change-incompressible-compressible.html

Limiting the pressure (p) typically applied for compressible solvers in OpenFoam. This is what Hrvoje Jason says himself: “What you are seeing is that compressible flows are much more sensitive to solver setting, initial field, time-step and boundary conditions.” See post Item 3 at https://www.cfd-online.com/Forums/openfoam-solving/60001-pressure-instability-rhosimplefoam.html Implementation in OpenFoam: rhoSimpleFoam uses “bool pLimited = pressureControl.limit(p);” in pEn.H. See line 90 of https://develop.openfoam.com/Development/openfoam/-/blob/OpenFOAM-v1812.200312/applications/solvers/compressible/rhoSimpleFoam/pEqn.H pressureControl is defined in OpenFOAM-v1906/src/finiteVolume/lnInclude/pressureControl.H and the corresponding *.C file. The online repositories do not show these source files. Are these files created in the installation process? The essence is that pMaxFactor (pMinFactor) bring a large (small) value for pressure within range. For future reference: we found this file using “for i in find . -name *.H ; do echo $i; grep -i pressureControl $i ; done” in the openFoam source tree; a sample usage is provided by the rhoSimpleFoam tutorial squareBendLiq. pMinFactor and pMaxFactor are specified as entries in the SIMPLE dictionary provided in the system/fvSolution file. the line “bool pLimited = pressureControl.limit(p);” is not used in the TonkomoLLC solver. We could simply add it and re-compile the code; https://github.com/TonkomoLLC

Limiting the density (rho) setting on rhoMin and rhoMax: see tutorial heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation in file system/bottomAir verify that bounds are effectively being used by the solver: see file pEqn.H in directory /multiRegionReactingFoam/OpenFOAM-v1812/multiRegionReactingFoam/multiRegionSimpleReactingFoam/fluid in line 53-54

Setting the density on the inlet patches Whether required or not, depends on the thermodynamics used?

SIMPLE { … rhoMin 0.2; rhoMax 2; }

Test Case for Compressible Turbulent Channel Flow

  1. G. N. Coleman, J. Kim andR. D. Moser, A numerical study of turbulent supersonic isothermal-wall channel flow
  2. Large eddy simulation of compressible channel flow

Section 19: Adding Chemistry

  • chemistry reader (CHEMKIN or CANTERA);
  • chemkinToFoam to convert chemkin files to the required format; chemkinToFoam sample run: chemkinToFoam twostepmethanemech.dat thermo30.dat transportProperties reactions thermo.compressibleGas

Section 20: Reactive Flow

solve for pressure (pEqn), velocity (UEqn), energy (EEqn) and species (YEqn): closure by equation of state;

energy equation in terms of total enthalpy that is sum of sensible enthalpy (integration over temperature) and formation enthalpy (extracted for enthalpy of formation from reaction mechaninism). See documentation of the thermodynamics.

convection-diffusion-reaction equation for each of the chemical species;

mass conservation is employed by not solving for the non-reactive species: https://www.cfd-online.com/Forums/openfoam-solving/229789-dumb-question-about-species-fractions-reactingfoam.html no wall function required (although not clear why);

expression for the heat source (Qdot) of the combustion;

question: how does wFuel_ enter the transport equation for the individual equation for the fuel and oxidizer;
unit Lewis number approximation and beyond (differential diffusion) (turbulent) Smith number, Damholker number, Karlowich number; the Chomiak and Karlsson conference paper this presentation gives a more complete discussion on what the PaSR model is conceived. Here is the DOI: https://doi.org/10.1016/S0082-0784(96)80088-9 rhoReactingFoam; https://www.cfd-online.com/Forums/openfoam-solving/226787-problem-nocombustion-rhoreactingfoam.html counterFlowFlame2D rhoReactingCentralFoam: an OpenFOAM solver that combines rhoCentralFoam and rhoReactingFoam for high-speed reactive flows. Tested extensively for ability to model detonations: https://github.com/duncanam/rhoReactingCentralFoam chemkintoFoam conversion: https://www.cfd-online.com/Forums/openfoam-pre-processing/224143-chemkintofoam-parsing-issues-v19.html#post757584 cfd-online: https://www.cfd-online.com/Wiki/Combustion#The_Damk.C3.B6hler_Number cad-online wiki on combustion: https://www.cfd-online.com/Wiki/Combustion

EDM Combustion Model

EDC Combustion Model

FlameletFoam

Flamelet Generated Model

Thickened Flame Model

Section 21: Conjugate Heat Transfer

case set-up option (1/2): mesh complete geometry and run splitMeshRegion: run “splitMeshRegions -cellZonesOnly -overwrite” to decompose the mesh into regions. Creates directories system/gas, system/solid, constant/gas/polyMesh, constant/solid/polyMesh, constant/cellToRegion, 0/gas, 0/solid, 0/cellToRegion. Files are system/gas/fvSchemes and system/gas/fvSolution.

case set-up option (2/2): mesh solid and gas domains seperately and perform a join

solver for conjugate heat transfer by Eric Daymo, see https://github.com/TonkomoLLC/multiRegionReactingFoam tutorial: chtMultiRegionSimpleFoam: http://openfoamwiki.net/index.php/Getting_started_with_chtMultiRegionSimpleFoam_-_planeWall2D definition of patches and boundary conditions: https://www.cfd-online.com/Forums/openfoam-solving/172293-coupling-patches-chtmultiregionsimplefoam.html

decoupling the time step in the solid and fluid domains: http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2016/YuzhuPearlLi/final_report_Jan2017.pdf

run checkMesh using checkMesh -allTopology -allGeometry

chtMultiRegionFoam does by default take into account source terms of radiative heat transfer and combustion. The thermodynamics is limited to incompressible; chtRhoMultiRegionFoam (developed by Franjo) takes thermodynamics for higher Mach numbers into account;

Creating mesh by extrusion of a single patch

Use extrudeMesh, mergeMesh and splitMeshRegions: three utilities require further documentation

Setting up the mesh file

in file constant/meshXX/boundary bottomAir_to_leftSolid { type mappedWall; nFaces 130; startFace 4680; sampleMode nearestPatchFace; sampleRegion leftSolid; samplePatch leftSolid_to_bottomAir; } sampleRegion and samplePatch are the neighbour region and patch

Setting up the boundary conditions

type compressible::turbulentTemperatureCoupledBaffleMixed; Tnbr T; kappa fluidThermo; kappaName none; value uniform 1873=[];

Tutorial Videos on Conjugate Heat Transfer

  1. DD Fluids Conjugate Heat Transfer in OpenFoam

Section 22: Pollutant Formation

NoxFOAM for thermal NOX using Zeldovic mechanism

Command-Line interface (CLI): https://cfd.direct/openfoam/user-guide/v6-post-processing-cli/ and scalartransport within this interface: https://www.openfoam.com/documentation/cpp-guide/html/guide-fos-solvers-scalar-transport.html

scalarTransportFoam: https://openfoamwiki.net/index.php/ScalarTransportFoam

NOx post-processing: https://github.com/nae9on/OpenFOAM_MSc_Thesis_Project (including an example); CFD Online post: https://www.cfd-online.com/Forums/openfoam-solving/60296-nox.html;

No vs. CO formation, residence time, returning of NOx

wiki on Zel’dovich: https://en.wikipedia.org/wiki/Zeldovich_mechanism

Section 23: Test Case

Motor Bike Tutorial

Sequential (using GAMG) and parallel (not using GAMG) runs;

Sandia Flame D

Sandia FlameD

Reverse Burner Test Case

Radiative and conjugate heat transfer;

DLR Confined Jet High Pressure Combustor

DLRCJH

Section 24: Code Constructs

fvOptions for Source Terms

E.g. ignition source

fvOptions for Scalar Transport

fvOptions using Object Registry

Object registry for the implementation of flamelet model for combustion.

Section 25: Criticisms against OpenFoam

Section 26: Conclusions

Section 27: References







  

    

      

  
About


      

        No description, website, or topics provided.
      



    
Resources

    

      
        
        Readme
    


  





  
  

  

    

      
        Activity    



    
Stars

    

      
        6
        stars    


    
Watchers

    

      
        1
        watching    


    
Forks

    

      
        2
        forks    



    

      
          Report repository
    




    

  


  
      

        

          

  Releases


    
No releases published


        

      


  
  
  
      

        

          

  
  

  Packages
      



      

        

        
 

      

      

        

        
 

      

      

        

        
 

      




  

        

      


  
  
    

      

        
  
  

    Contributors  


  
    
      
  • 
        
    
      
    • 
      
  • 
        
    
      
    • 
      
  • 
        
    
      
    • 
  


  

      

    


  
      

        

          
Languages