### Navigation map

- Formulation
- 1st step: choosing the solver
- 2nd step: creating the mesh
- 3rd step: defining properties and boundary conditions
- 4th step: defining the heat source
- 5th step: getting the solver ready
- 6th step: running the solver
- Question 1: solution
- Question 2: solution

### Formulation

The temperature distribution of the central cross section of a long semicylindrical component is to be studied. This part is fixed to a metal block and on an insulating base. At the start (t=0) the temperature of the whole component is and suddenly the part is quickly immersed in a forced water flow at . An electrical current circulates through the region 4 so that a heat power generation per length meter, (W), must be considered (Joule effect) in this region.

Determine:

- The temperature distribution when the component reaches steady state conditions.
- The power balance in the system.

**Data**:

# 1st step: choosing the solver

The first thing one has to do in order to solve a a case with OpenFOAM is to decide what solver is to be used. To do that it is necessary to know what kind of phenomena are involved in the problem, because OpenFOAM has a vast variety of solvers to choose, each one for a different kind of problem. The case here presented consists of a typical problem of **thermal conduction in a solid in 2D**. The governing equation of this phenomenon is, as it is known, the following one:

where , also known as , is the thermal diffusivity, is the laplacian of the temperature and is the heat power generated per unit volume. The term on the RHS of the equation is the time derivative of the temperature. Now a solver that solves this equation has to be found. One can find all kinds of solvers available with OpenFOAM by typing on the terminal the following command:

ls $FOAM_SOLVERSbasic compressible DNS financial incompressible multiphase combustion discreteMethods electromagnetics heatTransfer lagrangian stressAnalysis

The list shown above shows all kinds of problem that OpenFOAM is capable to solve. Each directory contains a set of solvers related to its field. Regarding the problem we are trying to solve one could think that the appropriate solver may be found within the *heatTransfer* directory, however, as this is a basic problem that only involves heat conduction in a solid, **the solver we need is located within the basic directory**, which only contains three solvers. If the contents within a solver directory are listed one can see a set of files ended with “.C” or “.H” and a directory called “Make”. The most interesting file and the one that contains the main code of the solver is the one called “nameOfTheSolver.C”. This file also contains a brief explanation of what kind of phenomena this solver is capable to solve. For instance, if one wants to know about the sort of problems that

*laplacianFoam*can solve, first its main file must be opened with a text editor, such as, emacs, gedit, kate, or another one by typing in the terminal

gedit $FOAM_SOLVERS/basic/laplacianFoam/laplacianFoam.C &

Once opened, a brief description of the physical phenomena that the solver is able to solve can be seen in the lines 24-28.

Application laplacianFoam Description Solves a simple Laplace equation, e.g. for thermal diffusion in a solid.

It seems that we are on the right track! However, to get a deeper understanding of what the solver does and what it doesn’t we have to scroll down a little bit more, until the lines 57-60, here we can read

solve ( fvm::ddt(T) - fvm::laplacian(DT, T) );

where DT is , the thermal diffusivity. The expression shown above translated into a mathematical notation would be something like

which is a very similar equation to the one that we need to solve, however, this one does not take into account the heat source, so **this solver is not suitable to solve our problem**. *But don’t fall apart! It was not a bad start at all! If you want to become an OpenFOAM master you will have to get used to miss some shots before you hit the target!*

After the first failed attempt, let’s see what can do *scalarTransportFoam* for us! First of all one has to open the main file of this solver like it has been done before

gedit $FOAM_SOLVERS/basic/scalarTransportFoam/scalarTransportFoam.C &

and take a look at the lines 24-28 to see the brief description of the solver purpose.

Application scalarTransportFoam Description Solves a transport equation for a passive scalar

This time it has a more general description that might cause some confusion and might make you reject this solver. However, before making a hasty decision let’s look into some lines below, namely the lines 60-67.

solve ( fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(DT, T) == fvOptions(T) );

Now let’s translate them into mathematical notation!

We can consider this equation as an extended version of the previous one, or let’s call it a more general equation. It has a couple of terms more than *laplacianFoam* has: one of them will help us a lot; but the other one, at a first glance, may be a little annoying. The divergence term (the annoying one) takes into account the flux of material () that transport our passive scalar (temperature) through the domain. In a solid domain, however, there are no flux of material, so this term must be neglected. We will do it later by setting the velocity field to 0 (or (0, 0, 0), velocity is a vector field…). So, eventually, this term is not that annoying at all! The mysterious term on the RHS of the equation is, as we’ll see some steps below, the helpful one. This term will allow us to define the heat source that we need so as to solve this case. So, finally **we have found the proper solver to solve our problem!**

Congratulations! 1st step completed…

* * **Tip**: The most common way to proceed from this point on in order to set up your case is to copy an original case from the $FOAM_TUTORIALS directory to your working directory, i.e., $FOAM_RUN. In our case

cp -r $FOAM_TUTORIALS/basic/scalarTransportFoam/pitzDaily $FOAM_RUN

*Afterwards this case is renamed *

mv $FOAM_RUN/pitzDaily $FOAM_RUN/semicylinder_withGeneration

*and modified so as to suit our needs. This is a very common practice known as CPA.*

* * **Warning**: From this point on it is assumed that you are located within the case directory ($FOAM_RUN/semicylinder_withGeneration).

# 2nd step: creating the mesh

Time to show your creativity! Now, the mesh of the problem has to be built and, since the geometry of it is quite simple, it can be done by using *blockMesh*, a tool provided with OpenFOAM. This tool takes the values stored in a file we have to create or modify in order to generate the geometry and the mesh of the problem. This file is called *blockMeshDict* and it is placed within the *constant* directory (in current versions of OpenFOAM this file is placed within the *system* directory instead, so if you are using a newer version maybe you have to look for it in this directory, not in the *constant* directory).

First we have to open the *blockMeshDict* file with the following command

gedit constant/polyMesh/blockMeshDict

*blockMesh* generates the mesh over a block-based geometry, this means we have to define a set of blocks in the *blockMeshDict* file and, to do this, first a set of points is required. To understand how *blockMeshDict* is structured a deep read to the User’s Guide (UG) is required, specifically the chapter 5.3. After having enjoyed the reading, I’m sure you can understand that the geometry shown in the animation below corresponds with the *blockMeshDict* file also attached below.

An interesting point to mention is the block definition used.

blocks ( hex (0 3 5 0 11 14 16 11) base (6 8 1) simpleGrading (1 1 1) // block 0 ... );

First of all, some vertices appear doubled in the list of vertices that make up the block (vertices 0 and 11 in the example above), this is because this block is not a hexaedron since **it has an edge collapsed to form a wedge-shaped block**. This procedure is also comented in the UG but I think it deserves a special mention because this kind of blocks may cause some trouble to the beginners.

In second place, it can be noticed that between the two first parentheses there is a word that has not been comented in the UG, the presence of this word is a little trick that can speed up a little bit the set up of some simple cases like this one. This word is “the block name”! *And why would anyone want to name a block? You are wondering…* Well, recall that we have to define a heat source in a region of our domain, so here’s the answer! We need to give it a name so that later we can tell OpenFOAM which region is the one that generates heat. **By naming a block we are creating a cellZone with this name** (this concept will be discussed further later).

Once the *blockMeshDict* file is correctly written it has to be executed on the terminal like that

blockMesh /*---------------------------------------------------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.3.x | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ Build : 2.3.x-13910e4c7ad3 Exec : blockMesh Date : Jan 10 2016 Time : 14:13:11 Host : "linux-cfd" PID : 6924 Case : /home/cfd/OpenFOAM/cfd-2.3.x/run/semicylinder_withGeneration nProcs : 1 sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE). fileModificationChecking : Monitoring run-time modified files using timeStampMaster allowSystemOperations : Allowing user-supplied system call operations // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Create time Creating block mesh from ... [lots of lines] ---------------- Patches ---------------- patch 0 (start: 586 size: 16) name: metalBlock patch 1 (start: 602 size: 12) name: insulation patch 2 (start: 614 size: 28) name: waterConvection patch 3 (start: 642 size: 630) name: frontAndBack End

A message like the one above means that the execution has ended with success, so **the mesh has been generated successfully!**

* * * Tip: Notice that I usually draw a sketch of the geometry at the top of the *blockMeshDict

*file in order to be able to check quickly any error that may appear after the*blockMesh

*execution (and believe me, more than one will show up when you build your own geometry from scratch!). Obviously it is only possible when the geometry is simple enough. However, I strongly recomend to take some annotations within the*blockMeshDict

*file using commenting tags (single-line comments starting with*

`//`

; multi-line comments between `/*`

and `*/`

). I also recomend to always draw a sketch of the geometry on paper!* * * Warning: Watch out! The links given above to the UG take you to the documentation of the last version of the software. If you are using an older version some features discussed in those links may not be availabe. In this case, to avoid confusions, better find the official documentation of your own version with the following command*:

evince $FOAM_INST_DIR/OpenFOAM-x.x.x/doc/Guides-a4/UserGuide.pdf &

# 3rd step: defining properties and boundary conditions

Now it’s time to define the material properties, the initial conditions and the BC’s for the variables to be solved in accordance with the problem formulation. The **variables to solve** are the **temperature** and the **velocity** (this last one is required to compute the flux, which is 0 in a solid), on the other hand, the **material properties** needed to solve the case are the **density**, , the **thermal conductivity**, , and the **heat capacity**, . However, recall that all of these properties are related with the concept of *thermal diffusivity*, , called DT in the solver code. So, this is the only property needed by the solver. The files where all these values are stored are found within the *constant* and the *0* directory and are opened with the following command lines

gedit 0/T &gedit 0/U &gedit constant/transportProperties &

The files with the proper values to solve our case can be found on the links below.

Open T fileOpen U fileOpen transportProperties

Comments about the boundary and initial conditions:

**U:**

Notice that both initial and boundary field values are set to `(0 0 0)`

. This is to prevent from flux existence in the domain. However, it is only like that in the internal field. The **values at the boundaries are set to (0 0 0) because of the no-slip condition!**

dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { metalBlock { type fixedValue; value uniform (0 0 0); } ... }

**T:**

First of all, I have to clarify that, although the** temperature should be expressed in Kelvin**, in these simple cases **it has been specified in Celsius**. It is possible to proceed in this manner because no heat transfer by radiation is taken into account in this problem and because **none of the temperatures are below 0ºC**.

Besides that, special attention is required when it comes to specify the convective boundary condition. This BC has been specified with the aid of one of the most useful packages one must have to extend OpenFOAM capabilites without the need of programming. This package is **swak4foam** and it provides a lot of useful tools to enlarge the benefits of OpenFOAM. **groovyBC**, one of the extended capabilities provided with swak4foam, allows you to create your custom BC’s by using user-defined expressions. This expressions are used in order to define the values of the internal variables required by a BC, these are *valueExpression*, *gradientExpression* and *fractionExpression*. These ones substitute the common variables *refValue*, *refGrad* and *f* used by the regular BC’s that come out of the box with OpenFOAM.

dimensions [0 0 0 1 0 0 0]; internalField uniform 200; boundaryField { ... convectionWater { type groovyBC; value uniform 200; valueExpression "Twater"; gradientExpression "0"; fractionExpression "1.0/(1.0 + k/(mag(delta())*hc))"; variables ( "hc=50.0;" "Twater=20.0;" "k=0.5;" ); } ... }

If you want to understand the theory behind the definition of the convection BC take a look at the document CONVECTION BC.

* * **Tip**: T can be specified in Celsius when the case doesn’t include radiation and if no boundary field has a negative temperature. However, this trick has to be used with a lot of care since, although no field is specified with a temperature below 0ºC, some cell can get a negative value for the temperature during the solver run and it will make the solver crash!

* * **Warning**: swak4foam needs to be installed in order to be able to define the convection BC in the manner it has been done in this tutorial!

# 4th step: defining the heat source

Now it’s the time to solve the main issue of this case, at least for a beginner-level user. This is the definition of the heat source! As it has been said before, for this purpose,** fvOptions framework has to be used**. So, the question is how is

*fvOptions*framework used in order to create a heat source? Well, first of all we have to find the proper type of it, since many types of

*fvOptions*have been developed to be able to create different types of source terms for different governing equations. We will use a basic type of source called

**semiImplicitSource**, whose header file (file appended with “.H”) explains briefly how it is used. The file can be accessed with the command

gedit $FOAM_SRC/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H &

And its description states

Description Semi-implicit source, described using an input dictionary. The injection rate coefficients are specified as pairs of Su-Sp coefficients, i.e. \f[ S(x) = S_u + S_p x \f] where \vartable S(x) | net source for field 'x' S_u | explicit source contribution S_p | linearised implicit contribution \endvartable Example of the source specification: \verbatim SemiImplicitSourceCoeffs { volumeMode absolute; // specific injectionRateSuSp { k (30.7 0); epsilon (1.5 0); } } \endverbatim Valid options for the \c volumeMode entry include: - absolute: values are given as \ - specific: values are given as \/m3

The description is quite self-explanatory, we have to give a couple of values for the variable to be solved (T), one of them is a constant value and the other one is a temperature-dependent value. But, how can we find these values in order to create the heat source we need? Well, first of all, we have to keep in mind what governing equation is to be solved, which is the heat conduction equation presented in the 1st step, also shown below

which is written in the code as

solve ( fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(DT, T) == fvOptions(T) );

Thus, we get

Which gives us the values of and . Knowing all the information required to write correctly the *fvOptions* file, one can open it by executing on the terminal the following command (it may happen that this file does not exist in the original case you have CPAed so, in this case, you will have to copy this file from another tutorial case which has it)

gedit system/fvOptions &

and make it look like the *fvOptions* file that can be found on the following link.

A few comments about the use of **semiImplicitSource**:

: This keyword is used in order to let OpenFOAM know what kind of entity (**selectionMode**`cellZone`

,`cellSet`

,…) is to be selected to apply the source. The entity type determines the next keyword that must be defined in the file, that is, the field where the name of the entity is stored (*source*in this example).

: It allows you to define the quantity of the source in**volumeMode**`absoulte`

terms (watts in the example) or volumetric/`specific`

terms (W/m³).

: Subdictionary where the variables that will be affected by the source term (**injectionRateSuSp**`T`

) are stored along with their corresponding Su-Sp values.- Note that 3W of power generated have been considered, this is because, as stated in the problem formulation, per length meter and the length considered is 10 mm.
- Also note that in the attached file a value has been commented out (0.6668…), this value corresponds with the volumetric generation ratio so, if

is set to**volumeMode**`specific`

the results obtained will be the same. To the reader concerns to prove it!

# 5th step: getting the solver ready

This is usually the most annoying part when you are a beginner since the files to check and modify have a very bad look when you are not used to work with OpenFOAM. These files are called *fvSchemes*, *fvSolution* and *controlDict* and are placed within the *system* directory. The first one contains the discretization schemes to be used for each one of the governing equation terms; the second one is to specify the solver to be used for each of the variables to solve and the last one is to adjust the parameters that control the solver run. I strongly recomend to take a look at the UG, sections 4.3 (controlDict), 4.4 (fvSchemes) and 4.5 (fvSolution), in order to get a better understanding of the way these files are implemented. Now, let’s open the files!

gedit system/{fvSolution,fvSchemes,controlDict} &

However, I have a good new! Since it is a very simple case CPAed from another very simple case there are just a few changes to be made. The most important one is to be applied in the *fvSchemes* in the `ddtSchemes`

subdictionary. Since this is a steady state case the time-dependent component must be set accordingly by applying the following value to the `ddtSchemes`

subdictionary

ddtSchemes { default steadyState; }

By applying this modification the case can already be solved. However, it is important to check first the *controlDict* file in order to check the time parameters of the solver, although being a steady state case the time step (`deltaT`

) is not a key point. The next change I propose is only an esthetic change since it will make no difference

endTime 100.0; deltaT 2.0; writeControl timeStep; writeInterval 5;

I like these values more than the original ones used in the *pitzDaily* case only because they are *integers* and I think it looks better this way if these values will not affect the final result (because it is a steady state problem). The whole files can be found on the following links.

Open fvSchemes Open fvSolution Open controlDict

# 6th step: running the solver

Finally, it’s time to have fun. Its time to finish this tough process and reap the rewards of your hard work. *It’s time to solve your first OpenFOAM case!* And certainly this is the easiest one you will ever solve, for sure… To run a solver is as easy as typing in the terminal its name and press enter. In this case

scalarTransportFoam /*---------------------------------------------------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.3.x | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ Build : 2.3.x-13910e4c7ad3 Exec : scalarTransportFoam Date : Jan 15 2016 Time : 15:20:19 Host : "linux-cfd" PID : 6352 Case : /home/cfd/OpenFOAM/cfd-2.3.x/run/semicylinder_withGeneration nProcs : 1 sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE). fileModificationChecking : Monitoring run-time modified files using timeStampMaster allowSystemOperations : Allowing user-supplied system call operations // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Create time Create mesh for time = 0 ... [lots of lines] Time = 98 DILUPBiCG: Solving for T, Initial residual = 3.67827e-08, Final residual = 3.67827e-08, No Iterations 0 Time = 100 DILUPBiCG: Solving for T, Initial residual = 3.67827e-08, Final residual = 3.67827e-08, No Iterations 0 End

The final word *End* indicates that the solver ended the run without problems, which does not necessarily mean that the solver has solved the problem correctly… But don’t worry, if you have taken all the steps above as you should, everything will be alright!

And finally you have solved your first OpenFOAM case successfully! Time to watch your results!

* * **Tip**: If you want to store the output of an application (solver or any other utility) into a file instead of printing it on the screen you have to type this order instead

scalarTransportFoam > nameYouWantForTheFile 2>&1

* * **Tip**: In the case you have followed all the steps in order to solve this case, you must have gotten an output a little more dirty when running the solver, namely, something like what it is shown below must appear at each time step

Time = X --> FOAM Warning : From function gaussConvectionScheme in file finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H at line 123 Reading "/home/cfd/OpenFOAM/cfd-2.3.x/run/semicylinder_withGeneration/system/fvSchemes.divSchemes.div(phi,T)" at line 30 Unbounded 'Gauss' div scheme used in steady-state solver, use 'bounded Gauss' to ensure boundedness. To remove this warning switch off 'warnUnboundedGauss' in "/home/cfd/OpenFOAM/OpenFOAM-2.3.x/etc/controlDict" DILUPBiCG: Solving for T, Initial residual = XXXXXXX, Final residual = XXXXXXX, No Iterations XX

*It is just a warning message that has no effect on the results, at least in this case. If you want to get rid of it, take a look at the fvSchemes file I proposed above, there you will find the answer!*