Open-source CFD toolbox: OpenFOAM is now a mature open-source CFD program with reliability matching that of commercial products. This page describes summary of utilities and dictionaries used in OpenFOAM meshing and visualization such as (blockMesh, snappyHexMesh, ParaView) and OpenFOAM CFD codes and pre-processors such as (simpleFoam, pimpleFoam, engineFoam...).
How many of these names have you heard of? Click the names to know more about them!
This open-source CFD code is released under GPL (GNU Public License), the product is free, powerful but comes with its own 'inconvenience'. Though the authors and programmers try their best to document everything, it still requires additional things to be learnt other than CFD itself, such as:
OpenFOAM, by its design is meant for Linux platform. However, there are few companies which provide installers for Microsoft Windows free of cost (other than providing other paid services). Some of them are:
A set of tutorials explained in a step-by-step manner can be found in this pdf file.
Few troubleshooting and best practices can be found out at:The workflow that require interaction with other programs such as data exchange for structural simulations shall require data interpolation and conversion as per the format required by them. For example, data export to ABAQUS or ANSYS format may require conversion of OpenFOAM data directly to these programs or use of an intermediate format (such as FLUENT or I-DEAS) which ABAQUS or ANSYS can read.
One of the complexity a new user faces is the files (with no pre-defined extensions, most of us are so used to now) and folder-structure used by OpenFOAM. A bare minimum requirement for geometry, assuming that OpenFOAM case is stored and run in a folder named "BaseDir", then the following items must exist, all with read and write access:
There are few additional files which is just to make simulations easier:
"(left|right)" { type zeroGradient; }
Step-1 | Commercial Program | OpenFOAM |
Pre-processing | Geometry creation in CAD, import into CFD simulation program as IGES, STEP, PARASOLID | Create geometry inside OpenFOAM using blockMesh utility, create geometry in any 3D CAD environment, import geometry in STL for meshing using snappyHexMesh (sHM) utility. sHM requires other utilities such as surfaceCheck which can check the input STL geometry and prints overall size (bounding box) of the geometry. |
Mesh generation in a meshing program such as ICEM CFD, ANSA, GridPro, PointWise and import into a solver dependent format such as *.msh, *.cfx5, *.ccm | Generate mesh using blockMesh, generate mesh in any other program (GMSH) and use converter utilities such as gmsh2foam, fluentMeshToFoam | |
Naming of zone for easy selection and application of boundary condition is done in mesh generation environment | The method is equally applicable to OpenFOAM as boundary condition needs to be applied on a set of contiguous faces. |
Reading surface from "constant/triSurface/Fan_in_Duct.stl" ... Statistics: Triangles : 23474 Vertices : 11749 Bounding Box : (-550 -185 -185) (540 185 185) Region Size ------ ---- patch0 23474 Surface has no illegal triangles. Triangle quality (equilateral=1, collapsed=0): 0.00 .. 0.05 : 0.1744 0.05 .. 0.10 : 0.1612 ... min 4.76e-06 for triangle 19519 max 0.999899 for triangle 11050 Edges: min 0.00290 for edge 352 points (-22.92 108.37 143.69)(-22.9 108.37 143.7) max 500.073 for edge 101 points (-550.0 152.62 104.55)(-50.0 157.29 97.39) Checking for points less than 1e-6 of bounding box ((1090 370 370) metre) apart. Found 0 nearby points. Surface is not closed since not all edges connected to two faces: connected to one face : 24 connected to >2 faces : 8 Conflicting face labels:48 Dumping conflicting face labels to "problemFaces" Paste this into the input for surfaceSubset Number of unconnected parts : 2 Splitting surface into parts ... Writing zoning to "zone_Fan_in_Duct.vtk"... writing part 0 size 1084 to "Fan_in_Duct_0.obj" writing part 1 size 2239 to "Fan_in_Duct_1.obj" Number of zones (connected area with consistent normal) : 10 More than one normal orientation.
Step-2: Thermophysical Models, Transport (Viscous) Model and Boundary Conditions - the definition of material properties, and selection of turbulence model are done through the files stored in folder 'constant' namely 'transportProperties' and 'turbulenceProperties'. The application of boundary conditions for velocity is done through file located in folder '0' as explained above. Thermophysical models models are concerned with energy, heat and physical properties. The transport modelling evaluates dynamic viscosity, thermal conductivity and thermal diffusivity (for internal energy and enthalpy equations). The thermodynamic models are concerned with evaluating the specific heat capacity from which other properties are derived. All models are located under $FOAM_SRC/thermophysicalModels folder.
hePsiThermo: General thermophysical model calculation based on compressibility: ψ = 1/(RT) and applicable only to gases. hRhoThermo: General thermophysical model calculation based on density (applicable to gases, liquids, solids). hSolidThermo: applicable only to solids.
mapFieldDict:
patchMap ( w-x1 w-x2 ); cuttingPatches ( fixedWalls );
The first line specifiy patches that have different names in the cases. The second line implies that the 'fixedWalls' patch is cutting through the source case.
Non-Newtonian models: Herschel-Bulkley type of material poses an initial resistance to flow. Once the level of yield stress τ_{0} has been exceeded, the material flows according to the stress-shear relation: τ = τ_{0} + μ_{PL}γ^{n} where μ_{PL} = plastic viscosity and n = shear rate of the material. Division of shear stress by shear rate and density of the material renders the dynamic viscosity, μ = τ / ρ γ^{n}.transportModel CrossPowerLaw; //Newtonian CrossPowerLawCoeffs { nu0 0.01; nuInf 10; m 0.4; n 3; } BirdCarreauCoeffs { nu0 1e-06; nuInf 1e-06; k 0; n 1; } transportModel HerschelBulkley; HerschelBulkleyCoeffs { tau0 0.005; k 0.001; n 0.85; nu0 1.0e5; }
Step-3: Solver setting - files stored in folder 'system' can be used to select a solver, specify convergence criteria, chose relaxation factor, time steps for psuedo-transient or physically transient simulations. Rotating domains can be set using dictionary 'MRFproperties' or 'SRFProperties' under 'constant' directory. There is an attempt to consolidate the solvers and the dynamic-mesh versions such as pimpleDyMFoam, interDyMFoam are being merged to their parent solvers such as pimpleFoam and interFoam in this case.
functions( avgT { functionObjectLibs ("libutilityFunctionObjects.so"); type coded; redirectType average; outputControl outputTime; code #{ const volScalarField& T = mesh().lookupObject("T"); Info<<"T avg:" << average(T) << endl; #}; } );
startTime 0; endTime 100; n 5; ... writeInterval #codeStream { code #{ label dt = ($endTime - $startTime); label i = $n; os << (dt / i); #}; };
foamMonitor: Monitor data with Gnuplot from time-value(s) graphs written by OpenFOAM e.g. by functionObjects. This utility requires gnuplot. Example: foamMonitor -l postProcessing/residuals/0/residuals.dat
functions { #include "minMax" // Print stats }where the minMax file (dictionary) contains:
minMax { type fieldMinMax; // Type of functionObject // Where to load it from (if not already in solver) libs ("libfieldFunctionObjects.so"); log true; // Log to output (default: false) fields ( // Fields to be monitored - runTime modifiable U p ); }
Step-4: Post-processing - ParaView is a native post-processor for OpenFOAM result. Create a file foam.foam in the case folder and read it in ParaView. All the necessary result files will get loaded. Alternatively, use conversion utilities such as foamMeshToFluent and foamDataToFluent to post-process the result in other commercial programs. Use utility postProcess to generate intermediate outputs from OpenFOAM for subsequent post-processing in GNUPlot or ParaView: e.g. postProcess -func sampleDict -latestTime.
Refer in the later sections of this page on dictionaries requires for postProcess utility. All the available utilities can be found in applications/ utilities/ postProcessing/. Just type the command 'app' which will take you to the folder 'applications'.
The list can also be generated by the command 'postProcess -list'. Use the -help flag for more information, as usual for any other application. Examples can be found in folder $WM_PROJECT_DIR/ etc/ caseDicts/ postProcessing.
Sometimes command "postProcess -func wallShearStress -latestTime" shall give following error.--> FOAM FATAL ERROR: Unable to find turbulence model in the database From function virtual bool Foam::functionObjects::wallShearStress::execute() in file wallShearStress/wallShearStress.C at line 200.Command "postProcess -func yPlus -latestTime" shall give following error.
--> FOAM FATAL ERROR: Unable to find turbulence model in the database From function virtual bool Foam::functionObjects::yPlus::execute() in file yPlus/yPlus.C at line 176.
Time = 50 Reading thermophysical properties Selecting thermodynamics package { type hePsiThermo; mixture pureMixture; transport const; thermo hConst; equationOfState perfectGas; specie specie; energy sensibleInternalEnergy; } Reading field U Reading/calculating face flux field phi pressureControl pMax 205740.4 pMin 10000 Creating turbulence model Selecting turbulence model type RAS Selecting RAS turbulence model kOmegaSST Selecting patchDistMethod meshWave RAS { RASModel kOmegaSST; turbulence on; printCoeffs on; alphaK1 0.85; alphaK2 1; alphaOmega1 0.5; alphaOmega2 0.856; gamma1 0.55555556; gamma2 0.44; beta1 0.075; beta2 0.0828; betaStar 0.09; a1 0.31; b1 1; c1 10; F3 false; } No MRF models present Creating finite volume options from "system/fvOptions" Selecting finite volume options model type limitTemperature Source: limitT - selecting all cells - selected 16000 cell(s) with volume 2.0996568 wallShearStress wallShearStress: processing all wall patches wallShearStress wallShearStress write: writing object wallShearStress min/max(aerofoil) = (-8.64 -0.0 -3.32), (-0.751 0.0 3.32)"rhoSimpleFoam -postProcess -func wallShearStress -latestTime" shall calculate the wall shear stress only at the latestTime step.
Function objects to derive additional data both during and after calculations, typically in the form of additional logging to the screen, or generating text, image and field files. This eliminates the need to store all run-time generated data, hence saving considerable resources. Furthermore, function objects are readily applied to batch-driven processes, improving reliability by standardising the sequence of operations and reducing the amount of manual interaction. Data derived from function objects are written to the case 'postProcessing' directory, typically in a subdirectory with the same name as the object. E.g. $FOAM_CASE/ postProcessing/ 'functionObjectName'/ 'time'/ 'data'. The function objects are options that are run time selectable. The results generated by function objects can be both obtained in the due course of simulation and after the calculation. Source codes can be found in $FOAM_SRC/ OpenFOAM/ db/ functionObjects and $FOAM_SRC/ functionObjects.
Few examples from "Evoking existing function objects and creating new user-defined function objects for Post-Processing" by Sankar Raju. N
Calculate y+ of the specific walls.
yPls { type yPlus; libs ("libfieldFunctionObjects.so"); patches (leftWall rightWall stepWall); }
Compute stress tensor parameter 'R'. One can additionally compute eddy kinematic viscosity (nut), effective kinematic viscosity (nuEff), deviatoric stress tensor (devReff), eddy dynamic viscosity (mut), effective dynamic viscosity (muEff), thermal eddy diffusivity (alphat), effective eddy thermal diffusivity (alphaEff), deviatoric stress tensor (devRhoReff).
turbFieldsR { type turbulenceFields; libs ("libfieldFunctionObjects.so"); field R; }
Function object Q-criterion calculates and outputs the second invariant of the velocity gradient tensor [1/s^{2}]. Q iso-surfaces are good indicators of turbulent flow structures. In the example below, result is displayed with the name usrQ.
QCrit { type Q; libs ("libfieldFunctionObjects.so"); field U; result usrQ; }To find all of the examples of Function Objects being used in the standard OpenFOAM tutorials, run the following command to get a list of tutorial cases that used them: find $FOAM_TUTORIALS -name controlDict | xargs grep 'functions' -sl
Note: All the incompressible solvers implemented in OpenFOAM namely icoFoam, simpleFoam, pisoFoam and pimpleFoam use "specific engergy [J/kg] or modified pressure P = p/ρ [J/kg] where p is actual static pressure in [Pa]. Hence, during post-processing the results need to be multiplied by the density to get the pressure in [Pa].
For specific applications such as internal combustion engine, constant/engineGeometry dictionary is used which contains settings and geometry data such as bore, stroke, engine speed and minimum valve lift. It also contains information about names of patches or boundary surfaces and cell zones required to be set moving during mesh motion.
Similarly, constant/combustionProperties dictionary contains settings for the combustion calculations such as combustion model parameters, fuel properties and ignition parameters (ignition time, location, duration and strength).
Dictionary and keywords are two important concepts in OpenFOAM. For example:
fvSolution <-- Dictionary file solvers { <-- Dictionaries p { <-- Sub-Dictionaries or Keywords solver PCG; preconditioner DIC; <-- Values assigned to Keywords tolerance 1e-06; <-- Values assigned to Keywords relTol 0; }Most of the applications in OpenFOAM requires instruction from a plain text file called dictionaries. These are equivalent command line operation to the GUI features available in commercial softwares. Though, most of the dictionaries are commented well, the learning curve might still be steep for most of the beginners.
Headers of the dictionaries convey information about class and object (note that OpenFOAM has been written in C++ which is an object-oriented programming language). For example:
FoamFile { version 2.0; format ascii; class volVectorField; //dictionary, volScalarField, volTensorField object U; //p, T, gradU, sigma, k, epsilon, nut, omega, caseProperties, caseSummary //controlDict, sampleDict, fvSchemes, fvSolution, transportProperties ... }
In the example above, object 'U is of class type 'volVectorField' and should always be specified / initialized as vector such as "internalField uniform (0 0 0);". Similarly, an object of class type 'volTensorField' will need to be specified as 3×3 matrix.
A dictionary in OpenFOAM can contain multiple data entries as well as many sub-dictionaries. For a dictionary entry, the name is follow by the keyword entries in curly braces '{'. E.g.
solvers { //Dictionary 'solver' p { //Sub-dictionary 'p' as it is nested inside dictionary 'solver' solver PCG; preconditioner DIC; tolerance 1e-06; relTol 0; } }Macro expansion can be used to duplicate and access variables in dictionaries. For example:
$p; //Create a copy of dictionary 'p' $p.tolerance; //Access variable 'tolerance' in the dictionary 'p'Another example of dictionary nested inside a dictionary:
boundaryField { //Dictionary 'boundaryField' inlet { //Sub-dictionary 'inlet' as it is nested inside dictionary 'boundaryField' type zeroGradient; } outlet { //Sub-dictionary 'outlet' as it is nested inside dictionary 'boundaryField' type fixedValue; value uniform 0; } }Lists or list entries on the other hand are following by parentheses '('. For example:
vertices ( //Name of the the 'list' (0 0 5) (0 5 5) )
The blockMeshDict file is here . The mesh generated with simple biasing over cylinder body is shown below.
The blockMeshDict file is here. The mesh generated for a pipe flow, pipe diameter of 50 [mm] and length of 1000 [mm] is shown below.
The blockMeshDict file is here. A recommended blocking for a quater-symmetric cylindrical domain is shown below.
Note that all the mathematical operations need to be inside double quotes such as "$Lu + $Lb" is correct, $Lu + $Lb is not and it will result in fatal error. As explaine in the link https://bugs.openfoam.org/view.php?id=3034: The #calc functionality needs to be able to compile some code which requires a working build environment. By default, Windows Subsystem for Linux (WSL) does not include such packages. Use command "sudo apt-get install build-essential" to get what is needed for make.
Use of include feature: OpenFOAM dictionaries can read B.C. information stored in other files using directive 'include' as demonstrated below. This is applicable not just blockMeshDict but any other dictionary like U, p, T... For example: initialization of velocity filed, the information stored in a file called 'initCondition' is
flowVelocity (10 0 0); pressure 10000; turbulentKE 0.1; turbulentOmega 0.3; #inputMode mergeThis file has been included in 'U' dictionary as
#include "initConditions" dimensions [0 1 -1 0 0 0 0]; internalField uniform $flowVelocity;Other example of 'include' directive for definition of variables in blockMeshDict is explained below.
Note that -$L1 is not valid, a new variable L2 #calc "0 -$L1" or L2 #calc "-1.0*$L1" needs to be defined for negative values of a variable. This is because C/C++ pre-processor expands the macro variable before compiling '-' reads the minus sign and expects a number to come up and not the $ sign. The sin and cos function needs argument which is to be converted into radians using degToRad(θ). Thus, cos(0.7854) will not work but cos(degToRad(45.0)) will work even though 45° = 0.7854 [rad].
Variables can be accessed within different levels of sub-dictionaries, or scope. Scoping is performed using a ‘. ’ (dot) syntax, illustrated by the following example, where b is set to the value of a, specified in a sub-dictionary called subdict.
subdict { a 10; } b $subdict.a;There are further syntax rules for macro expansions: to traverse up one level of sub-dictionary, use the '.. ' (double-dot) prefix, to traverse up two levels use '... ' (triple-dot) prefix, to traverse to the top level dictionary use the ':' (colon) prefix (most useful), for multiple levels of macro substitution, each specified with the '$' dollar syntax, '{}' brackets are required to protect the expansion.
a 10; b a; c ${${b}}; // returns 10, since $b returns "a", and $a returns 10
dict1 { b $..a; // double-dot takes scope up 1 level, then "a" is available subdict2 { c $:a; // colon takes scope to top level, then "a" is available } }
φ_{B} = f × φ_{REF} + (1 - f) × (φ_{i} + refGrad × δ)
The recommended use of variables to define the boundary conditions can be found in folder $FOAM_SRC/ finiteVolume/ fields/ fvPatchField/ derived.Inlet Boundary Condition from Volumetric / Mass Flow Rate: This method specifies uniform velocity at inlet face. For volumetric heat flow, use the following settings.
inletZone { type flowRateInletVelocity; //Calculate velocity field from volumetric flow rate volumetricFlowRate 1.0; //volume flow rate in [m^{3}/s] value uniform (1 0 0); //Initial value: does not affect the B.C. }To calculate velocity from mass flow rate, the solver needs to know the density. Hence, DataEntry type 'rho' needs to be specified in this case.
inletZone { type flowRateInletVelocity; //Calculate velocity field from volumetric flow rate massFlowRate 1.0; //volume flow rate in [m^{3}/s] rho rho; //Name of density field rhoInlet 1.185; //Density at inlet in [kg/m^{3}] }
inletZone { type cylindricalInletVelocity; //specify that it is in cyl CS. axis (0 1 0); //z-axis of cylindrical CS centre (1 0 0); //centre of CS w.r.t global CS axialVelocity constant 10; //velocity along axial direction radialVelocity constant -5; rpm constant 1000; }
outletZone { type totalPressure; U U; phi phi; //Name of flux field rho rho; //Name of density field psi none; //1/R/T: R = Gas constant gamma 1.4; //Adiabatic exponent = Cp / Cv p0 uniform 1e5; //Total pressure value }
Excerpt from a post on cfd-online.com: "codedFixedValue boundary condition method attempts to load a shared library using case-supplied C++ code, so you need to allow your FOAM to execute user written code. You should switch "allowSystemOperations" in .../etc/controlDict from 0 to 1." The implementation of codeStream for diction 0/U is described below.
Boundary conditions codedFixedValue and codedMixed are derived from codeStream and let you access more information of the simulation database (e.g. time). Another feature of these boundary conditions, is that the code section can be read from an external dictionary (system/codeDict), which is run-time modifiable. The boundary condition codedMixed works in similar way. This boundary condition gives you access to fixed values (Dirichlet BC) and gradients (Neumann BC).
dimensions [0 1 -1 0 0 0 0]; internalField uniform 10; inletZone { //options for type: codedMixed, fixedValue type codedFixedValue; //B.C. - designated by term 'Fixed' value $internalField; //Initial value or use #codeStream redirectType parabolicInlet; //Unique name for this boundary patch // Note about initial value - it is specified so the utility always has // a starting value to work with (e.g. when plotting boundary values using // ParaView at time 0). In this case the "uniform 10" value is substituted // for time '0'. At later times, boundary values are evaluated by solver. // codeStream reads 3 entries: 'code', 'codeInclude' (optional), 'codeOptions' // (optional) and uses those to generate library sources inside codeStream/ // Note the #{ ... #} syntax is a 'verbatim' input mode that allows inputting // strings with embedded newlines. Method to read mesh if different for the // internal field and boundary faces. code #{ const IOdictionary& d = static_cast<const IOdictionary&> ( dict.parent().parent() // access current dictionary ); //Access the mesh database. const fvMesh& mesh = refCast<const fvMesh>(d.db()); //Specify the actual name of the patch: get the label ID const label zoneId = mesh.boundary().findPatchID("inletZone"); //Access the boundary mesh information using label ID 'zoneId' const fvPatch& patch = mesh.boundary()[zoneId]; //Initialize velocity components - patch.size() = number of faces vectorField U(patch.size(), vector(0, 0, 0)); //codeStream recognises regular macro substitutions using the '$' syntax //Create boundary condition without loop vector dir = vector(0, 1, 0); //Get unit vector along x-direction scalarField var = patch().Cf()&dir; //& = dot product, Cf = face centre scalarField value = 10.0*(1.0 - pow(var, 2) / pow(0.5, 2)); vector(1, 0, 0) * value; //Velocity components as vector field const scalar pi = constant::mathematical::pi; V0 = 10.0; //Create boundary condition with loop forAll(U, i) { //for (int i=0; patch.size()<i; i++) const scalar y = patch.Cf()[i][1]; U[i] = vector($V0 * (1 - pow(var,2) / pow(0.5, 2)); } U.writeEntry("", os); //Write values of U to the dictionary #} //Optional - Files needed for compilation codeInclude #{ #include "fvCFD.H" #include <cmath> #include <iostream> #}; //Optional - Compilation options codeOptions #{ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude #}; //Optional - Libraries needed for compilation, required only to visualize //the output of the boundary condition at time zero codeLibs #{ -lmeshTools \ -lfiniteVolume #}; }
Either or "(lt|rt|bt)Wall"{ "*Wall"{ type fixedValue; type fixedValue; value uniform (0 0 0); value uniform (0 0 0); } }
Fan at Inlet
inletZone { type fanPressure; //Name of fan boundary condition fileName "fanPQdata"; //Name of file containing P-Q data outOfBounds clamp; //Instruction for handling out-of-bound values direction in; //Inflow direction p0 uniform 0; //Initial value only. Gets updated as per fanPQdata value uniform 0; }
Fan at Outlet
outletZone { type fanPressure; //Name of fan boundary condition fileName "fanPQdata"; //Name of file containing P-Q data outOfBounds clamp; //Instruction for handling out-of-bound values direction out; //Inflow direction p0 uniform 0; //Initial value only. Gets updated as per fanPQdata value uniform 0; }
Example: Specify a ramped velocity at inlet - U = MIN(5.0 [m/s], 0.01 * t). timeVaryingMappedFixedValue can be used to specify test data (e.g. measured values of U, k and ε). Boundary condition type 'codedFixedValue' constructs on-the-fly a new boundary condition (derived from fixedValueFvPatchField) which is then used to evaluate. A special form is if the 'code' section is not supplied. In this case the code is read from a (runTimeModifiable) dictionary system/codeDict which would have a corresponding entry:
patchName { code #{ operator==(min(10, 0.1*this->db().time().value())); #};When code' section is supplied in the field dictionary 'U':
inletZone { type codedFixedValue; //Description of B.C. type value uniform 0; //Initial value redirectType rampedVelcity; //Name of new boundary condition type code #{ const scalar t = this->db().time().value(); //Get time operator == (min(5.0, 0.01 * t)); #}; }Similarly boundary condition type 'codedMixed' can be used to generate Robin boundary conditions on-the-fly. This is a combination of fixed value and gradient value defined as φ_{FACE} = f * φ_{REF} + (1-f) * [φ_{CELL} + Δ(∇φ_{REF})] where f = value fraction and Δ = distance between face and the centroid of the cell. 'codedFixedValue' and 'mixed' can be combined to get a 'codedMixed' boundary condition type on-the-fly.
inletZone { type mixed; //It tells the solver it is Robin-type B.C. refValue 10.0; //Fixed component refGradient 2.0; //Gradient component valueFraction 0.5; //Fraction f in the above equation }
wallConvective { type externalWallHeatFluxTemperature; //to model only convection Ta uniform 300.0; // Reference temperature q''= h*(T-Ta) h uniform 10.0; // HTC value [W/m^2] value uniform 300.0; // Boundary Initialization Temperature kappa solidThermo; // -kappa x dT/dx = h *(T - Ta) //fluidThermo, solidThermo, lookup, directionalSolidThermo kappaName none; }
Porous Domain - reference tutorial case: incompressible\ simpleFoam\ simpleCar\ system - porous media is modeled by adding a source term Si to the momentum equations. A porosity value 0 < γ < 1 is added to the time derivative. The source term is represented in the Darcy-Forchheimer law as a viscous (linear) loss term and an inertial (parabolic) loss term. The source term can also be modeled as a power-law of the velocity magnitudein the form Si = C_{0}.|u|^{C1}.
porousZoneX { type explicitPorositySource; active true; explicitPorositySourceCoeffs { type DarcyForchheimer; selectionMode cellZone; cellZone porousZone; DarcyForchheimerCoeffs { d d [0 -2 0 0 0 0 0] (5e7 -1000 -1000); f f [0 -1 0 0 0 0 0] (0 0 0); //A local coordinate system belonging to the porous media is created from //the vectors e1 and e2. The pressure loss in the local coordinate system is //then defined with the vectors d and f. coordinateSystem { type cartesian; origin (0 0 0); coordinateRotation { type axesRotation; e1 (1 0 0); e2 (0 1 0); } } } } } //The local cylindrical coordinateSystem is specified by a coordinateRotation //of type localAxesRotation. The user specifies: e3, the axis of rotation (z) //of the local cylindrical coordinate system; and, the origin, from which the //radial direction (r) is constructed by the vector from the origin to local //cell centre. porosityCylindrical { type DarcyForchheimer; cellZone porosity; DarcyForchheimerCoeffs { d d [0 -2 0 0 0 0 0] (0 1e5 0); f f [0 -1 0 0 0 0 0] (0 0 0); coordinateSystem { type cartesian; origin (0 0 0); coordinateRotation { type localAxesRotation; e3 (0 0 1); } } } } // In this case, the viscous momentum sink (0 1e5 0) denotes a value of 1e5 // in the circumferential direction and zero in radial and axial directions.
volumeHeatSource { type scalarSemiImplicitSource; active true; scalarSemiImplicitSourceCoeffs { selectionMode cellZone; //cellZone, cellSet, points, all //Domain where the source is applied cellZone solidDomain; volumeMode specific; //specific | absolute injectionRateSuSp { h (100000 0); /* S(T) = Su + Sp * T -> Su = q, Sp = 0 q''' = 'specific' volumetric generation [W/m^3] = 100,000 q = 'absolute' volumetric generation [W] */ } } }
defaultFieldValues ( volScalarFieldValue alpha.water 0 ); regions ( surfaceToCell { file "./geom/patchZone.stl"; outsidePoints ((0.50 0.75 0.00)); includeInside true; includeOutside false; includeCut false; fieldValues ( volScalarFieldValue alpha.water 1 ); } );
Solution initialization using codeStream
The code section of the codeStream IC in the body of the file alpha.water is as follows,code #{ //Access internal mesh information const IOdictionary& d = static_cast<const IOdictionary&>(dict); const fvMesh& mesh = refCast<const fvMesh>(d.db()); scalarField alpha(mesh.nCells(), 0.); //Initialize scalar field to zero forAll(alpha, i) { //Access cell centers coordinates const scalar x = mesh.C()[i][0]; const scalar y = mesh.C()[i][1]; //Assign values to alpha if (y >= 1.0*cos(0.5*constant::mathematical::pi*x)) { alpha[i] = 1.; } } alpha.writeEntry("", os); //Write output to input dictionary #};
solidZone_to_waterDomain { //type calculated; type compressible::turbulentTemperatureCoupledBaffleMixed; value uniform 300; //"value" entry is only an 'initialization' value at the interface. Tnbr T; kappa solidThermo; kappaName none; /*define thermal contact resistance between regions by giving the thickness and the conductivity of the inter-region layers. */ thicknessLayers (1e-3); // value in [m] - Rk = [K/W] = [m] / [W/m.K] kappaLayers (5e-4); // Value in [W/m.K] }
barotropicCompressibilityModel linear; psiv 2.5e-06; rholSat 830; psil 5e-07; pSat 4500; rhoMin 0.001;
Excerpts from "6-DOF VOF-solver without Damping in OpenFOAM" by Erik Ekedahl: OpenFOAM uses quaternion and septernion that are objects to store translation or rotation and when calling return a tensor for rotation or a vector for translation. Quaternion is the object for handling the rotation. A quaternion is a mathematical object, a vector in four-dimensional vector space, written as a linear combination of the basis elements {1, i, j, k} as (a.1 + b.i + c.j + d.k) where a, b,c and d are real numbers and i^{2} = j^{2} = k^{2} = ijk = -1 as found by Hamilton. Here a is called the scalar part and {b.i + c.j + d.k} is called the vector part. The vector part of the quaternion can be identified to be the same as an element of R^{3} vector space.
The tutorial multiphase\compressibleInterDyMFoam\laminar contains two cases (sloshingTank2D and sphereDrop) on dynamic mesh with multiphase flows.
Other utilities for moving mesh are:fieldFunctions
functionObjects are general libraries that can be attached run-time to the solver without having to re-compile the solver. A functionObject is added to a solver by adding a functions entry in system/controlDict dictionary. E.g. The fieldMinMax functionObject can be found in damBreakLeft/system/controlDict. Output in stored in time directory according to the time when it was initialized. Method to develop your own fieldFunction is described in this tutorial.functions { minMaxU { type fieldMinMax; libs ("libfieldFunctionObjects.so"); enabled true; log false; write true; fields ( U ); mode magnitude; writeControl timeStep; writeInterval 1; } }
postProcessing Utilities and functionObjects
While the results from OpenFOAM can be post-processed both quantitatively and qualitatively in ParaView or many other post-proccessors, there are some postProcessing utility available inside the solver. This will help avoid the interpolation error as values inside OpenFOAM are calculated as same discretization and interpolation scheme as the solution of field variables. This feature can be accessed by adding object 'function' in 'controlDict' dictionary. Available options are - value of attribute 'type' in function examples described below. Note that this utility supersedes foamCalc utility available in older version.functionObject for Pressure Drop between two zones
deltaP { type fieldValueDelta; functionObjectLibs ("libfieldFunctionObjects.so"); enabled true; region1 { writeFields off; type surfaceRegion; regionType patch; name inlet; operation sum; fields ( phi ); } region2 { writeFields off; type surfaceRegion; regionType patch; name outlet; operation sum; fields ( phi ); } operation subtract; }
functionObject for mass flow rate at a boundary patch
outletMdot { type surfaceRegion; functionObjectLibs ("libfieldFunctionObjects.so"); enabled true; //writeControl outputTime; writeControl timeStep; writeInterval 1; log true; writeFields false; regionType patch; name outlet; operation sum; fields ( phi ); }
The fieldAverage functionObject calculates the time-average of specified fields and writes the results in the time directories. It needs to be added to the functions dictionary in controlDict. Example can be found in pisoFoam /les /pitzDaily case.
fieldAvgU { type fieldAverage; libs ("libfieldFunctionObjects.so"); enabled true; writeControl outputTime; fields ( U { mean on; prime2Mean on; //RMS base time; } ); }
#include "ExtFunctionObjects" can be used to include to call an external dictionary with the functionObjects definition.
postProcess -func surfaces is used to extract surfaces in VTK (ParaView) format. Dictionary file example can be located in folder $WM_PROJECT_DIR /etc /caseDicts /postProcessing /visualization /surfaces
postProcess can also be used to calculate new fields from existing ones. E.g. postProcess -func 'div(U)' -case pitzDaily
Example to calculate lift and drag coefficients are given below:
functions { forces { type forceCoeffs; //postProcessing -list: available variables libs ("libforces.so"); writeControl timeStep; //As multiple of timesteps writeInterval 1; //Every timestep patches (cylWall); //Walls to be used for force calc log true; rhoName rhoInf; //Indicates incompressible rhoInf 1.185; //Ref. density - not used in incompressible flows CofR (0 0 0); //Centre of Rotation liftDir (0 1 0); //Direction of list:Y-axis in this case dragDir (1 0 0); //Direction of drag: X-xis in this case pitchAxis (0 0 1); //Pitching axis: Z-axis in this case magUInf 50.0; //Reference velocity magnitude lRef 0.50; //Reference length scale Aref 1.00; //Reference area } }C_{D} = F_{D} / [2 * A_{REF} * ρ_{INF} * V_{INF}^{2}]. The reference length scale can be 'diameter' in case of flow over a cylinder, chord length in case of flow over an airfoil, length of the car along flow direction in case of external aerodynamics over a car. Similarly, A_{REF} is usually D×L in case of flow over a cylinder, maximum height * maximum width in case of a car.
Similarly, to get residuals for any field variable, following 'function' object needs to be added inside 'controlDict' dictionary.
functions { residuals { type residuals; //postProcessing -list: available options functionObjectLibs ("libutilityFunctionObjects.so"); enabled true; outputControl timeStep; outputInterval 5; //Every 5 timestep fields (p U k epsilon T); } }
Two other examples to write value of pressure at two 'probe' locations and average pressure of a wall surface can be found in examples multiphase / interDyMFoam / ras / sloshingTank2D / system / controlDict and incompressible / pisoFoam / pitzDaily / system / controlDict. The values can be used to plot as monitor points during and/or after the simulation.
The probes functionObject probes the development of the results during a simulation, writing to a file in the directory postProcessing/probes.
The probes function can be copied to the local system directory using "foamGet probes" and configured as per user's need. The 'singleGraph' function samples data for graph plotting. The singleGraph file should be copied using 'foamGet' into the system directory to be configured (example in pitzDaily case). It creates a folder 'singleGraph' inside folder postProcessing and saves x-y data at each time directory. The files are named as line_U.xy, line_p.xy, line_T.xy ...
functions { probeSolid { //Name of 'probe' block type probes; libs ("libsampling.so"); writeControl writeTime; regions solid; //Only if multiple zones such as in CHT probeLocations ( (0 9.95 19.77) (0 -9.95 19.77) ); fixedLocations false; fields( p Ux ); } //When solver is run, time-value data is written into p files in postProcessing/probes/0. probeFluid { type probes; libs ("libsampling.so"); writeControl writeTime; //enabled true; //writeControl timeStep; //writeInterval 1; region fluid; //Only if multiple zones such as in CHT probeLocations ( (0 5.55 5.555) (0 -5.55 5.555) ); fixedLocations false; fields ( p U T Ux Uy Uz ); } wallPressure { type surfaces; //On surfaces defined below libs ("libsampling.so"); writeControl writeTime; surfaceFormat raw; fields (p); interpolationScheme cellPoint; surfaces ( walls { type patch; patches (walls); triangulate false; } ); } }
Content inside the dictionary 'system/singleGraph'.
singleGraph { start (0.0 0.0 0.0); end (0.1 0.0 0.0); fields (p U T); #includeEtc "caseDicts/postProcessing/graphs/sampleDict.cfg" //Some installation does not recognize the caseDicts path //Use "${WM_PROJECT_DIR}/etc/caseDicts/postProcessing/graphs/sampleDict.cfg" setConfig { axis distance; //axis y; } // Must be last entry #includeEtc "caseDicts/postProcessing/graphs/graph.cfg" //"${WM_PROJECT_DIR}/etc/caseDicts/postProcessing/graphs/graph.cfg" }The variables p, U and T are extracted along a horizontal line at 100 points and the results are written in <case folder> /postProcessing /singleGraph.
The formatting of the graph is specified in configuration files in $FOAM_ETC /caseDicts /postProcessing /graphs. The sampleDict.cfg file in that directory contains a setConfig sub-dictionary.
'axis' specifies how to write point coordinates. The options are:sampleDict
The sampleDict file contains a sets subdictionary which defines where the fields should be sampled:
sets ( line { // name of the set type uniform; // type of the set axis x; start (-0.05 0 0); end ( 0.10 0 0); nPoints 1000; //nPoints number of sampling points for uniform - default is 100 } );
sets functionObject: This functionObject is used to sample field data on a line. The output of this functionObject is saved in ASCII format in the files line1_p.xy and line1_U.xy located in the time directories inside the folder postProcessing/setsDataOnLine. Here line1 in line1_p.xy is the name of the sets.
setsDataOnLine { type sets; functionObjectLibs ("libfieldFunctionObjects.so"); enabled true; writeControl outputTime; interpolationScheme cellPointFace; setFormat raw; sets ( line1 { type midPointAndFace; axis x; start (-0.002 -0.002 0.20); end ( 0.002 0.002 0.20); } ); fields ( U p ); }
The sampleDict contains entries to be set according to the user needs. The choice for the interpolationScheme can be one of the following. Vertex values are determined from neighboring cell-centre values whereas face values determined using the current face interpolation scheme for the field (linear, gamma ...).
interpolationScheme cellPoint;#includeEtc "caseDicts/postProcessing/visualization/surfaces.cfg" fields (p U); surfaces ( zNormal { $cuttingPlane; pointAndNormalDict { basePoint (0.05 0.05 0.005); // Overrides default basePoint (0 0 0) normalVector $z; // $y: macro for (0 0 1) } } p0 { $isosurface; isoField p; isoValue 0; } movingWall { $patchSurface; patches (movingWall); } );
More examples are:
functions { cuttingPlane { type surfaces; functionObjectLibs ("libsampling.so"); outputControl outputTime; surfaceFormat vtk; fields (U p T); interpolationScheme cellPoint; surfaces ( yNormal { type cuttingPlane; planeType pointAndNormal; pointAndNormalDict { basePoint (0.1 0.0 0.0); normalVector (1.0 0.0 0.0); } interpolate true; } ); } }
Streamlines
functions { streamLines { type streamLine; outputControl outputTime; setFormat vtk; //Options: gnuplot, xmgr, raw, jplot, csv, ensight UName U; trackForward true; fields (p U); //Streamlines can be coloured by scale of pressure lifeTime 1000; nSubCycle 5; cloudName particleTracks; seedSampleSet uniform; //Options: cloud, triSurfaceMeshPointSet uniformCoeffs { type uniform; axis x; //Option: distance start (0 0 0); end (5 0 0); nPoints 20; } } }
sampling
The runtime sampling of data on a line or surface ... can be done using sample utility alongwith sampleDict dictionary. An example of sampleDict dictionary with all relevant comment can be found here.Plots can also be generated with Python and Matplotlib. Example code (reference: Examples of how to use some utilities and functionObjects by Hakan Nilsson, Chalmers / Applied Mechanics / Fluid Dynamics) is given below. Type python plotPressure.py at the command prompt.
#!/usr/bin/env python description = """ Plot the pressure samples""" import os, sys import math from pylab import * from numpy import loadtxt def addToPlots( timeName ): fileName = "cavity/postProcessing/singleGraph/" + timeName + "/line_p.xy" i=[] time=[] abc =loadtxt(fileName, skiprows=4) for z in abc: time.append(z[0]) i.append(z[1]) legend = "Pressure at " + timeName plot(time,i,label="Time " + timeName ) figure(1); ylabel(" p/rho "); xlabel(" Distance (m) "); title(" Pressure along sample line ") grid() hold(True) for dirStr in os.listdir("cavity/postProcessing/singleGraph/"): addToPlots( dirStr ) legend(loc="upper left") savefig("myPlot.jpeg") show() #Problems with ssh
planeType pointAndNormal; pointAndNormalDict { basePoint (0 -1.5 0); normalVector (0 1.0 0); } planeTolerance 1e-3;This resulted in a new mesh with patch name same as original mesh. A topoSetDict was created for topoSet utility to generate a face set for upper cylinder.
actions ( // Make a set of all face elements of patch 'cylinder' { name cylinderFace; //name of the set type faceSet; //pointSet/faceSet/cellSet/faceZoneSet/cellZoneSet action new; //create new, add to existing, delete or subset source patchToFace; //What is type of 'source' for points/faces/cells sourceInfo { //method to specify the source name "cylinder"; } } { name cylinderFace; type faceSet; action subset; source boxToFace; sourceInfo { box (-1.01 -1.01 -1.01)(1.01 1.01 1.01); } } );
The topoSet utility created a faceSet as shown below.
Finally, the createPatch utility was used to create a new patch for upper cylinder.
pointSync false; patches ( { // Patches to create. name cylinder01; // Name of new patch patchInfo { // Type of new patch type patch; } constructFrom set; //How to construct: either from 'patches' or 'set' set cylinderFace; // name of faceSet created by topoSet } );
Other related utilities are:
autoPatch: by default it works on all the boundary (external) faces, no particular faceSet or patch name can be specified.
Step-by-step guide to use FLUENT Mesh in OpenFOAM
The tutorial file incompressible/icoFoam/elbow contains a sample FLUENT *.msh file. It cannot be viewed in ParaView as it is.Step-01: use utility "FluentMeshToFoam elbow.msh" and convert the mesh data into OpenFOAM format. Create a blank file foam.foam in the main folder which will allow reading the mesh converted into OpenFOAM format. Identify the names of the zones that were defined while generating mesh for FLUENT. e.g. frontAnfBackPlanes, pressure-outlet-7, velocity-inlet-5, velocity-inlet-6, wall-4 and wall-8 in this case.
Step-02: Once (FLUENT) zone names are known, define appropriate B.C. on those patches as described below.
In order to define fluid regions (such as in conjugate heat transfer or multiple reference frame cases), mesh manipulation utilities can be exploited. One can also use topoSet utility in conjunction with topoSetDict to create new cell zones.
Step-03: Make the simulation run in OpenFOAM. In order to convert the result files back into FLUENT format, use the utilities foamMeshToFluent and foamDataToFluent. Alterntively, the post-processing activities can be carried out directly in ParaView.
One can use utilities (i)foamToEnsight, (ii)foamToFireMesh, (iii)foamToStarMesh, (iv)foamToTetDualMesh, (v)foamToEnsightParts, (vi)foamToGMV, (vii)foamToSurface and (viii)foamToVTK to convert the OpenFOAM result in any other post-processing application. foamToEnsightParts - an Ensight part is created for each cellZone and patch. VTK - Visualization Toolkit is native format for ParaView.
angle 45.0; laplacianSchemes { #if #calc "${angle} < 75" default Gauss linear corrected; #else default Gauss linear limited corrected 0.5; #endif }The #ifEq compares a word or string and executes based on a match.
ddtSchemes { #ifeq ${FOAM_APPLICATION} simpleFoam default steadyState; #else default Euler; #endif }
List of Boundary Conditions
The recommended use of variables to define the boundary conditions can be found in folder $FOAM_SRC/ finiteVolume/ fields/ fvPatchField/ derived. E.g. the usage of 'flowRateInletVelocity' boundary type as per header file flowRateInletVelocityFvPatchVectorField.H is as follows.patchName { type flowRateInletVelocity; volumetricFlowRate 0.2; extrapolateProfile yes; value uniform (0 0 0); }
cylindricalInletVelocity | This boundary condition describes an inlet vector boundary condition in cylindrical co-ordinates given a central axis, central point, rpm, axial and radial velocity |
fanPressure | This boundary condition can be applied to assign either a pressure inlet or outlet total pressure condition for a fan |
flowRateInletVelocity | This boundary condition provides a velocity boundary condition, derived from the flux (volumetric or mass-based), whose direction is assumed to be normal to the patch freestream This boundary condition provides a free-stream condition. It is a 'mixed' condition derived from the inletOutlet condition, whereby the mode of operation switches between fixed (free stream) value and zero gradient based on the sign of the flux |
mappedVelocityFluxFixedValue | This boundary condition maps the velocity and flux from a neighbour patch to this patch |
outletMappedUniformInlet | This boundary conditon averages the field over the 'outlet' patch specified by name 'outletPatchName' and applies this as the uniform value of the field over this patch |
pressureDirectedInletOutletVelocity | This velocity inlet/outlet boundary condition is applied to pressure boundaries where the pressure is specified. A zero-gradient condition is applied for outflow (as defined by the flux). For inflow, the velocity is obtained from the flux with the specified inlet direction |
pressureDirectedInletVelocity | This velocity inlet boundary condition is applied to patches where the pressure is specified. The inflow velocity is obtained from the flux with the specified inlet direction |
pressureInletOutletVelocity | This velocity inlet/outlet boundary condition is applied to pressure boundaries where the pressure is specified. A zero-gradient condition is applied for outflow (as defined by the flux). For inflow, the velocity is obtained from the patch-face normal component of the internal-cell value |
pressureInletVelocity | This velocity inlet boundary condition is applied to patches where the pressure is specified. The inflow velocity is obtained from the flux with a direction normal to the patch faces |
pressureNormalInletOutletVelocity | This velocity inlet/outlet boundary condition is applied to patches where the pressure is specified. A zero-gradient condition is applied for outflow (as defined by the flux); for inflow, the velocity is obtained from the flux with a direction normal to the patch faces |
pressurePIDControlInletVelocity | This boundary condition tries to generate an inlet velocity that maintains a specified pressure drop between two face zones downstream. The zones should fully span a duct through which all the inlet flow passes |
rotatingPressureInletOutletVelocity | This velocity inlet/outlet boundary condition is applied to patches in a rotating frame where the pressure is specified. A zero-gradient is applied for outflow (as defined by the flux). For inflow, the velocity is obtained from the flux with a direction normal to the patch faces |
rotatingTotalPressure | This boundary condition provides a total pressure condition for patches in a rotating frame |
supersonicFreestream | This boundary condition provides a supersonic free-stream condition |
surfaceNormalFixedValue | This boundary condition provides a surface-normal vector boundary condition by its magnitude |
swirlFlowRateInletVelocity | This boundary condition provides a volumetric- OR mass-flow normal vector boundary condition by its magnitude as an integral over its area with a swirl component determined by the angular speed, given in revolutions per minute (RPM) |
syringePressure | This boundary condition provides a pressure condition, obtained from a zero-D model of the cylinder of a syringe |
timeVaryingMappedFixedValue | This boundary conditions interpolates the values from a set of supplied points in space and time. Supplied data should be specified in constant/boundaryData/patchname/0 and constant/boundaryData/patchname/points where points contains X, Y, Z coordiantes. The default mode of operation (mapMethod with planarInterpolation) is to project the points onto a plane (constructed from the first 3 points), construct a 2D triangulation, find for the face centres the triangle it is in and the weights to the 3 vertices. |
totalPressure | This boundary condition provides a total pressure condition |
totalTemperature | This boundary condition provides a total temperature condition |
turbulentInlet | This boundary condition generates a fluctuating inlet condition by adding a random component to a reference (mean) field |
turbulentIntensityKineticEnergyInlet | This boundary condition provides a turbulent kinetic energy condition, based on user-supplied turbulence intensity, defined as a fraction of the mean velocity. |
uniformTotalPressure | This boundary condition provides a time-varying form of the uniform total pressure boundary condition |
variableHeightFlowRateInletVelocity | This boundary condition provides a velocity boundary condition for multphase flow based on a user-specified volumetric flow rate |
variableHeightFlowRate | This boundary condition provides a phase fraction conditionbased on the local flow conditions, whereby the values are constrained to lay between user-specified upper and lower bounds |
waveSurfacePressure | This is a pressure boundary condition, whose value is calculated as the hydrostatic pressure based on a given displacement |
fixedFluxPressure | This boundary condition sets the pressure gradient to the provided value such that the flux on the boundary is that specified by the velocity boundary condition |
fixedPressureCompressibleDensity | This boundary condition calculates a (liquid) compressible density as a function of pressure and fluid properties |
freestream | This boundary condition provides a free-stream condition. It is a ’mixed’ condition derived from the i¸nletOutlet condition, whereby the mode of operation switches between fixed (free stream) value and zero gradient based on the sign of the flux |
freestreamPressure | This boundary condition provides a free-stream condition for pressure. It is a zero-gradient condition that constrains the flux across the patch based on the free-stream velocity |
mappedFlowRate | Describes a volumetric/mass flow normal vector boundary condition by its magnitude as an integral over its area |
outletInlet | This boundary condition provides a generic inflow condition, with specified outflow for the case of reverse flow |
pressureInletOutletParSlipVelocity | This velocity inlet/outlet boundary condition for pressure boundary where the pressure is specified. A zero-gradient is applied for outflow (as defined by the flux). For inflow, the velocity is obtained from the flux with the specified inlet direction |
pressureInletUniformVelocity | This velocity inlet boundary condition is applied to patches where the pressure is specified. The uniform inflow velocity is obtained by averaging the flux over the patch and then applying it in the direction normal to the patch faces |
advective / waveTransmissive | This boundary condition provides an advective / wave transmissive outflow condition, based on solving DDt(psi, U) = 0 at the boundary. Note these two B.C. are applicable on patches defined as outlet only. |
internalField #codeStream { //Use codeStream to set the value of the initial conditions { codeInclude #{ #include "fvCFD.H" #}; codeOptions #{ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude #}; codeLibs #{ -lmeshTools \ -lfiniteVolume #}; code #{ //Access internal mesh information const IOdictionary& d = static_cast<const IOdictionary&>(dict); const fvMesh& mesh = refCast<const fvMesh>(d.db()); //Initialize scalar field to zero scalarField alpha(mesh.nCells(), 0.); scalar he = 0.5; scalar ke = 0.5; scalar ae = 0.3; scalar be = 0.15; forAll(alpha, i) { //Access cell centers coordinates const scalar x = mesh.C()[i][0]; const scalar y = mesh.C()[i][1]; const scalar z = mesh.C()[i][2]; //Assign value to alpha according to conditional structure if ( pow(y-ke,2) <= ((1 - pow(x-he,2)/pow(ae,2) )*pow(be,2)) ) { alpha[i] = 1.0; } } writeEntry(os, "", alpha); //Write output to input dictionary #}; };
Radiation Modeling
OpenFOAM provides three radiation models: P1 model, fvDOM (finite volume discrete ordinates model) and viewFactor model.If a x L >> 1 then P1 model. If a × L < 1, use of fvDOM is recommended. Since fvDOM also captures the large optical length scales, it is the most accurate model though it is computationally more intensive since it solves the transport equation for each direction. fvDOM can handle non gray surfaces (dependence of radiation on the the solid angle is included). viewFactor method is recommended if non-participating mediums are present (such as spacecraft and solar radiation).
The tutorial heatTransfer\ chtMultiRegionFoam\ externalSolarLoadcase has 3 dictionaries dealing with radiation settings: viewFactorsDict, boundaryRadiationProperties and radiationProperties.
Other tutorial cases dealing with radiation are: heatTransfer\ buoyantSimpleFoam\ hotRadiationRoom, heatTransfer\ buoyantSimpleFoam\ hotRadiationRoomFvDOM and heatTransfer\ chtMultiRegionFoam\ multiRegionHeater.Universal fvSolution dictionary: the solver settings in various tutorial cases were compared and compiled to cover more than 90% of cases. The GAMG solver, short for Geometric Agglomerated Algebraic Multigrid solver, is a linear solver used for the pressure. SnGradSchemes is a user defined variable in OpenFOAM that allows the user to chose what surface normal gradient scheme to use. 'Limited' option means that a limited nonorthogonal correction is to be used. Gaussian integration is a second order discretization scheme and is defined in the fvSchemes dict together with the choice of interpolation scheme. In the Gaussian integration values from cell centers need to be interpolated to face centers.
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // /* The solution is considered converged when the residual has reached the tolerance 1.0e-6 or if it has been reduced by relTol 0.05 at each time step. relTol 0; disables relTol - e.g. in case of PISO as this algorithm solves each equation only once per time step. PCG - Preconditioned Conjugate Gradient solver for symmetric matrices GAMG - generalised geometric-algebraic multi-grid solver PBiCG - Preconditioned biConjugate Gradient solver for asymmetric matrices smoothSolver - solver using a smoother for symmetric / asymmetric matrices PCG - Preconditioned Conjugate Gradient solver for symmetric matrices DILU - Diagonal Incomplete Lower Upper decomposition DIC - Diagonal Incomplete Cholesky - decomposition. */ solvers { rho { //"rho.*" = (rho|rhoFinal) i.e. rho and rhoFinal solver PCG; preconditioner DIC; tolerance 1e-05; relTol 0.1; } rhoFinal { $rho; //copies setting from 'rho' tolerance 1e-05; relTol 0; } "(U|k|epsilon|f|v2|Urel)" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-06; relTol 0.1; //nSweeps 1.0; // used in simpleFoam/airfoil2D } //lagrangian\DPMFoam\Goldschmidt\system "(U.air)" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-05; relTol 0.1; } "(U.air)Final" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-05; relTol 0.1; } "(p|p_rgh)" { solver GAMG; tolerance 1e-7; relTol 0.1; smoother GaussSeidel; //DICGaussSeidel; nCellsInCoarsestLevel 20; nPreSweeps 0; nPostSweeps 2; cacheAgglomeration true; agglomerator faceAreaPair; //mergeLevels 1; } "(p|p_rgh)Final" { $p; tolerance 1e-06; relTol 0; maxIter 20; } "(U|k|epsilon|f|v2)Final" { $U; tolerance 1e-06; relTol 0; } nuTilda { solver smoothSolver; //PBiCGStab; smoother symGaussSeidel; //DILU; tolerance 0; relTol 0.1; } nuTildaFinal { $U; tolerance 1e-6; relTol 0; } "(omega|R)" { solver smoothSolver; smoother GaussSeidel; tolerance 1e-05; relTol 0; } "(omega|R)Final" { $U; tolerance 1e-06; relTol 0; } T { //scalarTransportFoam, laplacianFoam solver PBiCGStab; //PCG preconditioner DILU; //DIC tolerance 1e-06; relTol 0; } Phi { //potentialFoam solver GAMG; smoother DIC; tolerance 1e-06; relTol 0.01; } // -----------------------------Multiphase----------------------------------- //"alpha.*" { nAlphaSubCycles 3; } //multiphase\ interCondensatingEvaporatingFoam\ condensatingVessel /* ".*(rho|rhoFinal)" { solver diagonal; } */ "alpha.water.*" { //Update for alpha.air, alpha.fuel - each phase nAlphaCorr 2; nAlphaSubCycles 1; cAlpha 1; MULESCorr yes; nLimiterIter 3; solver smoothSolver; smoother symGaussSeidel; tolerance 1e-8; relTol 0; } "pcorr.*" { solver PCG; //preconditioner DIC; preconditioner { preconditioner GAMG; tolerance 1e-05; relTol 0; smoother DICGaussSeidel; } tolerance 1e-5; relTol 0; } "e.*" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-7; relTol 0; minIter 1; } "(Theta).*" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-7; relTol 0; minIter 1; } "(h).*" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-6; relTol 0; minIter 1; } "(kinematicCloud:theta)" { solver GAMG; tolerance 1e-06; relTol 0.01; smoother GaussSeidel; } "(kinematicCloud:theta)Final" { solver GAMG; tolerance 1e-06; relTol 0; smoother GaussSeidel; } // -----------------------------Combustion----------------------------------- "(Yi|CO2|O2|N2|CH4|H2|H2O|CO)" { solver PBiCGStab; preconditioner DILU; tolerance 1e-6; relTol 0; } "(Yi|CO2|O2|N2|CH4|H2|H2O|CO)Final" { solver PBiCGStab; preconditioner DILU; tolerance 1e-6; relTol 0; } h { $Yi; relTol 0.1; } hFinal { $Yi; } "(b|Su|Xi|ft|ha|hau)" { solver PBiCGStab; preconditioner DILU; tolerance 1e-05; relTol 0.1; } "(b|Su|Xi|ft|ha|hau)Final" { solver PBiCGStab; preconditioner DILU; tolerance 1e-05; relTol 0; } Ii { solver GAMG; tolerance 1e-4; relTol 0; smoother symGaussSeidel; nPostSweeps 1; cacheAgglomeration true; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; maxIter 5; nPreSweeps 0; nPostSweeps 1; } G { solver PCG; preconditioner DIC; tolerance 1e-04; relTol 0; } "(G)Final" { solver PCG; preconditioner DIC; tolerance 1e-04; relTol 0; } } // -----------------------------PV-Couplings------------------------------------ SIMPLE { nUCorrectors 2; nNonOrthogonalCorrectors 0; residualControl { p 1e-5; U 1e-5; nuTilda 1e-5; "(k|epsilon|omega|f|v2)" 1e-3; // found in simpleFoam/pitzDaily } } PISO { //p equation solved n = 2 time, thus pressure-velocity coupling is done twice nCorrectors 2; nNonOrthogonalCorrectors 0; //The pressure is set to pRefValue 0 in cell number pRefCell 0. //This is over-ridden if constant pressure boundary condition is used pRefCell 0; pRefValue 0; } PIMPLE { transonic no; nCorrectors 2; //p equation solved n = 2 time, thus pressure-velocity coupling is done twice nOuterCorrectors 5; nNonOrthogonalCorrectors 0; consistent no; momentumPredictor yes; hydrostaticInitialization yes; nHydrostaticCorrectors 5; correctPhi yes; moveMeshOuterCorrectors yes; //Compressible version pMaxFactor 1.5; pMinFactor 0.9; turbOnFinalIterOnly no; residualControl { "(U|k|epsilon)" { relTol 0; tolerance 0.0001; } } } // -----------------------------URFs------------------------------------------- relaxationFactors { // 0.9 is more stable but 0.95 more convergent fields { nuTilda 0.75; U 0.75; "p.*" 0.40; "rho.*" 0.90; T 0.90 } equations { k 0.50; epsilon 0.50; //"(U|k).*" 0.75; //Multiphase flows //".*" 0.90; //Same settings for all equations } equations { ".*Final" 0.95; "(CH4|O2|H2O|CO|CO2|h).*" 0.90; } } // ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 1.5 | | \\ / A nd | Web: http://www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default steadyState; //Euler, backward, CrankNicolson 0.5; } gradSchemes { // 'gradSchemes' is a <dictionary>, 'default' is a 'class' default Gauss linear; //'Gauss' is 'object' of class 'default' //default leastSquares; // electrostaticFoam, rhoSimpleFoam grad(p) Gauss linear; grad(U) Gauss linear; // cellLimited Gauss linear 1 } divSchemes { default none; //Gauss upwind div(U) Gauss linear; //Gauss vanLeerV, Gauss vanLeer div(phi,U) Gauss upwind; //Gauss LUST grad(U) div(phi,Urel) Gauss limitedLinear 1; div(phi,k) Gauss upwind; div(phi,epsilon) Gauss upwind; //Gauss linearUpwind limitedGrad div(phi,omega) bounded Gauss limitedLinear 0.2; div(phi,T) Gauss linearUpwind grad(T); div(phi,R) Gauss upwind; div(phi,Ekp) Gauss limitedLinear 1; div(phi,T) Gauss limitedLinear 1; div(R) Gauss linear; div(phi,e) bounded Gauss upwind; //Gauss limitedLinear 1 div(phi,nuTilda) Gauss upwind; //Gauss limitedLinear 1 div(phid,p) Gauss limitedLinear 1; div(phi,K) Gauss limitedLinear 1; div(phiv,p) Gauss limitedLinear 1; div((nuEff*dev(grad(U).T()))) Gauss linear; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; //Multiphase flow: for each phase div(rhoPhi,U) Gauss linear; div(phi,alpha) Gauss vanLeer; div(phi,s) Gauss vanLeer; div(phirb,s) Gauss linear; div(phirb,alpha) Gauss linear; div(phi,k) Gauss limitedLinear 1; div(phi,B) Gauss limitedLinear 1; div(B) Gauss linear; div(alphaPhi.air,U.air) Gauss linearUpwindV unlimited; div(((alpha.air*nuEff.air)*dev2(T(grad(U.air))))) Gauss linear; "div\(phi,alpha.*\)" Gauss vanLeer; "div\(phir,alpha.*,alpha.*\)" Gauss vanLeer; "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1; div(Rc) Gauss linear; "div\(phi.*,U.*\)" Gauss limitedLinearV 1; "div\(phir,alpha.*\)" Gauss vanLeer; "div\(alphaRhoPhi.*,U.*\)" Gauss limitedLinearV 1; "div\(alphaRhoPhi.*,(h|e).*\)" Gauss limitedLinear 1; "div\(alphaRhoPhi.*,K.*\)" Gauss limitedLinear 1; "div\(alphaPhi.*,p\)" Gauss limitedLinear 1; //reactingTwoPhaseEulerFoam div(alphaRhoPhi.particles,Theta.particles) Gauss limitedLinear 1; "div\(alphaRhoPhi.*,(k|epsilon).*\)" Gauss limitedLinear 1; div((((alpha.air*thermo:rho.air)*nuEff.air)*dev2(T(grad(U.air))))) Gauss linear; "div\(\(\(\(alpha.*\*thermo:rho.*\)\*nuEff.*\)\*dev2\(T\(grad\(U.*\)\)\)\)\)" Gauss linear; div((((thermo:rho.particles*nut.particles)*dev2(T(grad(U.particles)))) + (((thermo:rho.particles*lambda.particles)*div(phi.particles))*I))) Gauss linear; //Combustion models div(Ji,Ii_h) Gauss upwind; //div(phi,Yi_h) Gauss upwind; div(phi,Yi_h) Gauss multivariateSelection { O2 limitedLinear01 1; CH4 limitedLinear01 1; N2 limitedLinear01 1; H2O limitedLinear01 1; CO2 limitedLinear01 1; h limitedLinear 1; }; div((Su*n)) Gauss linear; div((U+((Su*Xi)*n))) Gauss linear; div(phi,ft_b_ha_hau) Gauss multivariateSelection { fu limitedLinear01 1; ft limitedLinear01 1; b limitedLinear01 1; ha limitedLinear 1; hau limitedLinear 1; }; div(phiXi,Xi) Gauss limitedLinear 1; div(phiSt,b) Gauss limitedLinear01 1; } laplacianSchemes { default Gauss linear orthogonal; //Gauss linear uncorrected laplacian(nuEff,U) Gauss linear corrected; //Gauss linear limited corrected 0.5 laplacian((1|A(U)),p) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected; laplacian(DepsilonEff,epsilon) Gauss linear corrected; laplacian(DREff,R) Gauss linear corrected; laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; laplacian(DT,T) Gauss linear corrected; //laplacianFoam laplacian(diffusivity,cellMotionU) Gauss linear uncorrected; } interpolationSchemes { default linear; interpolate(U) linear; } snGradSchemes { default orthogonal; //corrected, limited corrected 0.333 } fluxRequired { default no; p; } wallDist { method meshWave; } // ************************************************************************* //
The content on CFDyna.com is being constantly refined and improvised with on-the-job experience, testing, and training. Examples might be simplified to improve insight into the physics and basic understanding. Linked pages, articles, references, and examples are constantly reviewed to reduce errors, but we cannot warrant full correctness of all content.
Template by OS Templates