[subhojitkadiacfd]

Hi all,
I want to set a boundary condition for epsilon on topSurface in my channel which is a function of k and coordinates y and z. It shall also take the k data from the topSurface for each iteration (time dependent) and calculate epsilon accordingly. I am using v2012 and v 8.0. What would be the best way of doing it?
With the help of one member I tried to create a new boundary condition. The member function is as follows:

// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

void Foam::epsilonFreeSurfaceFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }


    const IOdictionary& d = static_cast<const IOdictionary&>(dict.parent().parent());
    const fvMesh& mesh = refCast<const fvMesh>(d.db());
    const label id = mesh.boundary().findPatchID("topSurface");
    const fvPatch& patch = mesh.boundary()[id];
	const scalar t = this->db().time().timeOutputValue();//Get the current time step
	const volScalarField& k = mesh.lookupObject<volScalarField>("k");
	scalarField kp = k.boundaryField()[id];
	scalarField epsp(patch.size(), 0.0);
				
    if (t >= t) 		//Perform calculation at each time step
		//Info << "t = " << t << endl;
	{
		forAll(epsp, facei)
			{
				scalar y_ = patch.Cf()[facei].y();
				scalar z_ = patch.Cf()[facei].z();
				scalar k_ = kp[facei];
				epsp[facei] = k_*y_*z_; //say
			}
			
	}			
}

The error I am getting is
error: ‘dict’ was not declared in this scope

Do I need to do anything about the write entry for epsp also?

Is the time calling is okay?

Any help in improving/modifying/correcting the code would be highly appreciated.

Thank you.

[subhojitkadiacfd]

Hello,
I have successfully compiled the boundary condition

............................../mnt/c/Windows/System32/epsilonFreeSurface$ wmake libso
wmake libso (epsilonFreeSurface)
Making dependency list for source file epsilonFreeSurfaceFvPatchScalarField.C
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2012 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas  -O3  -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2012/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2012/src/meshTools/lnInclude -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2012/src/OSspecific/POSIX/lnInclude   -fPIC -c epsilonFreeSurfaceFvPatchScalarField.C -o Make/linux64Gcc63DPInt32Opt/epsilonFreeSurfaceFvPatchScalarField.o
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2012 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas  -O3  -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2012/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2012/src/meshTools/lnInclude -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2012/src/OSspecific/POSIX/lnInclude   -fPIC -shared -Xlinker --add-needed -Xlinker --no-as-needed  Make/linux64Gcc63DPInt32Opt/epsilonFreeSurfaceFvPatchScalarField.o -L/opt/OpenFOAM/OpenFOAM-v2012/platforms/linux64Gcc63DPInt32Opt/lib \
    -lfiniteVolume -lmeshTools  -o /home/subhojik/OpenFOAM/subhojik-v2012/platforms/linux64Gcc63DPInt32Opt/lib/libepsilonFreeSurface.so

And added the library in the controlDict

But when running the simulation I am getting the following error

--> FOAM Warning :
    From void* Foam::dlLibraryTable::openLibrary(const Foam::fileName&, bool)
    in file db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C at line 188
    Could not load "libepsilonFreeSurface.so"
/..........................v2012/platforms/linux64Gcc63DPInt32Opt/lib/libepsilonFreeSurface.so: undefined symbol: _ZTIN4Foam36epsilonFreeSurfaceFvPatchScalarFieldE
Create mesh for time = 0

[1] --> FOAM FATAL IO ERROR: (openfoam-2012)
[1] Unknown patchField type epsilonFreeSurface for patch type patch

Any help/suggestion would be highly appreciated.

[Tobi]

Use the codedFixedValue boundary condition.

[subhojitkadiacfd]

Many thanks Tobi for you suggestion.
How can I access the k values for every iteration (time step)?

Subhojit

[Tobi]

Use google to get more information regarding the codedFixedValue condition.
You can get easy access to the "k" field inside the "epsilon" using something link that:

const volScalarField& k = this->db().lookupObject<volScalarField>("k");


// Then you need to get access to the patch fields and to the corresponding boundary

The idea of getting the data of k at the given boundary is given here: https://bitbucket.org/shor-ty/laserconvectionbc/src/3a413984d3cc96d3b8b6f8da82ec40d349a7be65/laserConvectionFvPatchField.C#lines-981
Keep in mind, «k» does not represent the turbulentKineticEnergy here and I worked with templates here.
You don't have to take care about time-steps etc.

Whenever we solve the matrix equations, we update the boundaries, hence we always do have actual values here.

[HPE]

In addition to Tobi's remarks, I would also give a try to use PatchFunction1 which allows you to use a myriad of ways to manipulate boundary conditions in a much generic way.

You can then use a coded or expression-based boundary condition or any other Function1 input types in your input dictionary by using PatchFunction1 in your code.

You can find some examples in the boundary conditions under src/atmosphericModels.

Hope helps.

[Tobi]

Hi,


good hint but at the end it is equal to the codedFixedValue if he uses the code Function1 option. As he wants to have access to other fields in the boundary condition, one needs to use a coded one - don't matter if it is the uniformFixedValue or codedFixedValue. Al least as far as I know

Nevertheless, I agree, one can also use the uniformFixedValue condition.

[HPE]

Hi,

If both options would require the same effort, I would use the generic one. Though, PathFunction1 would be a bit more involved. Depends what the inquirer wants. Both ways are good, anyhow.

[subhojitkadiacfd]

Many thanks Tobi and HPE for your valuable suggestions. As per my understanding I have tried the following two options:

Case1

topSurface
    {
		type            codedFixedValue;//changing with k
		value           uniform 0.00346;
		name		    myEpsilonFS;
		code
		#{
			const scalar b = 0.3; //channel width b (m) 
			const scalar h = 0.106; // flow depth h (m)
			const fvMesh& mesh = patch().boundaryMesh().mesh();
			const label id = mesh.boundary().findPatchID("topSurface");
			const fvPatch& boundaryPatch = mesh.boundary()[id];;
			const volScalarField& k = this->db().lookupObject<volScalarField>("k");
			scalarField kp = k.boundaryField()[id];
			scalarField epsilon(boundaryPatch.size(), 0.00346);
			
			forAll(epsilon, facei)
			{
				const scalar y_ = boundaryPatch.Cf()[facei][1];// 1 for y coordinate, 2 for z coordinate
				const scalar z_ = boundaryPatch.Cf()[facei][2];
				const scalar y1_ = 0.07*h;
				const scalar y2_ = min((0.5*b - mag(z_)), mag(y_));
				const scalar k_ = kp[facei];
				epsilon[facei] = 0.391*pow(k_, 1.5)*(1.0/y1_ + 1.0/y2_); 
			}
			// epsilon (free surface) = 0.391*k(free surface)^3/2*(1/y*+1/y')
		#};

		codeOptions
			#{
				-I$(LIB_SRC)/finiteVolume/lnInclude\
				-I$(LIB_SRC)/meshTools/lnInclude\
				-I$(LIB_SRC)/TurbulenceModels/lnInclude
			#};
			
		codeInclude
			#{
				#include "fvCFD.H"
				#include <cmath>
				#include <iostream>
			#};
			
		codeLibs 
			#{ -lfiniteVolume -lmeshTools 
			#};
	}

Case2:

...........
code
		#{
			const scalar b = 0.3; //channel width b (m) 
			const scalar h = 0.106; // flow depth h (m)
			const fvPatch& boundaryPatch = this->patch();
			const vectorField& Cf = boundaryPatch.Cf();
			const label& patchLabel = boundaryPatch.index();
			const volScalarField& k = this->db().lookupObject<volScalarField>("k");
			fvPatchScalarField kp = k.boundaryField()[patchLabel];
			scalarField& field = *this;
						
			forAll(Cf, facei)
			{
				const scalar y_ = Cf[facei].y();
				const scalar z_ = Cf[facei].z();
				const scalar y1_ = 0.07*h;
				const scalar y2_ = min((0.5*b - mag(z_)), mag(y_));
				const scalar k_ = kp[facei];
				field[facei] = 0.391*pow(k_, 1.5)*(1.0/y1_ + 1.0/y2_); 
			}
			// epsilon (free surface) = 0.391*k(free surface)^3/2*(1/y*+1/y')
		#};
............

However, the results are not satisfactory.

Could you please check the code and let me know if I need any corrections? Thanks.

Subhojit

[Tobi]

In the codedFixedValue you need the operator==

operator==(epsilon);

However, the second one should work. If the values do not fit, you have to check it. We cannot give you an answer on your calculation.

[subhojitkadiacfd]

Hi Tobi,
Many thanks for your suggestions. The code is working properly. Hoping to get the desired result soon.
Subhojit

[subhojitkadiacfd]

Further improvement was required, so I tried to calculate the average distance from the wall. I added the following code within the codedFixedValue condition

forAll(Cf, facei)
			{
				const scalar z_ = Cf[facei].z();// lateral coordinate of the face centre
				const scalar y_1 = 0.07*h;
				const scalar z_1 = 0.5*b - mag(z_);// horizontal distance from the left wall 
				const scalar k_ = kp[facei];
if (z_1 <= h)
				{
					const int n_1 = 100; // taking 100 parts along the flow depth
					const scalar delta_h = h/n_1; 
					const scalar sum_1 = 0.0;// initial value
					for (int i = 0; i <= n_1; ++i)
					{
						const scalar h_1 = delta_h*i;
						const scalar d_1 = sqrt(sqr(h_1) + sqr(z_1));// distance between the points
						const scalar sum_1 = sum_1 + d_1;
					}
					const scalar y_2 = sum_1/(n_1 + 1);// average distance
					epsilon[facei] = 0.401*pow(k_, 1.5)*(1.0/y_1 + 1.0/y_2);
				}
.......
.......

However, the simulation could not calculate epsilon in iteration 1. I think the for loop is not returning sum_1 correctly and needs some improvements.
Thanks in advance.

Subhojit

[Tobi]

Not sure what you calculate but in a programming point of view there are some mistakes. I just grab one snippet:

const scalar sum_1 = 0.0;
for (int i = 0; i <= n_1; ++i)
{
    const scalar h_1 = delta_h*i;
    const scalar d_1 = sqrt(sqr(h_1) + sqr(z_1));
    const scalar sum_1 = sum_1 + d_1;
}

const scalar y_2 = sum_1/(n_1 + 1);

You define a variable named sum_1 to be initialized with the double zero. However, you declare this variable to be in read-only-mode (const). Hence, you are not allowed to change that variable anymore (okay - some dirt stuff works but that is not the aim here).

After that you jump into the for-loop and you define another variable that takes the same name as before. Note, as you are inside a sub-frame (no idea about the correct name), you can create the same variable again but it is not the same as before - it is a new one. The validity of this second sum_1 ends after the bracket and it does not exist anymore. Example:

const scalar foobar = 1;      // Variable 1
const scalar foobar = 2;      // Variable 2 - compiler error as foobar already exists

const scalar foobar = 1;     // Variable 1
{
    const scalar foobar = 2; // Variable 2 - no compiler error
}

here you have two variables (better to say objects) with the same name. Lets go a bit futher:

const scalar foobar = 1;     // Variable 1
{
    Info << foobar << endl;  // Will return 1 as we check Variable 1

    // This definition of the object "foobar" will hide the previous defined one
    // No possible chance to access Variable 1
    const scalar foobar = 2;   // Variable 2 - no compiler error

    Info << foobar << endl;   // Will return 2 as we check Variable 2
}  // After that bracket, the Variable 2 is not valid anymore

Info << foobar << endl;      // Will return 1  as we check Variable 1

Now the semi-correct on:

const scalar foobar = 0;
foobar = 2;    // Compiler error as foobar is in read-only mode (const)

So what you need is the following:

scalar sum_1 = 0.0;
for (int i = 0; i <= n_1; ++i)
{
    const scalar h_1 = delta_h*i;
    const scalar d_1 = sqrt(sqr(h_1) + sqr(z_1));
    sum_1 += d_1;
}

Simply c++.

[subhojitkadiacfd]

Many thanks Tobi.