Quantum computing offers a promising avenue for advancing computational methods in science and engineering. In this work, we introduce the quantum asymptotic numerical method (qANM), a framework for solving nonlinear path-following problems using quantum computing. Based on the principle of high-order perturbation techniques, the proposed method uses Taylor series expansions to transform complex nonlinear systems into sequences of linear equations, which are then solved using quantum linear solvers. The central objective of this study is to demonstrate nonlinear path-following with the required linear systems solved on real quantum hardware. To realize this objective, we develop the quantum-enhanced Jacobi method (q-Jacobi), an iterative quantum linear solver used in the hardware experiment. Numerical simulations on a quantum simulator validate the convergence of the method. A highlight of this work is a proof-of-principle experiment on a superconducting quantum processor. Despite the noise inherent in near-term quantum hardware, the experiment achieves 98
In this study, an algebraically derived consistent tangent matrix and a pseudo-stress for the local residual are applied to elastoplastic problems with damage for shells. The authors define a coupled problem of the global equation, which denotes the weak form of the equilibrium equation for the overall structure, and the local equation, which represents the constraint condition satisfied at each material point. To solve the coupled boundary value problems, the consistent tangent matrix and the pseudo-stress are constructed algebraically by block Newton method. The proposed scheme does not involve any local iterative calculations and allows us to decrease the residuals of the coupled problem simultaneously, since the pseudo-stress for the local residual is incorporated into the linearized global equation. Some numerical examples using shell elements demonstrate the validity and effectiveness of the proposed approach in complex situations with material and geometrical nonlinearities, including the Lemaitre ductile damage model.
The discovery of constitutive laws from experimentally accessible measurements is a central problem in nonlinear computational mechanics. Many data-driven constitutive identification approaches rely either on paired strain–stress data or on full-field displacement measurements, both of which are difficult to obtain in realistic three-dimensional settings. We present a differentiable finite element framework for the discovery of hyperelastic material laws from partial observations, including boundary-only displacement measurements and global reaction forces. The method embeds the nonlinear finite element equilibrium problem directly into the learning loop, so that candidate strain-energy densities are assessed through the deformation fields and reactions they induce. This formulation enforces mechanical equilibrium as a constraint and allows the loss function to be evaluated only at observed locations. To ensure physical admissibility and promote numerical solvability throughout training, the constitutive response is represented by Hyperelastic Neural Networks, a structure-preserving neural class that enforces residual energy and stress-free conditions, frame indifference, isotropic material symmetry, polyconvexity, coercivity, and controlled volumetric growth by construction. The resulting PDE-constrained learning problem is solved using a quasi-Newton strategy combined with continuation and solver-aware backtracking. Numerical experiments in two- and three-dimensional finite elasticity demonstrate accurate recovery of hyperelastic isotropic responses from boundary-only data, robustness to measurement noise, and generalization across geometries, loading conditions, and boundary conditions.
Contact problems in mechanical system dynamics remain challenging due to the combined properties of coupling, nonlinearity, and phase transition in contact laws. In this work, we examine several representative models, including the second-order cone complementarity problem (SOCCP) model, the linear complementarity problem (LCP) polyhedral model, and the LCP cuboid model, focusing on how their complementarity structures represent the underlying physical contact law. Based on this analysis, we propose a nonlinear adaptive pyramid model that preserves the structural simplicity of the LCP cuboid model while recovering the accuracy of the SOCCP model. A two-stage solution procedure is developed, consisting of an outer stage that adaptively aligns a fixed friction pyramid with the ideal friction cone using a set of new state-dependent rescaling rules, and an inner stage that employs a contact-point pivoting solver for the fixed friction pyramid. Numerical results demonstrate that the proposed model and solution procedure accurately capture the combined contact properties while maintaining robustness and efficiency comparable to conventional pivoting solvers for LCP.
Most commercial structural mechanics software restrict access to system matrices, making it difficult to leave the simulation environment or apply external model order reduction (MOR) techniques. At the same time, extensive simulation data from contact-free structural models are already available. In this paper, we propose a novel operator inference-based MOR approach tailored for dynamic contact problems. Our method leverages snapshots from simulations without active contact to construct reduced-order models (ROMs) that capture contact mechanical behavior, without requiring access to full-order system matrices. This enables the reuse of existing data from non-contact simulations to efficiently model and solve dynamical contact problems. To enable the construction of the reduced system matrices and the coupling between contact and interior nodes, a Craig–Bampton-like substructuring method is integrated into the inference process. Characteristic properties of the inferred matrices, such as symmetry and positive definiteness, are enforced by appending additional constraints to the optimization problem. By applying adjoint methods, the dual system can be derived, providing contact forces. Due to the inferred matrix properties, the dual formulation yields a well-defined linear complementarity problem, which can be efficiently solved using pivoting algorithms. The performance of the proposed approach is demonstrated on three-dimensional finite element models involving dynamic contact.
Concurrent optimization of structural topology and material orientation is of great significance for improving the mechanical performance of composite structures. This paper presents a novel discrete-continuous parameterization model for structural topology optimization considering continuous material orientation. Building upon the previous discrete-continuous parameterization (DCP) framework, an improved sigmoid function-based weight interpolation scheme is introduced into the discrete material optimization model, enabling the selection among multiple candidate subintervals using only a single design variable. As a result, the number of design variables is significantly reduced. Furthermore, the convexity of the original optimization problem remains unchanged, because the constructed weight function is nearly linear with respect to the design variable during the initial optimization stage. The new interpolation model is further integrated into the DCP framework to establish a concurrent topology optimization method for structural layout and continuous material orientation. Several numerical examples are provided to verify the effectiveness of the proposed method. The results show that, despite the reduction in design variables, the S-DCP method still achieves favorable optimization solutions for stiffness problems using only two subintervals, which is identical to the original DCP method. This strongly demonstrates the effectiveness of the proposed approach, as it significantly reduces the number of design variables while ensuring that the convexity of the optimization model does not deteriorate.
The numerical simulation of nonlinear structural dynamics remains a challenging task in engineering practice. Conventional A-stable time integration schemes, though widely used, may diverge in highly nonlinear problems and yield unusable results. To overcome this limitation, BN-stable three-substep composite time integration schemes are investigated to ensure robust stability in nonlinear regimes. A comprehensive analysis is conducted for three-stage, third-order singly diagonally implicit Runge–Kutta schemes, with the three-substep composite formulations cast into the SDIRK3 class via identical displacement and velocity iterative forms, leading to the identification of all algebraically stable members and their extension into a family with flexible, controllable numerical dissipation tailored to structural dynamics. It is proved that any third-order BN-stable SDIRK3 scheme is necessarily dissipative, and a non-dissipative member does not exist within this class. Two optimal BN-stable schemes are selected from the resulting family to facilitate computations for large-scale structural systems. Strategies for acceleration output are also developed, addressing a classical shortcoming of Runge–Kutta formulations. The proposed family further includes a fourth-order member and a parameter setting that minimizes period errors. Nonlinear benchmarks confirm unconditional stability and tunable dissipation, and both a one-dimensional wave-propagation problem and simulations of a complex folding rudder with joint clearance demonstrate the applicability of the new methods to real structures.
Nonlinear partial differential equations (PDEs) arise ubiquitously in computational mechanics, requiring efficient and robust solvers for forward simulation and inverse identification. Physics-informed neural networks (PINNs) provide a flexible mesh-free framework, but their practical performance is often limited by slow and unstable training. The random feature method (RFM) offers an efficient alternative by recasting PDE solving as a convex linear least-squares problem, thereby avoiding gradient-based iterative optimization. 4existing RFMs for nonlinear PDEs commonly rely on stabilized Newton-type iterations, introducing substantial computational overhead through repeated Jacobian evaluations and line-search backtracking. In this work, we propose a robust RFM-Picard framework for forward and inverse problems governed by nonlinear PDEs. The method combines a space–time partition of unity with Jacobian-free Picard iteration, efficiently treating nonlinear terms without costly Jacobian-related matrix operations. To improve physical fidelity, an adaptive residual balancing strategy is introduced to strengthen the satisfaction of conservation laws and boundary conditions. The framework is assessed on representative nonlinear PDEs, including the nonlinear Helmholtz equation, the viscous Burgers’ equation, and the incompressible Navier–Stokes equations. Numerical results show that the method outperforms standard PINNs in both accuracy and computational efficiency, and achieves substantially higher efficiency than FEM at comparable accuracy. Compared with stabilized Newton-type RFM solvers, it achieves improved efficiency by eliminating Jacobian assembly and line-search overhead. For inverse problems, it remains robust under extremely sparse observations, enabling accurate reconstruction at computational costs comparable to forward simulations. These results demonstrate that RFM-Picard provides an efficient, robust, and accurate mesh-free RFM method for nonlinear PDEs.
In this paper, the transient topology optimization problem for non-viscously damped system is comparatively investigated. The damping forces of the non-viscously damped system depend on the past history of motion through convolution integrals with suitable kernel functions. Although previous studies confirmed that the order of discretization and differentiation significantly affects the transient topology optimization of viscously damped systems, its effect on non-viscously damped systems is still unclear. Motivated by this, two topology optimization methods for transient response problems of non-viscously damped systems are developed, namely the differentiate-then-discretize adjoint variable method (AVM) and the discretize-then-differentiate AVM. The discretize-then-differentiate approach discretizes via the trapezoidal rule and the Newmark- β method, then differentiates the equations using the AVM. Conversely, the differentiate-then-discretize approach follows the reverse order. The optimized designs, computational efficiency, sensitivity consistency and implementation effort for topology optimization are discussed. Three numerical examples are illustrated to compare the performances of the optimized structures obtained by the two methods. The results indicate that, unlike the viscously damped systems, both methods introduce consistency errors for non-viscously damped systems. However, the discretize-then-differentiate AVM produces smaller consistency errors, leading to better optimized results. Therefore, the discretize-then-differentiate AVM is more recommended for topology optimization of non-viscously damped systems.
This paper presents an improved edge-based smoothed finite element formulation for three-node triangular shell elements to reduce the excessive stiffness commonly observed in conventional triangular elements. Based on the Reissner–Mindlin shell theory, the formulation is developed within an updated Lagrangian framework, and can be applied to both static and explicit dynamic analyses of shell structures involving geometric, material, and contact nonlinearities. In this formulation, each node possesses six degrees of freedom, including a drilling rotation. The membrane strain contribution induced by the drilling rotation is consistently incorporated into the strain field, and the membrane and bending strains are subsequently smoothed over edge-based smoothing domains. To alleviate transverse shear locking, the central point–based discrete shear gap (CP-DSG) technique is adopted. Fundamental tests are conducted to verify the theoretical consistency and spatial isotropy of the formulation. A series of classical benchmark problems are investigated, showing improved accuracy on coarse meshes and good agreement with reference solutions in nonlinear analyses. A vehicle frontal impact experiment is performed to further demonstrate its applicability in engineering analyses. The predicted deformations and acceleration histories are consistent with the experimental measurements, confirming the effectiveness of the proposed formulation for nonlinear analysis of shell structures.
Heterogeneous piezoelectric structures exhibit coupled electromechanical responses that are sensitive to meso-scale material variability. Predicting how such uncertainty propagates to macro-scale quantities of interest remains computationally demanding, especially when direct numerical simulation (DNS) is used to explicitly resolve the underlying mesostructure. This study develops an energy-consistent stochastic DFE^2 –PCE framework for multiscale uncertainty analysis of heterogeneous piezoelectric structures. The proposed method extends the piezoelectric Direct FE^2 ( DFE^2 ) formulation to stochastic settings by introducing randomness among the computational unit cells associated with macro-scale quadrature points. For each realization, the governing problem remains a deterministic two-scale electromechanical analysis, and the macro–meso coupling is preserved through an electric-enthalpy-based Hill–Mandel-consistent formulation, which ensures realization-wise energetic consistency across scales. On this basis, repeated deterministic DFE^2 simulations are used to train a non-intrusive polynomial chaos expansion (PCE) surrogate for efficient approximation of the stochastic response surface. The framework is assessed using porous and multiphase piezoelectric panels, with DNS adopted as the reference model for numerical validation. The results demonstrate that the stochastic DFE^2 model reproduces the DNS responses accurately while requiring substantially fewer degrees of freedom and significantly lower computational cost. The trained PCE surrogates further enable efficient univariate and bivariate parametric analyses, revealing distinct response sensitivities of mechanical and electrical quantities to uncertain elastic and piezoelectric parameters. Under matching sampling assumptions, the proposed stochastic DFE^2 formulation provides a statistically consistent approximation of the corresponding DNS ensemble response. The developed framework thus offers an efficient and physically consistent tool for stochastic multiscale analysis and uncertainty-aware design of heterogeneous piezoelectric structures.
This paper introduces the Enhanced Minkowski Portal Refinement (EMPR) algorithm, a novel contact detection method for convex polygons in discrete element modeling (DEM). The key contributions are: (1) Extending the Minkowski Portal Refinement (MPR) algorithm, originally developed for collision detection in computer graphics, to accurately compute contact features between convex shapes. (2) Identifying and addressing an interior point dependency issue in MPR, which can lead to inconsistent final portals and inaccuracies in contact features and energy conservation. (3) Proposing the EMPR algorithm, which includes a criterion to detect potentially inaccurate final portals and a corrective scheme to ensure precise results. (4) Demonstrating that EMPR outperforms the combined GJK and EPA methods for contact detection in DEM by providing accurate contact states and features in a unified, simple, efficient, and easy-to-implement manner. Numerical examples with various polygonal shapes and particle numbers validate the algorithm's correctness, robustness, and effectiveness.
Classical mesh-based numerical methods, such as the finite element methods via the Galerkin weak formulation, suffer from complex and computationally expensive meshing procedures, especially for irregular geometries with holes, trimmed boundary surfaces and singularities. In contrast, immersed boundary methods use a background regular mesh that does not fit the boundary representation, eliminating classical conforming meshing. A problem arises for cut finite elements due to the demand for numerical integration, the imposition of essential Dirichlet boundary conditions, and the stabilization of the small cut elements or the supports of the used basis functions. These drawbacks are removed with the application of the shifted boundary finite element method where the cut elements are discarded from the analysis, but imposition of boundary conditions is replaced with internal surrogate boundaries. In this paper, we present a novel form of the shifted boundary method applying the procedure based on strong formulation. A system of equations is formed by the collocation method where differential equation is satisfied in internal Greville collocation points on regular background grid, and all the boundary conditions are satisfied only in the collocation points of the internal surrogate boundaries via a Taylor series expansion or an equivalent cut element polynomial up to the order of the Fup (spline) basis functions used. The methodology is demonstrated on two-dimensional Poisson examples, attaining the same convergence rate as other collocation procedures, which proves that a shifted imposition of boundary conditions does not reduce the accuracy or efficiency of the proposed immersogeometric method.
In this study, the thermo-mechanical behavior of a Shape Memory Alloy (SMA) smart beam is investigated using the nonlocal peridynamics (PD) theory. The beam is composed of two thin SMA layers and a substrate. The equations of motion for the nonlocal Timoshenko beam are derived, and axial, transverse, and bending displacements are analyzed. The PD model accurately captures the key characteristics of the SMA material. Using the Brinson constitutive equation, the shape memory effect and pseudo-elasticity of the SMA layer are simulated within the nonlocal PD framework. The stress–strain response of the SMA material is analyzed at various temperatures and validated against experimental results. The beam’s behavior under axial force and bending moment is demonstrated through implicit numerical solutions and compared with numerical and experimental data. The proposed PD model is verified effectively to track SMA characteristics under different stress and temperature conditions.
We present a cable isogeometric analysis method based on the Kirchhoff–Love beam theory coupled with Coulomb friction. In simulation of cable entanglement, twist plays a key role. Isogeometric rotation-free bending-stabilized cable formulations are attractive in that they can incorporate bending effects without introducing additional rotational degrees of freedom. However, since the cross-sectional orientation is entirely determined by the centerline geometry, such formulations inherently lack torsional degrees of freedom. As a result, they are unable to represent twisting behavior and cannot properly account for cross-sectional orientation, and that limits their applicability to problems involving pre-bent configurations or geometries with non-axisymmetric cross sections. This limitation is particularly significant in the presence of frictional contact, where tangential forces acting on the cable surface naturally induce torsional effects. Therefore, a formulation that can account for torsion is essential for accurately capturing the resulting mechanical response. One possible approach in incorporating torsional effects is to employ beam models with rotational degrees of freedom, such as the Kirchhoff–Love beam model, which has been widely used in classical finite element formulations. Isogeometric implementations of such beam models have also been proposed in the literature. In this work, rather than introducing additional complexity into the beam formulation itself, we revisit the classical Kirchhoff–Love theory and provide a concise and consistent formulation within the isogeometric framework. Our focus is on the coupling with Coulomb friction, enabling an accurate representation of contact behavior with torsional effects. We present test computations with contact between a rigid body and an inclined rigid plane and between a flexible cable and an inclined rigid cable. In both problems, the tests cover the frictionless, stick, and slip cases. The results show the accuracy and robustness of the cable isogeometric analysis method presented.
Metal components operating in high-temperature environments are prone to creep damage due to the combined effects of continuous constant stress and thermal loading, which may ultimately lead to structural failure. Conventional creep life prediction models, such as the Larson–Miller parameter and θ‑projection method, are largely empirical, ignoring energy dissipation and micro‑damage mechanisms, which makes them inherently limited for accurate creep failure prediction. To overcome the above limitations, this work develops two models which take thermodynamics as the theoretical basis and regard entropy generation as the core characterization indicator. Specifically, the first‑type model couples energy dissipation with temperature to quantify entropy generation from irreversible processes. The second‑type model uses creep strain as the dominant variable, links macroscopic stress to mesoscopic plastic deformation, and the thermodynamic state index (TSI) is introduced as a failure criterion. Finally, the proposed method is validated using material strain data obtained from high-temperature creep experiments. The results reveal that both models exhibit good predictive performance, and the second-type model achieves better accuracy under high-temperature conditions.
Structural mechanical simulations of short-fiber reinforced plastics are significant for evaluation of noise, vibration and harshness performance. In order to get an overall picture of the structural-dynamic behaviour of different conditioning states the dynamic elasticity constants to describe linear-viscoelasticity must be known. Whilst the dynamic modulus of elasticity can be determined experimentally as a function of frequency and fiber orientation for different temperatures and humidities, the challenge is significantly greater for dynamic, direction-dependent shear moduli and Poisson’s ratios. Based on a two-step homogenisation approach, a dataset was generated by matching homogenised stiffness predictions with experimentally determined directional elastic moduli to derive effective matrix properties and orthotropic elastic constants. Artificial neural networks were trained on the generated material dataset comprising three thermoplastic matrix materials to predict orthotropic elastic and damping properties as a function of fiber orientation, frequency, and environmental conditions, and were integrated into the finite element simulation as a surrogate model of the material behaviour. The prediction quality was evaluated by comparing simulated frequency response functions of bending specimens with experimental data. The results show sufficient correlation for 0° and 90° cutout orientations. For 30° and 45° orientations, the deviations were greater due to shear coupling. A method with constant Poisson’s ratios also showed sufficient prediction quality for samples with shear coupling. The two-step homogenisation approach tended to overestimate stiffness and underestimate damping.
Predicting rare failures in high-dimensional engineering systems remains challenging: existing reliability methods either incur prohibitive computational costs or fail to capture irregular, multimodal failure boundaries. Here we introduce an automated importance sampling framework using mixture von Mises–Fisher–Nakagami distributions that naturally captures directional statistics in high-dimensional spaces. The method employs subset simulation to map failure regions, constructs locally adaptive probability densities without geometric assumptions, and integrates Bayesian machine learning to handle expensive computational models. Across benchmark problems and practical applications in aerospace and semiconductor packaging, the approach achieves 10- to 100-fold efficiency gains over current techniques while eliminating manual parameter tuning. By enabling robust quantification of rare-event probabilities in previously intractable problems, this framework advances reliability-informed design for next-generation engineering systems operating under uncertainty.
Convexification-based approaches to damage mechanics emerged recently as powerful alternatives to non-local or gradient-augmented damage models. Convexifying the incremental potentials leads to a damage-transport equation at material-point level, and a simple modification based on entropic regularization was recently shown to give rise to a continuum damage model which shows softening yet whose incremental potentials are strongly convex, ensuring well-posedness of the boundary-value problem and efficiency of computational resolution. These salient properties come at a price: The damage-transport equation, an advection-type equation, must be solved at each material point. The latter involves a discretization of the attainable damage states and seeks a large number of transition coefficients. The work at hand introduces a novel reformulation of the damage-transport equation as a one-dimensional root-finding problem. Once this problem is solved, the transition coefficients may be obtained in post-processing. We discuss the monotonicity and uniqueness properties of the one-dimensional function whose root is sought, provide a simple yet robust computational strategy for its resolution and study the efficiency of the approach in dedicated computational experiments.
This paper presents a practical guide to building lightweight neural-network surrogates for the plastic return-mapping process in concurrent multiscale (FE2) simulations. Rather than proposing a new architecture, we show how a deliberately simple feed-forward network, structured to mirror the classical return-mapping update, can replace the prohibitively expensive nested fine-scale solves that dominate the cost of conventional FE2 schemes based on FFT homogenization at the meso-scale. We walk through the full workflow: generating training data from incremental homogenization analyses, constructing a compact yet sufficient dataset, embedding material symmetries directly into the mapping, and deploying the trained network as a user-defined material subroutine (UMAT) in a standard finite-element solver—enabling widespread use. A sensitivity study examines the model’s robustness to data density, increment size, and mesh refinement, and we characterize the regimes in which the surrogate holds and where it breaks down. For the macroscopically isotropic, two-dimensional plane-stress setting considered here, the surrogate reproduces the reference response while reducing the per-analysis cost from hours to seconds with speed-ups up to 30,000 over standard FE2. Moreover, by eliminating the need to store and update the extensive microscale state variables, it drastically reduces memory requirements and enables substantially finer finite element discretizations that would otherwise be impossible. The approach extends naturally to three dimensions and to weaker symmetry assumptions, given an appropriate sampling strategy and dataset.
