[Cyp]

Try to impose a delta P (fixedValue for inlet and fixedValue = 0 for outlet).

If you want a immobile system, I think a zeroGradient for inlet and fixedValue = 0 for outlet should work...

Note that since we solve only the pressure and the temperature in the solver, the boundary condition for the velocity field are not accounted for.

[nep101]

The problem we try to solve is this onset of thermal convection in a porous cylinder with conducting wall. So there is no inlet or outlet.
I also thougt that there is maybe a problem with the boundaries. In most cases when i got a floiting point exception, it was a problem with the bopundaries.

But the mainproblem is (i think), that the pressureequation is not properly solved or coded.

Cheers

[TomB]

Gday Cyp and Nep

Thanks so much for your input and collaboration, this has been a huge help to me.

I re-wrote my solver, based on the laplacianFoam. I wrote up some comments along the way.

I'd like to look into it to 'see' how pressure and velocity are linked. That is what I will do next, to make sure its working as expected at this stage.

I'm going to try to set up a pressure gradient with the boundary conditions as you two were discussing.

my method is
break the problem into a number of smaller steps on Cyps guidance
1. solve simple Darcy's Law, based on laplacianFoam (exclude the Temp term)
2. include a scalar transport equation
3. then look at adding temperature
4. then look at correcting T,U
5. add the (Ra*T*k) term to the U eqn
Cheers

Tom

[nep101]

Hi Tom,

wich equation you want to solve?
Is it necessary to handle the equations in a dimensionless form.
What i try to say is, that MAYBE (im no Pro in Openfoam) it is not so easy to handle nondimensionless equations in Openfoam.

Btw how is it going on??

Cheers
Werner

[TomB]

Hi Werner

Im solving for convection in a box of porous media. Sounds very similar to your problem otherwise. I do want to do it in terms of the Rayleigh number and non-dimensional temperature.

I will fudge the units if i have to ! it seems to work ok with the dimensions switched off in the /opt location. I just dont like fiddling with those settings bcoz im an amatuer.

At the moment, i am trying to make the darcyLaplacian solver with some temperature terms added. I seem to be missing a curly bracket, but I'm not sure where.

Did you get a pressure gradient working?

Cheers

Tom
edit

I added my solver and the TEqn it needs. Im having issues now with boundaryConditions. It doesnt like taking the gradient of a calculated internalField.

[nep101]

HI Tom,

no i dont got my continuity equation working correctly.

Send me your solver, maybe i find your problem.

Cheers
Werner

[TomB]

HI Werner

my solver is in the post above. I added another curly bracket and it works, but it doesnt make sense. i use bracket matching and indenting so i should be able to see where its missin..

I think my set up needs a setFields, so that the initial pressure is set before it starts to solve the fvm::laplacian(p)

[nep101]

Hi Tom,

which vers. of OF you use?

[TomB]

Hi Nep

Im using OpenFoam 2.0.1

Im having an issue with the boundary. I get the error message,

gradientInternalCoeffs cannot be called for a calculatedFvPatchField.

it suggests I migh be using a default boundary condition. I copied the boundary conditions from buoyantBoussinesqSimpleFoam

I'll work on this some more tmw ( about 8 hours from now zzzzzzzzzzzzzzzz)

[letzel]

I am just wondering whether this old thread has ever been solved. I think my thread may need a similar solution.

Cheers,
Marcus

[manuc]

Hello All
I am trying to include heat transfer in DarcyFoam solver available in porousMultiphaseFoam (https://github.com/phorgue/porousMul...ties/darcyFoam.

The solver is intended to simulate Natural convective flow in porous media filled cavity. So I would have zero velocity at all walls. Two of the walls would have a temperature (different from each other) with remaining being adiabatic. In the presence of gravitational force a flow is expected to be induced.

Intended equation to be solved:

u= K/mu (-grad(p)+ rho*g)= -K/mu(grad(p-rho_0*g*z)) -(K/mu)*(rho_0*beta*(T-Tref))*g

where

• u=velocity (m/s)
• K=permeability (m^2)
• mu=viscosity (kg/ms)
• p=pressure (kg/ms^2)
• g=gravity vector^(m/s^2)
• rho_0=reference density (kg/m3)
• beta=coeffcient of volume expansion (1/K)

In order to have such a solver, I modified the Main .C file as:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation, either version 3 of the License, or (at your
    option) any later version.

    foam-extend is distributed in the hope that it will be useful, but
    WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Application
    laplacianFoam

Description
    Solves a simple Laplace equation, e.g. for thermal diffusion in a solid.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "simpleControl.H"

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

int main(int argc, char *argv[])
{

#   include "setRootCase.H"

#   include "createTime.H"
#   include "createMesh.H"

    simpleControl simple(mesh);
    
#   include "readGravitationalAcceleration.H"
#   include "createFields.H"

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

    Info<< "\nCalculating flow and temperature distribution\n" << endl;

    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;
        
       
        while (simple.correctNonOrthogonal())
        {
        
        fvScalarMatrix TEqn
        (
         (eps*rhoCpf+(1.-eps)*rhoCps)*fvm::ddt(T) 
			+ rhoCpf*fvm::div(phi,T) 
			    == 
			  fvm::laplacian(DT, T)
        );
        
        TEqn.solve();
        TEqn.relax();

      

        rhok = rhof*(1.0- beta*(T - TRef));
        
        
        
      fvScalarMatrix pEqn
        (
            fvm::laplacian(-Mf,p)+ fvc::div(phiG)
        );
       
        pEqn.solve();
                phi = pEqn.flux() + phiG;
         
         pEqn.relax();
              
              U = fvc::reconstruct(phi);
              U.correctBoundaryConditions();
             
              Info<< "Vel done" << endl;   
              
      
        
     Info<< "Temp done" << endl;                 
     }   
          
        
#       include "write.H"

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

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

    return 0;
}


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

The createFields.H is modified as:

    Info<< "Reading thermophysical properties\n" << endl;

    Info<< "Reading field T\n" << endl;
    volScalarField T
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    Info<< "Calculating field g.h\n" << endl;
    volScalarField gh("gh", g & mesh.C());
    surfaceScalarField ghf("ghf", g & mesh.Cf());
      
    
    
    Info<< "Reading field p\n" << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
    
  

    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh, dimensionedVector("zero", dimensionSet(0,1,-1,0,0,0,0), vector::zero)
    );
      // mesh,    dimensionedVector("U", dimensionSet(0,1,-1,0,0,0,0), vector::zero)
    
     //permeability
    Info << "Reading field K\n" << endl;
  volScalarField K
   (
    IOobject
    (
        "K",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
  );

//#   include "createPhi.H"

surfaceScalarField phi("phi", fvc::interpolate(U) & mesh.Sf());

mesh.schemesDict().setFluxRequired(p.name());

  Info<< "Reading transportProperties\n" << endl;

    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    );


   Info<< "Reading diffusivity DT\n" << endl;

    dimensionedScalar DT
    (
        transportProperties.lookup("DT")
    );


     // Thermal expansion coefficient [1/K]
    dimensionedScalar beta(transportProperties.lookup("beta"));

    // Reference temperature [K]
    dimensionedScalar TRef(transportProperties.lookup("TRef"));


  // Reference density of fluid [kg/m^3]
   dimensionedScalar rhof
    (
        transportProperties.lookup("rhof")
    );
    
    
    //poroisity
    dimensionedScalar eps
    (
        transportProperties.lookup("eps")
    );


   //heat capacity solid
    dimensionedScalar rhoCps
    (
        transportProperties.lookup("rhoCps")
    );


    //heat capacity fluid
    dimensionedScalar rhoCpf
    (
        transportProperties.lookup("rhoCpf")
    );

   
    //viscosity fluid
    dimensionedScalar mu
    (
        transportProperties.lookup("mu")
    );


// Kinematic density for buoyancy force

volScalarField rhok
    (
        IOobject
        (
            "rhok",
            runTime.timeName(),
            mesh
        ),
        rhof*(1.0 - (beta*(T - TRef)))
    );
    
    surfaceScalarField buoyancyPhi(ghf*fvc::snGrad(rhok)*mesh.magSf());
    
    
     Info<< "Calculating field beta*(g.h)\n" << endl;
    surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf()));
    
    surfaceScalarField Kf=fvc::interpolate(K,"K");
    surfaceScalarField Mf ("Mf",Kf/mu);
    
    
    
    surfaceScalarField rhokf=fvc::interpolate(rhok);
    surfaceScalarField phiG("phiG",(rhokf * Mf)* (g & mesh.Sf()));

When I run the solver on a test case (attached herewith), the pressure diverges.
Test conditions tried and outcome

1. gravity set to zero--code works and gives a conduction profile
2. Equal temp at both walls with non zero gravity--Expected uniform temperature profile--outcome--core dumped
3. Two walls at different temperature, other wall adiabatic, gravity set to a value, coefficient of expansion beta in transport properties set to a value--outcome--core dumped
   

   DICPCG:  Solving for p, Initial residual = 1, Final residual = 1.2601e+08, No Iterations 1000
   Vel done
   Temp done
   DILUPBiCG:  Solving for T, Initial residual = 0.565344, Final residual = 1.85342e-12, No Iterations 1
   DICPCG:  Solving for p, Initial residual = 1, Final residual = 5.48003e+10, No Iterations 1000
   Vel done
   Temp done
   ExecutionTime = 1.51 s  ClockTime = 2 s
   
   Time = 4
   
   DILUPBiCG:  Solving for T, Initial residual = 0.818783, Final residual = 8.10268e-08, No Iterations 7
   DICPCG:  Solving for p, Initial residual = 1, Final residual = 8.65643e+08, No Iterations 1000
   Vel done
   Temp done

Important parameters to be changed for testing

• constant/g
• constant/transportProperties/beta (beta is the volume expansion coeffcient, that induces a forcing term equivalent to density change induced effect.

Please note: to calculate rhok I tried rhof*(1.0 - (beta*(T - TRef))) and rhof*(0.0 - (beta*(T - TRef))) depending on whether p is treated as pressure or pressure-rho*g*h. But both of them gave similar error.


Can someone please suggest what could be the issue and how I should approach the same.

please find the solver and testcase attached.