CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)
-   -   Splitting Momentum Eqns and Adding Extra Terms to Each (https://www.cfd-online.com/Forums/openfoam/217188-splitting-momentum-eqns-adding-extra-terms-each.html)

 ErenC May 2, 2019 14:25

Splitting Momentum Eqns and Adding Extra Terms to Each

Greetings again!

I am trying to validate a case, but my implementation doesn't work properly so I wanna try something different, well, actually the original. I made my research but couldn't find answers to my questions.

1) Is it possible to add different terms to x and y momentum equation? I am thinking of call them UxEqn and UyEqn and make two loops for UEqn in total. Thats how I did in CFD class with MATLAB in past.

2) How can I add those terms? Can I just type:

- (BT)*(Uy*sine(y)*cosine(y)-Ux*sine^2(y))

BT=thonf*B^2/rhonf

I don't think so, doesn't look right.

I am adding corresponding momentum equations, units are matching no problem there.(note: all extra terms are scalar excluding u-v)

Any help would be great!

https://i.ibb.co/RgTJ02L/mhd.jpg

 alexeym May 3, 2019 11:59

Hi,

1. You can always turn two scalars into vector multiplying them by corresponding unit vectors. I.e. vector = scalarx*i + scalary*j + scalarz*k. So your terms with sines and cosines could be rewritten as single term.

2. What is gamma?

3. The first term in the second expression is just buoyancy term. Select appropriate solver.

 ErenC May 3, 2019 12:11

Quote:
 Originally Posted by alexeym (Post 732671) Hi, 1. You can always turn two scalars into vector multiplying them by corresponding unit vectors. I.e. vector = scalarx*i + scalary*j + scalarz*k. So your terms with sines and cosines could be rewritten as single term. 2. What is gamma? 3. The first term in the second expression is just buoyancy term. Select appropriate solver.

1) How? I clearly can't think a way since I am newbie to OpenFOAM.

2) Its angle between Bx By vectors(magnetic field)

3) I did, I modified and re-compiled buoyantBoussinesqPimpleFoam and its working but It solves maxwell equations. I want to do exactly like this and my problem is with sine-cosine terms.

 alexeym May 3, 2019 16:34

a. Bx and By are usually scalars (since they are components of B vector). What is gamma?

b. Is gamma a constant or a field? I.e. does it depend on coordinate?

c. Is B0 a constant or a field?

Your sine/cosine terms can be rewritten, for example, as: where U is your usual velocity vector, and second term is constructed with something like

Code:

sqr(sin(gamma))*vector(1, 0, 0) + sin(gamma)*cos(gamma)*vector(0, 1, 0)
Though maybe it is in fact a product of velocity vector and some kind of rotation matrix and can be coded as vector-matrix multiplication.

 ErenC May 4, 2019 12:09

Sorry about late reply, I got forbidden error, after that I got another error. I had to switch browser.

Sorry for not beign clear.
gamma= tan-1(By/Bx)
B0= sqrt(By^2+Bx^2)

When I am solving maxwell equations, B is a volVectorField. But in that(like in the picture) NS form I can identify them as a transportproperty. At least I think that will work, I might be wrong.

About the 2nd equation, but that will be:

sqr(cos(gamma))*vector(0, 1, 0) + sin(gamma)*cos(gamma)*vector(1, 0, 0)

Am I right? but when adding those terms how do I implement them into:

 alexeym May 6, 2019 04:52

Hi,

You have got the idea correctly, though, I lost minus sign before sqr in my code example. So, the first term is:

Code:

T1 = -sqr(sin(gamma))*vector(1, 0, 0) + sin(gamma)*cos(gamma)*vector(0, 1, 0)
the second term will be:

Code:

T2 = -sqr(cos(gamma))*vector(0, 1, 0) + sin(gamma)*cos(gamma)*vector(1, 0, 0)
and resulting source term in momentum equation is:

Code:

U & (T1 + T2)
It is not quite clear, why you need to put these terms into fvm::grad. Also there is no fvm::grad (which constructs FV matrix), only fvc::grad (which constructs vector field).

 ErenC May 6, 2019 09:07

Thank you sir!

But as I said I am totally newbie, So I need just a little more assistance.

So, I created volVectorFields V1 and V2, created corresponding variables and Ueqn as following:

Code:
Code:

 // Momentum Eqns   V1 = (- sqr(sin(gamma))*U(1, 0, 0) + sin(gamma)*cos(gamma)*U(0, 1, 0) );      V2 = (- sqr(cos(gamma))*U(0, 1, 0) + sin(gamma)*cos(gamma)*U(1, 0, 0) );     MRF.correctBoundaryVelocity(U);     fvVectorMatrix UEqn     (         fvm::ddt(U) + fvm::div(phi, U) - dvc::grad(phi, (V1+V2))       + MRF.DDt(U)       + turbulence->divDevReff(U)       ==         fvOptions(U)     );     UEqn.relax();     fvOptions.constrain(UEqn);

And I am getting these errors:

Code:

 UEqn.H:4:37: error: no match for call to ‘(Foam::volVectorField {aka Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>}) (int, int, int)’     V2 = (- sqr(cos(gamma))*U(0, 1, 0) + sin(gamma)*cos(gamma)*U(1, 0, 0) );
Code:

UEqn.H:4:37: error: no match for call to ‘(Foam::volVectorField {aka Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>}) (int, int, int)’     V2 = (- sqr(cos(gamma))*U(0, 1, 0) + sin(gamma)*cos(gamma)*U(1, 0, 0) );

Code:

 UEqn.H:11:9: error: ‘dvc’ has not been declared         - dvc::grad(U, (V1+V2))

It says no match for call to but I have these lines inside createFields,H
Code:

Info<< "Reading field V1\n" << endl; volVectorField V1 (     IOobject     (         "V1",         runTime.timeName(),         mesh,         IOobject::MUST_READ,         IOobject::AUTO_WRITE     ),     mesh ); Info<< "Reading field V2\n" << endl; volVectorField V2 (     IOobject     (         "V2",         runTime.timeName(),         mesh,         IOobject::MUST_READ,         IOobject::AUTO_WRITE     ),     mesh );
And no Idea about dvc::grad because I am seeing this first time,

 alexeym May 6, 2019 11:12

So why you have decided you write U(1, 0, 0)?

dvc was (obviously) my typo. Yet I can not understand why do you need to put those expressions under grad operator. There are no grads in your original equations.

 ErenC May 6, 2019 13:11

Because of my ignorance. I though you were saying corresponding vector by that.
Code:

// Solve the momentum equation   V1 = (- sqr(sin(gamma))*vector(1, 0, 0)+sin(gamma)*cos(gamma)*vector(0, 1, 0) );      V2 = (- sqr(cos(gamma))*vector(0, 1, 0)+sin(gamma)*cos(gamma)*vector(1, 0, 0) );     MRF.correctBoundaryVelocity(U);     fvVectorMatrix UEqn     (       fvm::ddt(U) + fvm::div(phi, U)       -coef*V1       -coef*V2       + MRF.DDt(U)       + turbulence->divDevReff(U)
and its compiles. Sorry about my absurd comments, I am doing many things, I just got confused.

And I really havent understood the U&(V1+V2) Term. When I implement it, it doesn't work. What does "&" stands for? And in this form(my code) it is obviously wrong, I have to give BC's to V1 and V2 which will be absurd, vectors corresponding to nothing(or do they because they are in Ueqn?).

when I try
Code:

 fvVectorMatrix UEqn     (       fvm::ddt(U) + fvm::div(phi, U)       -U & (V1+V2)       + MRF.DDt(U)       + turbulence->divDevReff(U)
I am getting this error:

Code:

UEqn.H:11:10: error: no match for ‘operator&’ (operand types are ‘Foam::tmp<Foam::fvMatrix<Foam::Vector<double> > >’ and ‘Foam::tmp<Foam::fvMatrix<Foam::Vector<double> > >’)

 alexeym May 6, 2019 15:58

Well, yes, it is a little bit more complicated. & is overloaded operator and it means scalar product of two vectors. So to write source terms from your first post it is necessary to have something like:

Code:

T1 = -sqr(sin(gamma))*vector(1, 0, 0) + sin(gamma)*cos(gamma)*vector(0, 1, 0) Su = U & T1 T2 = -sqr(cos(gamma))*vector(0, 1, 0) + sin(gamma)*cos(gamma)*vector(1, 0, 0) Sv = U & T2 SU = Su * vector(1, 0, 0) + Sv*vector(0, 1, 0)
where SU is source term in momentum equation.

But if you are trying to add Lorentz force and, in fact, your source term is You can use either mhdFoam, where this term is already in equations, and equations for magnetic field are also solved. Or you can write this term directly using OpenFOAM(R) operator notation. Now it is a little bit clearer why you wanted to put something into grad operator.

 ErenC May 6, 2019 16:50

Well, yes. I actually implemented Maxwell equations into buoyantBoussinesqPimpleFoam, it is working nicely and it really wasn't diffucult like this(at least for me). But in paper they are using this formulation so I wanted to try this one.

So, I tried that formulation and I am getting
error: no match for 'operator='

Sv = U & V2

This error. I think this error is happens when left hand side is not same type as right hand side.

I defined T1, T2, Su, Sv, SU as volVectorField(which will be problem because it will ask for boundary conditions, but that is another days problem I think...) What am I missing?

 alexeym May 7, 2019 05:51

You do not need to declare them as MUST_READ / AUTO_WRITE, since they are calculated fields, so you do not read them. And these fields are for calculation of volumetric force, so BCs are of calculated type.

T1 is vector field, U is vector field, so U & T1 (scalar product) is scalar field.

Assignment is a bit of pseudocode. In real solver, you should use constructor a la:

Code:

const vector iHat(1, 0, 0); const vector jHat(0, 1, 0); const volVectorField T1(-sqr(sin(gamma))*iHat + sin(gamma)*cos(gamma)*jHat) const volScalarField Su(U & T1) const volVectorField T2(-sqr(cos(gamma))*jHat + sin(gamma)*cos(gamma)*iHat) const volScalarField Sv(U & T2) const volVectorField SU(Su*iHat + Sv*jHat)
Also during simulation you have to respect your choice of x, y (direction of g) and z (empty direction).

 ErenC June 3, 2019 18:55

Thank you so much for your patience and answers Dr. Matheichev. I prematurely tried to include the Lorenz force, so I went back to study on my hydrodynamics.

And I found this http://www.tfd.chalmers.se/~hani/kur...desTassone.pdf
Anyone willing to add lorenz force can follow up here. And I must say even student projects are really good in this University. I'll probably apply to post doc after my PhD, I am amazed by their thesis studies.

 Originally Posted by alexeym (Post 732930) You do not need to declare them as MUST_READ / AUTO_WRITE, since they are calculated fields, so you do not read them. And these fields are for calculation of volumetric force, so BCs are of calculated type. T1 is vector field, U is vector field, so U & T1 (scalar product) is scalar field. Assignment is a bit of pseudocode. In real solver, you should use constructor a la: Code: const vector iHat(1, 0, 0); const vector jHat(0, 1, 0); const volVectorField T1(-sqr(sin(gamma))*iHat + sin(gamma)*cos(gamma)*jHat) const volScalarField Su(U & T1) const volVectorField T2(-sqr(cos(gamma))*jHat + sin(gamma)*cos(gamma)*iHat) const volScalarField Sv(U & T2) const volVectorField SU(Su*iHat + Sv*jHat) Also during simulation you have to respect your choice of x, y (direction of g) and z (empty direction).