CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Conservative form of Navier Stokes equation. (http://www.cfd-online.com/Forums/openfoam-solving/96427-conservative-form-navier-stokes-equation.html)

balkrishna January 23, 2012 01:43

Conservative form of Navier Stokes equation.
 
Hi,

I wish to modify the pimpleFoam solver for varying density. As a first step, I decided to change the momentum equation to the conservative form. To do so I changed the momentum equation to
Code:

    fvVectorMatrix UEqn
    (
    fvm::ddt(rho, U) +
    fvm::div(phi, U)
    - fvm::laplacian(mu, U)
    - fvc::div(mu*dev2(fvc::grad(U)().T()))
    );

here rho is declared as a volScalarField. The pressure equation is coded in the following way :
Code:

{

    Info<<" In the most dreaded pEqn "<<endl;
    volScalarField rUA = 1.0/UEqn.A();
    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));

    U = rUA*UEqn.H();

    phi = fvc::interpolate(rho)*
    (
        (fvc::interpolate(U) & mesh.Sf())
      + fvc::ddtPhiCorr(rUA, rho, U, phi)
    );

    surfaceScalarField buoyancyPhi = -rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
    phi += buoyancyPhi;


    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
    {
        fvScalarMatrix p_rghEqn
        (
            fvc::ddt(rho) +                                                                                                     
          fvc::div(phi)
          - fvm::laplacian(rhorUAf, p_rgh)
        );

        p_rghEqn.solve
        (
            mesh.solver
            (
                p_rgh.select
                (
                    (
                        finalIter
                    && corr == nCorr-1
                    && nonOrth == nNonOrthCorr
                    )
                )
            )
        );

        if (nonOrth == nNonOrthCorr)
        {
            // Calculate the conservative fluxes                                                                                                                                 
            phi += p_rghEqn.flux();

            // Explicitly relax pressure for momentum corrector                                                                                                                   
            p_rgh.relax();

            // Correct the momentum source with the pressure gradient flux                                                                                                       
            // calculated from the relaxed pressure                                                                                                                               
            U += rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf);
        }                                                                                                                                                                         
    }                                                                                                                                                                             
                                                                                                                                                                                 
    p = p_rgh + rho*gh;                                                                                                                                                           
                                                                                                                                                                                 
                                                                                                                                 

    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

As a first test I decided to simulate for a case of non changing density for a pipe flow.My rho is specified in the following manner.
Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  1.7.0                                |
|  \\  /    A nd          | Web:      www.OpenFOAM.com                      |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version    2.0;
    format      ascii;
    class      volScalarField;
    location    "0";
    object      rho;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -3 0 0 0 0 0];

internalField  uniform 1000;

boundaryField
{
    outlet
    {
        type            fixedValue;
        value          uniform 1000;
    }
    inlet
    {
        type            fixedValue;
        value          uniform 1000;
    }
    wall
    {
      type            fixedValue;
      value          uniform 1000;
    }
    frontBack
    {
        type            empty;
    }
}


// ************************************************************************* //
0/rho (END)

This gives highly unphysical results as far as the velocity profiles are concerned. Could someone tell me where I might be going wrong in the momentum equation or pressure equation ?

alberto January 24, 2012 12:22

Take a look at rhoPimpleFoam, and for a solver using p_rgh, you could take a look at buoyantPimpleFoam.

Best,

balkrishna January 25, 2012 09:33

Solved.
 
thanks for d help.


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