PaperlayerView exact revision on arXiv

Source-linked AI summary

A Discontinuous Galerkin Method for the Intrinsic Beam Model with Kelvin-Voigt Damping

Shivam Sundriyal, Christian Bleffert, Lukas Dreyer, Melven Röhrig-Zöllner, Gregor Gassner, Vadym Aizinger

arXiv:2609.08397v1math.NAcs.NA

TL;DR

Large-motion beam models with Kelvin–Voigt damping pose a numerical challenge because the damping law introduces mixed space–time derivatives. The paper removes them with an auxiliary-variable formulation and constructs an energy-stable DG discretization; its experiments support near-fourth-order convergence and show relaxation toward the analytic constant-spin branch with an actuated energy balance.

  • Problem

    Kelvin–Voigt damping introduces mixed space–time derivatives into the nonlinear intrinsic beam equations, for which the paper reports no prior DG method with a semidiscrete energy analysis.

    Source: §1 · p. 3

  • Method

    An auxiliary variable splits the rate-dependent sectional-resultant contribution, removing mixed derivatives and enabling an energy-stable DG spatial discretization.

    Source: §1 · p. 3, §6 · p. 25

  • Results

    The cubic scheme supports near-fourth-order convergence of all twelve primary variables, while the rotating-beam simulation relaxes toward the analytic constant-spin branch and closes the actuated energy ledger.

    Source: §6 · p. 25

  • Takeaways & Limitations

    The reformulation provides a convenient way to treat Kelvin–Voigt damping while preserving the energy structure, with the DG method offering a foundation for further beam and rotor-blade simulations.

    Source: §6 · p. 25

    The energy analyses cover homogeneous boundary data; extending them to inhomogeneous boundary data and practical fully discrete time integrators remains future work.

    Source: §4.3 · p. 19, §7 · p. 26

Abstract

from arXiv · show

We present an energy-stable discontinuous Galerkin (DG) spatial discretization of the intrinsic beam equations with Kelvin-Voigt damping. We reformulate the governing equations by splitting the total sectional resultants into elastic and viscous parts, thereby exposing a mixed system without mixed space-time derivatives. For homogeneous cantilever boundary data, the continuous system obeys an energy dissipation identity. With characteristic upwinding and alternating auxiliary traces the DG discretization is energy stable at the semidiscrete level for the stated homogeneous split boundary closure. The method is implemented in the Julia-based simulation framework Trixi.jl. A physically compatible manufactured solution with all state and viscous components active assesses convergence, while a rotating-beam test verifies the analytic constant-spin branch and the actuated energy balance.

1. Introduction

Large-motion beam simulations motivate one-dimensional geometrically exact models, but Kelvin–Voigt damping complicates their numerical treatment. The paper introduces an auxiliary-variable formulation and an energy-stable DG method, with continuous and semidiscrete energy analyses under specified homogeneous boundary conditions.

1. Introduction

  • Large deflections, rotations, and vibrations make one-dimensional beam representations attractive when fully three-dimensional finite-element simulations are expensive for repeated system-level analyses.Classical linear Euler–Bernoulli and Timoshenko models do not describe large motions.Source: §1 · p. 2
  • The Kelvin–Voigt term creates mixed space–time derivatives, and prior DG energy analyses do not address the nonlinear damped intrinsic-beam structure and its boundary closure.The authors identify the absence of a DG method with a semidiscrete energy analysis for this model as the gap.Source: §1 · p. 3
  • An auxiliary variable for the rate-dependent sectional resultants removes the mixed derivatives and yields a first-order-in-time system suited to energy analysis and DG discretization.This reformulation extends energy-stable DG ideas from the undamped intrinsic beam equations [18].Source: §1 · p. 3
  • For homogeneous cantilever data, the continuous reformulated problem has a mechanical energy–dissipation identity and an a priori stability bound under stated positive-definiteness assumptions.The identity retains the boundary-work term.Source: §1 · p. 3
  • Characteristic upwinding and compatible auxiliary traces give semidiscrete energy stability under the stated homogeneous split cantilever closure.The method is implemented in Trixi.jl; experiments assess convergence and Kelvin–Voigt damping effects.Source: §1 · p. 3

2. Damped Intrinsic Beam Equations

The section reformulates Kelvin–Voigt damping as a constitutive resultant split, avoiding mixed space–time derivatives and supporting energy analysis and DG discretization. It also describes boundary data compatible with the reformulated system.

2. Damped Intrinsic Beam Equations

  • The intrinsic beam model uses strain, curvature, velocity, and sectional-resultant variables as primary unknowns, and its governing equations comprise coupled PDEs with constitutive closure.The model represents initially curved and twisted beams of anisotropic material using a one-dimensional reference line.Source: §2 · pp. 3–4
  • The formulation assumes uniformly symmetric positive-definite mass and flexibility matrices and an effective Kelvin–Voigt operator; singular compliance limits are excluded.The stated condition on the effective damping operator supports coercive material dissipation and a well-defined inverse auxiliary relation.Source: §2 · pp. 5–6
  • The Kelvin–Voigt law introduces mixed space–time derivatives, which complicate direct discretization and stability analysis.The paper introduces rτ as the rate-dependent contribution to the internal force and moment resultants, yielding an equivalent first-order-in-time system without mixed derivatives.Source: §2 · pp. 5–6
  • Unlike gradient-based auxiliary variables in many mixed formulations, rτ is a constitutive quantity representing the damping resultant and combining derivative and zero-order contributions.The resulting formulation confines nonlinear contributions to lower-order source terms and the auxiliary relation, simplifying stability analysis and numerical-flux design.Source: §2 · p. 6

2.1. Boundary Conditions

  • For the hyperbolic characteristics, boundary data are needed only for incoming components, and the damped model retains the undamped intrinsic beam's characteristic structure.In the undamped case, the stated boundary data correspond to sectional velocities at the root and sectional forces and moments at the tip.Source: §2.1 · p. 7
  • At the tip, prescribed total traction determines elastic and viscous boundary targets that the DG method imposes separately, with their sum satisfying the total-resultant condition.The traction condition has a strong-form endpoint ODE interpretation.Source: §2.1 · pp. 7–8

3. Continuous Analysis

The continuous energy analysis establishes a boundary-work identity for the reformulated damped system and a stability bound under homogeneous cantilever conditions. Its energy calculation depends on fixed boundary prescriptions; in the unforced homogeneous case, mechanical energy is non-increasing.

3. Continuous Analysis

  • The continuous analysis derives an energy identity retaining boundary work, then obtains an a priori stability bound for homogeneous cantilever data.The proof adapts prior energy arguments for the undamped and Kelvin–Voigt damped intrinsic beam equations.Source: §3 · pp. 8, 10
  • The energy calculation applies to fixed zero-boundary-power partitions, but well-posedness for arbitrary partitions requires separate analysis.The paper uses the homogeneous cantilever partition.Source: §3 · p. 8
  • With homogeneous boundary data and no external load, the mechanical energy is non-increasing because the damping contribution is non-negative.With bounded external forces and moments, the energy remains bounded by the initial data and total energy input.Source: §3 · p. 11

4. Discretization

The discretization combines characteristic fluxes with compatible auxiliary traces to reproduce the continuous energy structure. Under the stated homogeneous boundary closure, it satisfies a semidiscrete energy–dissipation identity and stability bound.

4. Discretization

  • The elementwise formulation approximates u and rτ with polynomials of degree at most k and uses numerical fluxes for the first-order operator and auxiliary coupling.The elementwise relations are assembled into the global semidiscrete DG formulation.Source: §4 · p. 11, §4.1 · p. 12, §4.2 · pp. 13–14

4.1. Numerical Fluxes and Boundary Treatment

  • At interfaces, σ = 0 gives the central flux and σ = 1 the characteristic upwind flux; alternating traces align the left essential u1 datum with the right natural rτ datum.The flux penalty uses the full characteristic matrix, and homogeneous ghost states give dissipative boundary contributions.Source: §4.1 · pp. 12–13

4.3. Energy Stability

  • With continuous material coefficients, exact element inner products, homogeneous boundary states, compatible auxiliary traces, and fully upwind boundary flux, the DG scheme satisfies a semidiscrete energy–dissipation identity and bound.The bound controls energy and cumulative interface, material, and boundary dissipation in terms of the initial state and external forcing.Source: §4.3 · pp. 15, 18, 19
  • Compatible auxiliary traces cancel interior-interface contributions exactly for both alternating LDG and central BR1 choices, leaving only boundary terms from the auxiliary coupling.This cancellation does not imply identical accuracy or convergence of the reconstructed rτh.Source: §4.1 · p. 13, §4.3 · pp. 16–17
  • The energy estimate does not cover inhomogeneous boundary closure or discontinuous material interfaces, which require further analysis or a separate two-material flux.For inhomogeneous data, the added trace-coupled terms are sign-indefinite and cannot be absorbed by the available dissipation.Source: §4.1 · p. 12, §4.3 · p. 19

5. Numerical Results

The numerical studies assess convergence with a compatible manufactured solution and test rotating-beam dynamics against an analytic constant-spin branch. They report near-fourth-order cubic convergence with alternating traces, relaxation toward the branch, and closure of the actuated energy ledger.

5.1. Convergence Tests via the Method of Manufactured Solutions

  • The MMS prescribes an exact solution and forcing while keeping the lower residual block zero and using nonhomogeneous split boundary data, so it verifies code rather than the homogeneous-data energy theorem.All twelve primary and six viscous-resultant components are active, but this construction does not excite the two compatibility coefficients independently.Source: §5.1 · pp. 19–20
  • For k = 3, alternating traces give finest-mesh EOCs of 3.92 for u1 and 3.96 for u2; BR1 u2 falls to 3.09.For k = 2, alternating block EOCs are 2.88 and 2.86 versus BR1’s 3.02 and 3.00; the comparison does not show uniform accuracy gains across every degree and block.Source: §5.1 · p. 21, Figure 1 · p. 21
  • The reported convergence rates are experimental; a rigorous optimal-order a priori error analysis for the full coupled nonlinear system remains unavailable.The alternating-flux proof would require one-sided Gauss–Radau projections matched to the selected traces.Source: §5.1 · pp. 21–22

5.2. Rotational Beam Experiment

  • The rotating-beam experiment ramps root angular velocity from rest, then compares the full damped transient with a reduced straight, axially stretched terminal branch that is steady in the co-rotating frame.The full time-dependent simulation does not impose the reduced branch restriction.Source: §5.2 · p. 22, §5.2.1 · p. 22

5.2.2. Analytic Steady Solution and Critical Speed

  • The reduced steady branch becomes singular as aℓ approaches π/2, and the numerical runs use terminal speeds well below the first critical speed.The dimensionless rotation parameter ωterℓ/ce reaches π/2 at the threshold.Source: §5.2.2 · pp. 22–23

5.2.3. Numerical Simulation Protocol

  • Transient runs use k = 3 on eight uniform cells with characteristic upwinding, alternating auxiliary traces, and a five-stage fourth-order Carpenter–Kennedy method at CFL 0.01.The baseline, doubled-damping, and undamped comparison runs extend to T = 40.Source: §5.2.3 · p. 23

5.2.4. Beam Dynamics and Results

  • After spin-up, the damped beam approaches the analytic branch: baseline discrete energy relaxes toward E⋆ = 264.3321651, with final profiles compared against the branch.Over 35 ≤ t ≤ 40, relative L∞ errors remain below 1.00 × 10−4 in f1 and 6.27 × 10−5 in v2.Source: §5.2.4 · p. 23, Figure 2 · p. 24
  • The accumulated root work of 377.935824 partitions into energy increase 264.324932 and material dissipation 113.610851, with cumulative ledger closure residual 1.10 × 10−7.The net characteristic boundary correction is 4.02 × 10−5 and the free-tip term is 2.91 × 10−8; the characteristic boundary contribution is not additional physical actuator work.Source: §5.2.4 · pp. 24–25

6. Conclusion

The conclusion presents an auxiliary-variable reformulation and an energy-structured DG discretization for Kelvin–Voigt damped intrinsic beams. Numerical tests support near-fourth-order cubic convergence and rotating-branch validation within the stated settings.

6. Conclusion

  • The auxiliary variable represents the rate-dependent sectional force and moment resultants, yielding a first-order-in-time mixed system with damping in an algebraic relation.The hyperbolic part retains the structure of the undamped intrinsic beam equations.Source: §6 · p. 25
  • Characteristic upwinding and alternating LDG traces support a discrete energy estimate controlling energy and dissipation under the stated homogeneous split cantilever closure.The estimate includes Kelvin–Voigt, interior jump, and boundary dissipation in terms of initial data and external forcing.Source: §6 · p. 25
  • Experiments support near-fourth-order cubic convergence, relaxation toward the analytic constant-spin branch, and closure of the actuated energy ledger; undamped runs retain oscillations.Increasing damping accelerates relaxation, while alternating traces are not uniformly more accurate for every degree and block.Source: §6 · p. 25
  • The authors describe the method as a foundation for further beam and rotor-blade simulations, combining high-order accuracy, element locality, and a discrete energy estimate.This is a stated prospective application, not a validation of more realistic coupled models.Source: §6 · p. 25

7. Future Work

Future work is identified in analysis, algorithms, and modeling applications. The authors point to open energy and error analyses, implementation extensions, and broader rotor-blade validation and coupling.

Analysis

  • Analysis remains open for inhomogeneous boundary data, fully discrete energy behavior of practical time integrators, and a complete a priori error theory for the mixed formulation.The boundary-work term can be sign-indefinite; the alternating-flux error route calls for a flux-matched one-sided Gauss–Radau projection.Source: §7 · p. 26

Algorithms

  • Algorithmic extensions include hp-adaptivity, locally implicit or IMEX time integration, and more systematic treatment of stiffness from strong damping or high-order spatial discretizations.The implementation currently uses uniform meshes and a fixed explicit time-stepping strategy.Source: §7 · p. 26

Modelling and applications

  • Future application studies could add anisotropic composite sections, initial twist or pre-curvature, spatially varying materials, 3D finite-element comparisons, and aerodynamic or multi-component coupling.The present numerical experiments use representative configurations for convergence, the analytic steady branch, and energy balance.Source: §7 · p. 26

CRediT authorship contribution statement

The authors’ contribution statement assigns conceptual, methodological, software, validation, analysis, supervision, and writing roles. Appendix A explains the sectional material matrices and states assumptions used for the energy estimates.

CRediT authorship contribution statement

  • Sundriyal is credited with software, validation, visualization, investigation, and formal analysis; the other listed authors are credited across conceptualization, methodology, supervision, and writing.The contribution statement lists individual roles for each of the six authors.Source: CRediT authorship contribution statement · p. 26

Appendix A. Sectional Material Matrices

  • The mass matrix maps cross-sectional velocities to momenta, while flexibility maps elastic resultants to strains and curvatures; the paper assumes mass and flexibility are uniformly symmetric positive definite.The Kelvin–Voigt matrix governs rate-dependent sectional resultants through rτ = Cτ ˙u2, with the effective operator H assumed uniformly symmetric positive definite.Source: Appendix A · pp. 27–28
  • For a homogeneous prismatic isotropic beam aligned with its principal axes and reference centers, extension, shear, torsion, and bending decouple and the flexibility matrix is diagonal.The appendix identifies the area, elastic and shear moduli, principal area moments, and shear correction factors in that special-case form.Source: Appendix A · pp. 27–28

Appendix B. An Energy-Stable DGSEM in One Dimension

Appendix B derives the one-dimensional DGSEM and establishes its quadrature energy stability under the stated hypotheses. Appendices C and D describe the reformulated system and the flux-matrix computations underlying the analysis.

Appendix B. An Energy-Stable DGSEM in One Dimension

  • The DGSEM uses degree-k polynomials at Legendre–Gauss–Lobatto nodes and approximates inner products by LGL quadrature.Mapping elements to a reference element allows the polynomial basis and quadrature weights to be reused.Source: Appendix B · pp. 28–29
  • Under the hypotheses of Theorem 4.2, the LGL-collocated DGSEM satisfies quadrature analogues of the energy–dissipation identity and stability estimate with either alternating LDG or central BR1 interior traces.The result uses the specified physical-boundary flux.Source: Appendix B · p. 29

Appendix C. Reformulation as a Nonlinear Advection–Diffusion System

  • The reformulated equations have a nonlinear hyperbolic–parabolic form: D collects Kelvin–Voigt diffusion, A(u) transport-like gradient terms, and S(u) nonlinear sources.The text notes that this mathematical structure resembles advection–diffusion models despite the absence of physical advection.Source: Appendix C · p. 31

Appendix D. Computations of Flux Matrices

  • Appendix D collects the algebraic derivations for flux splitting, weighted flux matrices, and boundary stability terms used in the discrete energy analysis.Source: Appendix D · p. 32

Appendix D.1. Eigendecomposition

The appendix constructs characteristic flux matrices from a positive-definite matrix eigendecomposition, then derives the absolute flux matrix and its stability-relevant property.

Appendix D.1. Eigendecomposition

  • Positive definiteness of C and M makes Ψ positive definite, so Ψ has an orthogonal eigendecomposition with strictly positive eigenvalues.Source: Appendix D.1 · p. 32

Appendix D.2. Characteristic Flux Splitting

  • The spectral parts A+ and A− represent contributions associated with positive and negative characteristic speeds, respectively.Their sum decomposition recovers A as right- and left-propagating characteristic contributions.Source: Appendix D.2 · pp. 32–33

Appendix D.3. Absolute Flux Matrix

Appendix E. Boundary Contribution Analysis

The boundary analysis evaluates the energy contribution using fully upwind characteristic fluxes and homogeneous ghost states, while noting additional terms for inhomogeneous data.

Appendix E. Boundary Contribution Analysis

  • The boundary energy contribution is evaluated using fully upwind characteristic flux with σ = 1.Source: Appendix E · p. 33
  • The stated boundary analysis uses homogeneous ghost states.Source: Appendix E · p. 33
  • For inhomogeneous data, the computation produces additional terms recorded in Remark 4.3.Source: Appendix E · p. 33

Appendix E.1. Left Boundary (x = 0)

The boundary contractions combine characteristic and physical-flux contributions; together, the characteristic boundary closure contributes nonnegative dissipation to the semidiscrete energy identity.

Appendix E.1. Left Boundary (x = 0)

  • At each endpoint, the analysis combines positive and negative characteristic contributions with a subtracted physical-flux contribution.Source: Appendix E.1 · p. 34, Appendix E.2 · p. 34
  • The endpoint estimates are combined using the boundary operator's right-minus-left orientation.Source: Appendix E.2 · p. 35
  • The characteristic boundary closure contributes nonnegative dissipation to the semidiscrete energy identity.Source: Appendix E.2 · p. 35

Appendix F. Tabulated Convergence Data

Appendix F tabulates average L2 errors and experimental orders of convergence for polynomial degrees k = 1, 2, and 3.

Appendix F. Tabulated Convergence Data

  • Table F.2 reports average L2 errors and experimental orders of convergence for polynomial degree k = 1.The EOC values are computed from the tabulated data.Source: Table F.2 · p. 35, Appendix F · p. 35
  • Table F.3 reports average L2 errors and experimental orders of convergence for polynomial degree k = 2.The EOC values are computed from the tabulated data.Source: Table F.3 · p. 35, Appendix F · p. 35
  • Table F.4 reports average L2 errors and experimental orders of convergence for polynomial degree k = 3.The EOC values are computed from the tabulated data.Source: Appendix F · p. 35, Table F.4 · p. 36
Loading 2609.08397v1…