diff --git a/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_BBP_base.xml b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_BBP_base.xml new file mode 100644 index 00000000000..092d959197c --- /dev/null +++ b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_BBP_base.xml @@ -0,0 +1,158 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_BBP_smoke.xml b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_BBP_smoke.xml new file mode 100644 index 00000000000..198e1abdcd4 --- /dev/null +++ b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_BBP_smoke.xml @@ -0,0 +1,122 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 36dc61610ac..5f77946517e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -298,7 +298,6 @@ if (NOT DEFINED ENABLE_ATS) set( ENABLE_ATS true CACHE BOOL "") endif() - if ( ENABLE_ATS ) if (NOT DEFINED ATS_WORKING_DIR) message( WARNING "ATS_WORKING_DIR is not defined (required for integrated testing system)" ) diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index 6af439fe23f..4b9a2823fa1 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -152,6 +152,7 @@ set( constitutive_headers permeability/PermeabilityBase.hpp permeability/PermeabilityFields.hpp permeability/PressurePermeability.hpp + permeability/BartonBandisPermeability.hpp permeability/ProppantPermeability.hpp permeability/SlipDependentPermeability.hpp permeability/WillisRichardsPermeability.hpp @@ -301,6 +302,7 @@ set( constitutive_sources permeability/ParallelPlatesPermeability.cpp permeability/PermeabilityBase.cpp permeability/PressurePermeability.cpp + permeability/BartonBandisPermeability.cpp permeability/ProppantPermeability.cpp permeability/SlipDependentPermeability.cpp permeability/WillisRichardsPermeability.cpp diff --git a/src/coreComponents/constitutive/ConstitutivePassThru.hpp b/src/coreComponents/constitutive/ConstitutivePassThru.hpp index d78d7bb1bf8..ccad125f20d 100644 --- a/src/coreComponents/constitutive/ConstitutivePassThru.hpp +++ b/src/coreComponents/constitutive/ConstitutivePassThru.hpp @@ -49,6 +49,7 @@ #include "permeability/CarmanKozenyPermeability.hpp" #include "permeability/ExponentialDecayPermeability.hpp" #include "permeability/ParallelPlatesPermeability.hpp" +#include "permeability/BartonBandisPermeability.hpp" #include "permeability/PressurePermeability.hpp" #include "permeability/ProppantPermeability.hpp" #include "permeability/SlipDependentPermeability.hpp" @@ -375,6 +376,7 @@ struct ConstitutivePassThru< CompressibleSolidBase > CompressibleSolid< PressurePorosity, CarmanKozenyPermeability >, CompressibleSolid< PressurePorosity, ExponentialDecayPermeability >, CompressibleSolid< PressurePorosity, ParallelPlatesPermeability >, + CompressibleSolid< PressurePorosity, BartonBandisPermeability >, CompressibleSolid< PressurePorosity, PressurePermeability >, CompressibleSolid< PressurePorosity, SlipDependentPermeability >, CompressibleSolid< PressurePorosity, WillisRichardsPermeability > @@ -389,6 +391,7 @@ struct ConstitutivePassThru< CompressibleSolidBase > CompressibleSolid< PressurePorosity, CarmanKozenyPermeability >, CompressibleSolid< PressurePorosity, ExponentialDecayPermeability >, CompressibleSolid< PressurePorosity, ParallelPlatesPermeability >, + CompressibleSolid< PressurePorosity, BartonBandisPermeability >, CompressibleSolid< PressurePorosity, PressurePermeability >, CompressibleSolid< PressurePorosity, SlipDependentPermeability >, CompressibleSolid< PressurePorosity, WillisRichardsPermeability > @@ -462,6 +465,7 @@ struct ConstitutivePassThru< CoupledSolidBase > CompressibleSolid< PressurePorosity, CarmanKozenyPermeability >, CompressibleSolid< PressurePorosity, ExponentialDecayPermeability >, CompressibleSolid< PressurePorosity, ParallelPlatesPermeability >, + CompressibleSolid< PressurePorosity, BartonBandisPermeability >, CompressibleSolid< PressurePorosity, PressurePermeability >, CompressibleSolid< PressurePorosity, SlipDependentPermeability >, CompressibleSolid< PressurePorosity, WillisRichardsPermeability >, @@ -503,6 +507,7 @@ struct ConstitutivePassThru< CoupledSolidBase > CompressibleSolid< PressurePorosity, CarmanKozenyPermeability >, CompressibleSolid< PressurePorosity, ExponentialDecayPermeability >, CompressibleSolid< PressurePorosity, ParallelPlatesPermeability >, + CompressibleSolid< PressurePorosity, BartonBandisPermeability >, CompressibleSolid< PressurePorosity, PressurePermeability >, CompressibleSolid< PressurePorosity, SlipDependentPermeability >, CompressibleSolid< PressurePorosity, WillisRichardsPermeability >, diff --git a/src/coreComponents/constitutive/permeability/BartonBandisPermeability.cpp b/src/coreComponents/constitutive/permeability/BartonBandisPermeability.cpp new file mode 100644 index 00000000000..bc33c482200 --- /dev/null +++ b/src/coreComponents/constitutive/permeability/BartonBandisPermeability.cpp @@ -0,0 +1,102 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file BartonBandisPermeability.cpp + */ + +#include "BartonBandisPermeability.hpp" + +#include "constitutive/permeability/PermeabilityFields.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace constitutive +{ + + +BartonBandisPermeability::BartonBandisPermeability( string const & name, Group * const parent ): + PermeabilityBase( name, parent ), + m_updateTransversalComponent( true ) +{ + registerWrapper( viewKeyStruct::transversalPermeabilityString(), &m_transversalPermeability ). + setApplyDefaultValue( -1 ). + setInputFlag( InputFlags::OPTIONAL ). + setSizedFromParent( 0 ). + setDescription( "Default value of the permeability normal to the surface. If not specified the permeability is updated using the cubic law. " ); + + /// TODO: must become a required parameter. + registerWrapper( viewKeyStruct::apertureZeroString(), &m_aperture0 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setApplyDefaultValue( 1e-6 ). + setDescription( "Reference hydraulic aperture. It is the aperture at zero normal stress." ); + + registerWrapper( viewKeyStruct::biotString(), &m_biot ). + setApplyDefaultValue( 1.0 ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Biot coefficient." ); + registerWrapper( viewKeyStruct::poissonString(), &m_poisson ). + setApplyDefaultValue( 0.3 ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Poisson ratio." ); + registerWrapper( viewKeyStruct::normalStiffnessString(), &m_normalStiffness ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Normal stiffness: Kni." ); + registerWrapper( viewKeyStruct::referencePressureString(), &m_referencePressure ). + setApplyDefaultValue( 1.0e5 ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Reference pressure: p_0." ); + registerWrapper( viewKeyStruct::referenceTotalStressString(), &m_referenceTotalStress ). + setApplyDefaultValue( { 85.0e6, 85.0e6, 105.0e6 } ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Total stress at reference state: sigmaT_0." ); +} + + +void BartonBandisPermeability::postInputInitialization() +{ + PermeabilityBase::postInputInitialization(); + + if( m_transversalPermeability > -1 ) + { + m_updateTransversalComponent = false; + } +} + +void BartonBandisPermeability::initializeState() const +{ + localIndex const numE = m_permeability.size( 0 ); + integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point + + auto permView = m_permeability.toView(); + + real64 const transversalPerm = m_transversalPermeability; + + forAll< parallelDevicePolicy<> >( numE, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + for( localIndex q = 0; q < numQuad; ++q ) + { + permView[ei][q][2] = transversalPerm; + } + } ); +} + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, BartonBandisPermeability, string const &, Group * const ) + +} +} /* namespace geos */ diff --git a/src/coreComponents/constitutive/permeability/BartonBandisPermeability.hpp b/src/coreComponents/constitutive/permeability/BartonBandisPermeability.hpp new file mode 100644 index 00000000000..9818e310d58 --- /dev/null +++ b/src/coreComponents/constitutive/permeability/BartonBandisPermeability.hpp @@ -0,0 +1,323 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file BartonBandisPermeability.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_PERMEABILITY_BARTONBANDISPERMEABILITY_HPP_ +#define GEOS_CONSTITUTIVE_PERMEABILITY_BARTONBANDISPERMEABILITY_HPP_ + +#include "constitutive/permeability/PermeabilityBase.hpp" +#include "LvArray/src/genericTensorOps.hpp" + + +namespace geos +{ +namespace constitutive +{ + +class BartonBandisPermeabilityUpdate : public PermeabilityBaseUpdate +{ +public: + + BartonBandisPermeabilityUpdate( arrayView3d< real64 > const & permeability, + arrayView3d< real64 > const & dPerm_dPressure, + bool const updateTransversalComponent, + real64 const aperture0, + real64 const biot, + real64 const poisson, + real64 const normalStiffness, + real64 const referencePressure, + R1Tensor const &referenceTotalStress ) + : PermeabilityBaseUpdate( permeability, dPerm_dPressure ), + m_numDimensionsToUpdate( 3 ), + m_aperture0( aperture0 ), + m_biot( biot ), + m_poisson( poisson ), + m_normalStiffness( normalStiffness ), // Kni + m_referencePressure( referencePressure ), + m_referenceTotalStress( referenceTotalStress ) + { + m_numDimensionsToUpdate = updateTransversalComponent ? 3 : 2; + + m_referenceEffectiveStress[0] = m_referenceTotalStress[0] - m_biot*m_referencePressure; + m_referenceEffectiveStress[1] = m_referenceTotalStress[1] - m_biot*m_referencePressure; + m_referenceEffectiveStress[2] = m_referenceTotalStress[2] - m_biot*m_referencePressure; + } + + GEOS_HOST_DEVICE + void compute( real64 const & oldHydraulicAperture, + real64 const & newHydraulicAperture, + arraySlice1d< real64 > const & permeability, + real64 & dPerm_dHydraulicAperture ) const + { + GEOS_UNUSED_VAR( oldHydraulicAperture ); + + real64 const perm = newHydraulicAperture*newHydraulicAperture / 12.0; + dPerm_dHydraulicAperture = newHydraulicAperture / 6.0; + + for( int dim=0; dim < m_numDimensionsToUpdate; dim++ ) + { + permeability[dim] = perm; + } + } + + GEOS_HOST_DEVICE + virtual void updateFromPressureApertureAndNormal( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & oldHydraulicAperture, + real64 const & newHydraulicAperture, + arraySlice1d< real64 const > const & normal, + real64 const & dHydraulicAperture_dNormalJump ) const override final + { + GEOS_UNUSED_VAR( q, dHydraulicAperture_dNormalJump); + // compute effective normal stress on the fracture + real64 dStress_dPressure = -1.0; + real64 const fractureStress = computeFractureStress(pressure, normal, dStress_dPressure); + // compute new aperture using Barton Bandis model + real64 dAperture_dStress = -1.0; + real64 hydraulicAperture = computeHydraulicAperture(pressure, fractureStress, normal, dAperture_dStress, k); + + real64 dPerm_dHydraulicAperture = -1.0; + compute( oldHydraulicAperture, + hydraulicAperture, + m_permeability[k][0], + dPerm_dHydraulicAperture ); + + real64 const dPerm_dPressure = dPerm_dHydraulicAperture * dAperture_dStress * dStress_dPressure; + for( localIndex i=0; i < m_permeability[k][0].size(); i++ ) // size = 3 + { + m_dPerm_dPressure[k][0][i] = dPerm_dPressure; + } + } + + + + GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + real64 computeHydraulicAperture( real64 const pressure, real64 const normalComponentOfStressOnFracture, + arraySlice1d< real64 const > const & normal, real64 & dAperture_dStress, int k ) const + { + real64 const referenceTotalStress[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3 (m_referenceTotalStress); + real64 const biot_pressure = m_biot * m_referencePressure; // biot is alpha in the equations + + // Computation of maximum fracture closure (Barton-Bandis parameter) + // Fracture traction via Terzaghi's Principle + real64 sigma_c0[3] = {0.0}; + LvArray::tensorOps::hadamardProduct< 3 >( sigma_c0, referenceTotalStress, normal ); + LvArray::tensorOps::scaledAdd< 3 >(sigma_c0, normal, -biot_pressure); + + real64 const sigma_n0 = LvArray::tensorOps::AiBi< 3 >( sigma_c0, normal ); + + // \dfrac{-K_{ni}a_0 + \sqrt{(K_{ni}a_0)^2 + 4K_{ni}a_0\sigma_{n0}}}{2K_{ni}}. + real64 const normalStiffApertureProduct = m_normalStiffness * m_aperture0; + real64 const sqrtTerm = normalStiffApertureProduct * (normalStiffApertureProduct + 4.0*sigma_n0) ; + real64 const g0 = ( -normalStiffApertureProduct + std::sqrt(sqrtTerm) ) / (2.0 * m_normalStiffness); + + real64 const maximumFractureClosure = g0 + m_aperture0; // V_m or a_m -> aperture at stress-free state + + // Normal effective stress on the fracture + // g_n(\sigma_n) = \dfrac{\sigma_n * V_m}{K_{ni} * V_m + \sigma_n} + real64 const fractureClosure = (normalComponentOfStressOnFracture*maximumFractureClosure) / (m_normalStiffness*maximumFractureClosure + normalComponentOfStressOnFracture); + + real64 const newHydraulicAperture = maximumFractureClosure - fractureClosure; + + // derivative + // \frac{da}{d\sigma}(\sigma) = -\dfrac{K_{ni} V_m^2}{(K_{ni} V_m + \sigma)^2} + real64 const normalStiffMaxClosureProduct = m_normalStiffness*maximumFractureClosure; + real64 const denom = normalStiffMaxClosureProduct + normalComponentOfStressOnFracture; + dAperture_dStress = -(normalStiffMaxClosureProduct*maximumFractureClosure) / (denom * denom); + + return newHydraulicAperture; + } + /*GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + real64 computeHydraulicAperture( real64 const pressure, real64 const normalComponentOfStressOnFracture, + array1d< real64 > const & normal, real64 & dAperture_dStress, int k ) const + { + real64 const biot_pressure = m_biot * m_referencePressure; // biot is alpha in the equations + + // Computation of maximum fracture closure (Barton-Bandis parameter) + // Fracture traction via Terzaghi's Principle + real64 const sigma_c0[3] = { m_referenceTotalStress[0] * normal[0] - biot_pressure * normal[0], + m_referenceTotalStress[1] * normal[1] - biot_pressure * normal[1], + m_referenceTotalStress[2] * normal[2] - biot_pressure * normal[2] }; + real64 const sigma_n0 = sigma_c0[0]*normal[0] + + sigma_c0[1]*normal[1] + + sigma_c0[2]*normal[2]; + real64 const g0 = (-m_normalStiffness*m_aperture0 + + std::sqrt((m_normalStiffness*m_aperture0)* + (m_normalStiffness*m_aperture0) + + 4.0*m_normalStiffness*sigma_n0*m_aperture0)) / (2.0*m_normalStiffness); + real64 const maximumFractureClosure = g0 + m_aperture0; // Vm -> a_m -> aperture at free-stress state + + // Normal effective stress on the fracture + real64 const fractureClosure = normalComponentOfStressOnFracture*maximumFractureClosure/(m_normalStiffness*maximumFractureClosure + normalComponentOfStressOnFracture); // gn_BB + + // Compute the new aperture which is equal to the aperture at the free-stress state + // minus the closure from the free-stress state to the current state + real64 const newHydraulicAperture = maximumFractureClosure - fractureClosure; + + // derivative + real64 const Kni_apert = m_normalStiffness*maximumFractureClosure; + real64 const Kni_aper_stress = Kni_apert + normalComponentOfStressOnFracture; + dAperture_dStress = -(Kni_apert*maximumFractureClosure) / (Kni_aper_stress * Kni_aper_stress); + + return newHydraulicAperture; + }*/ + +private: + + int m_numDimensionsToUpdate; + + + real64 m_aperture0; + + real64 m_biot; + real64 m_poisson; + real64 m_normalStiffness; // Kni + real64 m_referencePressure; // p_0 + + R1Tensor m_referenceTotalStress; // sigmaT_0 computed analytically + R1Tensor m_referenceEffectiveStress; // sigma_0 + + GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + real64 computeFractureStress( real64 const pressure, arraySlice1d< real64 const > const & normal, real64 & dStress_dPressure ) const + { + //real64 const normal[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3 (normal_); + + real64 const deltaSigmaZ = m_biot * (pressure - m_referencePressure); + real64 const poisson_deltaSigma = deltaSigmaZ * m_poisson/(1.0 - m_poisson); + // sigma: matrix diagonal + real64 effectiveStress[3] = { m_referenceEffectiveStress[0] - poisson_deltaSigma, + m_referenceEffectiveStress[1] - poisson_deltaSigma, + m_referenceEffectiveStress[2] - deltaSigmaZ }; + real64 effectiveStressOnFracture[3] = {0.0}; // sigma_c + LvArray::tensorOps::hadamardProduct< 3 >( effectiveStressOnFracture, normal, effectiveStress ); + real64 const normalComponentOfStressOnFracture = LvArray::tensorOps::AiBi< 3 >(effectiveStressOnFracture, normal); // sigmaN_N + + // derivative + dStress_dPressure = -m_biot; + + return normalComponentOfStressOnFracture; + } + + + /*GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + real64 computeFractureStress( real64 const pressure, array1d< real64 > const & normal, real64 & dStress_dPressure ) const + { + // TODO: remove this lambda expression + auto matmul = [](real64 const (&u)[3], array1d< real64 > const &v, array1d< real64 > &r) -> void + { + r[0] = u[0]*v[0]; + r[1] = u[1]*v[1]; + r[2] = u[2]*v[2]; + }; + + real64 const deltaSigmaZ = m_biot * (pressure - m_referencePressure); + real64 const poisson_deltaSigma = deltaSigmaZ * m_poisson/(1.0 - m_poisson); + // sigma: matrix diagonal + real64 effectiveStress[3] = { m_referenceEffectiveStress[0] - poisson_deltaSigma, + m_referenceEffectiveStress[1] - poisson_deltaSigma, + m_referenceEffectiveStress[2] - deltaSigmaZ }; + array1d< real64 > effectiveStressOnFracture(3); // sigma_c + matmul(effectiveStress, normal, effectiveStressOnFracture); + real64 normalComponentOfStressOnFracture = dot(effectiveStressOnFracture, normal); // sigmaN_N + // derivative + dStress_dPressure = -m_biot; + + return normalComponentOfStressOnFracture; + }*/ + + +}; + + +class BartonBandisPermeability : public PermeabilityBase +{ +public: + + BartonBandisPermeability( string const & name, dataRepository::Group * const parent ); + + static string catalogName() { return "BartonBandisPermeability"; } + + virtual string getCatalogName() const override { return catalogName(); } + + virtual void initializeState() const override final; + + /// Type of kernel wrapper for in-kernel update + using KernelWrapper = BartonBandisPermeabilityUpdate; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelWrapper() const + { + return KernelWrapper( m_permeability, + m_dPerm_dPressure, + m_updateTransversalComponent, + m_aperture0, + m_biot, + m_poisson, + m_normalStiffness, + m_referencePressure, + m_referenceTotalStress); + } + + + struct viewKeyStruct : public PermeabilityBase::viewKeyStruct + { + static constexpr char const * transversalPermeabilityString() { return "transversalPermeability"; } + + /// string/key for aperture under zero normal stress + static constexpr char const * apertureZeroString() { return "referenceAperture"; } + static constexpr char const * biotString() { return "biot"; } + static constexpr char const * poissonString() { return "poisson"; } + static constexpr char const * normalStiffnessString() { return "normalStiffness"; } + static constexpr char const * referencePressureString() { return "referencePressure"; } + static constexpr char const * referenceTotalStressString() { return "referenceTotalStress"; } + }; + +protected: + + virtual void postInputInitialization() override; + +private: + + real64 m_transversalPermeability; + + bool m_updateTransversalComponent; + + /// Reference hydraulic aperture. Aperture at zero normal stress + real64 m_aperture0; /// TODO: this will replace what is currently called defaultAperture. + real64 m_biot; + real64 m_poisson; + real64 m_normalStiffness; // Kni + real64 m_referencePressure; // p_0 + + R1Tensor m_referenceTotalStress; // sigmaT_0 computed analytically +}; + +}/* namespace constitutive */ + +} /* namespace geos */ + + +#endif //GEOS_CONSTITUTIVE_PERMEABILITY_BARTONBANDISPERMEABILITY_HPP_ diff --git a/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.hpp b/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.hpp index 80b5901714b..63ad58c1c9c 100644 --- a/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.hpp +++ b/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.hpp @@ -60,26 +60,30 @@ class ParallelPlatesPermeabilityUpdate : public PermeabilityBaseUpdate permeability[dim] = perm; dPerm_dDispJump[dim][0] = dPerm_dHydraulicAperture * dHydraulicAperture_dNormalJump; dPerm_dDispJump[dim][1] = 0.0; - dPerm_dDispJump[dim][2] = 0.0; + dPerm_dDispJump[dim][2] = 0.0; } } GEOS_HOST_DEVICE - virtual void updateFromAperture( localIndex const k, - localIndex const q, - real64 const & oldHydraulicAperture, - real64 const & newHydraulicAperture, - real64 const & dHydraulicAperture_dNormalJump ) const override final + virtual void updateFromPressureApertureAndNormal( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & oldHydraulicAperture, + real64 const & newHydraulicAperture, + arraySlice1d< real64 const > const & normal, + real64 const & dHydraulicAperture_dNormalJump ) const override final { - GEOS_UNUSED_VAR( q ); + GEOS_UNUSED_VAR( q, pressure, normal ); compute( oldHydraulicAperture, newHydraulicAperture, dHydraulicAperture_dNormalJump, m_permeability[k][0], m_dPerm_dDispJump[k][0] ); + } + GEOS_HOST_DEVICE virtual void updateFromApertureAndShearDisplacement( localIndex const k, localIndex const q, @@ -91,8 +95,8 @@ class ParallelPlatesPermeabilityUpdate : public PermeabilityBaseUpdate real64 const ( &traction )[3] ) const override final { GEOS_UNUSED_VAR( dispJump, traction, pressure ); - - updateFromAperture( k, q, oldHydraulicAperture, newHydraulicAperture, dHydraulicAperture_dNormalJump ); + array1d< real64 > const normal; + updateFromPressureApertureAndNormal( k, q, pressure, oldHydraulicAperture, newHydraulicAperture, normal, dHydraulicAperture_dNormalJump ); } private: diff --git a/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp b/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp index 190b8f1b4f6..388564ba865 100644 --- a/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp +++ b/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp @@ -56,14 +56,27 @@ class PermeabilityBaseUpdate GEOS_UNUSED_VAR( k, q, pressure, porosity ); } - GEOS_HOST_DEVICE - virtual void updateFromAperture( localIndex const k, + //@jafranc + //for consistency maybe we want to add a new one but keep the previous + /*virtual void updateFromAperture( localIndex const k, localIndex const q, real64 const & oldHydraulicAperture, real64 const & newHydraulicAperture, real64 const & dHydraulicAperture_dNormalJump ) const { GEOS_UNUSED_VAR( k, q, oldHydraulicAperture, newHydraulicAperture, dHydraulicAperture_dNormalJump ); + }*/ + + GEOS_HOST_DEVICE + virtual void updateFromPressureApertureAndNormal( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & oldHydraulicAperture, + real64 const & newHydraulicAperture, + arraySlice1d< real64 const > const & normal, + real64 const & dHydraulicAperture_dNormalJump ) const + { + GEOS_UNUSED_VAR( k, q, pressure, oldHydraulicAperture, newHydraulicAperture, normal, dHydraulicAperture_dNormalJump ); } GEOS_HOST_DEVICE diff --git a/src/coreComponents/constitutive/permeability/PressurePermeability.hpp b/src/coreComponents/constitutive/permeability/PressurePermeability.hpp index 9d654048a0e..6e4ed374512 100644 --- a/src/coreComponents/constitutive/permeability/PressurePermeability.hpp +++ b/src/coreComponents/constitutive/permeability/PressurePermeability.hpp @@ -87,7 +87,7 @@ class PressurePermeabilityUpdate : public PermeabilityBaseUpdate referencePermeability[0] = m_referencePermeability[k][0][0]; referencePermeability[1] = m_referencePermeability[k][0][1]; referencePermeability[2] = m_referencePermeability[k][0][2]; - + switch( m_presModelType ) { case PressureModelType::Exponential: diff --git a/src/coreComponents/constitutive/solid/CompressibleSolid.cpp b/src/coreComponents/constitutive/solid/CompressibleSolid.cpp index 5519d77543a..04791b16c44 100644 --- a/src/coreComponents/constitutive/solid/CompressibleSolid.cpp +++ b/src/coreComponents/constitutive/solid/CompressibleSolid.cpp @@ -24,6 +24,7 @@ #include "constitutive/permeability/CarmanKozenyPermeability.hpp" #include "constitutive/permeability/ExponentialDecayPermeability.hpp" #include "constitutive/permeability/ParallelPlatesPermeability.hpp" +#include "constitutive/permeability/BartonBandisPermeability.hpp" #include "constitutive/permeability/PressurePermeability.hpp" #include "constitutive/permeability/SlipDependentPermeability.hpp" #include "constitutive/permeability/WillisRichardsPermeability.hpp" @@ -48,6 +49,7 @@ typedef CompressibleSolid< PressurePorosity, CarmanKozenyPermeability > Compress typedef CompressibleSolid< PressurePorosity, PressurePermeability > CompressibleRockPressurePerm; typedef CompressibleSolid< PressurePorosity, ExponentialDecayPermeability > FaultED; typedef CompressibleSolid< PressurePorosity, ParallelPlatesPermeability > FractureRock; +typedef CompressibleSolid< PressurePorosity, BartonBandisPermeability > FractureRockBB; typedef CompressibleSolid< PressurePorosity, SlipDependentPermeability > Fault; typedef CompressibleSolid< PressurePorosity, WillisRichardsPermeability > FaultWR; @@ -55,6 +57,7 @@ REGISTER_CATALOG_ENTRY( ConstitutiveBase, CompressibleRockConstant, string const REGISTER_CATALOG_ENTRY( ConstitutiveBase, CompressibleRockCK, string const &, Group * const ) REGISTER_CATALOG_ENTRY( ConstitutiveBase, CompressibleRockPressurePerm, string const &, Group * const ) REGISTER_CATALOG_ENTRY( ConstitutiveBase, FractureRock, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, FractureRockBB, string const &, Group * const ) REGISTER_CATALOG_ENTRY( ConstitutiveBase, FaultED, string const &, Group * const ) REGISTER_CATALOG_ENTRY( ConstitutiveBase, Fault, string const &, Group * const ) REGISTER_CATALOG_ENTRY( ConstitutiveBase, FaultWR, string const &, Group * const ) diff --git a/src/coreComponents/constitutive/solid/CompressibleSolid.hpp b/src/coreComponents/constitutive/solid/CompressibleSolid.hpp index fb1832b2804..33e1c186939 100644 --- a/src/coreComponents/constitutive/solid/CompressibleSolid.hpp +++ b/src/coreComponents/constitutive/solid/CompressibleSolid.hpp @@ -63,16 +63,17 @@ class CompressibleSolidUpdates : public CoupledSolidUpdates< NullModel, PORO_TYP } GEOS_HOST_DEVICE - void updateStateFromPressureAndAperture( localIndex const k, + void updateStateFromPressureApertureAndNormal( localIndex const k, localIndex const q, real64 const & pressure, real64 const & oldHydraulicAperture, - real64 const & newHydraulicAperture ) const + real64 const & newHydraulicAperture, + arraySlice1d< real64 const > const & normal ) const { real64 const temperature = 0; real64 const dHydraulicAperture_dNormalJump = 1.0; m_porosityUpdate.updateFromPressureAndTemperature( k, q, pressure, temperature ); - m_permUpdate.updateFromAperture( k, q, oldHydraulicAperture, newHydraulicAperture, dHydraulicAperture_dNormalJump ); + m_permUpdate.updateFromPressureApertureAndNormal( k, q, pressure, oldHydraulicAperture, newHydraulicAperture, normal, dHydraulicAperture_dNormalJump); } GEOS_HOST_DEVICE diff --git a/src/coreComponents/constitutive/unitTests/CMakeLists.txt b/src/coreComponents/constitutive/unitTests/CMakeLists.txt index 9dfbfc2c8a3..84edcf931fc 100644 --- a/src/coreComponents/constitutive/unitTests/CMakeLists.txt +++ b/src/coreComponents/constitutive/unitTests/CMakeLists.txt @@ -30,6 +30,7 @@ set( gtest_geosx_tests testNegativeTwoPhaseFlash9Comp.cpp testNegativeTwoPhaseFlashModel.cpp testParticleFluidEnums.cpp + testPrescribedStressPath.cpp testPressureTemperatureCoordinates.cpp testPropertyConversions.cpp testSoreideWhitsonEOSPhaseModel.cpp diff --git a/src/coreComponents/constitutive/unitTests/testPrescribedStressPath.cpp b/src/coreComponents/constitutive/unitTests/testPrescribedStressPath.cpp new file mode 100644 index 00000000000..bfe168f9f4c --- /dev/null +++ b/src/coreComponents/constitutive/unitTests/testPrescribedStressPath.cpp @@ -0,0 +1,87 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +// Source includes +#include "constitutive/ConstitutiveManager.hpp" +#include "constitutive/permeability/BartonBandisPermeability.hpp" +#include "constitutive/permeability/PermeabilityFields.hpp" +#include "dataRepository/xmlWrapper.hpp" + +// TPL includes +#include +#include + +using namespace geos; +using namespace ::geos::constitutive; + + +TEST( BartonBandisPermeabilityTests, testNewAperture ) +{ + conduit::Node node; + dataRepository::Group rootGroup( "root", node ); + ConstitutiveManager constitutiveManager( "constitutive", &rootGroup ); + + real64 constexpr referenceAperture = 1.2421e-4; // in-situ + real64 constexpr referencePressure = 1.0e5; + real64 constexpr fractureStress = 1.049e+08; // only if Z component of normal is 1.0 and referenceTotalStress in Z is 105.0e6 + std::stringstream ss_aperture; + ss_aperture << std::scientific << std::setprecision(4) << referenceAperture; + std::stringstream ss_pressure; + ss_pressure << std::scientific << std::setprecision(4) << referencePressure; + + string const inputStream = + "" + " " + ""; + + xmlWrapper::xmlDocument xmlDocument; + xmlWrapper::xmlResult xmlResult = xmlDocument.loadString( inputStream ); + if( !xmlResult ) + { + GEOS_LOG_RANK_0( "XML parsed with errors!" ); + GEOS_LOG_RANK_0( "Error description: " << xmlResult.description()); + GEOS_LOG_RANK_0( "Error offset: " << xmlResult.offset ); + } + + xmlWrapper::xmlNode xmlConstitutiveNode = xmlDocument.getChild( "Constitutive" ); + constitutiveManager.processInputFileRecursive( xmlDocument, xmlConstitutiveNode ); + constitutiveManager.postInputInitializationRecursive(); + + BartonBandisPermeability & cm = constitutiveManager.getConstitutiveRelation< BartonBandisPermeability >( "fracturePerm" ); + + BartonBandisPermeability::KernelWrapper cmw = cm.createKernelWrapper(); + + { + array1d < real64 > normal(3); + normal[0] = 0.0; + normal[1] = 0.0; + normal[2] = 1.0; + //real64 const normal[ 3 ] = {0.0, 0.0, 1.0}; + + + real64 dAperture_dStress = -1.0; + real64 newHydraulicAperture = cmw.computeHydraulicAperture(referencePressure, fractureStress, normal, dAperture_dStress, 0); + + EXPECT_DOUBLE_EQ( newHydraulicAperture, referenceAperture ); + } +} + diff --git a/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp b/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp index 5de67fe2f11..bf37eec16a4 100644 --- a/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp +++ b/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp @@ -333,7 +333,7 @@ SurfaceElementStencilWrapper:: localIndex const er1 = m_elementRegionIndices[iconn][k[1]]; localIndex const esr1 = m_elementSubRegionIndices[iconn][k[1]]; localIndex const ei1 = m_elementIndices[iconn][k[1]]; - + real64 const t0 = m_weights[iconn][k[0]] * coefficient[er0][esr0][ei0][0][0]; // this is a bit insane to access perm real64 const t1 = m_weights[iconn][k[1]] * coefficient[er1][esr1][ei1][0][0]; @@ -341,10 +341,10 @@ SurfaceElementStencilWrapper:: real64 const arithmeticWeight = 0.25 * (t0+t1); real64 const value = m_meanPermCoefficient * harmonicWeight + (1 - m_meanPermCoefficient) * arithmeticWeight; - + weight[connectionIndex][0] = value; weight[connectionIndex][1] = -value; - + real64 const dt0 = m_weights[iconn][k[0]] * dCoeff_dVar[er0][esr0][ei0][0][0]; real64 const dt1 = m_weights[iconn][k[1]] * dCoeff_dVar[er1][esr1][ei1][0][0]; @@ -355,7 +355,7 @@ SurfaceElementStencilWrapper:: real64 dArithmetic[2]; dArithmetic[0] = 0.25 * dt0; dArithmetic[1] = 0.25 * dt1; - + dWeight_dVar[connectionIndex][0] = m_meanPermCoefficient * dHarmonic[0] + (1 - m_meanPermCoefficient) * dArithmetic[0]; dWeight_dVar[connectionIndex][1] = -( m_meanPermCoefficient * dHarmonic[1] + (1 - m_meanPermCoefficient) * dArithmetic[1] ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index bf63134e49d..6873848aabf 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -36,6 +36,7 @@ #include "physicsSolvers/fluidFlow/kernels/MinPoreVolumeMaxPorosityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/StencilWeightsUpdateKernel.hpp" + namespace geos { @@ -85,21 +86,32 @@ void updatePorosityAndPermeabilityFixedStress( POROUSWRAPPER_TYPE porousWrapper, } ); } +void nada(arraySlice1d< real64 const > const & normal) +{ + GEOS_LOG_RANK_0("nada normal "<( sigma_c0, referenceTotalStress, normal ); +} + template< typename POROUSWRAPPER_TYPE > -void updatePorosityAndPermeabilityFromPressureAndAperture( POROUSWRAPPER_TYPE porousWrapper, +void updatePorosityAndPermeabilityFromPressureApertureAndNormal( POROUSWRAPPER_TYPE porousWrapper, SurfaceElementSubRegion & subRegion, arrayView1d< real64 const > const & pressure, arrayView1d< real64 const > const & oldHydraulicAperture, - arrayView1d< real64 const > const & newHydraulicAperture ) + arrayView1d< real64 const > const & newHydraulicAperture, + arrayView2d< real64 const > const & normalVector) { forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_DEVICE ( localIndex const k ) { for( localIndex q = 0; q < porousWrapper.numGauss(); ++q ) { - porousWrapper.updateStateFromPressureAndAperture( k, q, - pressure[k], - oldHydraulicAperture[k], - newHydraulicAperture[k] ); + porousWrapper.updateStateFromPressureApertureAndNormal( k, q, + pressure[k], + oldHydraulicAperture[k], + newHydraulicAperture[k], + normalVector[k] ); } } ); } @@ -118,7 +130,7 @@ FlowSolverBase::FlowSolverBase( string const & name, setApplyDefaultValue( 0 ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "Flag indicating whether the problem is thermal or not." ); - + this->registerWrapper( viewKeyStruct::allowNegativePressureString(), &m_allowNegativePressure ). setApplyDefaultValue( 0 ). // negative pressure is not allowed by default setInputFlag( InputFlags::OPTIONAL ). @@ -639,7 +651,7 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su GEOS_MARK_FUNCTION; arrayView1d< real64 const > const & pressure = subRegion.getField< flow::pressure >(); - + arrayView2d< real64 const > const & normalVector = subRegion.getField< fields::normalVector >(); // mesh/MeshFields.hpp arrayView1d< real64 const > const newHydraulicAperture = subRegion.getField< flow::hydraulicAperture >(); arrayView1d< real64 const > const oldHydraulicAperture = subRegion.getField< flow::aperture0 >(); @@ -649,7 +661,8 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su constitutive::ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) { typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousWrapper = castedPorousSolid.createKernelUpdates(); - updatePorosityAndPermeabilityFromPressureAndAperture( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture ); + + updatePorosityAndPermeabilityFromPressureApertureAndNormal( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture, normalVector ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 0e709372102..06e80977102 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -82,9 +82,13 @@ class FlowSolverBase : public PhysicsSolverBase static constexpr char const * fluidNamesString() { return "fluidNames"; } static constexpr char const * solidNamesString() { return "solidNames"; } + static constexpr char const * permeabilityNamesString() { return "permeabilityNames"; } static constexpr char const * solidInternalEnergyNamesString() { return "solidInternalEnergyNames"; } static constexpr char const * thermalConductivityNamesString() { return "thermalConductivityNames"; } + + static constexpr char const * computePrescribedStressPathString() { return "computePrescribedStressPath"; } + static constexpr char const * hydraulicApertureRelationNameString() { return "hydraulicApertureRelationName"; } }; /** diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 7c1521128b7..8b08e19e195 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -49,7 +49,6 @@ #include "physicsSolvers/fluidFlow/kernels/singlePhase/FluidUpdateKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolidInternalEnergyUpdateKernel.hpp" - namespace geos { @@ -133,6 +132,7 @@ void SinglePhaseBase::setConstitutiveNames( ElementSubRegionBase & subRegion ) c { setConstitutiveName< SinglePhaseThermalConductivityBase >( subRegion, viewKeyStruct::thermalConductivityNamesString(), "singlephase thermal conductivity" ); } + } void SinglePhaseBase::initializeAquiferBC() const @@ -674,7 +674,7 @@ void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_ applyDeltaVolume( subRegion ); - // This should fix NaN density in newly created fracture elements + // This should fix NaN density in newly created fracture elements updatePorosityAndPermeability( subRegion ); updateFluidState( subRegion ); // for thermal simulations, update solid internal energy @@ -1233,7 +1233,7 @@ void SinglePhaseBase::keepVariablesConstantDuringInitStep( real64 const time, void SinglePhaseBase::updateState( DomainPartition & domain ) { GEOS_MARK_FUNCTION; - + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, string_array const & regionNames ) diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 9a9bb51d843..63180f07fdc 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -222,18 +222,10 @@ - - - - - - - - @@ -667,6 +659,10 @@ + + + + @@ -747,6 +743,10 @@ + + + + @@ -1624,26 +1624,10 @@ stress - traction is applied to the faces as specified by the inner product of i - - - - - - - - - - - - - - - - @@ -1652,18 +1636,6 @@ stress - traction is applied to the faces as specified by the inner product of i - - - - - - - - - - - - @@ -6760,6 +6732,7 @@ Information output from lower logLevels is added with the desired log level + @@ -6780,6 +6753,7 @@ Information output from lower logLevels is added with the desired log level + @@ -6877,6 +6851,24 @@ Information output from lower logLevels is added with the desired log level + + + + + + + + + + + + + + + + + + @@ -7481,6 +7473,18 @@ if the table is requested to be output in the log, and it is too large, a CSV fi + + + + + + + + + + + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index c7434b7b99c..a2b5cee581d 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -351,15 +351,11 @@ A field can represent a physical variable. (pressure, temperature, global compos - - - - @@ -526,7 +522,7 @@ A field can represent a physical variable. (pressure, temperature, global compos - + @@ -1603,13 +1599,14 @@ A field can represent a physical variable. (pressure, temperature, global compos - + + @@ -1630,6 +1627,7 @@ A field can represent a physical variable. (pressure, temperature, global compos + @@ -1720,6 +1718,12 @@ A field can represent a physical variable. (pressure, temperature, global compos + + + + + + @@ -2462,6 +2466,7 @@ A field can represent a physical variable. (pressure, temperature, global compos + diff --git a/src/docs/sphinx/basicExamples/Index.rst b/src/docs/sphinx/basicExamples/Index.rst index af43ba213fe..0d00922eaaf 100644 --- a/src/docs/sphinx/basicExamples/Index.rst +++ b/src/docs/sphinx/basicExamples/Index.rst @@ -19,6 +19,8 @@ Basic Examples triaxialDriver/Example + stressPathDriven/Example + diff --git a/src/docs/sphinx/basicExamples/stressPathDriven/Example.rst b/src/docs/sphinx/basicExamples/stressPathDriven/Example.rst new file mode 100644 index 00000000000..ea515412b52 --- /dev/null +++ b/src/docs/sphinx/basicExamples/stressPathDriven/Example.rst @@ -0,0 +1,161 @@ +.. _stressPathDrivenExample: + + +############################################################# +Stress-Path Driven Simulation +############################################################# + + +**Context** + +In this example, we demonstrate how to set up a single-phase flow simulation driven by the oedometric stress path. The domain contains a single fracture, whose aperture is updated according to the Barton–Bandis model. + +**Objectives** + +At the end of this example, you will know: + + - how to set up the Barton-Bandis parameters for the oedometric stress path simulation + - how to update the fracture aperture without a mechanics solver + + +**Input file** + +The XML file for this test case is located at: + +.. code-block:: console + + inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_BBP_smoke.xml + + +------------------------------------------------------------------ +Description of the case +------------------------------------------------------------------ + +This example is part of the research about integrating geomechanical effects into the Embedded Discrete Fracture Model (EDFM). Here, the EDFM formulation is coupled with analytical poromechanics under an Oedometric Stress Path assumption (flat reservoir with low contrast between it and adjacent formation stiffness), enabling a direct relationship between pore-pressure variations and changes in effective normal stress. These stresses, in turn, govern aperture evolution via the Barton–Bandis constitutive law. Also, we assume: + + - In-situ stress data are available + - Andersonian stress regime (no lateral stress, constant vertical total stress) + - Permeability computed by Parallel Plates model + +The simulation presented below is a simplified case designed to validate the implementation against the analytical solution detailed in the Section ..... The simplification include the use of a single-phase flow solver with TPFA as discretization method, an impermeable matrix, and a through-going fracture. + +------------------------------ +Mesh and geometry +------------------------------ + +The hexahedral mesh is generated internally and consists of 11 cells aligned along the Y-axis, each measuring 1x1x1 m. The fracture is represented as a plane that crosses the entire domain, dividing it into two halves. A pressure gradient is applied at the boundary cells of the domain and directly at the fracture tips. + +.. figure:: stressPathDrivenExampleParaview.png + :align: center + :width: 500 + :figclass: align-center + + Visualization of the first simulation step. + +------------------------------------------------------------------ +Initial Conditions (Input file setup) +------------------------------------------------------------------ + +This example needs the following input values to work properly: + +- Surface Element Region's default aperture and Stress-path driven reference aperture must be the same. +- Initial pressure in the rock and fracture must be the same. +- Stress-path driven's reference pressure must equal the initial pressure. + + +------------------------------ +Constitutive law +------------------------------ + +The new class ``BartonBandisPermeability`` follows a similar implementation to the class ``ParallelPlatesPermeability``. Both are defined as permeability model managed by ``CompressibleSolid`` class, and each contains a private kernel responsible for computing the updated hydraulic aperture. However, because the proposed formulation of the Barton–Bandis model is stress-dependent, the implementation of the new class’s kernel differs and requires the following inputs: + +.. literalinclude:: ../../../../../inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_BBP_base.xml + :language: xml + :start-after: + :end-before: + +In summary, by assuming the oedometric conditions, it is possible to relate pore pressure and total stresses by + +.. math:: + \Delta \sigma_V = \alpha\, \Delta p, \qquad \Delta \sigma_H = \Delta \sigma_h = \frac{\nu}{1 - \nu}\, \Delta \sigma_V + +where :math:`\alpha` is the Biot coefficient and :math:`\nu` is Poisson's ratio. + +Given the in-situ effective stress, + +.. math:: + \sigma_0 = \sigma_{\text{ref}} - \alpha\, p_0\, I, + +where :math:`p_0` is the reference pore pressure, we compute the normal stress acting on the fracture plane by projecting the updated total stress tensor onto the fracture normal direction: + +.. math:: + \sigma_n = \left( \sigma_0 - \sigma_T \right) : (\vec{n} \otimes \vec{n}), + +where :math:`\vec{n}` is the unit normal vector to the fracture plane. + +Once the normal effective stress is obtained, the fracture closure is computed using the Barton--Bandis hyperbolic law, which relates normal stress to updated mechanical aperture. The fracture closure is given by + +.. math:: + g_n(\sigma_n) = \frac{\sigma_n\, V_m}{K_{ni}\, V_m + \sigma_n}, + +where :math:`K_{ni}` is the initial normal stiffness and :math:`V_m` is the maximum possible closure. +The resulting closure is then used to update the aperture (:math:`a`), from which the intrinsic fracture permeability (:math:`K_f`) is obtained using the parallel-plates law + +.. math:: + K_f = \frac{a^2}{12}. + +Both quantities — the updated aperture and the resulting permeability — affect the EDFM transmissibilities. + +------------------------------------------------------------------ +Single-phase flow solver with fracture aperture update +------------------------------------------------------------------ + +The single-phase flow solver definition in the input file remains the same, with the exception of one optional flag. The flag ``computePrescribedStressPath`` enables the computation of the updated fracture aperture using the new ``BartonBandisPermeability`` class. If this constitutive law is not defined in the input file while the flag is set to 1, an error will be raised. On the other hand, if this flag is set to 0, the fracture aperture will not be updated, even if the constitutive law is defined. + +.. literalinclude:: ../../../../../inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_BBP_base.xml + :language: xml + :start-after: + :end-before: + +The new class integration is illustrated in the diagram below. Since the proposed permeability class needs the current pressure and fracture normal vector for the effective stress computation, the function's signature in ``CompressibleSolid`` and in ``PermeabilityBase`` class were modified to reflect this requirement. + + +.. figure:: bartonBandisPermeabilityDiagram.png + :align: center + :width: 800 + :figclass: align-center + + Diagram of the new implementation to incorporate Barton-Bandis model to flow solvers. + +------------------------ +Events +------------------------ + +The code below is to save the pressure, aperture and permeability of the fracture at each time step for later validation against the analytical solution. + +.. literalinclude:: ../../../../../inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_BBP_smoke.xml + :language: xml + :start-after: + :end-before: + +------------------------------ +Unit test +------------------------------ + +The normal effective stress and aperture update computation can be tested by providing the in-situ (reference) values of pressure and aperture. The resulting aperture should equal the in-situ aperture. + + +------------------------------ +Integration test +------------------------------ + +... + +------------------------------------------------------------------ +To go further +------------------------------------------------------------------ + + +**Feedback on this example** + +For any feedback on this example, please submit a `GitHub issue on the project's GitHub page `_. diff --git a/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisPermeabilityDiagram.png b/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisPermeabilityDiagram.png new file mode 100644 index 00000000000..fb4d15b279e Binary files /dev/null and b/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisPermeabilityDiagram.png differ diff --git a/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisSPDDiagram.png b/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisSPDDiagram.png new file mode 100644 index 00000000000..435eb647117 Binary files /dev/null and b/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisSPDDiagram.png differ diff --git a/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaview.png b/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaview.png new file mode 100644 index 00000000000..5da55276b52 Binary files /dev/null and b/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaview.png differ diff --git a/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaviewShortFrac.png b/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaviewShortFrac.png new file mode 100644 index 00000000000..45037257521 Binary files /dev/null and b/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaviewShortFrac.png differ