Excercise on advection schemes
In this excercise we will compare the performance of the different advection schemes. We will use the general advection-diffusion from previsou sections:
Solver
OpenFoam has a simplified solver for the advection-diffusion equation called scalarTransportFoam, which uses a prescribed velocity field to move the scalar field around.
As always, a good starting point is find a suitable tutorial, copy it over to your own case and start modifying it. In this case, we will use the tutorial basic/scalarTransportFoam/pitzDaily as a starting point. Copy the case to your own case directory and give it a reasonable name.
cp -r $FOAM_TUTORIALS/basic/scalarTransportFoam/pitzDaily $HOME/HydrothermalFoam_runs
The solver is set up to solve the advection-diffusion equation for a scalar field called T. The velocity field is prescribed in the file 0/U and the diffusion coefficient is set in constant/transportProperties. The initial condition is set in 0/T.
Mesh and initial conditions
Mesh
We will use a 2D box with dimensions \(-0.5 < x < 0.5\) and \(-0.5 < y < 0.5\). You will need to modify the blockMeshDict file to set up the mesh. The mesh should be uniform with 100 cells in each direction. The mesh should be centered around the origin.
Temperature
The scalar field initial conditions are set in the file 0/T. We will use a Gaussian distribution with a standard deviation of \(\sigma = 0.1\) and a maximum value of \(T_0 = 2\). It should be centered at x. Let’s put the intitial gaussian at \(x_0 = 0\) and \(y_0 = 0.25\).
Velocity
We will fist look into solid body rotation, in which the initial gaussian will be rotated in clockwise direction without any shearing. The respective velocity field is given by:
Implementation
To implement the initial conditions and prescribed velocity field, you will need to modify the files 0/T and 0/U. We will again use codestream statements for this.
1/*--------------------------------*- C++ -*----------------------------------*\
2========= |
3\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration | Website: https://openfoam.org
5 \\ / A nd | Version: 7
6 \\/ M anipulation |
7\*---------------------------------------------------------------------------*/
8FoamFile
9{
10 version 2.0;
11 format ascii;
12 class volScalarField;
13 object T;
14}
15// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
16
17dimensions [0 0 0 1 0 0 0];
18
19//internalField uniform 278.15;
20internalField #codeStream
21{
22 codeInclude
23 #{
24 #include "fvCFD.H"
25 #};
26 codeOptions
27 #{
28 -I$(LIB_SRC)/finiteVolume/lnInclude \
29 -I$(LIB_SRC)/meshTools/lnInclude
30 #};
31 codeLibs
32 #{
33 -lmeshTools \
34 -lfiniteVolume
35 #};
36 localCode
37 #{
38 static double calTemperature(const scalar sigma, const scalar maxT, scalar x, scalar y)
39 {
40 return ???;
41
42 }
43 #};
44 code
45 #{
46 const IOdictionary& d = static_cast<const IOdictionary&>(dict);
47 const fvMesh& mesh = refCast<const fvMesh>(d.db());
48 scalarField T(mesh.nCells(), 0);
49 scalar sigma = .1, maxT = 2;
50
51 forAll(T, i)
52 {
53 const scalar x = mesh.C()[i][0];
54 const scalar y = mesh.C()[i][1];
55 T[i]=calTemperature(sigma, maxT,x, y);
56 }
57
58
59 writeEntry(os, "", T); //
60 #};
61};
62
63
64
65boundaryField
66{
67 left
68 {
69 type zeroGradient;
70 }
71
72 right
73 {
74 type zeroGradient;
75 }
76
77 top
78 {
79 type zeroGradient;
80 }
81
82 bottom
83 {
84 type zeroGradient;
85 }
86
87 frontAndBack
88 {
89 type empty;
90 }
91}
92
93// ************************************************************************* //
We can do something similar for the velocity field. The velocity field is set in the file 0/U. Again, we will use codestream statements to set the velocity field.
1/*--------------------------------*- C++ -*----------------------------------*\
2========= |
3\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration | Website: https://openfoam.org
5 \\ / A nd | Version: 7
6 \\/ M anipulation |
7\*---------------------------------------------------------------------------*/
8FoamFile
9{
10 version 2.0;
11 format ascii;
12 class volVectorField;
13 location "0";
14 object U;
15}
16// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
17
18dimensions [0 1 -1 0 0 0 0];
19
20internalField #codeStream
21{
22 codeInclude
23 #{
24 #include "fvCFD.H"
25 #};
26 codeOptions
27 #{
28 -I$(LIB_SRC)/finiteVolume/lnInclude \
29 -I$(LIB_SRC)/meshTools/lnInclude
30 #};
31 codeLibs
32 #{
33 -lmeshTools \
34 -lfiniteVolume
35 #};
36 code
37 #{
38 const double pi = 3.141592653589793;
39 const IOdictionary& d = static_cast<const IOdictionary&>(dict);
40 const fvMesh& mesh = refCast<const fvMesh>(d.db());
41 volVectorField& U = const_cast<volVectorField&>(
42 mesh.lookupObject<volVectorField>("U")
43 );
44
45 forAll(U, celli)
46 {
47 const point& coord = mesh.C()[celli];
48 U[celli].x() = ???;
49 U[celli].y() = ???;
50 U[celli].z() = 0;
51 }
52
53 writeEntry(os, "", U); //
54 #};
55};
56
57
58boundaryField
59{
60 left
61 {
62 type zeroGradient;
63 }
64 right
65 {
66 type zeroGradient;
67 }
68 top
69 {
70 type zeroGradient;
71 }
72 bottom
73 {
74 type zeroGradient;
75 }
76 frontAndBack
77 {
78 type empty;
79 }
80}
81
82
83// ************************************************************************* //
Setting up the case
Since we want to look at advection without any physical diffusion, we also need to set the diffusion coefficient to something small. This happens in constant/transportProperties:
1/*--------------------------------*- C++ -*----------------------------------*\
2========= |
3\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration | Website: https://openfoam.org
5 \\ / A nd | Version: 7
6 \\/ M anipulation |
7\*---------------------------------------------------------------------------*/
8FoamFile
9{
10 version 2.0;
11 format ascii;
12 class dictionary;
13 location "constant";
14 object transportProperties;
15}
16// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
17
18DT DT [0 2 -1 0 0 0 0] 1e-9;
19
20
21// ************************************************************************* //
Finally, we need to set the time step and simulation time. This happens in system/controlDict.
1/*--------------------------------*- C++ -*----------------------------------*\
2========= |
3\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration | Website: https://openfoam.org
5 \\ / A nd | Version: 7
6 \\/ M anipulation |
7\*---------------------------------------------------------------------------*/
8FoamFile
9{
10 version 2.0;
11 format ascii;
12 class dictionary;
13 location "system";
14 object controlDict;
15}
16// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
17
18application scalarTransportFoam;
19
20startFrom startTime;
21
22startTime 0;
23
24stopAt endTime;
25
26endTime 20;
27
28deltaT 0.0025;
29
30writeControl timeStep;
31
32writeInterval 20;
33
34purgeWrite 0;
35
36writeFormat ascii;
37
38writePrecision 6;
39
40writeCompression off;
41
42timeFormat general;
43
44timePrecision 6;
45
46runTimeModifiable true;
47
48
49// ************************************************************************* //
We here use a constant rund time of 20 and a time step length of \(\Delta t = 0.0025\); we save the output every 20 steps. To run the case, call blockMesh and then the solver scalarTransportFoam.
Excercise
Finally, change the advection scheme (e.g. vanLeer, upwind, MUSCL) and compare the results in paraview! A good way is to make multiple directions for the individual cases.
Now, change the velocity field to a shear shell and repeat the exercise.