diff --git a/src/coreComponents/constitutive/contact/BartonBandis.hpp b/src/coreComponents/constitutive/contact/BartonBandis.hpp index e8949781147..f6f4639cca6 100644 --- a/src/coreComponents/constitutive/contact/BartonBandis.hpp +++ b/src/coreComponents/constitutive/contact/BartonBandis.hpp @@ -140,10 +140,21 @@ real64 BartonBandisUpdates::computeHydraulicAperture( real64 const aperture, using namespace fields::contact; // Note: compressive - real64 const hydraulicAperture = ( fractureState == FractureState::Open ) ? (aperture + m_aperture0) : m_aperture0 / ( 1.0 - 9.0 * normalTraction /m_referenceNormalStress ); - dHydraulicAperture_dNormalTraction = - ( fractureState == FractureState::Open ) ? 0.0 : hydraulicAperture / ( 1.0 - 9.0 * normalTraction /m_referenceNormalStress ) * 9.0/m_referenceNormalStress; - dHydraulicAperture_aperture = ( fractureState == FractureState::Open ) ? 1.0 : 0.0; + real64 hydraulicAperture = -1e99; + if( fractureState == FractureState::Open ) + { + hydraulicAperture = (aperture + m_aperture0); + dHydraulicAperture_dNormalTraction = 0.0 ; + dHydraulicAperture_aperture = 1.0; + } + else + { + hydraulicAperture = m_aperture0 / ( 1.0 - 9.0 * normalTraction /m_referenceNormalStress ); + dHydraulicAperture_dNormalTraction = hydraulicAperture / ( 1.0 - 9.0 * normalTraction /m_referenceNormalStress ) * 9.0/m_referenceNormalStress; + dHydraulicAperture_aperture = 0.0; + } + + return hydraulicAperture; ///It would be nice to change this to return a tuple. } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index 3cba71973ca..ff4c5e8d754 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -312,6 +312,8 @@ void FlowSolverBase::setConstitutiveNamesCallSuper( ElementSubRegionBase & subRe setConstitutiveName< PermeabilityBase >( subRegion, viewKeyStruct::permeabilityNamesString(), "permeability" ); + setConstitutiveName< HydraulicApertureBase >( subRegion, viewKeyStruct::hydraulicApertureNamesString(), "hydraulicAperture" ); + if( m_isThermal ) { setConstitutiveName< SolidInternalEnergy >( subRegion, viewKeyStruct::solidInternalEnergyNamesString(), "solid internal energy" ); @@ -630,6 +632,26 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su arrayView1d< real64 const > const newHydraulicAperture = subRegion.getField< flow::hydraulicAperture >(); arrayView1d< real64 const > const oldHydraulicAperture = subRegion.getField< flow::aperture0 >(); + string const & hydraulicApertureName = subRegion.getReference< string >( viewKeyStruct::hydraulicApertureNamesString() ); + HydraulicApertureBase & hydraulicApertureModel = subRegion.getConstitutiveModel< HydraulicApertureBase >( hydraulicApertureName ); + + constitutive::ConstitutivePassThru< HydraulicApertureBase >::execute( hydraulicApertureModel, [=, &subRegion] ( auto & castedHydraulicApertureModel ) + { + typename TYPEOFREF( castedHydraulicApertureModel )::KernelWrapper hydraulicApertureWrapper = castedHydraulicApertureModel.createKernelUpdates(); + for( localIndex k = 0; k < subRegion.size(); ++k ) + { + real64 const normalStress = 0.0; + real64 dHydraulicAperture_aperture = 0.0; + real64 dHydraulicAperture_dNormalTraction = 0.0; + newHydraulicAperture[k] = hydraulicApertureWrapper.computeHydraulicAperture( oldHydraulicAperture[k], + normalStress, + 0, + dHydraulicAperture_aperture, + dHydraulicAperture_dNormalTraction ); + } + } ); + + string const & solidName = subRegion.getReference< string >( viewKeyStruct::solidNamesString() ); CoupledSolidBase & porousSolid = subRegion.getConstitutiveModel< CoupledSolidBase >( solidName ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index e67bae98e2e..c8c3e7c15c8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -82,6 +82,7 @@ class FlowSolverBase : public PhysicsSolverBase static constexpr char const * fluidNamesString() { return "fluidNames"; } static constexpr char const * solidNamesString() { return "solidNames"; } + static constexpr char const * hydraulicApertureNamesString() { return "hydraulicApertureNames"; } static constexpr char const * permeabilityNamesString() { return "permeabilityNames"; } static constexpr char const * solidInternalEnergyNamesString() { return "solidInternalEnergyNames"; } static constexpr char const * thermalConductivityNamesString() { return "thermalConductivityNames"; }