Dear DuMux Community,
I hope this email finds you well.
I am writing to kindly ask for your support in the following:
a) I am using porous medium flow
b) I am using injection problem
I need to inject gas between injectionDuration_ = 4.7304e8 seconds and time =
injectionEnd_ = 4.80816e8 seconds (3months) but it is not working. Once it
reaches the time it will abort.
I tired the following:
Method:1
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes bcTypes;
bcTypes.setAllNeumann();
return bcTypes;
}
NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
{
NumEqVector values(0.0);
if (time_ > injectionDuration_ && time_ < injectionEnd
&& globalPos[0] <3.15 + eps_ && globalPos[0] > 3.075 - eps_ &&
globalPos[1] < 0.9*this->gridGeometry().bBoxMax()[1])
{
values[Indices::conti0EqIdx + FluidSystem::CH4Idx] = -5e-5;
values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
}
return values;
}
Method 2
NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
{
// initialize values to zero, i.e. no-flow Neumann boundary conditions
NumEqVector values(0.0);
if (time_ < injectionDuration_
&& globalPos[0] <3.15 + eps_ && globalPos[0] > 3.075 - eps_ &&
globalPos[1] < 0.9*this->gridGeometry().bBoxMax()[1])
{
values[Indices::conti0EqIdx + FluidSystem::CH4Idx] = 0.0;
values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
}
else if (time_ > injectionDuration_
&& globalPos[0] <3.15 + eps_ && globalPos[0] > 3.075 - eps_ &&
globalPos[1] < 0.9*this->gridGeometry().bBoxMax()[1])
{
values[Indices::conti0EqIdx + FluidSystem::CH4Idx] = -5e-5;
values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
}
else if (time_ < injectionEnd_
&& globalPos[0] <3.15 + eps_ && globalPos[0] > 3.075 - eps_ &&
globalPos[1] < 0.9*this->gridGeometry().bBoxMax()[1])
{
values[Indices::conti0EqIdx + FluidSystem::CH4Idx] = -5e-5;
values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
}
else if (time_ > injectionEnd_
&& globalPos[0] <3.15 + eps_ && globalPos[0] > 3.075 - eps_ &&
globalPos[1] < 0.9*this->gridGeometry().bBoxMax()[1])
{
values[Indices::conti0EqIdx + FluidSystem::CH4Idx] = 0.0
values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
}
return values;
}
Looking forward to your help.
Best Regards,
Mohammad Hodroj
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program 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 2 of the License, or *
* (at your option) any later version. *
* *
* This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief The two-phase porousmediumflow problem for exercise-basic
*/
#ifndef DUMUX_EX_BASIC_PROBLEM_2P_HH
#define DUMUX_EX_BASIC_PROBLEM_2P_HH
#include <dumux/common/properties.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/porousmediumflow/problem.hh>
namespace Dumux {
/*!
* \ingroup TwoPModel
* \ingroup ImplicitTestProblems
* \brief Gas injection problem where a gas (here nitrogen) is injected into a fully
* water saturated medium. During buoyancy driven upward migration the gas
* passes a high temperature area.
*
* The domain is sized 60 m times 40 m.
*
* For the mass conservation equation neumann boundary conditions are used on
* the top, on the bottom and on the right of the domain, while dirichlet conditions
* apply on the left boundary.
*
* Gas is injected at the right boundary from 7 m to 15 m at a rate of
* 0.001 kg/(s m), the remaining neumann boundaries are no-flow
* boundaries.
*
* At the dirichlet boundaries a hydrostatic pressure and a gas saturation of zero a
*
* This problem uses the \ref TwoPModel model.
*/
template<class TypeTag>
class Injection2PProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
enum { dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
};
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
Injection2PProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
// initialize the tables of the fluid system
FluidSystem::init(/*tempMin=*/273.15,
/*tempMax=*/423.15,
/*numTemp=*/50,
/*pMin=*/0.0,
/*pMax=*/30e6,
/*numP=*/300);
// name of the problem and output file
// getParam<TYPE>("GROUPNAME.PARAMNAME") reads and sets parameter PARAMNAME
// of type TYPE given in the group GROUPNAME from the input file
name_ = getParam<std::string>("Problem.Name");
// depth of the aquifer, units: m
aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth");
// the duration of the injection, units: second
injectionDuration_ = getParam<Scalar>("Problem.InjectionDuration");
injectionEnd_ = getParam<Scalar>("Problem.InjectionEnd");
unsigned int codim = GetPropType<TypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethod::box? dim : 0;
permeability_.resize(fvGridGeometry->gridView().size(codim));
//permeability_.resize(fvGridGeometry->numDofs());
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief Returns the problem name
*
* This is used as a prefix for files generated by the simulation.
*/
std::string name() const
{ return name_+"-2p"; }
/*!
* \brief Returns the temperature \f$ K \f$
*/
Scalar temperature() const
{
return 273.15 +15 ; // [K]
}
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param globalPos The position for which the bc type should be evaluated
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes bcTypes;
// set the left of the domain (with the global position in "0 = x" direction as a Dirichlet boundary
bcTypes.setAllNeumann();
return bcTypes;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet
* boundary segment
*
* \param globalPos The global position
*/
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* \param globalPos The position of the integration point of the boundary segment.
*/
NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
{
// initialize values to zero, i.e. no-flow Neumann boundary conditions
NumEqVector values(0.0);
// if we are inside the injection zone set inflow Neumann boundary conditions
// using < boundary + eps_ or > boundary - eps_ is safer for floating point comparisons
// than using <= or >= as it is robust with regard to imprecision introduced by rounding errors.
if (time_ > injectionEnd_
&& globalPos[0] <3.15 + eps_ && globalPos[0] > 3.075 - eps_ && globalPos[1] < 0.9*this->gridGeometry().bBoxMax()[1])
{
// inject nitrogen. negative values mean injection
// units kg/(s*m^2)
values[Indices::conti0EqIdx + FluidSystem::CH4Idx] = -5e-5;
values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
}
return values;
}
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*
* \param globalPos The position for which the source term should be evaluated
*/
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{
return NumEqVector(0.0);
}
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param globalPos The position for which the initial condition should be evaluated
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
// get the water density at atmospheric conditions
const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5);
// assume an intially hydrostatic liquid pressure profile
// note: we subtract rho_w*g*h because g is defined negative
const Scalar pw = 1.0e5 - densityW*this->spatialParams().gravity(globalPos)[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]);
values[Indices::pressureIdx] = pw;
values[Indices::saturationIdx] = 0.0;
return values;
}
const std::vector<Scalar>& getpermeability()
{
return permeability_;
}
void updateVtkOutput(const SolutionVector& curSol)
{
auto fvGeometry = localView(this->gridGeometry());
for (const auto& element : elements(this->gridGeometry().gridView()))
{
const auto elemSol = elementSolution(element, curSol, this->gridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
VolumeVariables volVars;
volVars.update(elemSol, *this, element, scv);
const auto dofIdxGlobal = scv.dofIndex();
permeability_[dofIdxGlobal] = volVars.permeability();
}
}
}
// \}
//! set the time for the time dependent boundary conditions (called from main)
void setTime(Scalar time)
{ time_ = time;
//std::cout<<"injection problem\n";
this->spatialParams().setTime(time_);
//std::cout<<"injection problem after\n";
}
private:
static constexpr Scalar eps_ = 1e-6;
std::vector<double> permeability_;
std::string name_; //! Problem name
Scalar aquiferDepth_; //! Depth of the aquifer in m
Scalar injectionDuration_; //! Duration of the injection in seconds
Scalar injectionEnd_;
Scalar time_;
};
} //end namespace Dumux
#endif
_______________________________________________
DuMux mailing list
DuMux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
