## Semicylinder with generation

### 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 $T_0$ and suddenly the part is quickly immersed in a forced water flow at $T_{water}$. An electrical current circulates through the region 4 so that a heat power generation per length meter, $g_4$ (W),  must be considered (Joule effect) in this region.

Determine:

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

Data:

$|~a=b=0.010\:m;~|$ $|~\lambda=0,5\:W/(m \cdot K)~~(thermal ~ conductivity);~|$ $|~c_p=850\:J/(kg \cdot K)~~(specific ~ heat);~|$ $|~\rho=2700 \: kg/m^3~~(density);~|$ $|~g_4=300 \:W (per ~ length ~ meter);~|$ $|~h_{water}=50 \:W/(m^2 \cdot K) \quad (water ~ convection ~ coeficient) ~|$ $|~T_{water}=20 \,^\circ{C} ~|$ $|~ T_0=200\,^\circ{C} ~|$

# 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:

$\dfrac{\lambda}{\rho ~ c_p} \Delta T + \dfrac{g}{\rho~c_p} = \dfrac{\partial T}{\partial t}$

where $\frac{\lambda}{\rho ~ c_p}$, also known as $\alpha$, is the thermal diffusivity, $\Delta T$ is the laplacian of the temperature and $g$ 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_SOLVERS basic 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 $\alpha$, the thermal diffusivity. The expression shown above translated into a mathematical notation would be something like

$\dfrac{\partial T}{\partial t} - \alpha \cdot \Delta T=0$

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! $\dfrac{\partial T}{\partial t} + \nabla \cdot (\phi T)- \alpha \cdot \Delta T=fvOptions(T)$ 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 ($\phi$) 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  |                                                 |
\*---------------------------------------------------------------------------*/
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, $\rho$, the thermal conductivity, $\lambda$, and the heat capacity, $c_p$. However, recall that all of these properties are related with the concept of thermal diffusivity, $\alpha$, 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. 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

$\dfrac{\partial T}{\partial t} - \alpha \Delta T = \dfrac{g}{\rho~c_p}$

which is written in the code as

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

Thus, we get

$\dfrac{g}{\rho~c_p} = fvOptions(T) = S_u + S_p \cdot T$

Which gives us the values of  $S_u = \frac{g}{\rho~c_p}$  and  $S_p = 0$. 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.

• selectionMode: This keyword is used in order to let OpenFOAM know what kind of entity (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).
• volumeMode: It allows you to define the quantity of the source in absoulte terms (watts in the example) or volumetric/specific terms (W/m³).
• injectionRateSuSp: Subdictionary where the variables that will be affected by the source term (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, $g_4 = 300 \: W$ 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 volumeMode is set to 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
{
}


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.

# 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  |                                                 |
\*---------------------------------------------------------------------------*/
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