Commit b74b0382 authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

Adding the energy equation to the Navier-Stokes solver.

parent 43c75758

examples/navier-stokes/navier-stokes.cfg.desc

+1 −0
Original line number Diff line number Diff line
@@ -16,6 +16,7 @@ group Problem 
   real start-up(0.0)                   [Start-up time of the riser.];

   real mu(!)                           [Viscosity.];

   real p0(!)                           [Pressure for the initial condition.];

   real gamma(!)                         [Heat capacity ratio.];

   real T(!)                            [Temperature.];

   real R(!)                            [Gas constant.]; 

   real gravity(9.81)                   [Gravity force.];

examples/navier-stokes/navierStokesBoundaryConditions.h

+3 −2
Original line number Diff line number Diff line
@@ -40,13 +40,14 @@ class navierStokesBoundaryConditions 
   void apply( const RealType& time,

               Vector& rho,

               Vector& u1,

               Vector& u2 );

               Vector& u2,

               Vector& temperature );



   protected:



   const MeshType* mesh;



   RealType maxInflowVelocity, startUp;

   RealType maxInflowVelocity, startUp, T;

};



#include "navierStokesBoundaryConditions_impl.h"

examples/navier-stokes/navierStokesBoundaryConditions_impl.h

+8 −1
Original line number Diff line number Diff line
@@ -34,6 +34,7 @@ bool navierStokesBoundaryConditions< Mesh >::init( const tnlParameterContainer& 
   this -> maxInflowVelocity = parameters. GetParameter< double >( "max-inflow-velocity" );

   //this -> maxOutflowVelocity = parameters. GetParameter< double >( "max-outflow-velocity" );

   this -> startUp = parameters. GetParameter< double >( "start-up" );

   this -> T = parameters. GetParameter< double >( "T" );

   return true;

}
@@ -48,7 +49,8 @@ template< typename Mesh > 
void navierStokesBoundaryConditions< Mesh >::apply( const RealType& time,

                                                    Vector& rho,

                                                    Vector& u1,

                                                    Vector& u2 )

                                                    Vector& u2,

                                                    Vector& temperature )

{

   /****

     * Set the boundary conditions.
@@ -85,6 +87,8 @@ void navierStokesBoundaryConditions< Mesh >::apply( const RealType& time, 


         rho[ c1 ] = rho[ c2 ];

         rho[ c3 ] = rho[ c4 ];

         temperature[ c1 ] = this->T;

         temperature[ c3 ] = this->T;

          //rho[ c3 ] = this -> p_0 / ( this -> R * this -> T );

      }
@@ -114,6 +118,9 @@ void navierStokesBoundaryConditions< Mesh >::apply( const RealType& time, 


         rho[ c1 ] = rho[ c2 ];

         rho[ c3 ] = rho[ c4 ];

         temperature[ c1 ] = this->T;

         temperature[ c3 ] = this->T;



      }

      /*rho_u1[ c1 ] = rho[ c1 ] * this -> u1[ c1 ];

      rho_u2[ c1 ] = rho[ c1 ] * this -> u2[ c1 ];

examples/navier-stokes/navierStokesSolver.h

+1 −3
Original line number Diff line number Diff line
@@ -93,11 +93,9 @@ class navierStokesSolver 


   MeshType mesh;



   //tnlVector< RealType, DeviceType, IndexType > rho, u1, u2, p;



   DofVectorType dofVector, rhsDofVector;



   RealType p_0, gravity;

   RealType p_0, gravity, T;



   EulerScheme eulerScheme;

examples/navier-stokes/navierStokesSolver_impl.h

+13 −10
Original line number Diff line number Diff line
@@ -47,7 +47,8 @@ __device__ void computeVelocityFieldCuda( const Index size, 


template< typename Mesh, typename EulerScheme >

navierStokesSolver< Mesh, EulerScheme > :: navierStokesSolver()

: p_0( 0.0 )

: p_0( 0.0 ),

  T( 0.0 )

{



   this -> mesh. setName( "navier-stokes-mesh" );
@@ -133,8 +134,9 @@ bool navierStokesSolver< Mesh, EulerScheme > :: init( const tnlParameterContaine 
    * Set-up model coefficients

    */

   this->p_0 = parameters. GetParameter< double >( "p0" );

   this->T =  parameters. GetParameter< double >( "T");

   nsSolver.setHeatCapacityRatio( parameters. GetParameter< double >( "gamma" ) );

   nsSolver.setMu( parameters. GetParameter< double >( "mu") );

   nsSolver.setT( parameters. GetParameter< double >( "T") );

   nsSolver.setR( parameters. GetParameter< double >( "R") );

   nsSolver.setGravity( parameters. GetParameter< double >( "gravity") );

   if( ! this->boundaryConditions.init( parameters ) )
@@ -174,16 +176,16 @@ bool navierStokesSolver< Mesh, EulerScheme > :: init( const tnlParameterContaine 
template< typename Mesh, typename EulerScheme >

bool navierStokesSolver< Mesh, EulerScheme > :: setInitialCondition( const tnlParameterContainer& parameters )

{

   tnlSharedVector< RealType, DeviceType, IndexType > dofs_rho, rho_u1, rho_u2;

   tnlSharedVector< RealType, DeviceType, IndexType > dofs_rho, dofs_rho_u1, dofs_rho_u2, dofs_e;

   const IndexType& dofs = mesh. getDofs();

   dofs_rho.    bind( & dofVector. getData()[ 0        ], dofs );

   rho_u1.   bind( & dofVector. getData()[     dofs ], dofs );

   rho_u2.   bind( & dofVector. getData()[ 2 * dofs ], dofs );

   dofs_rho_u1. bind( & dofVector. getData()[     dofs ], dofs );

   dofs_rho_u2. bind( & dofVector. getData()[ 2 * dofs ], dofs );



   dofs_rho. setValue( p_0 / ( this->nsSolver.getR() *

                               this->nsSolver.getT() ) );

   rho_u1. setValue( 0.0 );

   rho_u2. setValue( 0.0 );

   dofs_rho_u1. setValue( 0.0 );

   dofs_rho_u2. setValue( 0.0 );



   const IndexType& xSize = mesh. getDimensions(). x();

   const IndexType& ySize = mesh. getDimensions(). y();
@@ -199,8 +201,9 @@ bool navierStokesSolver< Mesh, EulerScheme > :: setInitialCondition( const tnlPa 


         dofs_rho. setElement( c, p_0 / ( this->nsSolver.getR() *

                                          this->nsSolver.getT() ) );

         rho_u1. setElement( c, 0.0 );

         rho_u2. setElement( c, 0.0 );

         dofs_rho_u1. setElement( c, 0.0 );

         dofs_rho_u2. setElement( c, 0.0 );

         dofs_e. setElement( c, p_0 / ( this->nsSolver.getHeatCapacityRatio() - 1.0 ) );



      }

   return true;