Physics-informed neural networks may struggle to represent hyperbolic conservation laws with steep gradients and shocks, whereas purely data-driven latent reduced models can accumulate rollout errors. We propose a two-stage framework that combines Latent Dynamics Networks with a physics-informed correction in latent space. In Stage I, we learn a globally trained parametric latent-dynamics model and a coordinate-based decoder from full-order trajectories. In Stage II, the pretrained dynamics network and decoder are frozen, while a separate temporal corrector is optimized for each unseen parameter instance. The corrector is constrained by a latent-consistency term, which anchors it to the pretrained rollout, and by the governing PDE residual evaluated on the decoded physical field, together with an optional regularization term designed to discourage latent collapse. The resulting method is therefore a trajectory-wise, test-time refinement of a reusable parametric latent surrogate. Five benchmarks involving the Burgers, Buckley–Leverett, and compressible Euler equations, including a two-dimensional problem and a coupled hyperbolic system, show that the correction reduces the prediction error of the uncorrected latent model across held-out parameter values within the sampled parameter domains. A representative temporal-extrapolation test further indicates reduced rollout drift and improved shock tracking beyond the Stage-I training horizon. A physics-informed Stage-I baseline achieves comparable accuracy in the considered test, indicating that the two-stage formulation should be interpreted as an optional per-instance refinement rather than as a universally superior alternative. We finally discuss the limitations of pointwise residual enforcement near discontinuities and possible extensions based on weak, entropy-consistent, and globally parameterized correction strategies.
This work proposes a hybrid computational framework that integrates the Virtual Element Method (VEM) with neural network techniques to efficiently solve three-dimensional hyperelastic problems. The framework combines graph neural networks (GNNs), VEM, and physics-informed neural networks (PINNs), and incorporates a transfer learning strategy to improve the model generalization ability. Owing to its ability to operate on arbitrary polyhedral meshes, VEM provides a flexible and accurate discretization for hyperelastic problems involving complex geometries. Building on this property, the proposed framework employs VEM to construct high-fidelity physics constraints, utilizes GNNs to encode the topology of polyhedral meshes, and adopts PINNs to solve the governing equations under these constraints. Furthermore, the introduced transfer learning strategy enables rapid adaptation of the trained model to new material parameters and loading conditions. Numerical experiments demonstrate that the proposed framework achieves accuracy comparable to that of traditional high-fidelity numerical simulations while improving computational efficiency by one to two orders of magnitude. In addition, the method exhibits strong robustness and generalization capability under variations in parameters and loading scenarios. These results indicate that the proposed approach provides a promising and efficient solution for the rapid analysis of multi-scenario and multi-parameter three-dimensional hyperelastic problems.
This study proposes a novel discrete hygro-thermo-mechanical model for concrete under freeze-thaw cycles (FTCs) within the framework of the Multi-physics Lattice Discrete Particle Model. The evolution of ice content and its hysteretic behavior are simulated using a modified liquid-solid interfacial energy function. The non-uniform local eigenstrains induced by temperature variation and ice formation govern the FTC-induced cracking behavior; accordingly, these eigenstrains are incorporated into the mechanical constitutive model via a one-way coupling scheme. Notably, the model successfully captures cracking patterns, which initiate at the surface and propagate inwards over hundreds of FTCs. It also accurately predicts the associated degradation in both tensile and compressive strength. Therefore, this study finds that an external compressive load exerts a complex influence on FTC-induced degradation. Specifically, a moderate compressive load equal to 50% of the material’s compressive strength increases the residual strength by approximately 1.5 MPa after 50 FTCs. This beneficial effect arises because the applied external load restrains frost-heaving deformation and thereby suppresses the initiation of FTC-induced micro-cracks. In contrast, when the external load exceeds 70% of the compressive strength, degradation is significantly accelerated.
In response to external stimuli, hydrogels undergo large deformation coupled with solvent diffusion. The chemomechanical coupling produces transient, strongly nonlinear responses that can trigger instabilities and bifurcations. Modeling such responses using conventional Finite Element Methods (FEMs) often breaks down due to significant mesh distortion. Here, we present an implicit, mixed material point method (MPM) grounded in non-equilibrium thermodynamics that remains robust under large distortions while consistently coupling deformation and solvent diffusion. Cut-cell ill-conditioning near boundaries is removed using extended B-splines that interpolate boundary degrees of freedom from interior ones via Lagrange polynomials, and subdivision-stabilized interpolation of displacement and chemical potential provides an inf–sup stable discretization. Essential boundary conditions are imposed weakly through a symmetric Nitsche formulation on boundary material points. Together, these improvements yield a symmetric, well-conditioned tangent stiffness suitable for eigenvalue-based stability and bifurcation analysis. We demonstrate the stability and accuracy of the method on benchmark problems, including one-dimensional constrained swelling, swell-induced buckling of a hydrogel column, and bifurcation and surface creasing of a hydrogel with cylindrical pores. The present MPM methodology provides a robust computational foundation for studying extreme, diffusion-coupled deformations in soft-matter applications.
This paper presents a theoretical analysis of a fourth-order exponential time differencing Runge–Kutta rational approximation scheme for the Allen–Cahn equation based on real and distinct poles, aiming to establish a theoretical framework for the method combining exponential time differencing with rational approximation. Starting from the scalar linearized equation, we rigorously prove that the scheme is L-stable and that its matrix norm is stable with three types of boundary conditions, revealing that the L-acceptable property of the RDP rational approximation is the key to ensuring stability. Regarding the error analysis, we establish that the fully discrete scheme attains fourth-order temporal convergence, which theoretically confirms that the scheme achieves fourth-order accuracy. Numerical examples in both 2D and 3D verify the accuracy and effectiveness of the method.
High fidelity analysis of nonlinear flutter in supersonic panels with complex configurations generally requires high-dimensional nonlinear dynamical models, whose direct simulation is prohibitively expensive. To reduce computational cost, this paper develops an equation-driven model reduction framework tailored to nonlinear panel flutter. Because the proposed reduction directly operates on the system matrices and nonlinear coefficient tensors of the full-order models (FOMs), a numerically stable finite element discretization is essential. To this end, a nonlinear Hermite-interpolation-based Bogner–Fox–Schmit element (HBE) is developed. By constructing explicit shape functions, the HBE inherently eliminates the numerical instability that traditional bicubic-interpolation-based Bogner–Fox–Schmit elements typically suffer during mesh refinement. Based on the spectral submanifold (SSM) theory, a master spectral subspace selection strategy incorporating both unstable and coupled stable modes is then formulated. The dimensionless dynamic pressure is embedded into an augmented system as a parameter coordinate, yielding a parameter-dependent mixed-mode SSM. The SSM and its reduced dynamics are constructed directly from the finite element governing equations without using FOM trajectories as training data. This framework reduces finite element models with thousands of degrees of freedom (DOFs) to reduced-order models (ROMs) with only a few coordinates while retaining the dominant nonlinear dynamics. Furthermore, since the dynamic pressure parameter is explicitly embedded in the reduced dynamics, highly efficient simulations under varying dynamic pressure conditions can be achieved via a single offline reduction. Finally, numerical examples demonstrate the computational accuracy and efficiency of the proposed framework, and the effects of various stiffening schemes on the critical dynamic pressure and limit cycle oscillation (LCO) amplitude are investigated.
Surrogate models are central to scientific machine learning, where they enable fast prediction, simulation, inference, and control for complex physical systems. For time-dependent problems, however, accurate interpolation of training trajectories is not sufficient: reliable surrogates should also respect the conservation laws, invariants, admissibility conditions, and dissipative structures that give those trajectories physical meaning. We introduce Physics-conforming Latent Twins, a framework for learning latent surrogate solution operators whose dynamics satisfy selected physical principles by design. The method builds on the Latent Twin formulation by jointly learning an encoder, a decoder, and a latent flow map between arbitrary time-indexed states, while constraining the latent dynamics to preserve or dissipate prescribed structural quantities. We develop a constraint-transfer viewpoint that separates pullback compatibility from latent conformity, connecting physical structure in the original state space with enforceable constraints in latent space, and prove structure-preservation bounds showing how latent enforcement improves control of physical defects after decoding. We also derive algebraic conditions for latent flow maps that preserve linear and quadratic invariants or enforce dissipative inequalities. Numerical experiments on representative ODE and PDE benchmarks demonstrate improved constraint satisfaction, structural fidelity, and qualitative long-time behavior while maintaining accurate surrogate prediction.
Numerical simulations of solids undergoing dynamic fragmentation, a problem characterized by dynamic fracture and dense contacts, require accurately capturing the transition from a solid continuum to a collection of interacting fragments. We use the finite-element method with the extrinsic cohesive zone model for fracture. For contact, conventional penalty-based methods often exhibit numerical instabilities in dynamic collision-rich settings. To address this, we adapt and validate a novel semi-explicit time-integration scheme: the Nonsmooth Newmark-β (NSN) method for unilateral contact. Based on the Nonsmooth Contact Dynamics (NSCD) method, this formulation strongly enforces contact constraints at the velocity level. Within this scheme, the bulk dynamics are non-impulsive and integrated explicitly with second-order accuracy, and the fracture model allows displacement discontinuities and integrates the cohesive softening explicitly. Contact, treated rigorously as nonsmooth, is integrated implicitly with first-order accuracy. Benchmark tests demonstrate that the NSN scheme achieves accuracy comparable to established nonsmooth methods, such as the semi-explicit CD-Lagrange and implicit Moreau–Jean schemes. Moreover, it outperforms penalty-based approaches by orders of magnitude. Although the NSN method incurs a higher per-step computational cost, its enhanced stability allows for significantly larger time steps. Consequently, for 1D benchmarks, overall computational efficiency is comparable to or better than that of purely explicit approaches. We applied this framework to 1D fragmentation under free and confined expansion. Results reveal that confinement shifts the fracture energy budget from local fragment kinetic energy to larger-scale global system kinetic energy. Additionally, we found, counterintuitively, that, compared to fully elastic contact, adding contact dissipation reduces fracture energy yet increases the final fragment count. It occurs because such dissipation reduces the vibration within damaged fragments, allowing cleaner stress-wave propagation and better damage localization, driving cracks to full separation rather than distributing damage. These results establish the NSN scheme as a robust tool for generating high-fidelity fragmentation statistics.
Inspired by the vectorial lattice Boltzmann method for linear elastodynamics (Boolakee et al., 2025), we construct a total-Lagrangian vectorial lattice Boltzmann formulation for two-dimensional finite-strain hyperelastic dynamics. The governing equations are first written as a conservative first-order system for the material velocity and the full deformation gradient. This representation separates the kinematic part of the dynamics from the constitutive closure: the first Piola–Kirchhoff stress is evaluated locally from the current deformation gradient and enters the lattice only through nonlinear flux moments. A D2Q4 stencil with six-component vector populations is then used to match the state and the two material-coordinate fluxes. The formulation includes a second-order population initialization, trapezoidally centered body forcing, displacement reconstruction by velocity quadrature, and half-way reconstructions for velocity Dirichlet and Neumann traction boundaries on grid-aligned domains. The resulting method preserves the local collide-stream structure of standard lattice Boltzmann schemes while adapting the vectorial first-order strategy from linear elastodynamics to hyperelastic finite-strain dynamics.
We present an implicit, fully-coupled hydro-mechanical solver for the three-dimensional simulation of fluid-driven rupture propagation along pre-existing discontinuities. The solver simultaneously handles frictional slip and tensile failure along arbitrary intersecting fractures and faults in a linearly elastic and impermeable rock matrix. Spatial discretization combines a displacement discontinuity boundary element method with a Galerkin finite element method for pore-fluid pressure diffusion. Frictional and tensile failure are governed by a poro-elastoplastic interface law incorporating slip-weakening friction, dilatancy, and tensile strength degradation. Block preconditioning of the coupled tangent system ensures robustness across a wide range of fracture behaviors, including friction and tensile hydraulic failure. Solver accuracy is verified – for the first time – against a comprehensive suite of semi-analytical rupture propagation solutions of increasing complexity: fluid-driven frictional ruptures, dilatant ruptures with permeability changes, and penny-shaped hydraulic fractures spanning the viscosity-to-toughness transition. Two multi-fracture examples further demonstrate the solver’s capabilities: injection into intersecting fractures, and a hydraulic fracture intersecting a strike-slip fault. These highlight the ability of the algorithm to capture frictional slip, dilatancy, permeability evolution, and tensile opening within a unified framework, making it well suited for fluid-driven rupture simulation in faulted and fractured rocks.
Extrapolative prediction of complex nonlinear dynamics remains a central challenge in engineering. This study proposes a one-shot learning method to identify global frequency-response curves from a single excitation time history by learning governing equations. We introduce MEv-SINDy (Multi-frequency Evolutionary Sparse Identification of Nonlinear Dynamics) to infer the governing equations of non-autonomous and multi-frequency systems. The methodology leverages the Generalized Harmonic Balance (GHB) method to decompose complex forced responses into a set of slow-varying evolution equations. We validated the capabilities of MEv-SINDy on two critical Micro-Electro-Mechanical Systems (MEMS). These applications include a nonlinear beam resonator and a MEMS micromirror. Our results show that the model trained on a single point accurately predicts softening/hardening effects and jump phenomena across a wide range of excitation levels. This approach significantly reduces the data acquisition burden for the characterization and design of nonlinear microsystems.
We consider an elastic model for a circular arch that incorporates membrane, transverse shear, and bending effects. The central line of the arch is partitioned into elements, and an ultra-weak variational formulation is developed alongside a discontinuous Petrov–Galerkin (DPG) approximation procedure based on so-called optimal test functions. The formulation uses discontinuous stress and displacement interpolations on the element mesh, with corresponding interface variables defined at the nodes. Theoretical analysis establishes well-posedness and quasi-optimal convergence properties of the DPG approximation, while also revealing potential error amplification influenced by the curvature of the arch and the imposed boundary conditions. The method is tested on examples with different support configurations. The numerical experiments confirm the theoretical predictions and further demonstrate that the accuracy of the DPG method can be improved by employing a suitably scaled test space norm.
This work aims at presenting a novel Stabilization-free Virtual Element method for 3D problems. In particular, we focus on the discretization of the linear elastic equation, but the new projection operator introduced can be applied to derive a self-stabilized formulation also for general scalar elliptic equations. The method is introduced in its lowest order formulation and for the analysis we consider the class of polyhedra with triangular faces, typically called deltahedra. We provide a sufficient condition on the polynomial projection space that implies the well-posedness. Several numerical tests assess the robustness of the method and confirm the theoretical convergence rates. Furthermore, we test the proposed method on a nonlinear elasticity problem, to show the robustness of the proposed self-stabilized formulation for solving nonlinear problems.
We present a systematic method for exactly enforcing Dirichlet, Neumann, and Robin type conditions on general quadrilateral domains with arbitrary curved boundaries. Our method is built upon exact mappings between general quadrilateral domains and the standard domain, and employs a combination of TFC (theory of functional connections) constrained expressions and transfinite interpolations. When Neumann or Robin boundaries are present, especially when two Neumann (or Robin) boundaries meet at a vertex, it is critical to enforce exactly the induced compatibility constraints at the intersection, in order to enforce exactly the imposed conditions on the joining boundaries. We analyze in detail and present constructions for handling the imposed boundary conditions and the induced compatibility constraints for two types of situations: (i) when Neumann (or Robin) boundary only intersects with Dirichlet boundaries, and (ii) when two Neumann (or Robin) boundaries intersect with each other. We describe a four-step procedure to systematically formulate the general form of functions that exactly satisfy the imposed Dirichlet, Neumann, or Robin conditions on general quadrilateral domains. The method developed herein has been implemented together with the extreme learning machine (ELM) technique we have developed recently for scientific machine learning. Ample numerical experiments are presented with several linear/nonlinear stationary/dynamic problems on a variety of two-dimensional domains with complex boundary geometries. Simulation results demonstrate that the proposed method has enforced the Dirichlet, Neumann, and Robin conditions on curved domain boundaries exactly, with the numerical boundary-condition errors at the machine accuracy.
A stabilisation-free hybrid virtual triangular element for the analysis of two-dimensional elasticity problems is presented. The formulation is developed within a hybrid framework in which the kinematic field is defined along the element boundary, while the stress field is independently approximated within the element domain through an energy-consistent projection procedure. The stress approximation employs a divergence-free, iso-stable polynomial basis, i.e., one in which the number of stress parameters matches the number of deformational kinematic modes. The basis isspecifically constructed to avoid spurious energy modes, thereby eliminating the need for additional stabilisation terms.The proposed element possesses nine degrees of freedom, namely two translational components and one drilling rotation at each vertex. The introduction of drilling rotations as independent kinematic variables within a hybrid virtual element framework constitutes a novel aspect of the formulation.The resulting element, denoted as HVT3-6, is assessed through standard benchmark problems in two-dimensional elasticity. The numerical results demonstrate optimal convergence rates of order h2 for displacement, stress, and complementary energy measures. Moreover, the proposed formulation yields improved accuracy compared to a displacement-based triangular finite element with the same number of degrees of freedom on both structured and unstructured meshes, including coarse discretisations.Finally, while the present formulation is developed for triangular elements, the proposed drilling-based hybrid framework constitutes a basis for future extensions of hybrid virtual elements with independent drilling rotations to general polygonal meshes.
In this paper, we develop and analyze an unconditionally stable fully discrete method for the diffusive-viscous wave equation, combining an hp-version continuous Petrov–Galerkin time discretization with an hp-conforming finite element method in space. We first reformulate the problem as a first-order system to facilitate the time-stepping construction; the combination of the hp-version temporal and spatial discretizations then yields simultaneous high-order accuracy in both time and space. Unconditional stability is established via an energy argument that employs a special linear weight function, and the discrete energy is proved to be monotonically non-increasing at the time nodes. A rigorous hp-version a priori error analysis yields optimal convergence rates in the L2(L2) and L2(H1) norms, with constants explicitly independent of the spatio-temporal discretization parameters. Extensive numerical experiments confirm the theoretical results, demonstrating optimal algebraic convergence under h-refinement and exponential convergence under p-refinement. Further numerical examples validate the method’s long-time stability, energy dissipation properties, and its ability to handle both homogeneous and heterogeneous media with discontinuous coefficients.
The construction of analysis-suitable spline parameterizations for multiply-connected planar domains remains a challenging task in isogeometric analysis, particularly when a single-patch representation is desired for superior smoothness and analysis performance. Although prevalent multi-patch approaches address geometric complexity through domain decomposition, they inevitably introduce artificial interfaces, induce singularities at patch junctions and require additional continuity constraints that compromise smoothness and convergence. Moreover, the general applicability of existing decomposition techniques for multiply-connected planar domains remains unclear, and their extension to three-dimensional cases is even more challenging. This paper introduces a novel and coherent framework for generating high-quality single-patch parameterizations of multiply-connected planar domains. Given an input multiply-connected domain, our method first establishes a quasi-conformal map onto a punctured unit square—a canonical parametric domain containing several interior square holes. Both the quasi-conformal map and the geometric configuration (centers and side lengths) of these inner squares are treated as unknowns. We then develop an alternating optimization algorithm that concurrently determines the quasi-conformal map and the layout of the interior squares. The final parameterization in the physical space is obtained by approximating the inverse map using truncated hierarchical B-splines. Extensive numerical experiments, ranging from simple annular geometries to high-genus engineering components, confirm that our approach produces high-quality parameterization results and significantly outperforms the conventional multi-patch alternatives. However, we also acknowledge that the current framework does not explicitly enforce geometric symmetry in the generated parameterizations—a limitation that may compromise the accuracy and reliability of subsequent simulations. Potential remedies are discussed for future work.
Water droplet erosion (WDE) compromises the structural integrity and efficiency of power-generation components, driving extensive research to elucidate its underlying mechanisms. Although finite element (FE) simulations provide insight into impact-induced stress and strain fields, modeling WDE onset and progression remains challenging because conventional FE formulations rely on classical continuum governing equations whose spatial derivatives become invalid at discontinuities. Peridynamics (PD), formulated without spatial derivatives, has proven effective for predicting complex damage patterns and offers a promising framework for investigating WDE mechanisms. Nevertheless, PD research has largely focused on solids, and existing PD-based fluid models mainly target incompressible and weakly compressible flows. Accordingly, this study develops and validates a two-dimensional Peridynamic Differential Operator (PDDO)-based formulation for compressible high-speed water droplet impact on a rigid wall. The PDDO-based formulation employs a horizon-based meshfree discretization of the classical Lagrangian Navier–Stokes equations and incorporates Roe's approximate Riemann solver with a dissipation limiter to improve numerical stability near sharp flow discontinuities. A GPU-accelerated MATLAB implementation incorporates free-surface treatment, a fictitious-layer wall boundary, particle regularization, and a boundary-collision correction. The model captures the key physical features of high-speed droplet impact and agrees well with semi-analytical predictions of maximum and spatially averaged wall pressure over impact velocities of 100–500 m/s, while demonstrating droplet-size independence, accurate early-stage contact-radius evolution, and good agreement of predicted lateral-jet velocities with experimental measurements. The present PDDO framework provides a foundation for future coupling with PD solid formulations, enabling a unified horizon-based computational framework for fluid–structure interaction simulations of WDE.
