[behzad-cfd]

Hi dear foamers,

I want to add new terms in momentum equation. I need some help to do. I'm trying to add magneto static equations in interFoam solver. To do this:
1- magnetic scalar potential should be solved, to get B (magnetic induction) and H (magnetic field). this part is done and i get reasonable results.
2- main part is to add a term (divergence of magnetic stress tensor ) to momentum equation that takes care of deformations caused by magnetic field. Actually I don't have any idea how should i write this terms in momentum equation(see page 5 eqn (2.12) & Apendix).

Any guidance is appreciated.

This is the link of my reference paper http://individual.utoronto.ca/afkham...Ferro_Drop.pdf

All of details are depicted in paper.

Thanks in advance,
Behzad

[behzad-cfd]

Let me make my question more specific. want to add this terms in u and p equations.


Now because 1/2*mu*H^2 term is proportional to identity tensor so it will be absorbed to P field (is that right?). i wrote this terms in uEqn like this:

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

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      - fvc::div(interfacePsi.muM()*sqr(Hm))
      - fvm::laplacian(muEff, U)
      - (fvc::grad(U) & fvc::grad(muEff))
    );

    UEqn.relax();

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    fvc::snGrad(psi)*fvc::interpolate(interfacePsi.sigmaC())*
                    fvc::interpolate(interfacePsi.delta())
                  - double(0.5)*interfacePsi.mu0()*dimMu*fvc::snGrad(magSqr(Hm))
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
            )
        );
    }

Edited part of pEqn:

 phi = phiU +
    (
                fvc::snGrad(psi)*fvc::interpolate(interfacePsi.sigmaC())*fvc::interpolate(interfacePsi.delta()) 
              - double(0.5)*interfacePsi.mu0()*dimMu*fvc::snGrad(magSqr(Hm)) 
              - ghf*fvc::snGrad(rho)
    )*rAUf*mesh.magSf();

I get good results for small deformations but bad results for high deformation in comparison with reference paper results. Now is this implementation right?
whats the deference of writing terms in left or right hand side of Ueqn?
my OF version and base solver is 2.1.1, interFoam.

Any comment is appreciated

[tom_flint2012]

Hi,

I know this thread is old, but I have just come across it.

I don’t think you need to add any terms to the reconstruct expression. The term that is multiplied by the identity tensor can just go as grad(A) instead of div (A I) where A is your expression in the momentum equation.
How are you getting on with this?

I’m working on a similar problem with an induction formulation.