fix: enable BartonBandis model for contact solver #3556

jhuang2601 · 2025-02-18T20:03:28Z

Following #3197, this PR

enable BartonBandis model for contact solver and other fracture models
add smoke tests for this feature

frankfeifan · 2025-02-20T17:45:24Z

src/coreComponents/constitutive/contact/BartonBandis.hpp

-  real64 const hydraulicAperture = ( aperture >= 0.0 ) ? (aperture + m_aperture0) : m_aperture0 / ( 1 + 9*normalTraction/m_referenceNormalStress );
-  dHydraulicAperture_dNormalStress = ( aperture >= 0.0 ) ? 0.0 : -hydraulicAperture / ( 1 + 9*normalTraction/m_referenceNormalStress ) * 9/m_referenceNormalStress;
-  dHydraulicAperture_aperture = ( aperture >= 0.0 ) ? 1.0 : 0.0;
+  real64 const hydraulicAperture = ( fractureState >= 3.0 ) ? (aperture + m_aperture0) : m_aperture0 / ( 1 + 9*normalTraction/m_referenceNormalStress );


We may want to be cautious about the convention of the normalTraction here. Here normalTraction is defined as positive for compression.

Yes, we need to ensure that the sign convention for normalTraction and referenceNormalStress is consistent.

CusiniM · 2025-02-20T22:46:04Z

src/coreComponents/constitutive/contact/BartonBandis.hpp

-  dHydraulicAperture_dNormalStress = ( aperture >= 0.0 ) ? 0.0 : -hydraulicAperture / ( 1 + 9*normalTraction/m_referenceNormalStress ) * 9/m_referenceNormalStress;
-  dHydraulicAperture_aperture = ( aperture >= 0.0 ) ? 1.0 : 0.0;
+  real64 const hydraulicAperture = ( fractureState >= 3.0 ) ? (aperture + m_aperture0) : m_aperture0 / ( 1 + 9*LvArray::math::abs(normalTraction)/m_referenceNormalStress );
+  dHydraulicAperture_dNormalStress = ( fractureState >= 3.0 ) ? 0.0 : -hydraulicAperture / ( 1 + 9*LvArray::math::abs(normalTraction)/m_referenceNormalStress ) * 9/m_referenceNormalStress;


Suggested change

dHydraulicAperture_dNormalStress = ( fractureState >= 3.0 ) ? 0.0 : -hydraulicAperture / ( 1 + 9*LvArray::math::abs(normalTraction)/m_referenceNormalStress ) * 9/m_referenceNormalStress;

dHydraulicAperture_dNormalStress = ( fractureState == FractureState::Open ) ? 0.0 : -hydraulicAperture / ( 1 + 9*LvArray::math::abs(normalTraction)/m_referenceNormalStress ) * 9/m_referenceNormalStress;

CusiniM · 2025-02-20T22:46:22Z

src/coreComponents/constitutive/contact/BartonBandis.hpp

-  dHydraulicAperture_aperture = ( aperture >= 0.0 ) ? 1.0 : 0.0;
+  real64 const hydraulicAperture = ( fractureState >= 3.0 ) ? (aperture + m_aperture0) : m_aperture0 / ( 1 + 9*LvArray::math::abs(normalTraction)/m_referenceNormalStress );
+  dHydraulicAperture_dNormalStress = ( fractureState >= 3.0 ) ? 0.0 : -hydraulicAperture / ( 1 + 9*LvArray::math::abs(normalTraction)/m_referenceNormalStress ) * 9/m_referenceNormalStress;
+  dHydraulicAperture_aperture = ( fractureState >= 3.0 ) ? 1.0 : 0.0;


Suggested change

dHydraulicAperture_aperture = ( fractureState >= 3.0 ) ? 1.0 : 0.0;

dHydraulicAperture_aperture = ( fractureState >= fractureState == FractureState::Open ) ? 1.0 : 0.0;

you may just need to include FractureState.hpp

CusiniM · 2025-02-20T22:47:03Z

src/coreComponents/constitutive/contact/HydraulicApertureTable.hpp

@@ -156,10 +157,11 @@ GEOS_HOST_DEVICE
 GEOS_FORCE_INLINE
 real64 HydraulicApertureTableUpdates::computeHydraulicAperture( real64 const aperture,
                                                                real64 const normalTraction,
+                                                                integer const fractureState,


Suggested change

integer const fractureState,

FractureState const fractureState,

CusiniM · 2025-02-20T22:47:24Z

src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp

@@ -770,6 +770,7 @@ void SinglePhasePoromechanicsConformingFractures< FLOW_SOLVER >::updateHydraulic
      arrayView1d< real64 > const aperture                 = subRegion.getElementAperture();
      arrayView1d< real64 > const hydraulicAperture        = subRegion.getField< flow::hydraulicAperture >();
      arrayView1d< real64 > const deltaVolume              = subRegion.getField< flow::deltaVolume >();
+      arrayView1d< integer > const & fractureState         = subRegion.getField< contact::fractureState >();


Suggested change

arrayView1d< integer > const & fractureState = subRegion.getField< contact::fractureState >();

arrayView1d< FractureState > const & fractureState = subRegion.getField< contact::fractureState >();

CusiniM · 2025-02-20T22:47:40Z

src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp

@@ -505,6 +505,8 @@ void SinglePhasePoromechanicsEmbeddedFractures::updateState( DomainPartition & d

      arrayView2d< real64 > const & fractureContactTraction = subRegion.template getField< contact::traction >();

+      arrayView1d< integer > const & fractureState = subRegion.template getField< contact::fractureState >();


Suggested change

arrayView1d< integer > const & fractureState = subRegion.template getField< contact::fractureState >();

arrayView1d< FractureState > const & fractureState = subRegion.template getField< contact::fractureState >();

CusiniM · 2025-02-20T22:48:23Z

...Components/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp

@@ -341,7 +341,8 @@ struct StateUpdateKernel
          arrayView1d< real64 > const & aperture,
          arrayView1d< real64 const > const & oldHydraulicAperture,
          arrayView1d< real64 > const & hydraulicAperture,
-          arrayView2d< real64 > const & fractureEffectiveTraction )
+          arrayView2d< real64 > const & fractureEffectiveTraction,
+          arrayView1d< integer > const & fractureState )


Suggested change

arrayView1d< integer > const & fractureState )

arrayView1d< FractureState > const & fractureState )

CusiniM · 2025-02-21T00:17:30Z

src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolverKernels.hpp

@@ -150,8 +152,10 @@ struct FluidMassResidualDerivativeAssemblyKernel
    real64 dHydraulicAperture_dNormalJump = 0.0;
    real64 dHydraulicAperture_dTraction = 0.0;
    real64 fractureTraction = 0.0;
+    integer fractureState = 0;


why forcing this to be 0?

frankfeifan · 2025-02-21T20:04:52Z

src/coreComponents/constitutive/contact/BartonBandis.hpp

+  using namespace fields::contact;
+
+  real64 const hydraulicAperture = ( fractureState == FractureState::Open ) ? (aperture + m_aperture0) : m_aperture0 / ( 1 + 9*LvArray::math::abs(normalTraction)/m_referenceNormalStress );
+  dHydraulicAperture_dNormalStress = ( fractureState == FractureState::Open ) ? 0.0 : -hydraulicAperture / ( 1 + 9*LvArray::math::abs(normalTraction)/m_referenceNormalStress ) * 9/m_referenceNormalStress;


For fractureState != FractureState::Open, shouldn't it be dHydraulicAperture_dNormalStress = -hydraulicAperture / ( 1 + 9*abs(normalTraction)/m_referenceNormalStress ) * 9/m_referenceNormalStress * normalTraction/abs(normalTraction)?

frankfeifan · 2025-02-21T20:12:59Z

src/coreComponents/constitutive/contact/BartonBandis.hpp

+
+  real64 const hydraulicAperture = ( fractureState == FractureState::Open ) ? (aperture + m_aperture0) : m_aperture0 / ( 1 + 9*LvArray::math::abs(normalTraction)/m_referenceNormalStress );
+  dHydraulicAperture_dNormalStress = ( fractureState == FractureState::Open ) ? 0.0 : -hydraulicAperture / ( 1 + 9*LvArray::math::abs(normalTraction)/m_referenceNormalStress ) * 9/m_referenceNormalStress;
+  dHydraulicAperture_aperture = ( fractureState == FractureState::Open ) ? 1.0 : 0.0;


For current HydroFractureSolver, if the fracture is not open, this cannot be zero if we update HydraulicAperture implicitly using the current aperture. It should be something like dHydraulicAperture_aperture = dHydraulicAperture_dNormalStress * dNormalStress_dAperture where dNormalStress_dAperture = penaltyStiffness.

add a static catalog name for BartonBandis model

dca2053

jhuang2601 requested a review from sytuannguyen February 18, 2025 20:03

jhuang2601 self-assigned this Feb 18, 2025

jhuang2601 requested review from rrsettgast, CusiniM, matteofrigo5 and Guotong-Ren as code owners February 18, 2025 20:03

jhuang2601 and others added 3 commits February 18, 2025 15:17

fix compilation error

a00051f

add catalog registration.

4f9d113

add a smoke test

4fb97dc

jhuang2601 requested review from castelletto1 and paveltomin as code owners February 19, 2025 16:10

check the fracture state for aperture update

a0f92f9

jhuang2601 requested review from frankfeifan and ryar9534 as code owners February 20, 2025 16:41

update cpp files with fracture state

3236cc5

frankfeifan reviewed Feb 20, 2025

View reviewed changes

jhuang2601 added 4 commits February 20, 2025 12:09

more changes with fracture state

3323a64

fix compilation error

1ee2026

ensure consistent sign convention

8ac2fd1

update testing cases

8ccb202

jhuang2601 added type: feature New feature or request Theme - FF Faults & fractures modeling flag: ready for review labels Feb 20, 2025

CusiniM reviewed Feb 20, 2025

View reviewed changes

jhuang2601 added 3 commits February 20, 2025 17:14

address review comments

c2cd64e

fix compilation errors

32c4ae7

scope fields::contact::FractureState

8c11aec

CusiniM reviewed Feb 21, 2025

View reviewed changes

frankfeifan reviewed Feb 21, 2025

View reviewed changes

Fix all compilation erros.

e9257db

CusiniM requested review from wrtobin and untereiner as code owners February 21, 2025 20:21

CusiniM and others added 4 commits February 21, 2025 12:51

one more fix in hydrofracture solver.

9d92e55

use integer to solve compilation error.

d6d720a

Merge branch 'develop' into feature/jhuang/enableBartonBandis

7057be2

Merge branch 'develop' into feature/jhuang/enableBartonBandis

5880bfd

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

fix: enable BartonBandis model for contact solver #3556

fix: enable BartonBandis model for contact solver #3556

jhuang2601 commented Feb 18, 2025 •

edited

Loading

frankfeifan Feb 20, 2025

jhuang2601 Feb 20, 2025

CusiniM Feb 20, 2025

CusiniM Feb 20, 2025

CusiniM Feb 20, 2025

CusiniM Feb 20, 2025

CusiniM Feb 20, 2025

CusiniM Feb 20, 2025

CusiniM Feb 20, 2025

CusiniM Feb 21, 2025

frankfeifan Feb 21, 2025

frankfeifan Feb 21, 2025

	dHydraulicAperture_dNormalStress = ( fractureState >= 3.0 ) ? 0.0 : -hydraulicAperture / ( 1 + 9LvArray::math::abs(normalTraction)/m_referenceNormalStress ) 9/m_referenceNormalStress;
	dHydraulicAperture_dNormalStress = ( fractureState == FractureState::Open ) ? 0.0 : -hydraulicAperture / ( 1 + 9LvArray::math::abs(normalTraction)/m_referenceNormalStress ) 9/m_referenceNormalStress;

	dHydraulicAperture_aperture = ( fractureState >= 3.0 ) ? 1.0 : 0.0;
	dHydraulicAperture_aperture = ( fractureState >= fractureState == FractureState::Open ) ? 1.0 : 0.0;

	integer const fractureState,
	FractureState const fractureState,

	arrayView1d< integer > const & fractureState = subRegion.getField< contact::fractureState >();
	arrayView1d< FractureState > const & fractureState = subRegion.getField< contact::fractureState >();

		@@ -505,6 +505,8 @@ void SinglePhasePoromechanicsEmbeddedFractures::updateState( DomainPartition & d

		arrayView2d< real64 > const & fractureContactTraction = subRegion.template getField< contact::traction >();

		arrayView1d< integer > const & fractureState = subRegion.template getField< contact::fractureState >();

fix: enable BartonBandis model for contact solver #3556

Are you sure you want to change the base?

fix: enable BartonBandis model for contact solver #3556

Conversation

jhuang2601 commented Feb 18, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jhuang2601 commented Feb 18, 2025 •

edited

Loading