CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Boundary condition for epsilon as a function of k, y, z and time

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 12, 2021, 11:12
Default Boundary condition for epsilon as a function of k, y, z and time
  #1
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 22
Rep Power: 2
subhojitkadiacfd is on a distinguished road
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 * * * * * * * * * * * * * //

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

Last edited by subhojitkadiacfd; May 15, 2021 at 06:55.
subhojitkadiacfd is offline   Reply With Quote

Old   May 13, 2021, 10:43
Default
  #2
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 22
Rep Power: 2
subhojitkadiacfd is on a distinguished road
Hello,
I have successfully compiled the boundary condition

Code:
............................../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
Code:
--> 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.

Last edited by subhojitkadiacfd; May 15, 2021 at 06:57.
subhojitkadiacfd is offline   Reply With Quote

Old   May 13, 2021, 18:10
Default
  #3
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,596
Blog Entries: 6
Rep Power: 47
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
Use the codedFixedValue boundary condition.
__________________
Keep foaming,
Tobias Holzmann
Tobi is online now   Reply With Quote

Old   May 15, 2021, 09:50
Default
  #4
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 22
Rep Power: 2
subhojitkadiacfd is on a distinguished road
Many thanks Tobi for you suggestion.
How can I access the k values for every iteration (time step)?

Subhojit
subhojitkadiacfd is offline   Reply With Quote

Old   May 15, 2021, 12:23
Default
  #5
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,596
Blog Entries: 6
Rep Power: 47
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to 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:


Code:

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.
__________________
Keep foaming,
Tobias Holzmann
Tobi is online now   Reply With Quote

Old   May 15, 2021, 16:44
Default
  #6
HPE
Senior Member
 
HPE's Avatar
 
Herpes Free Engineer
Join Date: Sep 2019
Location: The Home Under The Ground with the Lost Boys
Posts: 903
Rep Power: 10
HPE is on a distinguished road
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.
HPE is offline   Reply With Quote

Old   May 16, 2021, 05:12
Default
  #7
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,596
Blog Entries: 6
Rep Power: 47
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to 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.
__________________
Keep foaming,
Tobias Holzmann
Tobi is online now   Reply With Quote

Old   May 16, 2021, 05:47
Default
  #8
HPE
Senior Member
 
HPE's Avatar
 
Herpes Free Engineer
Join Date: Sep 2019
Location: The Home Under The Ground with the Lost Boys
Posts: 903
Rep Power: 10
HPE is on a distinguished road
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.
HPE is offline   Reply With Quote

Old   May 16, 2021, 11:57
Default
  #9
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 22
Rep Power: 2
subhojitkadiacfd is on a distinguished road
Many thanks Tobi and HPE for your valuable suggestions. As per my understanding I have tried the following two options:

Case1
Code:
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:
...........
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
subhojitkadiacfd is offline   Reply With Quote

Old   May 16, 2021, 16:41
Default
  #10
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,596
Blog Entries: 6
Rep Power: 47
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
In the codedFixedValue you need the operator==


Code:
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.
__________________
Keep foaming,
Tobias Holzmann
Tobi is online now   Reply With Quote

Old   May 22, 2021, 16:45
Default
  #11
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 22
Rep Power: 2
subhojitkadiacfd is on a distinguished road
Hi Tobi,
Many thanks for your suggestions. The code is working properly. Hoping to get the desired result soon.
Subhojit
subhojitkadiacfd is offline   Reply With Quote

Old   June 7, 2021, 07:49
Default
  #12
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 22
Rep Power: 2
subhojitkadiacfd is on a distinguished road
Further improvement was required, so I tried to calculate the average distance from the wall. I added the following code within the codedFixedValue condition
Code:
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

Last edited by subhojitkadiacfd; June 8, 2021 at 07:04.
subhojitkadiacfd is offline   Reply With Quote

Old   June 8, 2021, 08:05
Default
  #13
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,596
Blog Entries: 6
Rep Power: 47
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
Not sure what you calculate but in a programming point of view there are some mistakes. I just grab one snippet:

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

Code:
const scalar foobar = 1;      // Variable 1
const scalar foobar = 2;      // Variable 2 - compiler error as foobar already exists
Code:
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:
Code:
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:
Code:
const scalar foobar = 0;
foobar = 2;    // Compiler error as foobar is in read-only mode (const)
So what you need is the following:
Code:
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++.
__________________
Keep foaming,
Tobias Holzmann

Last edited by Tobi; June 8, 2021 at 09:33.
Tobi is online now   Reply With Quote

Old   June 8, 2021, 08:41
Default
  #14
New Member
 
Join Date: Jan 2021
Location: Norway
Posts: 22
Rep Power: 2
subhojitkadiacfd is on a distinguished road
Many thanks Tobi.
subhojitkadiacfd is offline   Reply With Quote

Reply

Tags
boundary condition, codedfixedvalue

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
CFD analaysis of Pelton turbine amodpanthee CFX 31 April 19, 2018 19:02
Wrong flow in ratating domain problem Sanyo CFX 17 August 15, 2015 07:20
Time dependant pressure boundary condition yosuke1984 OpenFOAM Verification & Validation 3 May 6, 2015 07:16
Upgraded from Karmic Koala 9.10 to Lucid Lynx10.04.3 bookie56 OpenFOAM Installation 8 August 13, 2011 05:03
Convective Heat Transfer - Heat Exchanger Mark CFX 6 November 15, 2004 16:55


All times are GMT -4. The time now is 18:59.