[CharIie]

I am trying to create a solver for ferrofluids under the influence of a static magnetic field from buoyantPimpleFoam. I want to add another source term, the Kelvin Body force, to the momentum equation. In my case the Kelvin Body Force defines as:
with
If the current temperature of the fluid is closer to its curie temperature the KBF gets smaller. In green you see where I want to add the term.
I managed to add H and T_Curie in createFields.H.
But my question is how do I depict the current temperature in the momentum equation to calculate the Magnetization M each time step?

It does not work like this:

tmp<fvVectorMatrix> tUEqn
    (
        fvm::ddt(rho, U) + fvm::div(phi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevTau(U)
     + Foam::constant::electromagnetic::mu0*grad(H)*H*(Tc - T);
     ==
        fvOptions(rho, U)

I have read How to add a body force to a momentum equation but in my case my source changes over time according to the temperature.

[SHUBHAM9595]

Hi CharIie,

You are almost there. Just two comments that might help you

1) The expression for Magnetization

M = H (T - T)

seems incorrect...

As you already know both M & H have dimensions of Amp/meter....However with the above expression the RHS will have dimesnsions of A.Kelvin/meter



2) Regarding the implementation, the only thing you need to worry about is the product of with M (or eventually H).....as

is a tensor of rank 2 (9 components)

whereas M or H is a tensor of rank 1 (3 components )

and the UEqn in which u wish to enter this source term is also a tensor of rank 1 (3 components)....Hence, the product of should be an inner product....which is implemented in the FOAM's love language as shown below

volTensorField gradHH = fvc::grad(H);
volVectorField KBF = Foam::constant::electromagnetic::mu0*(M & gradHH);

[CharIie]

Thank you SHUBHAM for the answer.
The solver finally compiles when I use the inner product as you showed.
(For the moment I replaced M with H.)

But now I am trying to depict M as:
This cancels out the Kelvin units and M_0 is the initial magnetization.

volVectorField M = Mo*(Tc/T);

This also does not work. I think I am messing around with Vectors, tensors and scalars again.


This is the error message I get when compiling:

UEqn.H: In function ‘int main(int, char**)’:
UEqn.H:13:25: error: no match for ‘operator/’ (operand types are ‘Foam::dimensioned<double>’ and ‘<unresolved overloaded function type>’)
   13 | volVectorField M = Mo*(Tc/T);
      |                     ~~~~~^~

[SHUBHAM9595]

Indeed it is because of the different declarations of ....

I assume u r providing both and as constant parameters.....as shown below

     // Curie Temp
    dimensionedScalar    Tc(transportProperties.lookup("Tc"));
      //Equilibrium Magnetization
    dimensionedScalar    Mo(transportProperties.lookup("Mo"));

So, now if u see....as per ur code

volVectorField M = Mo*(Tc/T);

on LHS is volVectorField however on RHS all of them are either dimensionedScalar () or volScalarField ()

One of the many workarounds is you create volVectorField using the dimensionedScalar

volVectorField Mnot
( 
	IOobject 
      ( 
	       "Mnot", 
	        runTime.timeName(), 
		mesh, 
		IOobject::NO_READ, 
		IOobject::AUTO_WRITE 
	), 
	mesh, 
dimensionedVector("vU", dimensionSet(0,-1,0,0,0,1,0),vector(Mo.value(), Mo.value(), Mo.value()))
);

Finally,

volVectorField M = Mnot*(Tc/T);

should do the job...

[CharIie]

Thank you again!
I have added Tc and Mo in createFileds.H like this: ( i will later a file with the magnetic properties to constant folder)

 
Info<< "Reading magnetic properties\n" << endl;
IOdictionary magneticProperties
(
    IOobject
    (
        "magneticProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

dimensionedScalar Mo ( "Mo",    magneticProperties.lookup("Mo"));
dimensionedScalar Tc ( "Tc",    magneticProperties.lookup("Tc"));

volVectorField Mnot
( 
    IOobject 
      ( 
           "Mnot", 
            runTime.timeName(), 
        mesh, 
        IOobject::NO_READ, 
        IOobject::AUTO_WRITE 
    ), 
    mesh, 
dimensionedVector("Mnot", dimensionSet(0,-1,0,0,0,1,0),vector(Mo.value(), Mo.value(), Mo.value()))
 );

(I have changed "vU" to "Mnot" from your example.)
and my UEqn.H looks like this now:

// Curie Temp
    dimensionedScalar    Tc(magneticProperties.lookup("Tc"));
//Equilibrium Magnetization
    dimensionedScalar    Mo(magneticProperties.lookup("Mo"));


volTensorField gradHH = fvc::grad(H);
volVectorField M = Mnot*(Tc/T);
volVectorField KBF = Foam::constant::electromagnetic::mu0*(M & gradHH);


    tmp<fvVectorMatrix> tUEqn
    (
        fvm::ddt(rho, U) + fvm::div(phi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevTau(U)
     + KBF
     
     ==
        fvOptions(rho, U)
        
     );

But still I get this error message:

 In file included from fhdFoam.C:157:
UEqn.H: In function ‘int main(int, char**)’:
UEqn.H:17:28: error: no match for ‘operator/’ (operand types are ‘Foam::dimensionedScalar’ {aka ‘Foam::dimensioned<double>’} and ‘<unresolved overloaded function type>’)
       17 | volVectorField M = Mnot*(Tc/T);
          |                          ~~^~

I also tried to add Tnot as volScalarField in similar way to Mnot, but it also results in the same error. A vector multiplied with a scalar should stay a vector. I don't understand where this error is still coming from.

[SHUBHAM9595]

are you sure you didn't forgot to declare the Temperature field ?

UEqn.H:17:28: error: no match for ‘operator/’ .....and ’<unresolved overloaded function type>’)

[CharIie]

What do you mean by declaring the temperature field? Isn't the temperature already declared onboard with buoyantPimpleFoam?

[SHUBHAM9595]

buoyantPimpleFoam solve either enthalpy based or internal energy based equation (EEeqn) instead of directly solving for Temperature....(https://cfd.direct/openfoam/energy-equation/#x1-70004)

You can easily check out the contents of createFields.H....(https://develop.openfoam.com/Develop...createFields.H)......

there does not exist any field variable "T"....

[CharIie]

Ok, this brings me back to my initial question in this thread:

Originally Posted by CharIie

But my question is how do I depict the current temperature in the momentum equation to calculate the Magnetization M each time step?

If the temperture is not directly solved, but via enthalpy based or internal energy based equation in EEeqn; how do I make use of the temperature, since a physical value for T is shown to me in Paraview?

How can I depict the current temperature in the momentum equation for simple calculations?

[SHUBHAM9595]

There are two approaches....assuming that you want to stick with buoyancy driven flow solvers....

1) You can switch to buoyantBoussinesqPimpleFoam which inherently solves Temperature based energy equation.

2) In the current buoyantPimpleFoam version..... comment the EEqn.H in buoyantPimpleFoam.C.....and instead add TEqn.H....more details can be found here https://openfoamwiki.net/index.php/H...ure_to_icoFoam

[CharIie]

I decided to go with approach 2 since buoyantBoussinesqPimpleFoam isnot part of my OpenFOAM V8 .

I added the TEqn.H according to https://openfoamwiki.net/index.php/H...ure_to_icoFoam, I did not comment out old the EEqn.H, since it seems to do no harm to my solver.

My solver also does compile now, but when I try to solve my test case I get the following error:

--> FOAM FATAL ERROR: 
incompatible dimensions for operation 
    [T[0 0 -1 1 0 0 0] ] + [T[1 -3 -1 1 0 0 0] ]

This error comes from the new temperature calculation in the EEqn. buoyantPimpleFoam is compressible, icoFoam incompressible. What changes could I make to get rid of the dimensions error?

[SHUBHAM9595]

Hi CharIie,

As you mentioned the issue seems to come from extra density....kg/m3 in one of your equations.....

--> FOAM FATAL ERROR: 
incompatible dimensions for operation 
    [T[0 0 -1 1 0 0 0] ] + [T[1 -3 -1 1 0 0 0] ]

Now, if we follow the link https://openfoamwiki.net/index.php/H...ure_to_icoFoam about adding Temp...…we don't have any term that have dimensions of kg/m3....(T = K, phi = m3/s, DT = m2/s)

So there's high probability that the error might be coming from the EEqn.H.....Have u tried any simulation without EEqn.H ?


regards,
Shubham

[CharIie]

Hi I have added TEqn like this:

            #include "UEqn.H"
            #include "EEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())
            {
                turbulence->correct();
                thermophysicalTransport->correct();
            }
        }

        rho = thermo.rho();
       //add these lines...
Info<< "Temperature calculation\n" << endl;
        fvScalarMatrix TEqn
        (
            fvm::ddt(rho,T) //added rho
            + fvm::div(phi, T)
            //- fvm::laplacian(DT, T) 
        );

        TEqn.solve();
//done adding lines...
   
        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

     return 0;

In accordance with How to add temperature to cavitatingFoam solver I added rho in the temperature calculation. It works now without or together with the energy equation.

I finally get results with my solver, but the courant number is way to high:

Courant Number mean: 49.3032 max: 353.349

I need help with including the the molecular transport (- fvm::laplacian(DT, T)) or does someone have a better idea on how to better address the current temperature?

[SHUBHAM9595]

Glad to know that u r not getting any error. However, you should not solve the same governing equation TWICE.....as in ur modified version u have both EEqn.H & contents of TEqn.H........

Out of curiosity, can u specify exact steps and command u r follwoing after making the changes in .C file.....I'm afraid we might be missing out some intermediate steps here

[CharIie]

Okay, here we go:
My new solver is a copy of buyantPimpleFoam saved in /OpenFOAM-8/applications/solvers/heatTransfer/fhdFoam.
So I far I have made changes to fhdFoam.C, UEqn, createFields and the make files of course.

I changed the column including the EEqn like this:
// #include "EEqn.H"
Then I type wclean and wmake in the terminal. The solver compiles.

Then I go into my test case where I have added a H field with mapFields from another solver. I delete old solved time steps greater than 0. Then I start my solver with the command fhdFoam.


Also when I only calculate the TEqn only, I get weird looking results and this courant result in the last time step:
Courant Number mean: 49.3029 max: 353.347

[SHUBHAM9595]

The steps look correct. Regarding the high Courant number.....have u already tried using the following options in ur controlDict ?

maxCo         1;

adjustTimeStep  true;

[CharIie]

Yes, that is what I did and gives me first results.
But computational time is 10x longer.

Do you have an idea how I can include the the molecular transport.
Now the temperature is static in my results and does not change and is not influenced by the fluid flow.

[CharIie]

Originally Posted by CharIie

Ok, this brings me back to my initial question in this thread:

If the temperture is not directly solved, but via enthalpy based or internal energy based equation in EEeqn; how do I make use of the temperature, since a physical value for T is shown to me in Paraview?

How can I depict the current temperature in the momentum equation for simple calculations?

If someone might be interested some day. T can be used buoyantPimpleFoam:
Simply use:

const volScalarField& T = mesh.lookupObject<volScalarField>("T");