Output of velocity gradients

wersoe · July 11, 2011, 09:32

Dear Foamers,

I am quite new to OpenFOAM and try to modify one solver to write the velocity gradients as output. The solver used is interDynMFoam.

I started to add divergence of U as output as follows, and this is working good (green in code below)

The next step was to add the velocity gradients output, as you seen in red in code below.

First problem was, that I cant use muEff * (grad U+ gradU.T()) with an error, that it is not the same dimension. So I skipped it, as you can see in code below to think about later... any help is appreciated here.

The second problem is, that it is able to compile now, and even able to run after adding it fvSolution and to 0/gradU and 0/tau, however, no output is written at all.
The files 1/gradU and 1/tau are the same asa in 0/. Can anybody help with that?

Or is there any other, even more simple possibility to get the velocity gradients in a interFoam solver?

UEqn.H added:

Code:

surfaceScalarField muEff
    (
        "muEff",
        twoPhaseProperties.muf()
      + fvc::interpolate(rho*turbulence->nut())
    );

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      - fvm::laplacian(muEff, U)
      - (fvc::grad(U) & fvc::grad(muEff))
    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
    );

    UEqn.relax();

	volTensorField gradU 
		(fvc::grad(U));


if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
            )
        );

	solve
	(
		fvm::ddt(d)
		==
		- fvm::div(phi,d)
	);
	
    }

own_interDynMFoam.C

Code:

-------------------snip/snap-------------------------
runTime.write();

	d=fvc::div(phi);
              tau=(gradU + gradU.T());

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-------------------snip/snap-------------------------

createFields.H

Code:

-------------------snip/snap-------------------------
//Creating field for gradU
Info<< "Reading field gradU\n" << endl;
    volTensorField gradU
    (
        IOobject
        (
            "gradU",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

//Creating field for tau
Info<< "Reading field tau\n" << endl;
    volTensorField tau
    (
        IOobject
        (
            "tau",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

//Creating field for divergence d
Info<< "Reading field d\n" << endl;
    volScalarField d
    (
        IOobject
        (
            "d",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
-------------------snip/snap-------------------------

owayz · October 21, 2011, 10:49

Hallo Soeren,
I want to know why do you want to calculate the gradients of velocity and write them?
I actually have an idea, I think you might also be trying to do the same. I want to do adaptive mesh refinement on the basis of velocity gradients. It would be really helpful for me if I could calculate the velocity gradients some how in each timestep I would then be able to do the adaptive mesh refinement. Tell me if you have the same idea.

Regards,
Awais Ali

wersoe · October 21, 2011, 11:44

Dear Awais ,

I manage to write the gradients since some weeks...
We are working in the field of biotech, and using CFD to optimize and/or develop bioreactors in which cells growth, which are quite sensitive to any stress. So, we would like to use the gradients to caculate shear strain and normal strain...

However, we know, that these values (shear strain, vel gradients, ...) are quite much dependent on the size of the mesh... so we should be able to compare just same cases with different working parameters, e.g. stirrer speed...

Later on, we would like to compare between different reactor sizes, which means different meshes, too. Therefore, we would need to adapt the mesh until the change in the gradients is small enough...

If you could help me with adapting the mesh, I would much appreciate...

Here the changes in the solver, I use interDyMFoam:

file interDyMFoam.C

Code:

.
.
.
#include "UEqn.H"
#include "shear.H" // added
.
.
.

]

file shear.H

Code:

volScalarField mu
    (
        "mu",
        twoPhaseProperties.mu()
    );

volTensorField gradU = fvc::grad(U);   // is not written, and needs not to be declared in createFields.H if declaration (volTensorField) is done here

//--------------------------------------------------------------------------------------------------------------------
	D = alpha1*(gradU+gradU.T());				// needs to be declared in createFields.H
	Tau = alpha1 * mu * (gradU+gradU.T());			// 
	Strain = alpha1 * pow(2,0.5) * mag(symm(fvc::grad(U))); //
	TauStrain= mu * Strain;

file ../createFields.H

Code:

.
.                   //just add as much as you defined new
.
//Creating field for Strain
Info<< "Reading field Strain\n" << endl;
    volScalarField Strain
    (
        IOobject
        (
            "Strain",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,    //is read, needs to be in folder 0
            IOobject::AUTO_WRITE   //is written
        ),
        mesh
    );
.
.
.

Hope this helps, dont hesitate to ask again...

If anybody finds this is really not good, please inform me, since this is my own hacking, and I did it the first time with OF... thanks

Best, Sören

owayz · October 22, 2011, 20:53

Thanks for your reply Soeren,
Well if you are using interDyMFoam you must have seen the dambreak tutorial. My idea is that if you want to refine on the basis of strain and as you are able to write the strain at every step. You just need to redefine the dynamicMeshdict. Instead of alpha1 you need to refine on the basis of mag or gradients or strain. If hope this helps you and best of luck to you.
But I have one problem. The solver that I have modiefied according to the information that you have provided asks me to provide a file for magnitude of gradients of U. And even if I provide a file it somehow is writing the same file every time step again.
If you can suggest me something which could help me.
Regards,
Awais Ali

Linse · October 24, 2011, 10:57

Quote:

Originally Posted by owayz

Hallo Soeren,
I want to know why do you want to calculate the gradients of velocity and write them?
I actually have an idea, I think you might also be trying to do the same. I want to do adaptive mesh refinement on the basis of velocity gradients. It would be really helpful for me if I could calculate the velocity gradients some how in each timestep I would then be able to do the adaptive mesh refinement. Tell me if you have the same idea.

Regards,
Awais Ali

Seems we are rowing the same boat here. At least this is one of the issues I plan to do in close future. So this thread is on my "watchlist" by now. ;-)
I'll keep you posted on any developments from my side...

wersoe · October 24, 2011, 11:09

Quote:

Originally Posted by owayz

Well if you are using interDyMFoam you must have seen the dambreak tutorial. My idea is that if you want to refine on the basis of strain and as you are able to write the strain at every step. You just need to redefine the dynamicMeshdict. Instead of alpha1 you need to refine on the basis of mag or gradients or strain. If hope this helps you and best of luck to you.

Thanks, will try it someday... Rite now i didnt care about refinement at all. So, I know the dambreak, but didnt have a closer look to the refinement at the surface...

Quote:

Originally Posted by owayz

But I have one problem. The solver that I have modiefied according to the information that you have provided asks me to provide a file for magnitude of gradients of U. And even if I provide a file it somehow is writing the same file every time step again.

Well, its quite hard to say anything...
First, you have to provide the file because there is a

Code:

IOobject::MUST_READ,    //is read, needs to be in folder 0

in createFields.H
Try to change it to NO_READ or something similar, I do not know the right keyword, so I just provided a file...
And here the error may be caused: gradU is a volTensorField.
So, in createFields it must be defined as volTensorField - and so it must be defined in 0/gradU:

Code:

.
.
.
dimensions [ 0 0 -1 0 0 0 0 ]
internalField uniform ( 0 0 0 0 0 0 0 0 0 )
.
.
.

This is valid for 3D, change it accordingly to 2D...

Hope this helps, Sören

owayz · November 28, 2011, 11:10

Dear Sören,
//--------------------------------------------------------------------------------------------------------------------
D = alpha1*(gradU+gradU.T()); // needs to be declared in createFields.H
Tau = alpha1 * mu * (gradU+gradU.T()); //
Strain = alpha1 * pow(2,0.5) * mag(symm(fvc::grad(U))); //
TauStrain= mu * Strain;

I have a question. In the above code what does the "symm" function do or why are you using it. Does it make any difference if I calculate the gradients without using the symm function.

Regards,
Awais

romant · November 28, 2011, 11:40

Quote:

Originally Posted by owayz

Dear Sören,
//--------------------------------------------------------------------------------------------------------------------
D = alpha1*(gradU+gradU.T()); // needs to be declared in createFields.H
Tau = alpha1 * mu * (gradU+gradU.T()); //
Strain = alpha1 * pow(2,0.5) * mag(symm(fvc::grad(U))); //
TauStrain= mu * Strain;

I have a question. In the above code what does the "symm" function do or why are you using it. Does it make any difference if I calculate the gradients without using the symm function.

Regards,
Awais

"symm" gives you the symmetric tensor of fvc::grad(U) since it is symmetric. I think it might be reduce computational and memory usage.

JFM · September 14, 2013, 01:34

Good day everyone

I am using OpenFOAM to assess a hydraulic structure which is designed to operate at critical condition (Fr=1). Previously I had assessed this structure using two commercial 2D FVM model, however the results were not acceptable.

As a result I have created a 3D interFoam model of the structure and are aiming to prove that OpenFOAM is a suitable tool for designers. I need to determine the influence, if any, of the cross momentum terms from the convective velocity field, div (phi, U). The terms I am particularly interested in are:
x-direction: u.du/dx + v.du/dy + w.du/dz
y-direction: u.dv/dx + v.dv/dy + w.dv/dz
z-direction: u.dw/dx + v.dw/dy + w.dw/dz
(at this stage I am not interested in the RAS or LES stresses)

I have not found anyway (in the forum or other references) to write this tensor out from the model. However, it appears that there is a way to extract these components from the model and your post provides some directions on how to do this. Could someone provide some direction on how this could be achieved in my case, please include file / dictionary names that may need to be amended.

I have been using OpenFOAM for around a year but still have limited confidence when it comes to manipulation of source files therefore any assistance will be greatly appreciated - hopefully I can return the favour one day.

Regards
JFM

wersoe · September 14, 2013, 09:15

You should use something like this in a seperate file:

-----------grad.H--------------------------

volScalarField X=U.component(0);
volScalarField Y=U.component(1);
volScalarField Z=U.component(2);

volScalarField XX=gradU.component(0);
volScalarField YY=gradU.component(4);
volScalarField ZZ=gradU.component(8);

volScalarField XY=gradU.component(1);
volScalarField YZ=gradU.component(5);
volScalarField ZX=gradU.component(6);

volScalarField YX=gradU.component(3);
volScalarField ZY=gradU.component(7);
volScalarField XZ=gradU.component(2);

Variable = X*XX + ...

Include it in the solver, dont forget to create the fields in createFields.H...
Just see my previous post...

Then you will be able to write it during the calculation.
It would be much easier and more comfortable if you calculate this solution after the simulation is done...
There are a couple of tools which could be the basis, need to be changed and recompiled too...

Best, Sören

huangxianbei · March 31, 2014, 21:57

Hi:
I just use this code

Code:

volTensorField graU = fvc::grad(UMean);
    volScalarField graUy(graU.component(1));

in calculateFields.H to calculate dU/dy in postChannel utility, and in the collapse.H

Code:

scalarField graUyValues(channelIndexing.collapse(graUy));
makeGraph(y, graUyValues, "graUy", path, gFormat);

while when I check the results, I found the results are wrong for the value is of e-23 order, I'm wondering if this can be used in a post-process?

j-avdeev · December 19, 2015, 05:15

Quote:

Originally Posted by owayz

//--------------------------------------------------------------------------------------------------------------------
D = alpha1*(gradU+gradU.T());

Hello.
Can someone explain what gradU.T() means?
Why here gradU+gradU.T(), but not just gradU?
As I understand we want to get gradient field of U here.

romant · January 6, 2016, 09:17

Quote:

Originally Posted by j-avdeev

Hello.
Can someone explain what gradU.T() means?
Why here gradU+gradU.T(), but not just gradU?
As I understand we want to get gradient field of U here.

The .T() operator gives you the transpose of the matrix http://foam.sourceforge.net/docs/cpp...61ae7a972b4f96

Mowgli · June 24, 2016, 10:54

I use the same code to calculate velocity gradients of the field in my case, but I get 0 gradients, which seems like the 0 file for gradU is being duplicated. Even with a NO_READ condition, I get 0 gradients. ANy idea as to why fvc::grad(U) outputs 0 gradient ?

October 21, 2011, 11:44		#3
wersoe Member Soeren Werner Join Date: Mar 2009 Location: Wädenswil, Switzerland Posts: 32 Rep Power: 17	Dear Awais , I manage to write the gradients since some weeks... We are working in the field of biotech, and using CFD to optimize and/or develop bioreactors in which cells growth, which are quite sensitive to any stress. So, we would like to use the gradients to caculate shear strain and normal strain... However, we know, that these values (shear strain, vel gradients, ...) are quite much dependent on the size of the mesh... so we should be able to compare just same cases with different working parameters, e.g. stirrer speed... Later on, we would like to compare between different reactor sizes, which means different meshes, too. Therefore, we would need to adapt the mesh until the change in the gradients is small enough... If you could help me with adapting the mesh, I would much appreciate... Here the changes in the solver, I use interDyMFoam: file interDyMFoam.C Code: . . . #include "UEqn.H" #include "shear.H" // added . . . ] file shear.H Code: volScalarField mu ( "mu", twoPhaseProperties.mu() ); volTensorField gradU = fvc::grad(U); // is not written, and needs not to be declared in createFields.H if declaration (volTensorField) is done here //-------------------------------------------------------------------------------------------------------------------- D = alpha1(gradU+gradU.T()); // needs to be declared in createFields.H Tau = alpha1 mu * (gradU+gradU.T()); // Strain = alpha1 * pow(2,0.5) * mag(symm(fvc::grad(U))); // TauStrain= mu * Strain; file ../createFields.H Code: . . //just add as much as you defined new . //Creating field for Strain Info<< "Reading field Strain\n" << endl; volScalarField Strain ( IOobject ( "Strain", runTime.timeName(), mesh, IOobject::MUST_READ, //is read, needs to be in folder 0 IOobject::AUTO_WRITE //is written ), mesh ); . . . Hope this helps, dont hesitate to ask again... If anybody finds this is really not good, please inform me, since this is my own hacking, and I did it the first time with OF... thanks Best, Sören arvindpj, meth and wbywbywby6 like this.

September 14, 2013, 01:34	Extract div(phi, U) tensor	#9
JFM Member John Anonymous Join Date: Jan 2011 Location: Melbourne Australia Posts: 37 Rep Power: 15	Good day everyone I am using OpenFOAM to assess a hydraulic structure which is designed to operate at critical condition (Fr=1). Previously I had assessed this structure using two commercial 2D FVM model, however the results were not acceptable. As a result I have created a 3D interFoam model of the structure and are aiming to prove that OpenFOAM is a suitable tool for designers. I need to determine the influence, if any, of the cross momentum terms from the convective velocity field, div (phi, U). The terms I am particularly interested in are: x-direction: u.du/dx + v.du/dy + w.du/dz y-direction: u.dv/dx + v.dv/dy + w.dv/dz z-direction: u.dw/dx + v.dw/dy + w.dw/dz (at this stage I am not interested in the RAS or LES stresses) I have not found anyway (in the forum or other references) to write this tensor out from the model. However, it appears that there is a way to extract these components from the model and your post provides some directions on how to do this. Could someone provide some direction on how this could be achieved in my case, please include file / dictionary names that may need to be amended. I have been using OpenFOAM for around a year but still have limited confidence when it comes to manipulation of source files therefore any assistance will be greatly appreciated - hopefully I can return the favour one day. Regards JFM

September 14, 2013, 09:15		#10
wersoe Member Soeren Werner Join Date: Mar 2009 Location: Wädenswil, Switzerland Posts: 32 Rep Power: 17	You should use something like this in a seperate file: -----------grad.H-------------------------- volScalarField X=U.component(0); volScalarField Y=U.component(1); volScalarField Z=U.component(2); volScalarField XX=gradU.component(0); volScalarField YY=gradU.component(4); volScalarField ZZ=gradU.component(8); volScalarField XY=gradU.component(1); volScalarField YZ=gradU.component(5); volScalarField ZX=gradU.component(6); volScalarField YX=gradU.component(3); volScalarField ZY=gradU.component(7); volScalarField XZ=gradU.component(2); Variable = X*XX + ... Include it in the solver, dont forget to create the fields in createFields.H... Just see my previous post... Then you will be able to write it during the calculation. It would be much easier and more comfortable if you calculate this solution after the simulation is done... There are a couple of tools which could be the basis, need to be changed and recompiled too... Best, Sören owayz, RuiVO, Markus Kolano and 3 others like this.

March 31, 2014, 21:57		#11
huangxianbei Senior Member Huang Xianbei Join Date: Sep 2013 Location: Yangzhou,China Posts: 302 Rep Power: 14	Hi: I just use this code Code: volTensorField graU = fvc::grad(UMean); volScalarField graUy(graU.component(1)); in calculateFields.H to calculate dU/dy in postChannel utility, and in the collapse.H Code: scalarField graUyValues(channelIndexing.collapse(graUy)); makeGraph(y, graUyValues, "graUy", path, gFormat); while when I check the results, I found the results are wrong for the value is of e-23 order, I'm wondering if this can be used in a post-process?

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Problem with assigned inlet velocity profile as a boundary condition	Ozgur_	FLUENT	5	August 25, 2015 05:58
Emergency:UDF for a time dependent parabolic velocity	zumaqiong	Fluent UDF and Scheme Programming	12	March 25, 2010 13:00
Velocity gradients from previous timesteps	Ed Kay	FLUENT	0	March 4, 2010 09:04
How to output unsteady velocity profiles?	Jason	FLUENT	3	August 8, 2007 10:40
Velocity Gradients at Inviscid Flow	Lada	FLUENT	0	January 15, 2001 08:20

October 21, 2011, 10:49		#2
owayz Senior Member Awais Ali Join Date: Feb 2010 Location: Germany Posts: 128 Rep Power: 17	Hallo Soeren, I want to know why do you want to calculate the gradients of velocity and write them? I actually have an idea, I think you might also be trying to do the same. I want to do adaptive mesh refinement on the basis of velocity gradients. It would be really helpful for me if I could calculate the velocity gradients some how in each timestep I would then be able to do the adaptive mesh refinement. Tell me if you have the same idea. Regards, Awais Ali

October 22, 2011, 20:53		#4
owayz Senior Member Awais Ali Join Date: Feb 2010 Location: Germany Posts: 128 Rep Power: 17	Thanks for your reply Soeren, Well if you are using interDyMFoam you must have seen the dambreak tutorial. My idea is that if you want to refine on the basis of strain and as you are able to write the strain at every step. You just need to redefine the dynamicMeshdict. Instead of alpha1 you need to refine on the basis of mag or gradients or strain. If hope this helps you and best of luck to you. But I have one problem. The solver that I have modiefied according to the information that you have provided asks me to provide a file for magnitude of gradients of U. And even if I provide a file it somehow is writing the same file every time step again. If you can suggest me something which could help me. Regards, Awais Ali

November 28, 2011, 11:10		#7
owayz Senior Member Awais Ali Join Date: Feb 2010 Location: Germany Posts: 128 Rep Power: 17	Dear Sören, //-------------------------------------------------------------------------------------------------------------------- D = alpha1(gradU+gradU.T()); // needs to be declared in createFields.H Tau = alpha1 mu * (gradU+gradU.T()); // Strain = alpha1 * pow(2,0.5) * mag(symm(fvc::grad(U))); // TauStrain= mu * Strain; I have a question. In the above code what does the "symm" function do or why are you using it. Does it make any difference if I calculate the gradients without using the symm function. Regards, Awais

June 24, 2016, 10:54		#14
Mowgli New Member gnompher Join Date: Jun 2016 Posts: 6 Rep Power: 10	I use the same code to calculate velocity gradients of the field in my case, but I get 0 gradients, which seems like the 0 file for gradU is being duplicated. Even with a NO_READ condition, I get 0 gradients. ANy idea as to why fvc::grad(U) outputs 0 gradient ?