[adkar]

Hello foamers, I have a simple question regarding writing equations in OpenFoam. Could anybody please help me understand this?

I would like to write an equation like the following where I have to add div A and grad A to U equation. Here, A is set in boundary condition.

But when I solve this Equation, it gives me the same result as momentum equation without the addition of these two terms.
Here div (constant1, A) - grad (A) acts something like opposing force.

Do I have use fvc::Su () or fvc:: Sp ()?
like fvc::Su(fvc::div(constant1, A) - fvc::grad(A) )
Or am I missing something very basic?

I am not sure how to write this piece of code.

Please see the following code:

****************************************:
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

solve(UEqn == -fvc::grad(p) + fvc::div(constant1, A) - fvc::grad(A);

Then goes the Piso loop to couple p and U without A getting involved.

********************************************

Thank you in advance

[adkar]

Could anybody please help??

[adrin]

I don't know much about openfoam, neither is it obvious what you're trying to achieve. But, I have a more fundamental question/concern regarding what you're trying to do.

It's possible openfoam's objects have certain (unusual) overloading options, but from equation -fvc::grad(p) + fvc::div(constant1, A) - fvc::grad(A); I can only suggest that you consider the accuracy of your formulation first, because:

grad(p) is a vector, and if A is a vector then div(A) is a scalar and grad(A) is a tensor. If A is a tensor then div(A) would be a vector but grad(A) would still not be correct. So, again, unless openfoam's objects have certain optional capabilities that allow one to do poor math you have a serious rank mismatch, both physically and mathematically.

adrin

[adkar]

Originally Posted by adrin

I don't know much about openfoam, neither is it obvious what you're trying to achieve. But, I have a more fundamental question/concern regarding what you're trying to do.

It's possible openfoam's objects have certain (unusual) overloading options, but from equation -fvc::grad(p) + fvc::div(constant1, A) - fvc::grad(A); I can only suggest that you consider the accuracy of your formulation first, because:

grad(p) is a vector, and if A is a vector then div(A) is a scalar and grad(A) is a tensor. If A is a tensor then div(A) would be a vector but grad(A) would still not be correct. So, again, unless openfoam's objects have certain optional capabilities that allow one to do poor math you have a serious rank mismatch, both physically and mathematically.

adrin

First of all, I would like to thank you Adrin for taking your time in answering my question. I make a mistake in generalizing the terms as A when I pasted the piece of code thinking it would be harder to explain, but it ruined the meaning of my question.

The actual code is the following from mhdFoam solver.
Could you please be kind enough for having a look at it one more time?

B is a vector used in the boundary condition. Eg. (0 20 0)


My problem is that when I remove the B-Piso loop part, the solver will not calculate the effect of the extra terms in the momentum equation.


Am I missing something?
Please help. Any help is appreciated.
Thank you very much.

/////////////////////////////////////////////////////////////////////////////////////////////////////////////
#include "fvCFD.H"
#include "OSspecific.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"

    #include "createTime.H"
    #include "createMesh.H"
    #include "createFields.H"
    #include "initContinuityErrs.H"

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< nl << "Starting time loop" << endl;

    while (runTime.loop())
    {
        #include "readPISOControls.H"
        #include "readBPISOControls.H"

        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "CourantNo.H"

        {
            fvVectorMatrix UEqn
            (
                fvm::ddt(U)
              + fvm::div(phi, U)
              - fvc::div(phiB, 2.0*DBU*B)     //Extra term//
              - fvm::laplacian(nu, U)
              + fvc::grad(DBU*magSqr(B))  //Extra term//
            );

            solve(UEqn == -fvc::grad(p));


            // --- PISO loop

            for (int corr=0; corr<nCorr; corr++)
            {
                volScalarField rAU(1.0/UEqn.A());
                volVectorField HbyA("HbyA", U);
                HbyA = rAU*UEqn.H();

                surfaceScalarField phiHbyA
                (
                    "phiHbyA",
                    (fvc::interpolate(HbyA) & mesh.Sf())
                  + fvc::ddtPhiCorr(rAU, U, phi)
                );

                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
                {
                    fvScalarMatrix pEqn
                    (
                        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                    );

                    pEqn.setReference(pRefCell, pRefValue);
                    pEqn.solve();

                    if (nonOrth == nNonOrthCorr)
                    {
                        phi = phiHbyA - pEqn.flux();
                    }
                }

                #include "continuityErrs.H"

                U = HbyA - rAU*fvc::grad(p);
                U.correctBoundaryConditions();
            }
        }

/*
        // --- B-PISO loop                                       /////// I tried to remove this part of B-PISO loop since I would only need the first 
                                                                      /////// equation and the B-Piso loop part would take too much time for calculation 
                                                                      ///////and its effect  is very less. B is given by boundary condition.

        for (int Bcorr=0; Bcorr<nBcorr; Bcorr++)
        {
            fvVectorMatrix BEqn
            (
                fvm::ddt(B)
              + fvm::div(phi, B)
              - fvc::div(phiB, U)
              - fvm::laplacian(DB, B)
            );

            BEqn.solve();

            volScalarField rBA(1.0/BEqn.A());

            phiB = (fvc::interpolate(B) & mesh.Sf())
                + fvc::ddtPhiCorr(rBA, B, phiB);

            fvScalarMatrix pBEqn
            (
                fvm::laplacian(rBA, pB) == fvc::div(phiB)
            );
            pBEqn.solve();

            phiB -= pBEqn.flux();

            #include "magneticFieldErr.H"
        }

        runTime.write();
    }
    */

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

    return 0;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

[FMDenaro]

Originally Posted by adkar

Here, A is set in boundary condition.

what do you mean?

[adrin]

As I mentioned earlier, I'm not familiar enough with openfoam to be of help there. But there is actually an openfoam-specific forum on cfd-online, which I strongly recommend you post your question(s) to. I'm sure openfoam has its own user forum as well, so you're strongly recommended to try them as well

adrin