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 11:29.