Add mass-proportional Rayleigh damping for solid mechanics

src/assemble/solidMechanics/solidDisplacement/solidDisplacementAssemblerNodeTerms.cpp

@@ -1,7 +1,8 @@
 // File       : solidDisplacementAssemblerNodeTerms.cpp
 // Created    : Fri Jan 31 2026 10:00:00 (+0100)
 // Author     : Mhamad Mahdi Alloush
-// Description:
+// Modified   : Added mass-proportional Rayleigh damping (alpha * rho * d(D)/dt)
+//              matching solids4foam's "dampingCoeff" convention.
 // Copyright 2026 CCFNUM HSLU T&A. All Rights Reserved.
 
 #include "solidDisplacementAssembler.h"
@@ -45,13 +46,30 @@ void solidDisplacementAssembler::assembleNodeTermsFusedSecondOrderUnsteady_(
 #ifdef USE_CVFEM_SOLID_MECHANICS
     // ========================================================================
     // Second-order time discretization for structural dynamics:
-    // ρ * d²D/dt² ≈ ρ/dt² * (D^{n+1} - 2*D^n + D^{n-1})
     //
-    // For constant density (typical for solid mechanics):
-    // LHS diagonal: ρ*V/dt²
-    // RHS: ρ*V/dt² * (D - 2*D_old + D_old_old)
+    //   rho * d^2 D / dt^2 ~ rho / dt^2 * (D^{k+1} - 2*D^n + D^{n-1})
     //
-    // Note: Uses reference volume V₀ for total Lagrangian formulation
+    // where D^{k+1} is the new outer-iteration displacement, D^n is the
+    // converged displacement at the end of the previous time step, and
+    // D^{n-1} is the converged displacement two time steps back.
+    //
+    // ------------------------------------------------------------------------
+    // Mass-proportional Rayleigh damping (alpha)
+    // ------------------------------------------------------------------------
+    // When damping_coeff = alpha > 0, the additional term
+    //
+    //   alpha * rho * d(D)/dt ~ alpha * rho / dt * (D^{k+1} - D^n)
+    //
+    // is added to the equation. The discretization on the first time
+    // derivative is BDF1 (Euler implicit), which matches the default
+    // behaviour of solids4foam's fvm::ddt(D) in its solid models. The term
+    // contributes
+    //   LHS diagonal: + alpha * rho * V_0 / dt
+    //   RHS         : - alpha * rho * V_0 / dt * (D^{k+1} - D^n)
+    // in the residual form used elsewhere in this assembly; at convergence
+    // (D^{k+1} -> D^n+1) it correctly enforces the damped momentum balance.
+    // Reduces the under-damped oscillation seen in undamped FSI runs and
+    // matches the behaviour of s4f's "dampingCoeff" entry in solidProperties.
     // ========================================================================
 
     auto& mesh = field_broker_->meshRef();
@@ -74,12 +92,17 @@ void solidDisplacementAssembler::assembleNodeTermsFusedSecondOrderUnsteady_(
     scalar* p_lhs = &lhs[0];
     scalar* p_rhs = &rhs[0];
 
-    // Get displacement fields
+    // Get displacement fields (current iter, prev time, prev-prev time)
     const STKScalarField* DSTKFieldPtr = phi_->stkFieldPtr();
     const STKScalarField* DSTKFieldPtrOld = phi_->prevTimeRef().stkFieldPtr();
     const STKScalarField* DSTKFieldPtrOldOld =
         phi_->prevTimeRef().prevTimeRef().stkFieldPtr();
 
+    // Rayleigh mass-proportional damping coefficient (1/s). When > 0, the
+    // term alpha * rho * d(D)/dt is added to the equation (BDF1 stencil).
+    const scalar dampingCoeff = domain->solidMechanics_.dampingCoeff_;
+    const bool useDamping = (dampingCoeff > 0.0);
+
     // Get density field (constant for solid mechanics)
     const STKScalarField* rhoSTKFieldPtr = this->rhoRef().stkFieldPtr();
 
@@ -90,6 +113,7 @@ void solidDisplacementAssembler::assembleNodeTermsFusedSecondOrderUnsteady_(
     // Time step
     const scalar dt = mesh.controlsRef().getTimestep();
     const scalar rDeltaT2 = 1.0 / (dt * dt);
+    const scalar rDeltaT = 1.0 / dt; // needed for damping (1st time deriv)
 
     // Get interior parts
     const stk::mesh::PartVector& partVec = domain->zonePtr()->interiorParts();
@@ -137,9 +161,15 @@ void solidDisplacementAssembler::assembleNodeTermsFusedSecondOrderUnsteady_(
             const scalar rho = rhob[iNode];
             const scalar vol = volb[iNode];
 
-            // LHS coefficient: ρ*V/dt²
+            // LHS coefficient: rho * V / dt^2
             const scalar lhsfac = rho * vol * rDeltaT2;
 
+            // Damping LHS factor (alpha * rho * V_0 / dt). Computed once per
+            // node and zeroed when damping is disabled, so the inner-loop
+            // branch is predictable and cheap.
+            const scalar lhsfacDamp =
+                useDamping ? (dampingCoeff * rho * vol * rDeltaT) : 0.0;
+
             for (label i = 0; i < SPATIAL_DIM; ++i)
             {
                 const scalar Di = Db[SPATIAL_DIM * iNode + i];
@@ -149,8 +179,20 @@ void solidDisplacementAssembler::assembleNodeTermsFusedSecondOrderUnsteady_(
                 // LHS: add to diagonal
                 p_lhs[i * SPATIAL_DIM + i] += lhsfac;
 
-                // RHS: ρ*V/dt² * (D - 2*D_old + D_old_old)
+                // RHS: rho*V/dt^2 * (D - 2*D_old + D_old_old)
                 p_rhs[i] -= lhsfac * (Di - 2.0 * DiOld + DiOldOld);
+
+                // Mass-proportional Rayleigh damping
+                //   LHS diag : + alpha * rho * V / dt
+                //   RHS      : - alpha * rho * V / dt * (D - D_old)
+                // (residual form; at convergence D -> D^{n+1} this enforces
+                //  alpha*rho*(D^{n+1} - D^n)/dt, matching s4f's
+                //  dampingCoeff*rho*fvm::ddt(D))
+                if (useDamping)
+                {
+                    p_lhs[i * SPATIAL_DIM + i] += lhsfacDamp;
+                    p_rhs[i] -= lhsfacDamp * (Di - DiOld);
+                }
             }
 
             Base::applyCoeff_(

src/domain/domain.h

@@ -67,6 +67,13 @@ struct solidMechanics
     bool planeStress_ = false;
     bool lumpedMass_ =
         true; // true = lumped (diagonal), false = consistent (full mass matrix)
+
+    // Mass-proportional Rayleigh damping coefficient [1/s].
+    // Adds the term alpha * rho * d(D)/dt to the solid momentum equation,
+    // implicit in the new displacement (BDF1 stencil on the first time
+    // derivative). Matches the "dampingCoeff" setting in solids4foam's solid
+    // models. Default 0.0 reproduces the undamped solver exactly.
+    scalar dampingCoeff_ = 0.0;
 };
 
 struct buoyancy

src/domain/domainIO.cpp

@@ -1,4 +1,4 @@
-// File       : domain.cpp
+// File       : domainIO.cpp
 // Created    : Wed Jan 03 2024 13:38:51 (+0100)
 // Author     : Mhamad Mahdi Alloush
 // Description:
@@ -1020,6 +1020,27 @@ void domain::read_()
                              "not implemented");
                 }
             }
+
+            // ================================================================
+            // Mass-proportional Rayleigh damping coefficient (alpha, units
+            // 1/s). Adds the term alpha*rho*d(D)/dt (BDF1 implicit) to the
+            // solid momentum equation. Matches solids4foam's "dampingCoeff"
+            // entry in the solidProperties dictionary. Default 0 = no damping,
+            // in which case the assembly reproduces the original undamped
+            // solver exactly (no contributions added to the LHS or RHS at the
+            // node-terms stage).
+            // ================================================================
+            if (solidMechanicsBlock["damping_coeff"])
+            {
+                solidMechanics_.dampingCoeff_ =
+                    solidMechanicsBlock["damping_coeff"]
+                        .template as<scalar>();
+
+                if (solidMechanics_.dampingCoeff_ < 0.0)
+                {
+                    errorMsg("damping_coeff must be non-negative");
+                }
+            }
         }
     }
     else