CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

How to add Brinkman Term to porousSimpleFoam solver?

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By AmirBaqa1987

Reply
 
LinkBack Thread Tools Display Modes
Old   December 24, 2013, 04:34
Default How to add Brinkman Term to porousSimpleFoam solver?
  #1
Member
 
Arjang Behnoud
Join Date: Oct 2012
Posts: 57
Rep Power: 5
AmirBaqa1987 is on a distinguished road
Hi everyone.

the momentum equation of porousSimpleFoam is :



the last two terms are darcy and Forchheimer resistance in order.UEqn.h and PEqn.h in solver are:

Ueqn.H
Code:
 tmp<fvVectorMatrix> UEqn
    (
        fvm::div(phi, U)
      + turbulence->divDevReff(U)
    );

    UEqn().relax();

    // Include the porous media resistance and solve the momentum equation
    // either implicit in the tensorial resistance or transport using by
    // including the spherical part of the resistance in the momentum diagonal

    tmp<volScalarField> trAU;
    tmp<volTensorField> trTU;

    if (pressureImplicitPorosity)
    {
        tmp<volTensorField> tTU = tensor(I)*UEqn().A();
        pZones.addResistance(UEqn(), tTU());
        trTU = inv(tTU());
        trTU().rename("rAU");

        volVectorField gradp = fvc::grad(p);

        for (int UCorr=0; UCorr<nUCorr; UCorr++)
        {
            U = trTU() & (UEqn().H() - gradp);
        }
        U.correctBoundaryConditions();
    }
    else
    {
        pZones.addResistance(UEqn());

        eqnResidual = solve
        (
            UEqn() == -fvc::grad(p)
        ). initialResidual();

        maxResidual = max(eqnResidual, maxResidual);

        trAU = 1.0/UEqn().A();
        trAU().rename("rAU");
    }
PEqn.H
Code:
if (pressureImplicitPorosity)
{
    U = trTU()&UEqn().H();
}
else
{
    U = trAU()*UEqn().H();
}

UEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
    tmp<fvScalarMatrix> tpEqn;

    if (pressureImplicitPorosity)
    {
        tpEqn = (fvm::laplacian(trTU(), p) == fvc::div(phi));
    }
    else
    {
        tpEqn = (fvm::laplacian(trAU(), p) == fvc::div(phi));
    }

    tpEqn().setReference(pRefCell, pRefValue);
    // retain the residual from the first iteration
    if (nonOrth == 0)
    {
        eqnResidual = tpEqn().solve().initialResidual();
        maxResidual = max(eqnResidual, maxResidual);
    }
    else
    {
        tpEqn().solve();
    }

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

#include "continuityErrs.H"

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

if (pressureImplicitPorosity)
{
    U -= trTU()&fvc::grad(p);
}
else
{
    U -= trAU()*fvc::grad(p);
}

U.correctBoundaryConditions();
now I want to add Brinkman Term to momentum equation, I mean:



for this I've just define in createFields.H and modify UEqn.H like below:

UEqn.H

Code:
 tmp<fvVectorMatrix> UEqn
    (
            fvm::div(phi, U)
         -  fvm::laplacian(mue,U)
         + turbulence->divDevReff(U)
    );
am I right?
shouldn't I modify PEqn.H?

thanks a lot.
your hints are appreciated .
Best Regards,
Arjang
Caio Martins likes this.

Last edited by AmirBaqa1987; December 24, 2013 at 08:44.
AmirBaqa1987 is offline   Reply With Quote

Old   December 25, 2013, 02:39
Default
  #2
Member
 
Arjang Behnoud
Join Date: Oct 2012
Posts: 57
Rep Power: 5
AmirBaqa1987 is on a distinguished road
I'm still waiting for HELP.
AmirBaqa1987 is offline   Reply With Quote

Old   April 30, 2016, 14:59
Default
  #3
New Member
 
Caio Martins Ramos de Oliveira
Join Date: Apr 2016
Posts: 6
Rep Power: 2
Caio Martins is on a distinguished road
Were you able to find a solution? Were you right?
Caio Martins is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
pyFoam and paraview eelcovv OpenFOAM Paraview & paraFoam 28 May 30, 2016 09:23
A problem with porousSimpleFoam solver AmirBaqa1987 OpenFOAM Running, Solving & CFD 4 October 30, 2015 11:54
How to add Radiation to existing buoyantBoussinesqSimpleFoam solver MasAmi OpenFOAM Running, Solving & CFD 4 January 16, 2014 07:17
Add source term in alphaEqn.H of interFoam tayo OpenFOAM 1 October 23, 2013 03:40
add source term in energy equation chaolian OpenFOAM Programming & Development 4 November 8, 2012 23:22


All times are GMT -4. The time now is 19:20.