Implemented spring boundary conditions

- Spring bc is zero robin bc, following Boon et al., Computat. Geosci. 2026 and Nordbotten & Keilegavlen, Comput. Math. Appl. 2025

Author: svenn-t

opm/models/discretization/common/tpsalinearizer.hpp

@@ -244,17 +244,19 @@ class TpsaLinearizer
     {
         for (auto& bdyInfo : boundaryInfo_) {
             // Get boundary information from problem
-            const auto [type, displacementAD] = problem_().mechBoundaryCondition(bdyInfo.cell, bdyInfo.dir);
+            const auto& mechBC = problem_().mechBoundaryCondition(bdyInfo.cell, bdyInfo.dir);
 
             // Strip the unnecessary (and zero anyway) derivatives off displacement
             std::vector<double> displacement(3, 0.0);
             for (std::size_t ii = 0; ii < displacement.size(); ++ii) {
-                displacement[ii] = displacementAD[ii].value();
+                displacement[ii] = mechBC.displacement[ii].value();
             }
 
             // Update boundary information
-            bdyInfo.bcdata.type = type;
+            bdyInfo.bcdata.type = mechBC.type;
             bdyInfo.bcdata.displacement = displacement;
+            bdyInfo.bcdata.shearModulus = mechBC.shearModulus;
+            bdyInfo.bcdata.distance = mechBC.distance;
         }
     }
 
@@ -419,16 +421,21 @@ class TpsaLinearizer
                         }
 
                         // Get boundary information from problem()
-                        const auto [type, displacementAD] = problem_().mechBoundaryCondition(myIdx, dir_id);
+                        const auto& mechBC = problem_().mechBoundaryCondition(myIdx, dir_id);
 
                         // Strip the unnecessary (and zero anyway) derivatives off displacement
                         std::vector<double> displacement(3, 0.0);
                         for (std::size_t ii = 0; ii < displacement.size(); ++ii) {
-                            displacement[ii] = displacementAD[ii].value();
+                            displacement[ii] = mechBC.displacement[ii].value();
                         }
 
                         // Insert boundary condition data in container
-                        BoundaryConditionData bcdata { type, displacement, bfIndex, bfArea };
+                        BoundaryConditionData bcdata{mechBC.type,
+                                                     displacement,
+                                                     mechBC.shearModulus,
+                                                     mechBC.distance,
+                                                     bfIndex,
+                                                     bfArea};
                         boundaryInfo_.push_back( { myIdx, dir_id, bfIndex, bcdata } );
                         ++bfIndex;
                         continue;
@@ -782,6 +789,8 @@ class TpsaLinearizer
     {
         BCMECHType type;
         std::vector<double> displacement;
+        Scalar shearModulus;
+        Scalar distance;
         unsigned boundaryFaceIndex;
         double faceArea;
     };

opm/models/tpsa/elasticitylocalresidualtpsa.hpp

@@ -43,6 +43,19 @@
 
 namespace Opm {
 
+namespace detail {
+
+/*! \brief Computes modulo 3 of possibly negative integers to get indices in cross product.
+ *
+ * If i = x-dir(=0), we want y-dir(=1) and z-dir(=2), hence -1 mod 3 must equal 2 and not -1
+*/
+inline int modNeg(int i)
+{
+    return ((i % 3) + 3) % 3;
+}
+
+} // namespace detail
+
 /*!
 * \brief Calculation of (linear) elasticity model terms for the residual
 *
@@ -167,15 +180,11 @@ class ElasticityLocalResidual
         // ///
         // Displacement, rotation and solid pressure (directional-dependent) equations
         // ///
-        // Lambda function for computing modulo 3 of possibly negative integers to get indices in cross product.
-        // E.g. if i = x-dir(=0), we want y-dir(=1) and z-dir(=2), hence -1 mod 3 must equal 2 and not -1
-        auto modNeg = [](int i) { return ((i % 3) + 3) % 3; };
-
         // Loop over x-, y- and z-dir (corresponding to dirIdx = 0, 1, 2)
         for (int dirIdx = 0; dirIdx < 3; ++dirIdx) {
             // Direction indices in cross-product
-            unsigned dirIdxNeg = modNeg(dirIdx - 1);
-            unsigned dirIdxPos = modNeg(dirIdx + 1);
+            unsigned dirIdxNeg = detail::modNeg(dirIdx - 1);
+            unsigned dirIdxPos = detail::modNeg(dirIdx + 1);
 
             // Displacement equation
             const Scalar faceNormalDir = faceNormal[dirIdx];
@@ -247,6 +256,13 @@ class ElasticityLocalResidual
                                     problem,
                                     globalIndex);
             break;
+        case BCMECHType::SPRING:
+            computeBoundaryTermSpring(bndryTerm,
+                                    materialState,
+                                    bdyInfo,
+                                    problem,
+                                    globalIndex);
+            break;
         default:
             throw std::logic_error("Unknown boundary condition type " +
                                     std::to_string(static_cast<int>(bdyInfo.type)) +
@@ -295,15 +311,11 @@ class ElasticityLocalResidual
         // ///
         // Displacement equation
         // ///
-        // Lambda function for computing modulo 3 of possibly negative integers to get indices in cross product.
-        // E.g. if i = x-dir(=0), we want y-dir(=1) and z-dir(=2), hence -1 mod 3 must equal 2 and not -1
-        auto modNeg = [](int i) { return ((i % 3) + 3) % 3; };
-
         // Loop over x-, y- and z-dir (corresponding to dirIdx = 0, 1, 2)
         for (int dirIdx = 0; dirIdx < 3; ++dirIdx) {
             // Direction indices in cross-product
-            unsigned dirIdxNeg = modNeg(dirIdx - 1);
-            unsigned dirIdxPos = modNeg(dirIdx + 1);
+            unsigned dirIdxNeg = detail::modNeg(dirIdx - 1);
+            unsigned dirIdxPos = detail::modNeg(dirIdx + 1);
 
             // Displacement equation
             const Scalar faceNormalDir = faceNormal[dirIdx];
@@ -362,10 +374,6 @@ class ElasticityLocalResidual
         // ///
         // Rotation and solid pressure (directional-dependent) equations
         // ///
-        // Lambda function for computing modulo 3 of possibly negative integers to get indices in cross product.
-        // E.g. if i = x-dir(=0), we want y-dir(=1) and z-dir(=2), hence -1 mod 3 must equal 2 and not -1
-        auto modNeg = [](int i) { return ((i % 3) + 3) % 3; };
-
         // Pre-compute dot product for rotation equation
         Evaluation dotProd = 0;
         for (int dirIdx = 0; dirIdx < 3; ++dirIdx) {
@@ -375,8 +383,8 @@ class ElasticityLocalResidual
         // Loop over x-, y- and z-dir (corresponding to dirIdx = 0, 1, 2)
         for (int dirIdx = 0; dirIdx < 3; ++dirIdx) {
             // Direction indices in cross-product
-            unsigned dirIdxNeg = modNeg(dirIdx - 1);
-            unsigned dirIdxPos = modNeg(dirIdx + 1);
+            unsigned dirIdxNeg = detail::modNeg(dirIdx - 1);
+            unsigned dirIdxPos = detail::modNeg(dirIdx + 1);
 
             // Rotation equation
             const Scalar faceNormalDir = faceNormal[dirIdx];
@@ -400,6 +408,98 @@ class ElasticityLocalResidual
         }
     }
 
+    /*!
+    * \brief Calculate spring boundary condition in TPSA formulation
+    *
+    * \param bndryTerm Boundary term vector
+    * \param materialState Material state container
+    * \param bdyInfo Boundary condition info container
+    * \param problem Flow problem
+    * \param globalIndex Cell index
+    *
+    * \note Spring boundary condition sets a shear modulus at fictitious point outside a certain
+    * distance outside the boundary
+    */
+    template <class BoundaryConditionData>
+    static void computeBoundaryTermSpring(Dune::FieldVector<Evaluation, numEq>& bndryTerm,
+                                          const MaterialState& materialState,
+                                          const BoundaryConditionData& bdyInfo,
+                                          Problem& problem,
+                                          unsigned globalIndex)
+    {
+        // Reset bondary term
+        bndryTerm = 0.0;
+
+        // Extract cell and boundary information
+        const unsigned bfIdx = bdyInfo.boundaryFaceIndex;
+        const auto& faceNormal = problem.cellFaceNormalBoundary(globalIndex, bfIdx);
+        const Scalar distIn = problem.normalDistanceBoundary(globalIndex, bfIdx);
+        const Scalar distEx = bdyInfo.distance;
+        const Scalar sModulusIn = problem.shearModulus(globalIndex);
+        const Scalar sModulusEx = bdyInfo.shearModulus;
+
+        // Calculate face properties
+        const Scalar weightIn = distIn / sModulusIn;
+        const Scalar weightEx = distEx / sModulusEx;
+        const Scalar weightAvgIn = weightIn / (weightIn + weightEx);
+        const Scalar weightAvgEx = 1.0 - weightAvgIn;
+        const Scalar weightProd = weightIn * weightEx;
+        const Scalar normDist = distIn + distEx;
+
+        // Effective shear modulus
+        const Scalar eff_sModulus = weightAvgIn * sModulusIn + weightAvgEx * sModulusEx;
+
+        // Distance ratio
+        const Scalar distRatio = 0.5 * weightProd / normDist;
+
+        // Solid pressure equation (direction-independent equation)
+        const Evaluation& solidP = materialState.solidPressure();
+        bndryTerm[contiSolidPresEqIdx] +=
+            distRatio * eff_sModulus * solidP;
+
+        // Pre-compute dot product for rotation equation
+        Evaluation dotProd = 0;
+        for (int dirIdx = 0; dirIdx < 3; ++dirIdx) {
+            dotProd += faceNormal[dirIdx] * materialState.rotation(dirIdx);
+        }
+
+        // Loop over x-, y- and z-dir (corresponding to dirIdx = 0, 1, 2)
+        for (int dirIdx = 0; dirIdx < 3; ++dirIdx) {
+            // Direction indices in cross-product
+            unsigned dirIdxNeg = detail::modNeg(dirIdx - 1);
+            unsigned dirIdxPos = detail::modNeg(dirIdx + 1);
+
+            // Displacement equation
+            const Scalar faceNormalDir = faceNormal[dirIdx];
+            const Scalar faceNormalNeg = faceNormal[dirIdxNeg];
+            const Scalar faceNormalPos = faceNormal[dirIdxPos];
+
+            const Evaluation& disp = materialState.displacement(dirIdx);
+
+            const Evaluation& rotNeg = materialState.rotation(dirIdxNeg);
+            const Evaluation& rotPos = materialState.rotation(dirIdxPos);
+
+            bndryTerm[conti0EqIdx + dirIdx] +=
+                2.0 * (eff_sModulus / normDist) * disp
+                - weightAvgIn * (faceNormalNeg * rotPos - faceNormalPos * rotNeg)
+                - faceNormalDir * weightAvgIn * solidP;
+
+            // Rotation equation
+            const Evaluation& dispNeg = materialState.displacement(dirIdxNeg);
+            const Evaluation& dispPos = materialState.displacement(dirIdxPos);
+
+            const Evaluation& rot = materialState.rotation(dirIdx);
+
+            bndryTerm[contiRotEqIdx + dirIdx] +=
+                - weightAvgEx * (faceNormalNeg * dispPos - faceNormalPos * dispNeg)
+                + distRatio * eff_sModulus * (dotProd * faceNormalDir - rot);
+
+            // Solid pressure (directional-dependent) equation
+            bndryTerm[contiSolidPresEqIdx] +=
+                - faceNormalDir * weightAvgEx * disp;
+        }
+    }
+
     /*!
     * \brief Calculate source term in TPSA formulation
     *

opm/simulators/flow/FlowProblemTPSA.hpp

@@ -88,10 +88,20 @@ class FlowProblemTPSA : public FlowProblemBlackoil<TypeTag>
 
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
     using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+    using EvalDimVector = Dune::FieldVector<Evaluation, dimWorld>;
     using FaceProperties = FacePropertiesTPSA<Grid, GridView, ElementMapper, CartesianIndexMapper, Scalar>;
     using InitialMaterialState = MaterialStateTPSA<Scalar>;
     using Toolbox = MathToolbox<Evaluation>;
 
+    // Boundary condition helper struct
+    struct MechBCData
+    {
+        BCMECHType type;
+        EvalDimVector displacement;
+        Scalar distance;
+        Scalar shearModulus;
+    };
+
     // ///
     // Public functions
     // ///
@@ -277,31 +287,36 @@ class FlowProblemTPSA : public FlowProblemBlackoil<TypeTag>
     * \note Only BCMECHTYPE = FREE and NONE implemented. FIXED will/should throw an error when computed in local
     * residual!
     */
-    std::pair<BCMECHType, Dune::FieldVector<Evaluation, 3>>
+    MechBCData
     mechBoundaryCondition(const unsigned int globalSpaceIdx, const int directionId)
     {
         // Default boundary conditions if BCCON/BCMECH not defined
         if (!this->nonTrivialBoundaryConditions()) {
-            return { BCMECHType::NONE, Dune::FieldVector<Evaluation, 3>{0.0, 0.0, 0.0} };
+            return {BCMECHType::NONE, EvalDimVector{0.0, 0.0, 0.0}, 0.0, 0.0};
         }
 
         // Default for BCMECH index = 0 or no BCMECH defined at current episode
         FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
         const auto& schedule = this->simulator().vanguard().schedule();
         if (this->bcindex_(dir)[globalSpaceIdx] == 0
             || schedule[this->episodeIndex()].bcstate.size() == 0) {
-            return {BCMECHType::NONE, Dune::FieldVector<Evaluation, 3>{0.0, 0.0, 0.0} };
+            return {BCMECHType::NONE, EvalDimVector{0.0, 0.0, 0.0}, 0.0, 0.0};
         }
 
         // Get current BC
         const auto& bc =
             schedule[this->episodeIndex()].bcstate[this->bcindex_(dir)[globalSpaceIdx]];
         if (bc.bcmechtype == BCMECHType::FREE) {
-            return { BCMECHType::FREE, Dune::FieldVector<Evaluation, 3>{0.0, 0.0, 0.0} };
+            return {BCMECHType::FREE, EvalDimVector{0.0, 0.0, 0.0}, 0.0, 0.0};
         }
-        else {
-            return { bc.bcmechtype, Dune::FieldVector<Evaluation, 3>{0.0, 0.0, 0.0} };
+        if (bc.bcmechtype == BCMECHType::SPRING) {
+            const auto& mechbcval = bc.mechbcvalue;
+            return {BCMECHType::SPRING, EvalDimVector{0.0, 0.0, 0.0},
+                    mechbcval.distance, mechbcval.shearmodulus};
         }
+
+        // Default return
+        return {bc.bcmechtype, EvalDimVector{0.0, 0.0, 0.0}, 0.0, 0.0};
     }
 
     /*!