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

Splitting Momentum Eqns and Adding Extra Terms to Each

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 2, 2019, 14:25
Default Splitting Momentum Eqns and Adding Extra Terms to Each
  #1
Member
 
Eren
Join Date: Aug 2018
Posts: 41
Rep Power: 4
ErenC is on a distinguished road
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!



Last edited by ErenC; May 2, 2019 at 16:59.
ErenC is offline   Reply With Quote

Old   May 3, 2019, 11:59
Default
  #2
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,922
Rep Power: 34
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to 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.
alexeym is offline   Reply With Quote

Old   May 3, 2019, 12:11
Default
  #3
Member
 
Eren
Join Date: Aug 2018
Posts: 41
Rep Power: 4
ErenC is on a distinguished road
Quote:
Originally Posted by alexeym View Post
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.
ErenC is offline   Reply With Quote

Old   May 3, 2019, 16:34
Default
  #4
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,922
Rep Power: 34
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
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.
alexeym is offline   Reply With Quote

Old   May 4, 2019, 12:09
Default
  #5
Member
 
Eren
Join Date: Aug 2018
Posts: 41
Rep Power: 4
ErenC is on a distinguished road
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( ) ?
ErenC is offline   Reply With Quote

Old   May 6, 2019, 04:52
Default
  #6
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,922
Rep Power: 34
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
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
alexeym is offline   Reply With Quote

Old   May 6, 2019, 09:07
Default
  #7
Member
 
Eren
Join Date: Aug 2018
Posts: 41
Rep Power: 4
ErenC is on a distinguished road
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,
ErenC is offline   Reply With Quote

Old   May 6, 2019, 11:12
Default
  #8
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,922
Rep Power: 34
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
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.
alexeym is offline   Reply With Quote

Old   May 6, 2019, 13:11
Default
  #9
Member
 
Eren
Join Date: Aug 2018
Posts: 41
Rep Power: 4
ErenC is on a distinguished road
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> > >’)

Last edited by ErenC; May 6, 2019 at 15:39.
ErenC is offline   Reply With Quote

Old   May 6, 2019, 15:58
Default
  #10
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,922
Rep Power: 34
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
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.
alexeym is offline   Reply With Quote

Old   May 6, 2019, 16:50
Default
  #11
Member
 
Eren
Join Date: Aug 2018
Posts: 41
Rep Power: 4
ErenC is on a distinguished road
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?
ErenC is offline   Reply With Quote

Old   May 7, 2019, 05:51
Default
  #12
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,922
Rep Power: 34
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to 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).
alexeym is offline   Reply With Quote

Old   June 3, 2019, 18:55
Default
  #13
Member
 
Eren
Join Date: Aug 2018
Posts: 41
Rep Power: 4
ErenC is on a distinguished road
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.
ErenC is offline   Reply With Quote

Old   May 22, 2020, 23:30
Default
  #14
New Member
 
Join Date: Jan 2017
Posts: 20
Rep Power: 5
sadsid is on a distinguished road
Quote:
Originally Posted by alexeym View Post
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?
sadsid 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



All times are GMT -4. The time now is 08:53.