Splitting Momentum Eqns and Adding Extra Terms to Each

ErenC · May 2, 2019, 14:25

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!

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

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.

Thank you for your answer.

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

OK. Let's continue asking questions:

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:

$v\sin\left(\gamma\right)\cos\left(\gamma\right) - u\sin^2\left(\gamma\right) = \vec{U}\cdot \left\| \begin{array}{c} -\sin^2\left(\gamma\right) \\ \sin\left(\gamma\right)\cos\left(\gamma\right) \\ 0\end{array}\right\|$

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:
-fvm::grad( ) ?

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

$\frac{\left(\vec{B}\cdot{}\nabla\right)\vec{B}}{\mu_0} - \nabla\left(\frac{B^2}{2\mu_0}\right)$

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.

sadsid · May 22, 2020, 23:30

Quote:

Originally Posted by alexeym

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).

Is it possible to use fvOptions for the same source term?

May 2, 2019, 14:25	Splitting Momentum Eqns and Adding Extra Terms to Each	#1
ErenC Member Eren Join Date: Aug 2018 Posts: 86 Rep Power: 8	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)(Uysine(y)cosine(y)-Uxsine^2(y)) BT=thonfB^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! Last edited by ErenC; May 2, 2019 at 16:59.*

May 3, 2019, 16:34		#4
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 38	OK. Let's continue asking questions: 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: $v\sin\left(\gamma\right)\cos\left(\gamma\right) - u\sin^2\left(\gamma\right) = \vec{U}\cdot \left\\| \begin{array}{c} -\sin^2\left(\gamma\right) \\ \sin\left(\gamma\right)\cos\left(\gamma\right) \\ 0\end{array}\right\\|$ 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.

May 6, 2019, 04:52		#6
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 38	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). Last edited by alexeym; May 6, 2019 at 11:09. Reason: Typo

May 6, 2019, 13:11		#9
ErenC Member Eren Join Date: Aug 2018 Posts: 86 Rep Power: 8	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) -coefV1 -coefV2 + 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> > >’) Last edited by ErenC; May 6, 2019 at 15:39.

May 6, 2019, 15:58		#10
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 38	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 $\frac{\left(\vec{B}\cdot{}\nabla\right)\vec{B}}{\mu_0} - \nabla\left(\frac{B^2}{2\mu_0}\right)$ 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.

May 3, 2019, 11:59		#2
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 38	Hi, 1. You can always turn two scalars into vector multiplying them by corresponding unit vectors. I.e. vector = scalarxi + scalaryj + 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.

May 4, 2019, 12:09		#5
ErenC Member Eren Join Date: Aug 2018 Posts: 86 Rep Power: 8	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: -fvm::grad( ) ?

May 6, 2019, 11:12		#8
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 38	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.

May 6, 2019, 16:50		#11
ErenC Member Eren Join Date: Aug 2018 Posts: 86 Rep Power: 8	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?

May 7, 2019, 05:51		#12
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 38	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(SuiHat + SvjHat) Also during simulation you have to respect your choice of x, y (direction of g) and z (empty direction).

June 3, 2019, 18:55		#13
ErenC Member Eren Join Date: Aug 2018 Posts: 86 Rep Power: 8	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.