From a18bd5e634445982fd1adadfd8052caeb27b52ab Mon Sep 17 00:00:00 2001 From: Mohamad Nasreddine Date: Thu, 30 Jul 2026 05:34:38 +0300 Subject: [PATCH] Add mass-proportional Rayleigh damping for solid mechanics Adds an optional damping_coeff (alpha, units 1/s) to the solid_mechanics block. When alpha > 0, the term alpha*rho*d(D)/dt is added to the solid momentum equation, discretized BDF1 (Euler implicit) to match solids4foam's fvm::ddt(D) and its dampingCoeff convention. Contributes alpha*rho*V_0/dt to the LHS diagonal and -alpha*rho*V_0/dt*(D - D_old) to the RHS residual. Default 0 reproduces the undamped solver exactly. Validated on the cavityFlexibleBottom case. --- .../solidDisplacementAssemblerNodeTerms.cpp | 60 ++++++++++++++++--- src/domain/domain.h | 7 +++ src/domain/domainIO.cpp | 23 ++++++- 3 files changed, 80 insertions(+), 10 deletions(-) diff --git a/src/assemble/solidMechanics/solidDisplacement/solidDisplacementAssemblerNodeTerms.cpp b/src/assemble/solidMechanics/solidDisplacement/solidDisplacementAssemblerNodeTerms.cpp index 056dc8f..9780e58 100644 --- a/src/assemble/solidMechanics/solidDisplacement/solidDisplacementAssemblerNodeTerms.cpp +++ b/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_( diff --git a/src/domain/domain.h b/src/domain/domain.h index 54bb24d..a13295a 100644 --- a/src/domain/domain.h +++ b/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 diff --git a/src/domain/domainIO.cpp b/src/domain/domainIO.cpp index fca1c3e..d237dff 100644 --- a/src/domain/domainIO.cpp +++ b/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(); + + if (solidMechanics_.dampingCoeff_ < 0.0) + { + errorMsg("damping_coeff must be non-negative"); + } + } } } else