## Multiregion 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. Furthermore, this time a new material for the generating region is to be tested.

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;~|$ $|~h_{water}=50 \:W/(m^2 \cdot K) \quad (water ~ convection ~ coeficient) ~|$ $|~T_{water}=20 \,^\circ{C} ~|$ $|~ T_0=200\,^\circ{C} ~|$
$\underline{ \forall i \neq 4} :$
$|~\lambda_i=0.5\:W/(m \cdot K)~~(thermal ~ conductivity);~|$ $|~c_{p,i}=850\:J/(kg \cdot K)~~(specific ~ heat);~|$ $|~\rho_i=2700 \: kg/m^3~~(density);~|$
$\underline{i = 4} :$
$|~\lambda_4=10\:W/(m \cdot K)~~(thermal ~ conductivity);~|$ $|~c_{p,4}=850\:J/(kg \cdot K)~~(specific ~ heat);~|$ $|~\rho_4=4200 \: kg/m^3~~(density);~|$ $|~g_4=300 \:W (per ~ length ~ meter);~|$

# 1st step: choosing the solver(s)

As it has been said before, the first thing one has to do in order to solve a case with OpenFOAM is to decide what solver is to be used. This time it is easier to decide because, as it is a heat conduction problem through multiple solids in contact, we only have a couple of solvers to choose from. These solvers are found within the heatTransfer directory and their names are chtMultiRegionFoam and chtMultiRegionSimpleFoam, being the first one for transient problems and the second one for steady-state problems. Let’s recall the governing equation of this phenomenon, which is 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. This equation solves the energy equation for temperature (T), but it can also be solved for enthalpy (h) instead. To do so one only has to remember that $h = c_p \cdot T$. Then, if we multiply the equation above by $\rho \cdot c_p$ we obtain

$\dfrac{\lambda}{c_p} \, \Delta h + g = \rho \, \dfrac{\partial h}{\partial t}$

This process is necessary so as to understand how chtMultiRegion* solvers develop the energy equation in solid regions. Herafter, the parameter $\frac{\lambda}{c_p}$, which must not be confused with the thermal diffusivity, will be called $\alpha'$. In order to study the energy equation that OpenFOAM solves for the case of these two solvers we have to open the corresponding solveSolid.H file located within the solid directory for each of the solvers. The following commands will open both files

  gedit $FOAM_SOLVERS/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H & gedit$FOAM_SOLVERS/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H &

Note that chtMultiRegionSimpleFoam is located within chtMultiRegionFoam directory. Once opened one can see the lines below

• chtMultiRegionFoam/solid/solveSolid.H (lines 9-19)
    tmp hEqn
(
fvm::ddt(betav*rho, h)
- (
thermo.isotropic()
? fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)")
: fvm::laplacian(betav*taniAlpha(), h, "laplacian(alpha,h)")
)
==
fvOptions(rho, h)
);

• chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H (lines 4-13)
    fvScalarMatrix hEqn
(
(
thermo.isotropic()
? -fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)")
: -fvm::laplacian(betav*taniAlpha(), h, "laplacian(alpha,h)")
)
==
fvOptions(rho, h)
);


The first time one sees theses lines may get a shock. In order to be able to understand them, first a little bit of knowledge in C++ is required, at least to understand the construction sentence ? action 1 : action 2. This is a conditional sentence and must be read as if sentence=true then action 1; else action 2. In this case it only checks for the solid isotropy, so one has only to take a look at the corresponding action line. We will focus on the action 1 line, since our case is isotropic.

At this stage, another weird parameter shows up, this is betav. This parameter is a factor to account for porosity in the solid region and takes values between 0 and 1. In order to use it, you have to create a file called betavSolid within the 0/ directory and set it up according to your needs. However, if you don’t need to account for porosity, just do nothing (that’s what we are going to do in this tutorial), and consider betav as 1 in the solver code.

Now let’s see another expression that may cause some confusion, this is the laplacian term in quotation marks within the main laplacian term. Is this the laplacian of the laplacian? Someone may think… No! This is the name given to the whole laplacian term to be used in other parts of the code, that is, something like a nickname! Thereby, in the fvSchemes you will define the discretization schemes for the laplacian term by calling laplacian(alpha,h) instead of calling the whole expression of the laplacian term.

Now that we know what each of the terms mean we can transform the whole code into something like

    tmp hEqn
(
fvm::ddt(rho, h)
- "laplacian(alpha,h)"
==
fvOptions(rho, h)
);


which is more understandable and can be directly compared to the expression obtained before, also shown below,

$\rho \, \dfrac{\partial h}{\partial t} - \alpha' \, \Delta h = g$

where, as it has been said above, $\alpha' = \frac{\lambda}{c_p}$, which is named alpha in the solver code.

### Brief explanation about multiregion solvers

At this stage, a brief explanation on how multiregion solvers work will be given. First of all, multiregion cases are seen by OpenFOAM like a set of multiple simple cases running simultaneoulsy and linked all together by means of some particular BC’s that take the variables from one region boundary to the corresponding neighbour region boundary.

So the solver acts as if it was solving multiple cases (one per region) at a time. But, how does OpenFOAM handle this? Where are the sub-case directories located within the case directory? Easy! The directories corresponding to each region (sub-case) are located within each one of the main directories of the case, that is to say, within the 0 -and the rest of the time directories that will be created after the execution-, constant and system directories! Thus, if one lists the content within the system directory of the multiRegionHeater case the following directories (in bold) and files are shown

  ls $FOAM_TUTORIALS/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system bottomWater decomposeParDict fvSolution leftSolid rightSolid topoSetDict controlDict fvSchemes heater README topAir  where the directories correspond with the multiple regions that comprise the case. The region directories contain the region characteristic files so that each region can have its own parameters -which is, obviously, the main goal of this solver!-. One of the most important files to mention in multiregion cases is the regionProperties file, which is where the multiple regions of the case are defined and classified by its nature, that is to say, the distinction beteween fluid and solid regions is done within this file. Tip: Remember, creating a new case from scratch is not a common practice -I never did it…-, CPA it instead!  cp -r$FOAM_TUTORIALS/heatTransfer/chtMultiRegionFoam/multiRegionHeater $FOAM_RUN mv$FOAM_RUN/multiRegionHeater $FOAM_RUN/semicylinder_withGeneration_multiregion Tip: I also recomend to first rename the regions from the original case you copy-pasted with your current region names according to their nature -solids or fluids, check first the “constant/regionProperties” file in order to know quickly each region’s nature- in the “constant” and “system” directories. You can do so by applying any of the following orders.  cp -r constant/ constant/ & mv constant/ constant/ & Finally, remove the remaining original regions. I usually do it because if you run the preprocessing utilities in order to create all the geometry and meshes for all the regions, if they don’t exist previously OpenFOAM will create the corresponding directories for all the non-existent regions with no information contained in it, whereas if you previously create them from an original region they already have all the files needed and you will only have to adapt the information contained in them with your own information. Warning: From this point on it is assumed that you are located within the case directory ($FOAM_RUN/semicylinder_withGeneration_multiregion).

# 2nd step: creating the mesh

Time to create the mesh of the case geometry and, this time, it is not as straightforward -at least the first time you have to prepare a multiregion case, after you have completed your first attempt everything becomes easier and faster- as it was in the previous case. However, although the process is a little harder and, perhaps, cryptic in the eyes of a newbie, if you understand the logic behind each step you will finally understand the reason why so many utilites -not that many, actually, but a newbie may get scared the first time- are to be used to set up properly the mesh of the case. Well, let’s define the path to follow in order to set up a (simple) multiregion mesh.

1. First of all, the global mesh, that comprises all the regions, has to be defined.
2. Afterwards, we need to define the zones (cellZones) where the different regions will be situated.
3. Next, the mesh is split according to the cellZones defined in order to create the multiple meshes corresponding to each of the regions.
4. Finally, the BC’s for the interregion faces must be defined in order to “connect” all the regions (some other patches may need to be redefined so as to get the physical behaviour desired).
Utilities to use in both meshing modes
General Mode “Shortcut” Mode
blockMesh blockMesh (including cellZone names)
topoSet
splitMeshRegions -cellZones -overwrite
changeDictionary -region

In general, one utility has to be executed in order to take each of the steps above. However, there is a trick or shortcut that allows to avoid running one of the utilities. The utility that can be skipped -mainly in simple cases- is topoSet and is used in order to define the cellZones required which, as it has been shown in the previous case, it can also be done within the blockMeshDict dictionary.

The possibility to skip the execution of a utility can be seen as a simplification or acceleration of the process, but it depends on the geometry of the case you are solving. As this is a very simple case, it can simplify -or, at least, speed up- the set up of the geometry, for this reason, both ways to prepare the case geometry will be explained.

Furthermore, some variables will be used to define the geometrical parameters in the blockMeshDict dictionary so as to automate the mesh generation process. The use of these variables makes really easy and fast to change the dimensions of some parts of the geometry just by changing the value of the corresponding variable. In the following links you can find the text files used for the geometry creation following both modes mentioned so far.

General Mode “Shortcut” Mode

Notice that the blockMeshDict file used for the “Shortcut” Mode is exactly the same used in the previous case, the only difference is the use of variables that make the file larger since a subdictionary with the variables is created before the mesh is defined. Also notice that, although the General Mode uses an extra stage -with the corresponding extra file to be defined-, this extra stage can simplify the geometry creation, because in order to skip this step more blocks -with the corresponding boundary faces, edges…- need to be defined in the blockMeshDict file, which can complicate the geometry generation.

The corresponding topoSetDict dictionary is used in order to manipulate the base mesh to get the multiregion mesh. To do so, first a set of cellSets is created and, afterwards, they are converted into the cellZones that define the different regions. As a picture is worth a thousand words the following animation shows what topoSet utility does when the topoSetDict dictionary above is used. Note that the colors used in the animation are the same used in the first set of blocks defined in topoSetDict. Guess why…

Once the cellZones have been created it’s time to split the overall mesh so as to obtain the regions corresponding to each cellZone. The following command executes the split:

  splitMeshRegions -cellZones -overwrite

Warning: It is necessary to give some dummy values to the variables within the 0 directory before splitting the mesh to avoid some errors. Something like that would be enough

...
internalField   uniform 300;

boundaryField
{
".*"
{
}
}

In case that zeroGradient still gives you troubles, you can try with fixedValue or calculated intead. But don’t worry the values you give them now it’s kind of meaningless, just to avoid problems. The real values will be given in the following step.

# 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. However, first of all, it’s necessary to specify which regions are to be defined as fluids and which ones as solids. This is done by means of the regionProperties dictionary located within the constant diectory. In this case, this dictionary has to look like it is shown below

regions
(
fluid       ()
solid       (base source)
);

It is noticed that no fluid region is to be computed since this problem only contains solid regions (base and source). Next, the material properties of each region need to be defined, these are the density, $\rho$, the thermal conductivity, $\lambda$, and the heat capacity, $c_p$. The files where all these values are stored are found within the constant/ directory and are opened with the following command line

   gedit constant/*/thermophysicalProperties &

Below you can find an example of thermophysicalProperties file, in this case only for the region base (for the source region only kappa and rho need to be redefined).

The thermoType subdictionary is worth mentioning since its use can be a little tricky because of the large list of possible choices you can take for each entry. I strongly recomend to take a look at the section 7.1 of the UG in order to know all the possible thermophysical models available. Another trick commonly used by many users is the well-known “Banana Trick”, which consists in giving a dummy value to any of the entries (e.g., banana) so as to get an error message with all the possible combinations accepted by OpenFOAM.

Finally, the variables to solve are the temperature and the pressure (this last one is required by the solver although it’s kind of meaningless in the case of a solid). In this case, we are not going to define tha initial conditions and BC’s in the 0 directory like we did in the previous case since we also have to define the coupling patches between regions. We could do that by hand within the 0/ directories but OpenFOAM provides a utility in order to automate this process and update all the BC’s just after the mesh is split. In order do so, first the corresponding dictionaries must be properly set up for each region within the system/ directory. These files are called changeDictionaryDict and below you can find the ones used for the current case.

As you can see this file looks very similar to the files where the varaibles are usually stored in single region cases, that is, the files within the 0 directory. The only difference is that in this cases all the information is stored within a subdictionary called dictionaryReplacement. Besides that, if you want to replace information for more than one variable you have to do it in the same file in contrast to what it is done normally (1 file, 1 variable) following the pattern below

dictionaryReplacement
{
variable_1
{
internalField   uniform 200;
boundaryField
{
boundary_1
{
...
}
...
boundary_n
{
...
}
}
...
variable_n
{
...
}
}


The BC used in the waterConvection patch is worth a mention, that is, externalWallHeatFluxTemperature. This is a BC that allows you to define a condition of heat transfer by convection with the environment only by giving a value for the ambient temperature (Ta) and the heat transfer coefficient (h) or, otherwise, you can also specify a constant heat flux (q). Besides that, more entries need to be specified, such as the way thermal conductivity will be computed (kappa) according to the physics of the case. For a further description of this BC take a look at the corresponding code file by executing the line below

  gedit $FOAM_SRC/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.H &  This BC can also be defined by means of the use of groovyBC, as we did in the previous tutorial. However, it is better to use original BC’s provided with OpenFOAM instead of using our custom BC’s, since they are supposed to work in a more optimal way. Nevertheless, the groovyBC specification is also provided, as you should have noticed, and commented out just in case you want to test that the two BC’s do exactly the same. And finally, the most remarkable BC and the most important one when it comes to a multiregion case is the one used to couple regions, this is compressible::turbulentTemperatureCoupledBaffleMixed and it can couple regions no mather what type they are. An interesting option of this BC is that it allows you to define a thermal contact resistance between regions by giving the thikness and the conductivity of the interregion layers. For further information about this BC take a look at the corresponding code file by executing the line below  gedit$FOAM_SRC/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.H &


After having defined properly the BC’s for all the regions it’s time to update them by running the application as follows

  changeDictionary -region


# 4th step: defining the heat source

Now it’s time to define the heat source. The procedure will be almost the same as in the previous case so I will not expand this section too much. Again, we will use a source of the type semiImplicitSource and, by following the same procedure we developed in the previous case, in the end we get

$g = fvOptions(rho, h) = S_u + S_p \cdot h$

Which gives us the values of  $S_u = g$  and  $S_p = 0$, as you can see, this time the values to use in fvOptions are more straightforward. 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/source/fvOptions &

Notice that this time the fvOptions file is placed within the region directory under system, not within the system directory itself. Below you can find the specification used for the case here studied, with the main differences from the previous case bolded.

volumetricHeatSource
{
type            scalarSemiImplicitSource;
active          true;
selectionMode   all;  //cellSet // points //all  (other options)

scalarSemiImplicitSourceCoeffs
{
volumeMode      absolute;
injectionRateSuSp
{
h (3 0);   // S(T)=S_u+S_p*T  ===>  S_u=g    S_p=0
}              // g= specific volumetric generation (W/m3)
}                  // g= absolute volumetric generation  ===> in the example 3W are considered
}

# 5th step: getting the solver ready

As it has been stated in the previous tutorial, this can be a tedious task when you are a beginner. However, this step is almost the same for both tutorials, they only differ in the variables computed in both solvers since one computes the temperature (T) and the other one calculates the enthalpy (h). Besides that, as it is a multi region solver, it needs that each region has their own files –fvSolution and fvSchemes– properly defined, so they will be defined within the system/ directory. Let’s open these files!

  gedit system/{source,base}/{fvSolution,fvSchemes} &


The last file to define –controlDict– it is placed in the same directory –system– regardless of the solver to be considered. The whole files can be found on the following links –fvSolution and fvSchemes have been considered to be the same in both regions.

Tip: 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.

# 6th step: running the solver

If you have followed all the previous steps correctly now it’s time to solve the case by running the solver, that is, typing its name in the terminal.

  chtMultiRegionSimpleFoam
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.x                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Exec   : chtMultiRegionSimpleFoam
Date   : Jun 19 2016
Time   : 19:21:24
Host   : "linux-cfd"
PID    : 6937
Case   : /home/cfd/OpenFOAM/cfd-2.3.x/run/semicylinder_withGeneration_multiregion
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 solid mesh for region base for time = 0

Create solid mesh for region source for time = 0

...[lots of lines, 32832 to be more precise]

Solving for solid region base
DICPCG:  Solving for h, Initial residual = 9.90509e-07, Final residual = 9.90509e-07, No Iterations 0
Min/max T:min(T) [0 0 0 1 0 0 0] 36.6266 max(T) [0 0 0 1 0 0 0] 200

Solving for solid region source
DICPCG:  Solving for h, Initial residual = 9.38133e-07, Final residual = 9.38133e-07, No Iterations 0
Min/max T:min(T) [0 0 0 1 0 0 0] 180.549 max(T) [0 0 0 1 0 0 0] 191.665
ExecutionTime = 1.48 s  ClockTime = 2 s

End

Once the execution has been finished it’s time to check that the solver has run properly and convergence has been reached. To do so, it’s necessary to check the residuals. There are several ways to plot the residuals and if you are a gnuplot geek you can create your own script in order to get your own residual plot. However, as it is intended to be a beginner’s guide I will assume that you are not a gnuplot geek and I will introduce you to an application created with the aim of ease and speed up the workfow of a case solved with OpenFOAM. This software is called pyFoam and one of the many tasks that this application allows you to carry out is the residual plots generation. Now I will show you how to use it in order to plot the residuals for all the regions, which is kind of tricky…

The tool to be used in order to study the residuals after the solver has finished its calculations is called pyFoamPlotWatcher and it is very easy to use for a single region case, you only need to invoke the tool and pass the log file of the run as an argument. Something like the following would be enough

  pyFoamPlotWatcher.py --silent --progress log.solverName &


However, if you execute the same command in a multi region case it will only plot the residual plot (and some others) for the first region solved by OpenFOAM, the rest of the regions will be omited. In order to take all the regions (or the ones you want to analyze) into consideration a “magic” text file must be included in the case directory. This file must be named as customRegexp and the lines it has to contain should look like the following ones

base_residuals
{
expr "Solving for solid region base\nDICPCG:  Solving for h, Initial residual = (%f%), Final residual = .+, No Iterations .+";
name Custom_base_res;
theTitle "Residuals for region: base";
titles
(
"h (1st iter.)"
);
type regular;
logscale true;
}

source_residuals
{
expr "Solving for solid region source\nDICPCG:  Solving for h, Initial residual = (%f%), Final residual = .+, No Iterations .+";
name Custom_source_res;
theTitle "Residuals for region: source";
titles
(
"h (1st iter.)"
);
type regular;
logscale true;
}

As it can be seen, there is a set of lines for each of the regions that make up the case, that is, for each of the plots we want to plot. The first line is the name of the dictionary that contains the set of instructions to draw the plot.
The most important field within the dictionary is the keyword expr which contains the expression that pyFoamPlotWatcher will look for in the solver log file. This expression has been copy-pasted directly from the log file and afterwards some tricky modifications have been done. The first one is the substitution of the newline appearing on the log file for \n, which are the characters that the application need to read a newline. If you leave the newline, the application won’t be able to read properly the expression, so no information will be plotted. Also, the numerical values of the residuals and iterations have been replaced by some mysterious characters: the first ones ((%f%)) indicate the value we want to plot, whereas the others (.+) indicate the values the application will omit.
The keyword titles is a list of the labels to be shown for each of the plotted values. In this case only the residuals for the enthalpy (h) are plotted for the 1st iteration.
And last but not least, the keyword logscale which is required so that the plots will be shown using logarithmic scale, otherwise these plots won’t be much useful…

Once this file has been properly defined you can run the following command so as to analyze the residual plots for the whole set of regions that comprise the case

  pyFoamPlotWatcher.py --silent --progress --no-default log.chtMultiRegionSimpleFoam &


After the application has finished running (it may take a few seconds) you will be able to see the following plots.

Note that once the residual tolerance has been reached the solver stops iterating, so we can accept that the solver has converged to a solution for both regions.