Hi Mahmoud,
you also need to extract the correct energy carried by the fluid in the
production well. Otherwise your domain will heat up there.
For solution-dependent point sources you can use the split interface. First you
add your point sources (in addPointSources) only providing positions as
argument to the constructor.
Then you overload the function pointSource (also in the problem) with something
like
template<class ElementVolumeVariables>
void pointSource(PointSource& source,
const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
const auto& pos = source.position();
const auto& volVars = elemVolVars[scv];
if (pos[0] < 40.0) // injection
{
const Scalar volumeSource = injectionRate; // injectionRate is
positive and m^3/s
const Scalar massSource =
volumeSource*IAPWSH2O::density(injectionTemperature, injectionPressure);
const Scalar energySource =
massSource*IAPWSH2O::enthalpy(injectionTemperature, injectionPressure);
source = NumEqVector({ massSource, energySource });
}
else // production
{
const Scalar volumeSource = productionRate; // productionRate is
negative and m^3/s
const Scalar massSource = volumeSource*volVars.density(0); // using
current control volume's water density
const Scalar energySource = massSource*volVars.enthalpy(0); // using
current control volume's water enthalpy
source = NumEqVector({ massSource, energySource });
}
}
This function will be called called for every point source that you added
whenever the point source contribution is added to the equation residual.
Whether your values are realistic or not I cannot help you. Initial conditions
seem fine.
Hope this helps
Timo
> On 9. Apr 2021, at 12:51, Mahmoud Atef Mahmoud Mohamed Aboelseoud S277151
> <s277...@studenti.polito.it> wrote:
>
>
>
> Dear Timo,
>
> I cannot thank you enough for all your help and guidance. I wrote the
> pointsource function as shown below. could you please let me know if it is OK
> now or not ? I used 3 pointsources to represent each of the production well
> and the injection well because I wanted to represent fully penetrating wells.
>
> void addPointSources(std::vector<PointSource>& pointSources) const
> {
>
> const auto enthalpy = IapwsH2O::liquidEnthalpy(293, 13.0e5);
> const Scalar energyFlow = 1.93 * enthalpy; // 0.8 Kg/s * enthalpy
> J/Kg
> // The injection well (source term)
> pointSources.push_back(PointSource({35, 45, 15}, {1.93,
> energyFlow}));
> pointSources.push_back(PointSource({35, 45, 25}, {1.93,
> energyFlow}));
> pointSources.push_back(PointSource({35, 45, 35}, {1.93,
> energyFlow}));
> // The production well (sink term)
> pointSources.push_back(PointSource({65, 45, 15}, {-1.93, 0}));
> pointSources.push_back(PointSource({65, 45, 25}, {-1.93, 0}));
> pointSources.push_back(PointSource({65, 45, 35}, {-1.93, 0}));
> }
>
>
> However, the results I'm getting on paraview are a bit confusing to me
> because I imposed Dirichlet boundary conditions to all boundaries which are
> basically pressure and temperature gradients and they are the same as the
> initial conditions
>
> PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
> {
>
> return initialAtPos(globalPos);
> }
>
> PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
> {
> PrimaryVariables values(0.0);
> Scalar densityW = 990; // Kg/m^3
> Scalar depth = this->gridGeometry().bBoxMax()[2] - globalPos[2];
> // Hydrostatic Pressure
> values[pressureIdx] = 11.0e5 -
> densityW*this->spatialParams().gravity(globalPos)[2]*depth; //Pascal
> // Geothermal Gradient
> values[temperatureIdx] = 358.0 + depth*0.03; //Kelvin
>
> return values;
> }
>
>
> However I could see what looks like numerical errors affecting certain spots
> of the pressure boundaries (the low permeability layers on top and bottom)
> (picture attached) and I also could see high temperatures (400+ Kelvin) in
> the production well while the highest temperature initially is around 359.5
> Kelvin (50 m thick domain) (picture attached). I would really appreciate it
> if you could let me know if there's something wrong here. Another thing
> causing confusion is that the simulation proceeds really quickly and when I
> impose a mass rate of 0.38 Kg/s for each of the 3 pointsources of a well
> which is equivalent to 100 m^3/day for the well, it's like I am not injecting
> anything at all so I had to impose higher mass rates !
>
> PS: I used "using IapwsH2O = Components::H2O<Scalar>;" and
> "const auto enthalpy = IapwsH2O::liquidEnthalpy(293, 13.0e5);"
>
> because when I tried
>
> "using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;"
> "const auto enthalpy = FluidSystem::liquidEnthalpy(293, 13.0e5);" it gave me
> this error:
>
>
> error: ‘liquidEnthalpy’ is not a member of
> ‘Dumux::OnePNIProblem<Dumux::Properties::TTag::OnePNICCMpfa>::FluidSystem’
> {aka ‘Dumux::FluidSystems::OnePLiquid<double, Dumux::Components::H2O<double>
> >’}
> 239 | const auto enthalpy = FluidSystem::liquidEnthalpy(293,
> 13.0e5);
>
>
> Thanks again Timo for your help !
>
>
> Best Regards,
>
> Mahmoud Atef
>
>
>
>
>
>
>
>
>
>
>
>
> -------------------------------------------------------------------------------------------------
> Il 2021-04-08 15:47 Timo Koch ha scritto:
>>> On 8. Apr 2021, at 13:42, Mahmoud Atef Mahmoud Mohamed Aboelseoud S277151
>>> <s277...@studenti.polito.it> wrote:
>>> Dear Timo,
>>> Thanks a lot for your reply. I didn't use the pointSources because I didn't
>>> know how to assign the temperature inside such function for the injection
>>> well because I want to impose a rate and also a temperature.
>> Hi Mahmoud,
>> for the pointSource for the OnePNI model you need to specify mass
>> (kg/s) and energy source (J/s) (Indices::conti0EqIdx (=0) and
>> Indices::energyEqIdx (=1)).
>> “NumEqVector" has size 2 because you have two equations.
>> The energy flow can be computed as massFlow*enthalpy(p,T) (neglecting
>> heat conduction which is usually unimportant here). So the temperature
>> enters via the enthalpy.
>> You can evaluate the correct enthalpy for your fluid system with
>> using FluidSystem = GetPropType<TypeTag, Properties:: FluidSystem>;
>> const auto enthalpy = FluidSystem::liquidEnthalpy(T, p);
>> for a temperature T (in K) and pressure p (in Pa).
>> In your examples below you build a temporary object of type FluidState
>> and set the temperature on this object. However this fluid state
>> object is used nowhere, so it does nothing.
>> Consequently you seem to inject mass without accounting for the energy
>> the injected fluid carries.
>> Hope this helps
>> Timo
>>> I'm using the OnePNI model in my problem. I will try to be more specific
>>> this time as I'm afraid I could be making some mistake that I'm not aware
>>> of. To define both the injection and production wells inside the
>>> sourceAtPos function, I'm writing the following:
>>> NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
>>> {
>>> NumEqVector values;
>>> using FluidState = GetPropType<TypeTag, Properties::FluidState>;
>>> FluidState fs;
>>> // we define one source term and one sink term
>>> //Injection well (Source Term)
>>> if ((globalPos[2] > 15 - eps_ && globalPos[2] < 35 + eps_) &&
>>> (globalPos[1] > 45 - eps_ && globalPos[1] < 45 + eps_) && (globalPos[0] >
>>> 35 - eps_ && globalPos[0] < 35 + eps_))
>>> {
>>> fs.setTemperature(293);
>>> values[contiWEqIdx] = 1e-1; // kg/(s*m^3)
>>> }
>>> //Production well (Sink Term)
>>> else if ((globalPos[2] > 15 - eps_ && globalPos[2] < 35 + eps_) &&
>>> (globalPos[1] > 45 - eps_ && globalPos[1] < 45 + eps_) && (globalPos[0] >
>>> 65 - eps_ && globalPos[0] < 65 + eps_))
>>> {
>>> values[contiWEqIdx] = -1e-1; // kg/(s*m^3)
>>> }
>>> else {
>>> values[contiWEqIdx]=0;
>>> }
>>> return values;
>>> }
>>> where contiWEqIdx is defined in an enum as contiWEqIdx =
>>> Indices::conti0EqIdx
>>> and my model dimensions and cells are as follows:
>>> [Grid]
>>> Positions0 = 0 10 30 70 90 100
>>> Positions1 = 0 10 30 70 90 100
>>> Positions2 = 0 10 40 50
>>> Cells0 = 1 2 4 2 1
>>> Cells1 = 1 2 4 2 1
>>> Cells2 = 1 3 1
>>> where there are upper and lower boundary layers with very low permeability
>>> (order of e-25 m^2) each 10 m thick. I'm not using any refinement for the
>>> moment around the wells. Is there anything wrong I'm making here in the
>>> above function and which causes me to have a very limited range of mass
>>> rate below which the simulation proceeds very quickly and above which I
>>> have really high pressures and the simulator aborts ??!
>>> In fact I want to impose a rate of 100 m^3/day in each well as it is a
>>> common liquid rate and that's why I do the calculations mentioned in the
>>> previous e-mail. but in order to implement what you proposed Timo in your
>>> response which is using "elemVolVars[scv].density(phaseIdx)" for the
>>> density and "scv.volume()" for the volume, shouldn't I divide the
>>> volumetric rate of the well which is 100 m^3/day by (3 cells x 8 scv)
>>> before dividing by the volume and multiplying by the density ? Anyway, I
>>> tried to follow what you proposed by writing the following for the
>>> injection well for example
>>> NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
>>> {
>>> SubControlVolume scv;
>>> ElementVolumeVariables elemVolVars;
>>> NumEqVector values;
>>> using FluidState = GetPropType<TypeTag, Properties::FluidState>;
>>> FluidState fs;
>>> // we define one source term and one sink term
>>> //Injection well (Source Term)
>>> if ((globalPos[2] > 15 - eps_ && globalPos[2] < 35 + eps_) &&
>>> (globalPos[1] > 45 - eps_ && globalPos[1] < 45 + eps_) && (globalPos[0] >
>>> 35 - eps_ && globalPos[0] < 35 + eps_))
>>> {
>>> fs.setTemperature(293);
>>> values[contiWEqIdx] =
>>> (1.157e-3*elemVolVars[scv].density(LiquidIdx))/(24*scv.volume()); //
>>> kg/(s*m^3)
>>> }
>>> where LiquidIdx is defined in an enum as LiquidIdx = FluidSystem::comp0Idx
>>> but I'm getting the following error:
>>> /home/mahmoud/dumux/dumux/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh:288:7:
>>> note: candidate expects 1 argument, 0 provided
>>> I would really appreciate it if you let me know how to correctly accurate
>>> for the 100 m^3/day volumetric rate either using the sourceAtPos function
>>> with elaborating the mistakes or using the pointSource function showing how
>>> to impose the 293 kelvin temperature for the injection well.
>>> I'm extremely sorry for the long e-mail. I just wanted to clarify my
>>> problem a bit more. Thanks in advance !
>>> Best Regards,
>>> Mahmoud Atef
>>> ------------------------------------------------------------------------------------------------
>>> Il 2021-04-07 20:27 Timo Koch ha scritto:
>>>> Dear Mahmoud,
>>>> what you are doing seems correct to me. You can get the relevant cell
>>>> volume with “scv.volume()” so you don’t have to compute that
>>>> yourself.
>>>> You can also get the density with
>>>> "elemVolVars[scv].density(phaseIdx)". But that shouldn’t change
>>>> anything if you already used the correct volume in your computation.
>>>> You can also use pointSources (see e.g
>>>> https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/test/porousmediumflow/1p/pointsources/timeindependent/problem.hh)
>>>> which is simpler for the case you describe.
>>>> There you just specify the position of the well(s) and the rate(s)
>>>> (this time directly in kg/s without the need for considering cell
>>>> volumes (that will be done automatically)).
>>>> I guess you have to recheck your scenario to see if you would expect
>>>> an influence on pressure for the specified production rate.
>>>> Best wishes
>>>> Timo
>>>>> On 6. Apr 2021, at 18:05, Mahmoud Atef Mahmoud Mohamed Aboelseoud
>>>>> S277151 <s277...@studenti.polito.it> wrote:
>>>>> Hello Dumux team,
>>>>> I'm defining two wells inside a structured 3D Yasp grid, one is
>>>>> injection and the other is production and I'm using the MPFA
>>>>> discretization. The system is a geothermal system for which I'm
>>>>> defining a certain low temperature for injection. I'm defining both
>>>>> wells inside the sourceAtPos method and they are both defined inside
>>>>> 3 cells in the z direction. each cell is 10x10x10 m^3 Now I
>>>>> understood from the doxygen documentation that the unit inside that
>>>>> source function is Kg/(m^3.s) and I want to impose a rate of 100
>>>>> m^3/day == 1.157e-3 m^3/s in each well so I'm dividing by (3 x 1000)
>>>>> and then multiplying by 1000 kg/m^3 for water to get an estimation
>>>>> of the mass rate to use so I get a value around 4e-4 Kg/(m^3.s) to
>>>>> use. However whenever I use that value, the simulation proceeds
>>>>> really quickly like I'm not injecting or producing anything at all
>>>>> (static condition) and whenever I use a random higher mass rate like
>>>>> 0.1 for example, I get very high pressures and warnings during the
>>>>> simulation and the simulator aborts. Could somebody please let me
>>>>> know how to correctly account for the 100 m^3/day that I want to
>>>>> impose in each of the production and injection wells ? tried using
>>>>> some middle values and the simulation proceeded successfully but I
>>>>> don't know the basis. Your guidance is much appreciated.
>>>>> Best Regards,
>>>>> Mahmoud Atef
>>>>> _______________________________________________
>>>>> DuMux mailing list
>>>>> DuMux@listserv.uni-stuttgart.de
>>>>> https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
>>>> _______________________________________________
>>>> DuMux mailing list
>>>> DuMux@listserv.uni-stuttgart.de
>>>> https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
>>> _______________________________________________
>>> DuMux mailing list
>>> DuMux@listserv.uni-stuttgart.de
>>> https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
>> _______________________________________________
>> DuMux mailing list
>> DuMux@listserv.uni-stuttgart.de
>> https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
> <high temperature in the production well.png><pressure boundaries
> modified.png>_______________________________________________
> DuMux mailing list
> DuMux@listserv.uni-stuttgart.de
> https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
_______________________________________________
DuMux mailing list
DuMux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
