fix: EFEM bugfixes - effective traction + oldStress

.integrated_tests.yaml

@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3280-8463-4595d21
+  baseline: integratedTests/baseline_integratedTests-pr3408-8504-5db9ba0
 
 allow_fail:
   all: ''

BASELINE_NOTES.md

@@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
 
+PR #3408 (2024-11-06)
+====================
+EFEM bugfixes: effective traction + oldStress.
+
 PR #3280 (2024-11-05)
 ====================
 Added Sprig-slider test.

src/coreComponents/constitutive/solid/CoupledSolidBase.hpp

@@ -170,6 +170,14 @@
     return getBaseSolidModel().getDensity();
   }
 
+  /*
+   * @brief get the current solid effective stress
+   * return a constant arrayView3d to effective stress in Voigt form
+   */
+  arrayView3d< real64 const, solid::STRESS_USD > const getEffectiveStress() const
+  {
+    return getBaseSolidModel().getStress();
+  }
 
   /*
    * @brief get the current biot coefficient

src/coreComponents/finiteVolume/TwoPointFluxApproximation.cpp

@@ -591,7 +591,7 @@
       stencilCellsRegionIndex[1] = fractureRegionIndex;
       stencilCellsSubRegionIndex[1] = 0;
       stencilCellsIndex[1] = kes;
-      stencilWeights[1] = 4. * faceArea[fractureRegionIndex] / hydraulicAperture[fractureRegionIndex][0][kes];
+      stencilWeights[1] = 4. * faceArea[kes] / hydraulicAperture[fractureRegionIndex][0][kes];
 
       edfmStencil.add( 2,
                        stencilCellsRegionIndex.data(),

src/coreComponents/physicsSolvers/contact/SolidMechanicsEFEMKernels.hpp

@@ -220,8 +220,9 @@
 
     // Compute the local residuals
     LvArray::tensorOps::Ri_add_AijBj< 3, 3 >( stack.localRw, stack.localKww, stack.wLocal );
-    LvArray::tensorOps::Ri_add_AijBj< 3, nUdof >( stack.localRw, stack.localKwu, stack.uLocal );
     LvArray::tensorOps::Ri_add_AijBj< nUdof, 3 >( stack.localRu, stack.localKuw, stack.wLocal );
+    // add EqM * effStress into the residual of enrichment nodes
+    LvArray::tensorOps::add< 3 >( stack.localRw, stack.localEqMStress );
 
     // Add traction contribution
     LvArray::tensorOps::scaledAdd< 3 >( stack.localRw, stack.tractionVec, -1 );

src/coreComponents/physicsSolvers/contact/SolidMechanicsEFEMKernelsBase.hpp

@@ -107,6 +107,7 @@
           inputDt,
           inputGravityVector ),
     m_w( embeddedSurfSubRegion.getField< fields::contact::dispJump >().toView() ),
+    m_effStress( inputConstitutiveType.getStress()),
     m_tractionVec( embeddedSurfSubRegion.getField< fields::contact::traction >().toViewConst() ),
     m_dTraction_dJump( embeddedSurfSubRegion.getField< fields::contact::dTraction_dJump >().toViewConst() ),
     m_nVec( embeddedSurfSubRegion.getNormalVector().toViewConst() ),
@@ -144,6 +145,7 @@
       localKww{ { 0.0 } },
       localKwu{ { 0.0 } },
       localKuw{ { 0.0 } },
+      localEqMStress{ 0.0 },
       wLocal(),
       uLocal(),
       hInv(),
@@ -171,6 +173,9 @@
     /// C-array storage for the element local Kuw matrix.
     real64 localKuw[numUdofs][numWdofs];
 
+    /// C-array storage for the element local EqM*effStress vector.
+    real64 localEqMStress[numWdofs];
+
     /// Stack storage for the element local jump vector
     real64 wLocal[3];
 
@@ -249,6 +254,9 @@
     // Gauss contribution to Kww, Kwu and Kuw blocks
     real64 Kww_gauss[3][3], Kwu_gauss[3][nUdof], Kuw_gauss[nUdof][3];
 
+    // Gauss contirbution to eqMStress which is EqMatrix*effStress, all stresses are in Voigt notation
+    real64 eqMStress_gauss[3]{};
+
     //  Compatibility, equilibrium and strain operators. The compatibility operator is constructed as
     //  a 3 x 6 because it is more convenient for construction purposes (reduces number of local var).
     real64 compMatrix[3][6], strainMatrix[6][nUdof], eqMatrix[3][6];
@@ -292,17 +300,25 @@
     LvArray::tensorOps::Rij_eq_AikBkj< 3, nUdof, 6 >( Kwu_gauss, matED, strainMatrix );
     // transp(B)DB
     LvArray::tensorOps::Rij_eq_AikBjk< nUdof, 3, 6 >( Kuw_gauss, matBD, compMatrix );
+    // EqMatrix * effStress
+    LvArray::tensorOps::Ri_eq_AijBj< 3, 6 >( eqMStress_gauss, eqMatrix, m_effStress[k][q] );
+
+    /// FIX: add old Equilibrium operator times oldStress (in Voigt notation)
 
     // multiply by determinant and add to element matrix
     LvArray::tensorOps::scaledAdd< 3, 3 >( stack.localKww, Kww_gauss, -detJ );
     LvArray::tensorOps::scaledAdd< 3, nUdof >( stack.localKwu, Kwu_gauss, -detJ );
     LvArray::tensorOps::scaledAdd< nUdof, 3 >( stack.localKuw, Kuw_gauss, -detJ );
+    LvArray::tensorOps::scaledAdd< 3 >( stack.localEqMStress, eqMStress_gauss, -detJ );
   }
 
 protected:
 
   arrayView2d< real64 > const m_w;
 
+  /// The effective stress at the current time
+  arrayView3d< real64 const, solid::STRESS_USD > m_effStress;
+
   arrayView2d< real64 const > const m_tractionVec;
 
   arrayView3d< real64 const > const m_dTraction_dJump;

src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp

@@ -489,9 +489,7 @@
 
       arrayView1d< real64 const > const area = subRegion.getElementArea().toViewConst();
 
-      arrayView2d< real64 > const & fractureTraction = subRegion.template getField< fields::contact::traction >();
-
-      arrayView1d< real64 >  const & dTdpf = subRegion.template getField< fields::contact::dTraction_dPressure >();
+      arrayView2d< real64 > const & fractureContactTraction = subRegion.template getField< fields::contact::traction >();
 
       arrayView1d< real64 const > const & pressure =
         subRegion.template getField< fields::flow::pressure >();
@@ -524,8 +522,7 @@
                                               aperture,
                                               oldHydraulicAperture,
                                               hydraulicAperture,
-                                              fractureTraction,
-                                              dTdpf );
+                                              fractureContactTraction );
 
         } );
       } );

src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp

@@ -115,6 +115,7 @@
       localKww{ { 0.0 } },
       localKwu{ { 0.0 } },
       localKuw{ { 0.0 } },
+      localEqMStress { 0.0 },
       localKwpm{ 0.0 },
       localKwpf( 0.0 ),
       wLocal(),
@@ -154,6 +155,9 @@
     /// C-array storage for the element local Kuw matrix.
     real64 localKuw[numUdofs][numWdofs];
 
+    /// C-array storage for the element local EqM*effStress vector.
+    real64 localEqMStress[numWdofs];
+
     /// C-array storage for the element local Kwpm matrix.
     real64 localKwpm[numWdofs];
 
@@ -233,6 +237,9 @@
 
   arrayView2d< real64 const > const m_w;
 
+  /// The effective stress at the current time
+  arrayView3d< real64 const, solid::STRESS_USD > m_effStress;
+
   /// The global degree of freedom number
   arrayView1d< globalIndex const > const m_matrixPresDofNumber;
 
@@ -249,6 +256,9 @@
   /// The rank-global fluid pressure array.
   arrayView1d< real64 const > const m_matrixPressure;
 
+  /// The rank-global fluid pressure array.
+  arrayView1d< real64 const > const m_fracturePressure;
+
   /// The rank-global delta-fluid pressure array.
   arrayView2d< real64 const > const m_porosity_n;
 
@@ -314,8 +324,7 @@
    * @param[out] deltaVolume the change in volume
    * @param[out] aperture the aperture
    * @param[out] hydraulicAperture the effecture aperture
-   * @param[out] fractureTraction the fracture traction
-   * @param[out] dFractureTraction_dPressure the derivative of the fracture traction wrt pressure
+   * @param[out] fractureContactTraction the fracture contact traction
    */
   template< typename POLICY, typename POROUS_WRAPPER, typename CONTACT_WRAPPER >
   static void
@@ -330,8 +339,7 @@
           arrayView1d< real64 > const & aperture,
           arrayView1d< real64 const > const & oldHydraulicAperture,
           arrayView1d< real64 > const & hydraulicAperture,
-          arrayView2d< real64 > const & fractureTraction,
-          arrayView1d< real64 > const & dFractureTraction_dPressure )
+          arrayView2d< real64 > const & fractureEffectiveTraction )
   {
     forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
     {
@@ -341,23 +349,21 @@
       real64 dHydraulicAperture_dNormalJump = 0.0;
       real64 dHydraulicAperture_dNormalTraction = 0.0;
       hydraulicAperture[k] = contactWrapper.computeHydraulicAperture( aperture[k],
-                                                                      fractureTraction[k][0],
+                                                                      fractureEffectiveTraction[k][0],
                                                                       dHydraulicAperture_dNormalJump,
                                                                       dHydraulicAperture_dNormalTraction );
 
       deltaVolume[k] = hydraulicAperture[k] * area[k] - volume[k];
 
-      // traction on the fracture to include the pressure contribution
-      fractureTraction[k][0] -= pressure[k];
-      dFractureTraction_dPressure[k] = -1.0;
-
       real64 const jump[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( dispJump[k] );
-      real64 const traction[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( fractureTraction[k] );
+      real64 const effectiveTraction[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( fractureEffectiveTraction[k] );
 
+      // all perm update models below should need effective traction instead of total traction
+      // (total traction is combined forces of fluid pressure and effective traction)
       porousMaterialWrapper.updateStateFromPressureApertureJumpAndTraction( k, 0, pressure[k],
                                                                             oldHydraulicAperture[k], hydraulicAperture[k],
                                                                             dHydraulicAperture_dNormalJump,
-                                                                            jump, traction );
+                                                                            jump, effectiveTraction );
 
     } );
   }

src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp

@@ -71,6 +71,7 @@
   m_disp( nodeManager.getField< fields::solidMechanics::totalDisplacement >() ),
   m_deltaDisp( nodeManager.getField< fields::solidMechanics::incrementalDisplacement >() ),
   m_w( embeddedSurfSubRegion.getField< fields::contact::dispJump >() ),
+  m_effStress( inputConstitutiveType.getEffectiveStress()),
   m_matrixPresDofNumber( elementSubRegion.template getReference< array1d< globalIndex > >( inputFlowDofKey ) ),
   m_fracturePresDofNumber( embeddedSurfSubRegion.template getReference< array1d< globalIndex > >( inputFlowDofKey ) ),
   m_wDofNumber( jumpDofNumber ),
@@ -80,6 +81,7 @@
   m_dFluidDensity_dPressure( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >(
                                                                                                                      fluidModelKey ) ).dDensity_dPressure() ),
   m_matrixPressure( elementSubRegion.template getField< fields::flow::pressure >() ),
+  m_fracturePressure( embeddedSurfSubRegion.template getField< fields::flow::pressure >() ),
   m_porosity_n( inputConstitutiveType.getPorosity_n() ),
   m_tractionVec( embeddedSurfSubRegion.getField< fields::contact::traction >() ),
   m_dTraction_dJump( embeddedSurfSubRegion.getField< fields::contact::dTraction_dJump >() ),
@@ -188,6 +190,12 @@
                        FUNC && kernelOp ) const
 {
 
+  // The quarature kernal deals with the fracture force balance (eq. 29 in https://onlinelibrary.wiley.com/doi/epdf/10.1002/nag.3168)
+  // The total stress in matrix: sigma_tau = simga_eff_old + sigma_eff_incr - biot * p_m = sigma_eff_new - biot * p_m
+  // We can either use formulation: sigma_eff_old + stiffness_matrix * incremental_strain - biot * p_m or directly effective_stress_current
+  // - biot * p_m
+  // The latter one is adopted here.
+
   localIndex const embSurfIndex = m_cellsToEmbeddedSurfaces[k][0];
 
   // Get displacement: (i) basis functions (N), (ii) basis function
@@ -202,6 +210,9 @@
   // Gauss contribution to Kww, Kwu and Kuw blocks
   real64 Kww_gauss[3][3]{}, Kwu_gauss[3][nUdof]{}, Kuw_gauss[nUdof][3]{}, Kwpm_gauss[3]{};
 
+  // Gauss contirbution to eqMStress which is EqMatrix*effStress, all stresses are in Voigt notation
+  real64 eqMStress_gauss[3]{};
+
   //  Compatibility, equilibrium and strain operators. The compatibility operator is constructed as
   //  a 3 x 6 because it is more convenient for construction purposes (reduces number of local var).
   real64 compMatrix[3][6]{}, strainMatrix[6][nUdof]{}, eqMatrix[3][6]{};
@@ -247,6 +258,8 @@
   LvArray::tensorOps::Rij_eq_AikBkj< 3, nUdof, 6 >( Kwu_gauss, matED, strainMatrix );
   // transp(B)DB
   LvArray::tensorOps::Rij_eq_AikBjk< nUdof, 3, 6 >( Kuw_gauss, matBD, compMatrix );
+  // EqMatrix * effStress
+  LvArray::tensorOps::Ri_eq_AijBj< 3, 6 >( eqMStress_gauss, eqMatrix, m_effStress[k][q] );
 
   LvArray::tensorOps::fill< 3 >( Kwpm_gauss, 0 );
   for( int i=0; i < 3; ++i )
@@ -260,6 +273,7 @@
   LvArray::tensorOps::scaledAdd< 3, 3 >( stack.localKww, Kww_gauss, -detJ );
   LvArray::tensorOps::scaledAdd< 3, nUdof >( stack.localKwu, Kwu_gauss, -detJ );
   LvArray::tensorOps::scaledAdd< nUdof, 3 >( stack.localKuw, Kuw_gauss, -detJ );
+  LvArray::tensorOps::scaledAdd< 3 >( stack.localEqMStress, eqMStress_gauss, -detJ );
 
   /// TODO: should this be negative???
   // I had No neg coz the total stress = effective stress - porePressure
@@ -285,20 +299,25 @@
 
   // Compute the local residuals
   LvArray::tensorOps::Ri_add_AijBj< 3, 3 >( stack.localJumpResidual, stack.localKww, stack.wLocal );
-  LvArray::tensorOps::Ri_add_AijBj< 3, nUdof >( stack.localJumpResidual, stack.localKwu, stack.dispLocal );
   LvArray::tensorOps::Ri_add_AijBj< nUdof, 3 >( stack.localDispResidual, stack.localKuw, stack.wLocal );
+  // add EqM * effStress into the residual of enrichment nodes
+  LvArray::tensorOps::add< 3 >( stack.localJumpResidual, stack.localEqMStress );
 
   // add pore pressure contribution
   LvArray::tensorOps::scaledAdd< 3 >( stack.localJumpResidual, stack.localKwpm, m_matrixPressure[ k ] );
 
   localIndex const embSurfIndex = m_cellsToEmbeddedSurfaces[k][0];
 
-  // Add traction contribution tranction
+  // Add total traction contribution from penalty force and fracture pressure
+  // total traction is T_total = -k * dispJump + pf (where dispJump < 0)
+  // -1 is because k*dispJump was saved in tractionVec
   LvArray::tensorOps::scaledAdd< 3 >( stack.localJumpResidual, stack.tractionVec, -1 );
   LvArray::tensorOps::scaledAdd< 3, 3 >( stack.localKww, stack.dTractiondw, -1 );
 
-  // JumpFractureFlowJacobian
-  real64 const localJumpFracPressureJacobian = -m_dTraction_dPressure[embSurfIndex] * m_surfaceArea[embSurfIndex];
+  // fracture pressure only affects normal direction
+  stack.localJumpResidual[0] += m_fracturePressure[embSurfIndex] * m_surfaceArea[embSurfIndex];
+  // fracture force balance residual w.r.t. fracture pressure
+  real64 const localJumpFracPressureJacobian = m_surfaceArea[embSurfIndex];
 
   // Mass balance accumulation
   real64 const newVolume = m_elementVolume( embSurfIndex ) + m_deltaVolume( embSurfIndex );