CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Community Contributions

[swak4Foam] problem with a parabolic velocity profile

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

Like Tree2Likes
  • 2 Post By billie

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 4, 2013, 07:53
Default problem with a parabolic velocity profile
  #1
Member
 
Claudio
Join Date: Mar 2012
Location: Milano, Italy
Posts: 49
Rep Power: 10
Claudio87 is on a distinguished road
Hi all!

I'm trying to apply, as BC, a parabolic velocity profile at the inlet of my problem (a pipe) using groovyBC (from swak4foam).
It works, in the sense that I have a parabolic profile, but:

I would like to have velocity values (only in x direction) between 0 and 2
instead:
I obtain a profile going from -0.11924 (negative!!) to 1.993 , looking only Ux in paraview,
and from 1.929e-5 to 1.993 considering the Umagnitude. (see attached pictures)

I'm worried about the negative value overall. Do you have any idea why I obtain these results instead of what I expected (zero?!

Here how I applied the condition:

Code:
    {
        type			groovyBC;
	value			uniform (1 0 0);
        variables		"Vm=2;yp=pos().y;zp=pos().z+0.005;h=-max(zp);a=-Vm/(h*h);b=Vm;";
	fractionExpression	"1"; 
        valueExpression		"vector(a*((yp*yp)+(zp*zp))+b,0,0)";
    }
I used the expression for a poiseuille flow ( V=Vm*(1-(r/R)^2) ).
I defined zp=pos().z+0.005 because z interval is [-0.0075 : -0.0025].


Thank you in advance!
Claudio
Attached Images
File Type: jpg U_x.jpg (42.4 KB, 152 views)
File Type: jpg U_magnitude.jpg (42.8 KB, 117 views)
Claudio87 is offline   Reply With Quote

Old   June 4, 2013, 13:10
Default
  #2
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,139
Rep Power: 46
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by Claudio87 View Post
Hi all!

I'm trying to apply, as BC, a parabolic velocity profile at the inlet of my problem (a pipe) using groovyBC (from swak4foam).
It works, in the sense that I have a parabolic profile, but:

I would like to have velocity values (only in x direction) between 0 and 2
instead:
I obtain a profile going from -0.11924 (negative!!) to 1.993 , looking only Ux in paraview,
and from 1.929e-5 to 1.993 considering the Umagnitude. (see attached pictures)

I'm worried about the negative value overall. Do you have any idea why I obtain these results instead of what I expected (zero?!

Here how I applied the condition:

Code:
    {
        type			groovyBC;
	value			uniform (1 0 0);
        variables		"Vm=2;yp=pos().y;zp=pos().z+0.005;h=-max(zp);a=-Vm/(h*h);b=Vm;";
	fractionExpression	"1"; 
        valueExpression		"vector(a*((yp*yp)+(zp*zp))+b,0,0)";
    }
I used the expression for a poiseuille flow ( V=Vm*(1-(r/R)^2) ).
I defined zp=pos().z+0.005 because z interval is [-0.0075 : -0.0025].


Thank you in advance!
Claudio
At first: never use point-values when claiming "it is that way". The only way you see the real values are cell values. Point values are paraview interpolating them and what you see might (doesn't have to) be an artefact of the interpolation. For checking boundary conditions it is even better to not import the internalMesh but only the patch in question (both readers support this) and USE THE CELL VALUE

The max not hitting 2 is OK: if the peak of the parabola doesn't "hit" a face-center then you will get a slightly lower value. The "+0.005"-trick can be avoided. Use max(pts().z) to get the maximum of the vertices of the patch. Get the minimum the same way (see how it is done in the pulsating pitzDaily).
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request
gschaider is offline   Reply With Quote

Old   June 5, 2013, 05:41
Default
  #3
Member
 
Claudio
Join Date: Mar 2012
Location: Milano, Italy
Posts: 49
Rep Power: 10
Claudio87 is on a distinguished road
Hi Bernhard,
thank you for your quick answer!

I tried to do what you said, and now it works!
I mean:

Quote:
For checking boundary conditions it is even better to not import the internalMesh but only the patch in question (both readers support this) and USE THE CELL VALUE
I was yet looking to the patch and not to the internalMesh; and using the cell value, it seems the same (look the attached pictures called "old").


Quote:
The "+0.005"-trick can be avoided. Use max(pts().z) to get the maximum of the vertices of the patch. Get the minimum the same way (see how it is done in the pulsating pitzDaily).
I modified the "trick" in this way:

Code:
    {
        type			groovyBC;
	value			uniform (1 0 0);
        variables		"Vm=2;yp=pos().y;zp=pos().z+0.005;h=-max(pts().z);a=-Vm/(h*h);b=Vm;";
	fractionExpression	"1"; 
        valueExpression		"vector(a*((yp*yp)+(zp*zp))+b,0,0)";
    }
and now it works (see the "new" pict), even if I didn't understand exactly why!
If you could explain it to me someway it could be usefull!

Thank you very much!
Attached Images
File Type: jpg U_magnitudeCell(old).jpg (53.1 KB, 87 views)
File Type: jpg U_xCell(old).jpg (52.8 KB, 80 views)
File Type: jpg U_x(new).jpg (50.4 KB, 96 views)
File Type: jpg geometry.jpg (39.7 KB, 75 views)
Claudio87 is offline   Reply With Quote

Old   June 18, 2013, 09:34
Default
  #4
Member
 
Daniel Pielmeier
Join Date: Apr 2012
Posts: 97
Rep Power: 10
billie is on a distinguished road
Quote:
Originally Posted by Claudio87 View Post
I modified the "trick" in this way:

Code:
    {
        type            groovyBC;
    value            uniform (1 0 0);
        variables        "Vm=2;yp=pos().y;zp=pos().z+0.005;h=-max(pts().z);a=-Vm/(h*h);b=Vm;";
    fractionExpression    "1"; 
        valueExpression        "vector(a*((yp*yp)+(zp*zp))+b,0,0)";
    }
and now it works (see the "new" pict), even if I didn't understand exactly why!
If you could explain it to me someway it could be usefull!

Thank you very much!
I think what Bernhard meant by avoiding the "+0.005"-trick is that you can acquire the the min and max position of the patch and and use them to calculate the offsets from the origin. This way you don't have to worry where your patch is positioned but only in which plane (e. g. XZ) it is located.

I tried the same but wanted to achieve a fixed volumetric flow rate. I post this here, maybe someone can make use of it. If somebody has a better solution I would be glad to hear it.

First I obtain the list of face centres for the x and z coordinate (xp, zp) corrected by the offset between the centre of the pipe inlet and the origin. This way the maximum value is at the centre of the inlet patch. For the parabolic profile (para) I used the same function Claudio used. What I call normalizedFlux is the resulting flux for a parabola with a maximum velocity of 1m/s. The ratio between my desired flux and the normalizedFlux is then used in the valueExpression along with the parabolic profile. This results in the same flux as for the flowRateInletVelocity boundary condition with a uniform value.

Code:
    inflow
    {
        //type            flowRateInletVelocity;
        //volumetricFlowRate 0.000583333;
        //value           uniform (0 0 0);
        type            groovyBC;
        variables (
            "flux=-0.000583333;"
            "minX=min(pts().x);"
            "maxX=max(pts().x);"
            "minZ=min(pts().z);"
            "maxZ=max(pts().z);"
            "r=mag((maxZ-minZ)/2.0);"
            "xOffset=minX+r;"
            "zOffset=minZ+r;"
            "xp=pos().x-xOffset;"
            "zp=pos().z-zOffset;"
            "para=1-((xp*xp)+(zp*zp))/pow(r,2);"
            "normalizedFlux=sum(para*area());"
            "ratio=flux/normalizedFlux;"
        );
        valueExpression "vector(0,ratio*para,0)";
        value           uniform (0 1 0);
    }
arvindpj and Romulus like this.
billie is offline   Reply With Quote

Old   May 29, 2014, 05:12
Default
  #5
Senior Member
 
Samuele Z
Join Date: Oct 2009
Location: Mozzate - Co - Italy
Posts: 516
Rep Power: 14
samiam1000 is on a distinguished road
Dear All,

I am trying to do the same trick, using the coded BC instead of using swak4Foam, since I can't compile it on the machine where I have to launch my simulations.

Do you know how I can fix this?

I tried something like
Code:
    movingWall
	{
		type codedFixedValue; 
		value uniform (0 -2 0); 
		redirectType parabolicInlet; 
		code 
		#{
			scalar U_0 = 2;
			scalar a = 0;
			scalar yp = 0;
			scalar zp = 0;
			fixedValueFvPatchVectorField myPatch(*this);
			forAll(this->patch().Cf(),i)
			{
				a = U_0 / (0.5*0.5);
				yp = pos().y;
				zp = pos().z+0.005;
				myPatch=vector(0,-2*a,0);
			}
			operator==(myPatch);
		#};

	}
but it does not work properly, since it "does not like" the lines where I evaluate yp and zp.

Could you help, please?

Thanks,
Samuele
samiam1000 is offline   Reply With Quote

Old   May 29, 2014, 09:30
Default
  #6
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,139
Rep Power: 46
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by samiam1000 View Post
Dear All,

I am trying to do the same trick, using the coded BC instead of using swak4Foam, since I can't compile it on the machine where I have to launch my simulations.

Do you know how I can fix this?

I tried something like
Code:
    movingWall
	{
		type codedFixedValue; 
		value uniform (0 -2 0); 
		redirectType parabolicInlet; 
		code 
		#{
			scalar U_0 = 2;
			scalar a = 0;
			scalar yp = 0;
			scalar zp = 0;
			fixedValueFvPatchVectorField myPatch(*this);
			forAll(this->patch().Cf(),i)
			{
				a = U_0 / (0.5*0.5);
				yp = pos().y;
				zp = pos().z+0.005;
				myPatch=vector(0,-2*a,0);
			}
			operator==(myPatch);
		#};

	}
but it does not work properly, since it "does not like" the lines where I evaluate yp and zp.

Could you help, please?

Thanks,
Samuele
You know that you've succeeded when people start using swak-idioms in their C++-code
I think what you want is (I'm doing this from memory so it might be slightly off. Especially the .y()-part)
Code:
yp=this->patch().Cf()[i].y();
zp=this->patch().Cf()[i].z()+0.005;
Although I honestly don't understand why you'd want to do that. These values are not used afterwards.
Also replace the next line with
Code:
myPatch[i]=vector(0,-2*a,0);
because otherwise the whole myPatch would have the value from the last iteration (not that it matters here because it is the same value anyway for the present code)

Concerning "can't compile it on the machine where I have to launch my simulations": if it is about outdated bison&flex on that machine: the development-version (and next release) of swak have a script to download and build a more up-to-date version of bison and use that. Nothing more than a C++-compiler is needed. And you must have access to that on that machine (otherwise coded wouldn't work)
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request
gschaider is offline   Reply With Quote

Reply

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
3D UDF Paraboilc Velocity Profile (Can't Maintain) Sing FLUENT 12 August 7, 2017 06:25
atmBoundaryLayerInletVelocity - Velocity Profile not continuous through domain sdfij6354 OpenFOAM Running, Solving & CFD 3 July 26, 2017 16:16
pressure contour and velocity profile problem alee1293 Main CFD Forum 0 July 18, 2014 10:31
Multi step transient UDF velocity profile problem William177 FLUENT 1 February 3, 2008 06:47
Variables Definition in CFX Solver 5.6 R P CFX 2 October 26, 2004 02:13


All times are GMT -4. The time now is 01:10.