Search

In numerical simulations such as computational fluid dynamics simulations or finite element analyses, mesh quality affects simulation accuracy directly and significantly. Smoothing is one of the most widely adopted methods to improve unstructured mesh quality in mesh generation practices. Compared with the optimization-based smoothing method, heuristic smoothing methods are efficient but yield lower mesh quality. The balance between smoothing efficiency and mesh quality has been pursued in previous studies. In this paper, we propose a new smoothing method that combines the advantages of the heuristic Laplacian method and the optimization-based method based on the deep reinforcement learning method under the Deep Deterministic Policy Gradient framework. Within the framework, the actor artificial neural network predicts the optimal position of each interior free node with its surrounding ring nodes. At the same time, a critic-network is established and takes the mesh quality as input and outputs the reward of the action taken by the actor-network. Training of the networks will maximize the cumulative long-term reward, which ends up maximizing the mesh quality. Training and validation of the proposed method are presented both on 2-dimensional triangular meshes and 3-dimensional surface meshes, which demonstrates the efficiency and mesh quality of the proposed method. Finally, numerical simulations on perturbed meshes and smoothed meshes are carried out and compared which prove the influence of mesh quality on the simulation accuracy. [ABSTRACT FROM AUTHOR]

In response to the problem that the standard AT2 phase field model cannot effectively model quasi-brittle failure and the existence of length dependence, a new modified phase field model is presented in this paper. By introducing an additional energy, the competing relationship between elastic strain energy and dissipation energy during fracture is changed. A new crack dissipation functional is established using the energy equivalent approach. By introducing a novel rational degradation function, not only can the strength of material failure be effectively utilized, but the model can also reproduce the cohesive softening relationship. A multi-field finite element method is used to discretize the model governing equations, and the equations are solved by an efficient BFGS monolithic algorithm. Finally, some representative numerical examples are used to analyze the effects of parameters in degradation function, length scale and mesh size on the results. The presented numerical simulation results demonstrate length scale and mesh scale independence, and are in good agreement with the experimental results and previous numerical results. At the same time, the numerical results also exhibit cohesive softening properties similar to the current phase field cohesive zone model. These results verify the robustness and effectiveness of the modified phase field model presented in this paper for simulating quasi-brittle failure and brittle fracture. [ABSTRACT FROM AUTHOR]

This paper presents a general Convolution Hierarchical Deep-learning Neural Networks (C-HiDeNN) computational framework for solving partial differential equations. This is the first paper of a series of papers devoted to C-HiDeNN. We focus on the theoretical foundation and formulation of the method. The C-HiDeNN framework provides a flexible way to construct high-order C n approximation with arbitrary convergence rates and automatic mesh adaptivity. By constraining the C-HiDeNN to build certain functions, it can be degenerated to a specification, the so-called convolution finite element method (C-FEM). The C-FEM will be presented in detail and used to study the numerical performance of the convolution approximation. The C-FEM combines the standard C 0 FE shape function and the meshfree-type radial basis interpolation. It has been demonstrated that the C-FEM can achieve arbitrary orders of smoothness and convergence rates by adjusting the different controlling parameters, such as the patch function dilation parameter and polynomial order, without increasing the degrees of freedom of the discretized systems, compared to FEM. We will also present the convolution tensor decomposition method under the reduced-order modeling setup. The proposed methods are expected to provide highly efficient solutions for extra-large scale problems while maintaining superior accuracy. The applications to transient heat transfer problems in additive manufacturing, topology optimization, GPU-based parallelization, and convolution isogeometric analysis have been discussed. [ABSTRACT FROM AUTHOR]

This paper presents a novel plasticity-based formulation for three-dimensional (3D) topology optimization of continuum structures. The proposed formulation addresses the optimization problem by combining mixed rigid-plastic analysis with density-based topology optimization, resulting in a volume minimization approach. Unlike conventional stress-constrained topology optimization methods that rely on linear elastic structural analysis, our developed formulation focuses on enhancing the loading capacity of the designed structures based on the plastic limit theory, leading to more cost-effective designs. To improve computational efficiency, we employ the smoothed finite element technique in our proposed method, enabling the utilization of linear tetrahedral elements for 3D mesh refinement. Moreover, the final formulation of our developed method can be efficiently solved using the advanced primal–dual interior point method, eliminating the need for a separate nonlinear finite element structural analysis. Numerical examples are presented to demonstrate the effectiveness of the proposed approach in offering enhanced design possibilities for continuum structures. [ABSTRACT FROM AUTHOR]

Artificial intelligence (AI) methods have been shown to be effective in aiding topology optimization (TO). This paper proposes an artificial neural network (ANN)-based structural TO method, where the advantages of both the traditional gradient-based methods and the population-based approaches are combined to achieve the converged design or design diversity efficiently. In this method, we adopted the ANN as the structural descriptor for a population of structures that broadly samples the design space and the popular gradient-based solution framework (i.e., the solid isotropic material with penalization, SIMP) to integrate the ANN and finite element analysis of the structures. In such a solution framework, the weights and bias associated with the ANN become the design variables to re-parameterize the density fields of a set of designs used in SIMP but independent of the finite element mesh. These design variables are optimized via the ANN's built-in back-propagation, where a group of points in the design space synergistically evolve under the gradient-based solution framework. Novel loss functions encoding all the sampled structural performances and topology awareness are proposed, including a competitive mechanism-based loss function for converged design and a population diversity-preserving strategy-based loss function to obtain diverse and competitive designs. Correspondingly, a combination of the adjoint method and automatic differentiation algorithm is introduced for the sensitivity analysis. We refer to this method as a multi-point synergistic gradient evolution method. Its efficiency and applicability are demonstrated through 2D and 3D examples. The proposed method is expected to serve as a powerful design tool and further advance the use of AI in TO. [ABSTRACT FROM AUTHOR]

This paper presents an isoparametric tangled finite element method (i-TFEM) method for handling tangled high order/curvilinear meshes. Tangled elements, i.e. elements with negative Jacobian determinant, frequently occur during various stages of analysis and optimization, leading to erroneous results in standard finite element method (FEM). The proposed i-TFEM is an extension of standard FEM to allow for tangled elements. Specifically, a novel variational formulation is proposed that leads to a simple modification of the standard FEM stiffness matrix and additional piece-wise compatibility constraints. Moreover, i-TFEM reduces to the standard FEM for non-tangled (regular) meshes. The accuracy of the proposed i-TFEM is demonstrated for tangled 9-node quadrilateral (Q9) and 6-node triangular (T6) elements. Numerical experiments involving linear and nonlinear elasticity and Poisson problems illustrate that the accuracy and convergence rate of the proposed i-TFEM over a tangled mesh is comparable to that of the standard FEM over a non-tangled mesh. [ABSTRACT FROM AUTHOR]

A new computational framework is developed in this paper for investigating the time-dependent behaviour of concrete including creep, shrinkage and cracking. The developed model aims to explain certain aspects of the time-dependent cracking and creep of concrete that cannot be captured using homogeneous models. The model is based on the scaled boundary finite element method, and it is coupled with a quadtree decomposition algorithm which converts digital images of concrete meso-structures into meshes. Concrete is treated as a two-phase composite which consists of elastic aggregates and mortar that is subjected to time-dependent deformation. The basic creep behaviour is treated as viscoelastic, which is modelled based on a rate-type rheological model corresponding to a Kelvin chain. Drying creep is modelled using a viscous unit which depends on the stress level, and drying shrinkage is stress independent. Both drying creep and drying shrinkage are related to the internal humidity. The humidity distribution within concrete is determined using a diffusion analysis. The moisture movement within mortar is governed by a nonlinear diffusion equation, whereas the aggregates are assumed impermeable. The cracking of concrete is explicitly modelled on the meso-scale through coupling of the continuum damage model for cracking within the mortar phase, and the cohesive zone model for debonding between aggregates and mortar. The proposed model is verified by simulating well-documented experimental studies in the literature. The capability of the proposed model in simulating the time-dependent behaviour of concrete and capturing the crack patterns has also been demonstrated. [ABSTRACT FROM AUTHOR]

This paper presents a new triangular nonlinear shell finite element with a novel kinematic model suitable for simulation with large displacements and rotations, herein introduced as "T6-3iKL". This element has 6 nodes, a quadratic displacement field, and a linear rotation field based on Rodrigues incremental rotation parameters, giving in total 21 degrees of freedom. The novelty of this shell element concerns a new kinematic model with properties from Kirchhoff-Love shell theory, making it possible to eliminate the drilling DOF in the rotation field (compared to Mindlin-Reissner models), approximating the rotation continuity between adjacent elements by a single scalar and allowing multiple branch connections in the mesh, making this element very simple and interesting, with no artificial parameters imputed by the code (such as penalties or Lagrange multipliers). The element permits the implementation of different material constitutive equations, including elastic anisotropic models, and the thickness of the shell is optionally allowed to change through the simulation. The model developed in this article is numerically implemented and the results compared to different references in multiple examples, showing the consistency and reliability of the new formulation. It is believed that this new versatile triangular shell element, with no necessity of artificial penalty calibration, simple kinematics, a relatively small number of DOFs, geometric exactness, the possibility to use 3D material constitutive models, and easy connection with multiple branched shells and beams, implemented together with reliable mesh generation, may be an effective option for shell simulation in many engineering applications. [ABSTRACT FROM AUTHOR]

Studies have found that the widely-used one-pass node-to-segment (NTS) algorithm for contacts fails to pass the contact patch test. Solution errors at the contacting surfaces resulting from the classical NTS algorithm do not necessarily decrease with the mesh refinement. This paper presents a novel NTS with the explicit dynamic formulation in detail, which can pass the contact patch test. The proposed algorithm consists of the following technical ingredient: (1) A weighted smoothing technique is developed to accurately evaluate contact forces for contact problems; (2) The central difference time integration scheme is utilized for efficiency, to save storage space, and to avoid the use of tangent stiffness matrix; (3) the penalty method is used for ease in implementation; (4) Our NTS algorithm is coded in the smoothed finite element method (S-FEM). Two S-FEM models are adopted for the proposed NTS algorithm: the node-based smoothed finite element method (NS-FEM) and the edge-based smoothed finite element method (ES-FEM) so that automatically generable triangular elements can be used for excellent adaptation. Numerical examples of the proposed algorithm are presented to examine the performance. It is found that the present algorithm alleviates the nonuniformity of contact pressure on both flat and curve boundaries. ES-FEM and NS-FEM combined with the new NTS algorithm have higher convergence rates. In addition, the ES-FEM has more accurate displacement and stress solutions at the same grid level. NS-FEM can even provide an upper bound solution in the strain energy for contact problems. [ABSTRACT FROM AUTHOR]

This paper develops two-step methods for solving contact problems with uncertainties. In the first step, we propose stochastic Lagrangian multiplier/penalty methods to compute a set of reduced basis. In the stochastic Lagrangian multiplier method, the stochastic solution is represented as a sum of products of a set of random variables and deterministic vectors. In the stochastic penalty method, the problem is divided into the solutions of non-contact and possible contact nodes, which are represented as sums of the products of two different sets of random variables and deterministic vectors, respectively. The original problems are then transformed into deterministic finite element equations and one-dimensional (corresponding to stochastic Lagrangian multiplier method)/two-dimensional (corresponding to stochastic penalty method) stochastic algebraic equations. The deterministic finite element equations are solved by existing numerical techniques, and the one-/two-dimensional stochastic algebraic equations are solved by a sampling method. Since the computational cost for solving stochastic algebraic equations does not increase dramatically as the stochastic dimension increases, the proposed methods avoid the curse of dimensionality in high-dimensional problems. Based on the reduced basis, we propose semi-reduced order Lagrangian multiplier/penalty equations with two components in the second step. One component is a reduced order equation obtained by smooth solutions of the reduced basis and the other is the full order equation for the nonsmooth solutions. A significant amount of computational cost is saved since the sizes of the semi-reduced order equations are usually small. Numerical examples of up to 100 dimensions demonstrate the good performance of the proposed methods. [ABSTRACT FROM AUTHOR]

In this paper, a novel crack and void identification method is proposed based on the scaled boundary finite element method (SBFEM) with the pixel-based quadtree mesh and the parameter level set method (PLSM), in which the particle swarm optimization method (PSO) is applied to optimize the objective function. It fully exploits the advantages of fast image-based quadtree mesh generation and efficient assembly of global stiffness matrix to accomplish a fast analysis of the forward problem and accelerate the efficiency of the inversion. Furthermore, multiple level set functions are introduced to represent cracks and voids simultaneously. The application of PSO dexterously avoids the difficulty of solving sensitivity. Numerical examples show that the proposed method can identify the position and shape of cracks and voids accurately by discussing different influencing factors. In addition, the method has shown certain anti-noise performance. [ABSTRACT FROM AUTHOR]

In recent years, a lot of progress has been made in the field of material modeling with artificial neural networks (ANNs). However, the following drawbacks persist to this day: ANNs need a large amount of data for the training process. This is not realistic, if real world experiments are intended to be used as data basis. Additionally, the application of ANN material models in finite element (FE) calculations is challenging because local material instabilities can lead to divergence within the solution algorithm. In this paper, we extend the approach of constrained neural network training from [28] to elasto-plastic material behavior, modeled by an incrementally defined feedforward neural network. Purely stress and strain dependent equality and inequality constraints are introduced, including material stability, stationarity, normalization, symmetry and the prevention of energy production. In the Appendices, we provide a comprehensive framework on how to implement these constraints in a gradient based optimization algorithm. We show, that ANN material models with training enhanced by physical constraints leads to a broader capture of the material behavior that underlies the given training data. This is especially the case, if a limited amount of data is available, which is important for a practical application. Furthermore, we show that these ANN models are superior to classically trained ANNs in FE computations when it comes to convergence behavior, stability, and physical interpretation of the results. [ABSTRACT FROM AUTHOR]

Crystal plasticity finite element model (CPFEM) is a powerful numerical simulation in the integrated computational materials engineering toolboxes that relates microstructures to homogenized materials properties and establishes the structure–property linkages in computational materials science. However, to establish the predictive capability, one needs to calibrate the underlying constitutive model, verify the solution and validate the model prediction against experimental data. Bayesian optimization (BO) has stood out as a gradient-free efficient global optimization algorithm that is capable of calibrating constitutive models for CPFEM. In this paper, we apply a recently developed asynchronous parallel constrained BO algorithm to calibrate phenomenological constitutive models for stainless steel 304 L, Tantalum, and Cantor high-entropy alloy. [ABSTRACT FROM AUTHOR]

Physics Informed Neural Networks (PINNs) are shown to be a promising method for the approximation of partial differential equations (PDEs). PINNs approximate the PDE solution by minimizing physics-based loss functions over a given domain. Despite substantial progress in the application of PINNs to a range of problem classes, investigation of error estimation and convergence properties of PINNs, which is important for establishing the rationale behind their good empirical performance, has been lacking. This paper presents convergence analysis and error estimates of PINNs for a multi-physics problem of thermally coupled incompressible Navier–Stokes equations. Through a model problem of Beltrami flow it is shown that a small training error implies a small generalization error. Posteriori convergence rates of total error with respect to the training residual and collocation points are presented. This is of practical significance in determining appropriate number of training parameters and training residual thresholds to get good PINNs prediction of thermally coupled steady state laminar flows. These convergence rates are then generalized to different spatial geometries as well as to different flow parameters that lie in the laminar regime. A pressure stabilization term in the form of pressure Poisson equation is added to the PDE residuals for PINNs. This physics informed augmentation is shown to improve accuracy of the pressure field by an order of magnitude as compared to the case without augmentation. Results from PINNs are compared to the ones obtained from stabilized finite element method and good properties of PINNs are highlighted. [ABSTRACT FROM AUTHOR]

We propose the Convolution Hierarchical Deep-learning Neural Network (C-HiDeNN) that can be tuned to have superior accuracy, higher smoothness, and faster convergence rates like higher order finite element methods (FEM) while using only linear element's degrees of freedom. This is based on our newly developed convolution interpolation theory (Lu et al. in Comput Mech, 2023) and this article focuses on the deep-learning interpretation of C-HiDeNN with graphics processing unit (GPU) programming using JAX library in Python. Instead of increasing the degrees of freedom like higher order FEM, C-HiDeNN takes advantage of neighboring elements to construct the so-called convolution patch functions. The computational overhead of C-HiDeNN is reduced by GPU programming and the total solution time is brought down to the same order as commercial FEM software running on a CPU, however, with orders of magnitude better accuracy and faster convergence rates. C-HiDeNN is locking-free regardless of element types (even with 3-node triangular elements or 4-node tetrahedral elements). C-HiDeNN is also capable of r-h-p-mesh adaptivity like its predecessor HiDeNN (Zhang et al. in Comput Mech 67:207–230, 2021) with additional "a" (dilation parameter) adaptivity that stems from the convolution patch function and "p" adaptivity with higher accuracy and with the same degrees of freedom as that of the linear finite elements. C-HiDeNN potentially has myriad future applications in multiscale analysis, additive and advanced manufacturing process simulations, and high-resolution topology optimization. Details on these applications can be found in the companion papers (Lu et al. 2023; Saha et al. in Comput Mech, 2023; Li et al. in Comput Mech, 2023) published in this special issue. [ABSTRACT FROM AUTHOR]

The paper deals with high-order discontinuous Galerkin (DG) method with the approximation order that exceeds 20 and reaches 100 and even 1000 with respect to one-dimensional case. To achieve such a high order solution, the DG method with finite difference method has to be applied. The basis functions of this method are high-order orthogonal Legendre or Chebyshev polynomials. These polynomials are defined in one-dimensional space (1D), but they can be easily adapted to two-dimensional space (2D) by cross products. There are no nodes in the elements and the degrees of freedom are coefficients of linear combination of basis functions. In this sort of analysis the reference elements are needed, so the transformations of the reference element into the real one are needed as well as the transformations connected with the mesh skeleton. Due to orthogonality of the basis functions, the obtained matrices are sparse even for finite elements with more than thousands degrees of freedom. In consequence, the truncation errors are limited and very high-order analysis can be performed. The paper is illustrated with a set of benchmark examples of 1D and 2D for the elliptic problems. The example presents the great effectiveness of the method that can shorten the length of calculation over hundreds times. [ABSTRACT FROM AUTHOR]

The Material Point Method (MPM) is an established and powerful numerical method particularly useful for simulating large-scale, rapid soil deformations. Therefore, it is often used for the numerical investigation of mass movement hazards such as landslides, debris flows, or avalanches. It combines the benefits of both mesh-free and mesh-based continuum-based discretization techniques by discretizing the physical domain with Lagrangian moving particles carrying the history-dependent variables while the governing equations are solved at the Eulerian background grid, which brings many similarities with commonly used finite element methods. However, due to this hybrid nature, the material boundaries do not usually coincide with the nodes of the computational grid, which complicates the imposition of boundary conditions. Furthermore, the position of the boundary may change at each time step and, moreover, may be defined at arbitrary locations within the computational grid that do not necessarily coincide with the body contour, leading to different interactions between the material and the boundary. To cope with these challenges, this paper presents a novel element-wise formulation to weakly impose non-conforming Dirichlet conditions using Lagrange multipliers. The proposed formulation introduces a constant Lagrange multiplier approximation within the constrained elements in combination with a methodology to eliminate superfluous constraints. Therefore, in combination with simple element-wise interpolation functions classically utilized in MPM (and FEM) to approximate the unknown field, a suitable Lagrange multiplier discretization is obtained. In this way, we obtain a robust, efficient, and user-friendly boundary imposition method for immersed methods specified herein for implicit MPM. Furthermore, the extension to frictionless slip conditions is derived. The proposed methodologies are assessed by comparing the numerical results with both analytical and experimental data to demonstrate their accuracy and wide range of applications. [ABSTRACT FROM AUTHOR]

The hierarchical deep-learning neural network (HiDeNN) (Zhang et al. Computational Mechanics, 67:207–230) provides a systematic approach to constructing numerical approximations that can be incorporated into a wide variety of Partial differential equations (PDE) and/or Ordinary differential equations (ODE) solvers. This paper presents a framework of the nonlinear finite element based on HiDeNN approximation (nonlinear HiDeNN-FEM). This is enabled by three basic building blocks employing structured deep neural networks: (1) A partial derivative operator block that performs the differentiation of the shape functions with respect to the element coordinates, (2) An r-adaptivity block that improves the local and global convergence properties and (3) A materials derivative block that evaluates the material derivatives of the shape function. While these building blocks can be applied to any element, specific implementations are presented in 1D and 2D to illustrate the application of the deep learning neural network. Two-step optimization schemes are further developed to allow for the capabilities of r-adaptivity and easy integration with any existing FE solver. Numerical examples of 2D and 3D demonstrate that the proposed nonlinear HiDeNN-FEM with r-adaptivity provides much higher accuracy than regular FEM. It also significantly reduces element distortion and suppresses the hourglass mode. [ABSTRACT FROM AUTHOR]

This paper presents an immersed, isogeometric finite element framework to predict the response of multi-material, multi-physics problems with complex geometries using locally refined discretizations. To circumvent the need to generate conformal meshes, this work uses an extended finite element method (XFEM) to discretize the governing equations on non-conforming, embedding meshes. A flexible approach to create truncated hierarchical B-splines discretizations is presented. This approach enables the refinement of each state variable field individually to meet field-specific accuracy requirements. To obtain an immersed geometry representation that is consistent across all hierarchically refined B-spline discretizations, the geometry is immersed into a single mesh, the XFEM background mesh, which is constructed from the union of all hierarchical B-spline meshes. An extraction operator is introduced to represent the truncated hierarchical B-spline bases in terms of Lagrange shape functions on the XFEM background mesh without loss of accuracy. The truncated hierarchical B-spline bases are enriched using a generalized Heaviside enrichment strategy to accommodate small geometric features and multi-material problems. The governing equations are augmented by a formulation of the face-oriented ghost stabilization enhanced for locally refined B-spline bases. We present examples for two- and three-dimensional linear elastic and thermo-elastic problems. The numerical results validate the accuracy of our framework. The results also demonstrate the applicability of the proposed framework to large, geometrically complex problems. [ABSTRACT FROM AUTHOR]

In this paper, the field-enriched finite element method is proposed for simulating three-dimensional crack propagation, in which the influences of three-dimensional crack on the physical fields of computational model are described by enriching the field variable on the nodes. Then, through a benchmark example, it is found that under different mesh densities, the convergence rate of the stress intensity factor and energy norm error obtained by the field-enriched finite element method is faster than that of the traditional finite element method. Moreover, the field-enriched finite element method can calculate the stress intensity factor which is basically close to the precision of the extended finite element method and the analytical solution. Finally, comparing with the other numerical results or experimental results on three numerical examples, it is found that the proposed method can effectively and accurately simulate the three-dimensional crack propagation process under different loading conditions. [ABSTRACT FROM AUTHOR]

Computational homogenization is commonly used to predict the responses of composite materials. However, it poses practical issues due to large computational cost especially in the material-by-design setting when various design parameters are to be examined. This paper presents the development of a parametric model order reduction strategy for the micromechanical analysis of composites when fibre distribution is the parameter of interest. The reduced order model is obtained by applying Galerkin projection in combination with proper orthogonal decomposition. The presented framework enables a significantly reduced computational load during parametric studies as the model dimension of the microscale analyses is significantly smaller. The results show that the proposed approach can reproduce the homogenized properties of material and local stress distributions in the microstructures very well. [ABSTRACT FROM AUTHOR]

This paper addresses the need of rapid thermomechanical simulation of metal additive manufacturing by presenting a fully vectorized implementation of predicting the displacement field and residual stress for computation on graphical processing units. The formulation is based on implicit time discretization and the finite element method, where the incremental elastoplastic problem is solved using the conjugate gradient method at each Newton iteration. Sparse representation of algorithmic (tangent) stiffness matrix and the strain–displacement operator, are used in this formulation. A combined hardening plastic model, along with temperature-dependent material properties, is utilized. The temperature field is obtained by conducting detailed part-level thermal simulation using the explicit finite element method and then used in the mechanical simulation to calculate residual stresses. The details of the implementation of the proposed method are provided. Three simulation examples are performed to validate the thermomechanical model, compare the computational efficiency on GPU and CPU, and study the influence of toolpath strategy on residual stresses, respectively. In the example cases, the developed GPU implementation is 10–25 × faster than the CPU version. The success of this development enables fast prediction of residual stress in additive manufacturing to improve the effectiveness of process design and to avoid process defects such as distortion and residual-stress induced fracture. [ABSTRACT FROM AUTHOR]

This paper presents an extended/generalized finite element method for bridging scales in linear elastodynamics in the absence of scale separation. More precisely, the GFEMgl framework is expanded to enable the numerical solution of multiscale problems through the automated construction of specially-tailored shape functions, thereby enabling high-fidelity finite element modeling on simple, fixed finite element meshes. This introduces time-dependencies in the shape functions in that they are subject to continuous adaptation with time. The temporal aspects of the formulation are investigated by considering the Newmark- β time integration scheme, and the efficacy of mass lumping strategies is explored in an explicit time-stepping scheme. This method is demonstrated on representative wave propagation examples as well as a dynamic fracture problem to assess its accuracy and flexibility. [ABSTRACT FROM AUTHOR]

3D Concrete Printing (3DCP) is a rapidly evolving technology that allows for the efficient and accurate construction of complex concrete objects. In this paper, a numerical modelling approach is presented for the simulation of the printing process of cementitious materials, based on the homogeneous fluid assumption. To cope with the large deformations of the domain and the nonlinearity resulting from the use of a non-Newtonian rheological law, the Navier–Stokes equations are solved in the framework of the Particle Finite Element Method (PFEM). Furthermore, tailored solutions have been formulated and implemented for the time-dependent moving boundary conditions at the nozzle outlet and for the efficient handling of the inter-layer contact in the same PFEM framework. The overall computational cost is decreased by the implementation of an adaptive de-refinement technique, which drastically reduces the number of degrees of freedom in time. The proposed modelling approach is finally validated by simulating the printing process of six rectilinear layers and one multi-layer "wall". The results show good agreement with the experimental data and provide valuable insights into the printing process, paving the way for the use of numerical modelling tools for the optimization of materials and processes in the field of 3D Concrete Printing. [ABSTRACT FROM AUTHOR]

We present a new computational framework to simulate the multi-phase convective conjugate heat transfer (CHT) problems emanating from realistic manufacturing processes. The paper aims to address the challenges of boundary-fitted and immersed boundary approaches, which cannot simultaneously achieve fluid-solid interface accuracy and geometry-flexibility in simulating this class of multi-physics systems. The method development is built on a stabilized Arbitrary Lagrangian-Eulerian (ALE)-based finite element thermal multi-phase formulation, which is discretized by overlapping one boundary-fitted mesh and non-boundary-fitted mesh with a quasi-direct coupling approach via Schwarz alternating method. The framework utilizes a volume-of-fluid (VoF)-based multi-phase flow model coupled with a thermodynamics model with phase transitions to capture the conjugate heat transfer between the solid and multi-phase flows and the multi-stage boiling and condensation phenomena. The quasi-direct coupling approach allows the exact and automatic enforcement of temperature and heat-flux compatibility at the fluid-solid interface with large property discontinuities. From the perspective of method development, the proposed framework fully exploits boundary-fitted approach's strength in resolving fluid-solid interface and boundary layers and immersed boundary approach's geometry flexibility in handling moving objects while circumventing each individual's limitations. From the perspective of industry applications, such as water quenching processes, the resulting model can enable accurate temperature prediction directly from process parameters without invoking the conventional empirical heat transfer coefficient (HTC)-based approach that requires intensive calibration. We present the mathematical formulation and numerical implementation in detail and demonstrate the claimed features of the proposed framework through a set of benchmark problems and real-world water quenching processes. The accuracy of the proposed framework is carefully assessed by comparing the prediction with other computational results and experimental measurements. [ABSTRACT FROM AUTHOR]

Dual analysis, requiring a pair of compatible and equilibrated solutions, has been well used for strict error estimation of many problems, including contact problems. Compatible solution is easily obtained by conventional displacement-based finite element method (FEM), while getting the equilibrated solution is not straightforward. To this end, an equilibrium finite element method (EFEM) is developed in this paper for equilibrated solution of contact problems. In doing so, there are two key steps. The first is the establishment of complementary energy principle for contact problems, along with inequality pressure constraint on contact force. In the second step, a traction-based equilibrium element is developed to construct the discrete equilibrated stress field, with which the complementary energy formulation turns into a quadratic programming problem. How to get strict bounds for the discretization error in compatible or equilibrated solution by dual analysis is also addressed. Numerical examples are finally conducted to see the performance and the availability to strict error estimation of the proposed EFEM. [ABSTRACT FROM AUTHOR]

This paper proposes a new shell element formulation using the scaled boundary finite element (SBFE) method. A shell element is treated as a three-dimensional continuum. Its bottom surface is approximated with a quadrilateral spectral element and the shell geometry is represented through normal scaling of the bottom surface. Neumann expansion is applied to approximate the inversions of the matrix polynomials of the thickness coordinate ξ, including the Jacobian matrix and the coefficient of the second-order term in the SBFE equation. This permits the solution along the thickness to be expressed as a matrix exponential function whose exponent is a high-order matrix polynomial of ξ. After introducing the boundary conditions on the top and bottom surfaces and evaluating the resulting matrix exponential via Padé expansion, we derive the element stiffness and mass matrices. Poisson thickness locking is avoided fundamentally. Numerical examples demonstrate the applicability and efficiency of the formulation. [ABSTRACT FROM AUTHOR]

In this paper an enhanced version of the semi-explicit Particle Finite Element Method for incompressible flow problems is presented. This goal is achieved by improving the solution of the advective sub-problem that results of applying the Strang operator splitting to the Navier–Stokes equations. An acceleration term is taken into account in the solution of the advective step and the Stokes problem. The solution of the advetive step is perfomed using a SPH kernel. Two test cases are solved for validating the methodology and estimating its accuracy. The numerical results demonstrate that the proposed scheme improves the accuracy of the semi-explicit PFEM scheme. [ABSTRACT FROM AUTHOR]

The main motivation of this work is to develop a numerical strategy for the shape optimization of a rolling elastic structure under contact with respect to a uniform rolling criterion. A first objective is to highlight the influence on the treatment of the contact terms. To do so, we present a numerical comparison between a penalty-based approach and the use of Nitsche's method which is known to have good consistency properties. A second task concerns the construction of an objective functional to force the uniform rolling criterion. Here, we present and compare two different strategies that will lead to quite similar results. All the numerical experiments proposed in this paper were performed using a fictitious domain approach coupled with a level set representation of the shape and the use of a cut finite element method to approximate the elastic equation. [ABSTRACT FROM AUTHOR]

In this paper, the use of the node-dependent kinematics concept for the geometrical nonlinear analysis of composite one-dimensional structures is proposed With the present approach, the kinematics can be independent in each element node. Therefore the theory of structures changes continuously over the structural domain, describing remarkable cross-section deformation with higher-order kinematics and giving a lower-order kinematic to those portion of the structure which does not require a refinement. In this way, the reliability of the simulation is ensured, keeping a reasonable computational cost. This is possible by Carrera unified formulation, which allows writing finite element nonlinear equilibrium and incremental equations in compact and recursive form. Compact and thin-walled composite structures are analyzed, with symmetric and unsymmetric loading conditions, to test the present approach when dealing with warping and torsion phenomena. Results show how finite element models with node-dependent behave as well as ones with uniform highly refined kinematic. In particular, zones which undergo remarkable deformations demand high-order theories of structures, whereas a lower-order theory can be employed if no local phenomena occur: this is easily accomplished by node-dependent kinematics analysis. [ABSTRACT FROM AUTHOR]

In this two-part paper we begin the development of a new class of methods for modeling fluid-structure interaction (FSI) phenomena for air blast. We aim to develop accurate, robust, and practical computational methodology, which is capable of modeling the dynamics of air blast coupled with the structure response, where the latter involves large, inelastic deformations and disintegration into fragments. An immersed approach is adopted, which leads to an a-priori monolithic FSI formulation with intrinsic contact detection between solid objects, and without formal restrictions on the solid motions. In Part I of this paper, the core air-blast FSI methodology suitable for a variety of discretizations is presented and tested using standard finite elements. Part II of this paper focuses on a particular instantiation of the proposed framework, which couples isogeometric analysis (IGA) based on non-uniform rational B-splines and a reproducing-kernel particle method (RKPM), which is a meshfree technique. The combination of IGA and RKPM is felt to be particularly attractive for the problem class of interest due to the higher-order accuracy and smoothness of both discretizations, and relative simplicity of RKPM in handling fragmentation scenarios. A collection of mostly 2D numerical examples is presented in each of the parts to illustrate the good performance of the proposed air-blast FSI framework. [ABSTRACT FROM AUTHOR]

The present paper explores the solution of a heat conduction problem considering discontinuities embedded within the mesh and aligned at arbitrary angles with respect to the mesh edges. Three alternative approaches are proposed as solutions to the problem. The difference between these approaches compared to alternatives, such as the eXtended Finite Element Method (X-FEM), is that the current proposal attempts to preserve the global matrix graph in order to improve performance. The first two alternatives comprise an enrichment of the Finite Element (FE) space obtained through the addition of some new local degrees of freedom to allow capturing discontinuities within the element. The new degrees of freedom are statically condensed prior to assembly, so that the graph of the final system is not changed. The third approach is based on the use of modified FE-shape functions that substitute the standard ones on the cut elements. The imposition of both Neumann and Dirichlet boundary conditions is considered at the embedded interface. The results of all the proposed methods are then compared with a reference solution obtained using the standard FE on a mesh containing the actual discontinuity. [ABSTRACT FROM AUTHOR]

In this paper, a new finite element method, termed the generalized modal element method (GMEM), is proposed. In GMEM, the element stiffness is derived by decomposing element deformation patterns into individual element generalized modes, where different methods are used to construct the generalized modes. Specifically, three different modal construction methods, including analytical method, assumed displacement method and traditional finite element technique, are proposed for developing the element generalized modes. The concept of modal local coordinate systems is also proposed to ensure the element frame invariance when using polynomial displacement functions, which successfully enables one to use the analytical solutions derived from governing differential equations to develop high accuracy element formulations. An asymmetric hexahedral solid element and a symmetric hexahedral solid element are subsequently derived by using GMEM. The displacement functions of the elemental 24 generalized modes are expressed in terms of Cartesian coordinates so that the element behavior is independent of mesh distortions. Furthermore, the first 21 generalized modes are derived from analytical method making the element capable of avoiding common locking phenomena. Several benchmark problems are performed to demonstrate the accuracy and performance of the new element formulations in linear static and frequency analysis. [ABSTRACT FROM AUTHOR]

Meshfree methods such as the Reproducing Kernel Particle Method and the Element Free Galerkin method have proven to be excellent choices for problems involving complex geometry, evolving topology, and large deformation, owing to their ability to model the problem domain without the constraints imposed on the Finite Element Method (FEM) meshes. However, meshfree methods have an added computational cost over FEM that come from at least two sources: increased cost of shape function evaluation and the determination of adjacency or connectivity. The focus of this paper is to formally address the types of adjacency information that arises in various uses of meshfree methods; a discussion of available techniques for computing the various adjacency graphs; propose a new search algorithm and data structure; and finally compare the memory and run time performance of the methods. [ABSTRACT FROM AUTHOR]

In this paper, we propose a new homogenization algorithm, virtual clustering analysis (VCA), as well as provide a  mathematical framework for the recently proposed self-consistent clustering analysis (SCA) (Liu et al. in Comput Methods Appl Mech Eng 306:319-341, 2016). In the mathematical theory, we clarify the key assumptions and ideas of VCA and SCA, and derive the continuous and discrete Lippmann-Schwinger equations. Based on a key postulation of “once response similarly, always response similarly”, clustering is performed in an offline stage by machine learning techniques (k-means and SOM), and facilitates substantial reduction of computational complexity in an online predictive stage. The clear mathematical setup allows for the first time a convergence study of clustering refinement in one space dimension. Convergence is proved rigorously, and found to be of second order from numerical investigations. Furthermore, we propose to suitably enlarge the domain in VCA, such that the boundary terms may be neglected in the Lippmann-Schwinger equation, by virtue of the Saint-Venant’s principle. In contrast, they were not obtained in the original SCA paper, and we discover these terms may well be responsible for the numerical dependency on the choice of reference material property. Since VCA enhances the accuracy by overcoming the modeling error, and reduce the numerical cost by avoiding an outer loop iteration for attaining the material property consistency in SCA, its efficiency is expected even higher than the recently proposed SCA algorithm. [ABSTRACT FROM AUTHOR]

In the stochastic finite element analysis of irregular wall structures considering material uncertainties, the random fields simulation and deterministic finite element analysis (FEA) are the two focuses. In this paper, we present an efficient stochastic finite element analysis procedure for irregular wall structures with inelastic random field properties. In the procedure, the geometry domains of irregular wall structures are deemed as two-dimensional (2D) manifolds. Then Isometric feature mapping (Isomap) is used to map the manifold to a 2D Euclidean domain, over which the well-developed stochastic harmonic function representation is applied to simulate the random field. Meanwhile, to accurately reflect the nonlinear behaviors of the irregular wall structures, we adopt an enhanced deterministic FEA method, which combines the multi-layered shell element, the softened damage-plasticity concrete model and the quasi-Newton solution with a two-level secant stiffness in a framework. Finally, the proposed approach is applied to the stochastic analysis of a U-shaped reinforced concrete shear wall to illustrate its feasibility and applicability. The results demonstrate that the proposed approach can effectively simulate the random fields over irregular geometry domains and can reproduce the representative stochastic inelastic behaviors of irregular wall structures. The approach can be further used in reliability or safety evaluation of structures. [ABSTRACT FROM AUTHOR]

This paper introduces a new concept called self-updated finite element (SUFE). The finite element (FE) is activated through an iterative procedure to improve the solution accuracy without mesh refinement. A mode-based finite element formulation is devised for a four-node finite element and the assumed modal strain is employed for bending modes. A search procedure for optimal bending directions is implemented through deep learning for a given element deformation to minimize shear locking. The proposed element is called a self-updated four-node finite element, for which an iterative solution procedure is developed. The element passes the patch and zero-energy mode tests. As the number of iterations increases, the finite element solutions become more and more accurate, resulting in significantly accurate solutions with a few iterations. The SUFE concept is very effective, especially when the meshes are coarse and severely distorted. Its excellent performance is demonstrated through various numerical examples. [ABSTRACT FROM AUTHOR]

Neural networks (NN) have been studied and used widely in the field of computational mechanics, especially to approximate material behavior. One of their disadvantages is the large amount of data needed for the training process. In this paper, a new approach to enhance NN training with physical knowledge using constraint optimization techniques is presented. Specific constraints for hyperelastic materials are introduced, which include energy conservation, normalization and material symmetries. We show, that the introduced enhancements lead to better learning behavior with respect to well known issues like a small number of training samples or noisy data. The NN is used as a material law within a finite element analysis and its convergence behavior is discussed with regard to the newly introduced training enhancements. The feasibility of NNs trained with physical constraints is shown for data based on real world experiments. We show, that the enhanced training outperforms state-of-the-art techniques with respect to stability and convergence behavior within FE simulations. [ABSTRACT FROM AUTHOR]

In this paper, we consider parallelization of the loop over elements using OpenMP in FEAP (Taylor, ), which is a research FE code, very popular at universities. Even for a serial version of FEAP (a cluster version also exists) such a parallelization is a non-trivial task due to the existing architecture of this code, which complicates efficient parallelization. First, we compare the serial version of FEAP to the parallel code Warp3D (Dodds et al., ), considering the usage of time and memory. As we found, Warp3D is much faster but uses more memory than FEAP. An analysis of Warp3D helps us to devise our method of parallelization of the loop over elements. Next, we describe several changes in FEAP, which were necessary to parallelize the loop over elements using OpenMP. In particular, the subroutine assembling elemental matrices is identified as crucial to good performance, and several directives for the mutual exclusion synchronization of OpenMP are implemented and tested. Finally, we demonstrate the performance of the parallelized FEAP, designated as ompFEAP, on numerical examples involving 3D and shell elements of FEAP as well as user's elements. We conclude that ompFEAP, using the directive ATOMIC for synchronization of the assembling, provides a very good speedup and efficiency. [ABSTRACT FROM AUTHOR]

This paper proposes a method of solving 3D large deformation frictional contact problems with the Cartesian Grid Finite Element Method. A stabilized augmented Lagrangian contact formulation is developed using a smooth stress field as stabilizing term, calculated by Zienckiewicz and Zhu Superconvergent Patch Recovery. The parametric definition of the CAD surfaces (usually NURBS) is considered in the definition of the contact kinematics in order to obtain an enhanced measure of the contact gap. The numerical examples show the performance of the method. [ABSTRACT FROM AUTHOR]

Stresses in gear contact simulations, performed using elastic multibody systems, are recovered. A single gear pair is used for stress investigations and an impact is chosen as simulation case representing an extremely dynamical situation. The gears are modeled as a reduced elastic multibody system allowing a fast computation of the dynamical problem. Depending on the projection matrix which is used for model order reduction, stresses can sometimes not be recovered accurately throughout the whole gear. Thus, the main focus in this paper lies on the selection of the functions which make up the projection matrix and, therefore, determine the elastic deformations and the quality of recovered stresses. However, the chosen set of modes does not only affect stress calculation, it also strongly affects the computation of the dynamics of the gear system and, thus, the computational effort, and may lead to serious drawbacks. This issue is discussed, too. Upon that, several different mode sets are analyzed trying to minimize the computational effort of the elastic multibody system for the given problem while still being able to recover accurate stress values on distinct geometric areas. The stress values are compared with a finite element reference computation. The novel contribution of this paper is the determination of a minimal set of modes including ones assigned to the nodes of the gear contact surface, which are able to accurately recover stresses but minimize the numerical drawbacks. [ABSTRACT FROM AUTHOR]

This paper deals with the modelling and the prediction of the dynamic instabilities for a rubbing system. Two hybrid approaches are introduced for dynamic instability analysis and applied to a reduced disc brake system. The methods are based on stochastic algorithms coupled with the finite element method (FEM) using the complex eigenvalue analysis technique. By considering the input parameters as random variables, the uncertainty analysis is performed through two approaches to predict the unstable frequencies of a braking system (1) Monte-Carlo (MC) using Mersenne-Twister (MT19937) algorithm and (2) periodic sampling technique. Since the mechanism of brake squeal involves many design parameters, stochastic finite element approaches will be coupled with sensitivity algorithms, e.g. Variance-Based Sensitivity Analysis and Fourier Sensitivity Analysis Test, to analyze the contribution of each random variable on the dynamic instabilities. First, a comparison between the two stochastic algorithms is performed on standard analytical models. The objective is to validate the accuracy and to assess the numerical efficiency that FAST presents to (1) propagate the uncertainties upstream of the model and (2) to compute the partial variances of the model output. Secondly, the coupling of the previous stochastic algorithms with FEM is carried out and tested through a reduced brake system consisting of a rotating disc with two flat pads. Results show that the hybrid approach FAST-FE is more robust and computationally more efficient compared to the widely used MC-FE for these types of problems. FAST-FE solver converges, within a reasonable computing time, either to approximate the probability density function of the random variables or to compute the partial variances of the dynamic instabilities. Hence, it can be considered as an efficient numerical method for squeal instability analysis in order to reduce squeal noise of such a mechanical system. [ABSTRACT FROM AUTHOR]

This paper proposes a description of a granular medium as a stochastic heterogeneous continuum medium. The heterogeneity of the material properties field recreates the heterogeneous stress field in a granular medium. The stochastic approach means that only statistical information, easily available, is required to construct the model. The heterogeneous continuum model is Calibrated with respect to discrete simulations of a set of railway ballast samples. As they are continuum-based, the equilibrium equations can be solved on a large scale using a parallel implementation of an explicit time discretization scheme for the Finite Element Method. Simulations representative of the influence on the environment of the passage of a train on a ballasted railway track clearly show the influence of the heterogeneity. These simulations seem to correlate well with previously unexplained overly damped measurements in the free field. [ABSTRACT FROM AUTHOR]

This paper presents an adaptive strategy for phase-field simulations with transition to fracture. The phase-field equations are solved only in small subdomains around crack tips to determine propagation, while an extended finite element method (XFEM) discretization is used in the rest of the domain to represent sharp cracks, enabling to use a coarser discretization and therefore reducing the computational cost. Crack-tip subdomains move as cracks propagate in a fully automatic process. The same mesh is used during all the simulation, with an h-refined approximation in the elements in the crack-tip subdomains. Continuity of the displacement between the refined subdomains and the XFEM region is imposed in weak form via Nitsche's method. The robustness of the strategy is shown for some numerical examples in 2D and 3D, including branching and coalescence tests. [ABSTRACT FROM AUTHOR]

Fully coupled global equations are proposed for enhancing the performance of Finite Element analysis of unsaturated soils. The governing equation describing mechanical equilibrium is formulated in terms of net stress, and in the mass conservation equation the contribution of this net stress in determining the change of degree of saturation is also included. The novelty of this paper is the development of new global finite element equations that can be used to find an approximate solution to these governing equations. The new equations have a mechanical term appearing in the flow matrix that is additional to the usual hydraulic term. This is in contrast to previous studies in which the coupling matrices ignore this effect. A performance study has been conducted for undrained footing problems, which shows that the additional mechanical term appearing in the flow matrix has a large influence on the accuracy of the numerical results. [ABSTRACT FROM AUTHOR]

In this paper we alleviate limitations of the conventional embedded reinforcement formulation for applications with dense fiber contents. We demonstrate that by condensing the fiber degrees of freedom during the assembly stage, the condition number of the resultant stiffness matrix is effectively reduced and the use of iterative solvers is facilitated even without preconditioning. Numerical benchmarks consisting of very large numbers of high aspect ratio discrete fibers are performed, and major advantages are reported in terms of the computational efficiency of the solution method. We apply the solution method to a set of examples in the modeling of the fiber reinforced composites, specifically the estimation of the homogenized mechanical properties of the discontinuous fiber composites and modeling of the compression molding process. In the latter case, a specimen containing more than 2 million discrete fibers is analyzed. [ABSTRACT FROM AUTHOR]

The paper is devoted to non-linear vibrations of plates, made of the Zener viscoelastic material modelled with the Caputo fractional derivative, and in particular to their response to harmonic excitation. The plate geometric non-linearity is of the von Kármán type. In the formulation shear effects and rotary inertia are considered, too. The problem is solved in the frequency domain. A one-harmonic form of the solution for plate displacements corresponding to the plate formulation is assumed. The amplitude equation is obtained from the time averaged principle of virtual work. The time averaging precedes the use of the harmonic balance method. In the space discretization the finite element method is used involving 8-noded rectangular plate elements with selective-reduced integration. Several numerical examples are analyzed and the response curves are found using a path-following method. The purpose of these analyses is to identify material features of the adopted model of viscoelasticity with the fractional derivative. [ABSTRACT FROM AUTHOR]

In this two-part paper we begin the development of a new class of methods for modeling fluid-structure interaction (FSI) phenomena for air blast. We aim to develop accurate, robust, and practical computational methodology, which is capable of modeling the dynamics of air blast coupled with the structure response, where the latter involves large, inelastic deformations and disintegration into fragments. An immersed approach is adopted, which leads to an a-priori monolithic FSI formulation with intrinsic contact detection between solid objects, and without formal restrictions on the solid motions. In Part I of this paper, the core air-blast FSI methodology suitable for a variety of discretizations is presented and tested using standard finite elements. Part II of this paper focuses on a particular instantiation of the proposed framework, which couples isogeometric analysis (IGA) based on non-uniform rational B-splines and a reproducing-kernel particle method (RKPM), which is a Meshfree technique. The combination of IGA and RKPM is felt to be particularly attractive for the problem class of interest due to the higher-order accuracy and smoothness of both discretizations, and relative simplicity of RKPM in handling fragmentation scenarios. A collection of mostly 2D numerical examples is presented in each of the parts to illustrate the good performance of the proposed air-blast FSI framework. [ABSTRACT FROM AUTHOR]

This paper develops an accurate and efficient finite element model for simulating coupled transient electromagnetic and dynamic mechanical fields that differ widely in the frequency ranges. This coupled modeling framework is necessary for effective modeling and simulation of structures such as antennae that are governed by multi-physics problems operating in different frequency and temporal regimes. A key development is the wavelet transformation induced multi-time scaling or WATMUS method that is designed to overcome shortcomings of modeling coupled multi-physics problems that are governed by disparate frequencies. The WATMUS-based FE model is enhanced in this paper with a scaled and preconditioned Newton-GMRES solver for efficient solution. Results from the WATMUS-based FE model show the accuracy and highly improved computational efficiency in comparison with single time-scale methods. The coupled FE model is used to solve two different antenna problems with large electromagnetic to mechanical frequency ratios. The examples considered are a monopole antenna and a microstrip patch antenna. Comparing the electromagnetic fields with the progression of mechanical cycles demonstrate complex multi-physics relations in these applications. [ABSTRACT FROM AUTHOR]

This paper focuses on the development and engineering applications of a new resolvent sampling based Rayleigh-Ritz method (RSRR) for solving large-scale nonlinear eigenvalue problems (NEPs) in finite element analysis. There are three contributions. First, to generate reliable eigenspaces the resolvent sampling scheme is derived from Keldysh's theorem for holomorphic matrix functions following a more concise and insightful algebraic framework. Second, based on the new derivation a two-stage solution strategy is proposed for solving large-scale NEPs, which can greatly enhance the computational cost and accuracy of the RSRR. The effects of the user-defined parameters are studied, which provides a useful guide for real applications. Finally, the RSRR and the two-stage scheme is applied to solve two NEPs in the FE analysis of viscoelastic damping structures with up to 1 million degrees of freedom. The method is versatile, robust and suitable for parallelization, and can be easily implemented into other packages. [ABSTRACT FROM AUTHOR]