This study presents a gradient-based shape optimization over a fixed mesh using a NURBS-based Interface-enriched Generalized Finite Element Method (NIGFEM), applicable to multi-material structures. In the proposed method, NURBS are used to parameterize the design geometry precisely and compactly by a small number of design variables. An analytical shape sensitivity analysis is developed to compute derivatives of the objective and constraint functions with respect to the design variables. Subtle, but important new terms involve the sensitivity of shape functions and their spatial derivatives. Verification and illustrative problems are solved to demonstrate the precision and capability of the method. This article is protected by copyright. All rights reserved.

When computing the homogenized response of a Representative Volume Element (RVE), a popular choice is to impose periodic boundary conditions on the RVE. Despite their popularity, it is well known that standard periodic boundary conditions lead to inaccurate results if cracks or localization bands in the RVE are not aligned with the periodicity directions. A previously proposed remedy is to use modified strong periodic boundary conditions, that are aligned with the dominating localization direction in the RVE. In the present work, we show that alignment of periodic boundary conditions can also conveniently be performed on weak form. Starting from a previously proposed format for weak micro-periodicity that does not require a periodic mesh, we show that aligned weakly periodic boundary conditions may be constructed by only modifying the mapping (mirror function) between the associated parts of the RVE boundary. In particular, we propose a modified mirror function, that allows alignment with an identified localization direction. This modified mirror function corresponds to a shifted stacking of RVEs, and thereby ensures compatibility of the dominating discontinuity over the RVE boundaries. The proposed method leads to more accurate results compared to using unaligned periodic boundary conditions, as demonstrated by the numerical examples. This article is protected by copyright. All rights reserved.

Fourier solvers have become efficient tools to establish structure-property relations in heterogeneous materials. Introduced as an alternative to the Finite Element (FE) method, they are based on fixed-point solutions of the Lippmann-Schwinger type integral equation. Their computational efficiency results from handling the kernel of this equation by the Fast Fourier Transform (FFT). However, the kernel is derived from an auxiliary homogeneous linear problem, which renders the extension of FFT-based schemes to non-linear problems conceptually difficult. This paper aims to establish a link between FE- and FFT-based methods, in order to develop a solver applicable to general history- and time-dependent material models. For this purpose, we follow the standard steps of the FE method, starting from the weak form, proceeding to the Galerkin discretization and the numerical quadrature, up to the solution of non-linear equilibrium equations by an iterative Newton-Krylov solver. No auxiliary linear problem is thus needed. By analyzing a two-phase laminate with non-linear elastic, elasto-plastic, and visco-plastic phases, and by elasto-plastic simulations of a dual-phase steel microstructure, we demonstrate that the solver exhibits robust convergence. These results are achieved by re-using the non-linear FE technology, with the potential of further extensions beyond small-strain inelasticity considered in this paper. This article is protected by copyright. All rights reserved.

A computationally efficient numerical model that describes carbon sequestration in deep saline aquifers is presented. The model is based on the multiphase flow and vertically averaged mass balance equations, requiring the solution of two partial differential equations - a pressure equation and a saturation equation. The saturation equation is a non-linear advective equation for which the application of Galerkin Finite Element Method (FEM) can lead to non-physical oscillations in the solution. In this article we extend three stabilized FEM formulations, which were developed for uncoupled systems, to the governing non-linear coupled PDEs. The methods developed are based on the Streamline Upwind (SU), the Streamline Upwind / Petrov-Galerkin (SUPG), and the Least Squares Finite Element Method (LSFEM). Two sequential solution schemes are developed: a single step (SS) and a predictor-corrector (PC). The range of Courant numbers yielding smooth and oscillation-free solutions is investigated for each method. The useful range of Courant numbers found depends upon both the sequential scheme (SS vs PC) and also the time integration method used (Forward Euler, Backward Euler, or Crank-Nicolson). For complex problems such as when two plumes meet, only the SU stabilization with an amplified stabilization parameter gives satisfactory results when large time steps are used. This article is protected by copyright. All rights reserved.

Response sensitivity is an essential component to understanding the complexity of material and geometric nonlinear finite element formulations of structural response. The direct differentiation method (DDM), a versatile approach to computing response sensitivity, requires differentiation of the equations that govern the state determination of an element and produces accurate and efficient results. The DDM is applied to a force-based element formulation that utilizes curvature-shear based displacement interpolation (CSBDI) in its state determination for material and geometric nonlinearity in the basic system of the element. The response sensitivity equations are verified against finite difference computations and a detailed example shows the effect of parameters that control flexure-shear interaction for a stress resultant plasticity model. The developed equations make the CSBDI force-based element available for gradient-based applications such as reliability and optimization where efficient computation of response sensitivities are necessary for convergence of gradient-based search algorithms. This article is protected by copyright. All rights reserved.

A mixture theory based model for multi-constituent solids is presented where each constituent is governed by its own balance laws and constitutive equations. Interactive forces between constituents that emanate from maximization of entropy production inequality provide the coupling between constituent specific balance laws and constitutive models. The deformation of multi-constituent mixtures at the Neumann boundaries requires imposing inter-constituent coupling constraints such that the constituents deform in a self-consistent fashion. A set of boundary conditions is presented that accounts for the non-zero applied tractions, and a variationally consistent method is developed to enforce inter constituent constraints at Neumann boundaries in the finite deformation context. The new method finds roots in a local multiscale decomposition of the deformation map at the Neumann boundary. Locally satisfying the Lagrange multiplier field and subsequent modeling of the fine scales via edge bubble functions results in closed-form expressions for a generalized penalty tensor and a weighted numerical flux that are free from tunable parameters. The key novelty is that the consistently derived constituent coupling parameters evolve with material and geometric nonlinearity, thereby resulting in optimal enforcement of inter-constituent constraints. Various benchmark problems are presented to validate the method and show its range of application.

The singular boundary method (SBM) is a recent strong-form boundary collocation method which uses a linear combination of the fundamental solution of the governing equation to approximate the field variables. Due to its full interpolation matrix, the SBM solution encounters the high computational complexity and storage requirement that limit its applications to large-scale engineering problems. This paper presents a way to overcome this drawback by introducing the diagonal form fast multipole method (FMM). A diagonal form fast multipole singular boundary method (DF-FMSBM) is then developed to reduce the computational operations of the SBM with direct solvers from *O*(*N ^{3}*) to

As a typical high dimensional nonlinear dynamic problem, the difficulties of the dynamic analysis on the spatial flexible damping beam mostly result from the coupling between the spatial motion and the transverse vibration. Considering the coupling effect and the weak structure damping, the dynamic behaviors of the spatial flexible beam are investigated by a complex structure-preserving method in this paper. Based on the variational principle, the dynamic model of the spatial flexible damping beam is established which can be decoupled into two parts approximately, one controls the spatial motion and another mainly controls the transverse vibration. For the first part, the classic fourth order Runge-Kutta method can be used to discrete it expediently. For the latter part, based on the generalized multi-symplectic idea, the approximate symmetric form as well as the generalized multi-symplectic conservation law is formulated and a fifteen-point scheme equivalent to the Preissmann scheme is constructed. Numerical iterations are performed for six typical initial cases between the two parts to study the dynamic behaviors of the spatial flexible beam. In the numerical experiments, the effects of the damping factor and the initial conditions (including the initial radial velocity and the initial attitude angle) on the dynamic behaviors of the spatial flexible beam are investigated in detail. From the numerical results, it can be concluded that: the damping effect on the long-time dynamic behaviors of the spatial flexible beam could not be neglected even if it is weak; the numerical method proposed in this paper owns the tiny numerical dissipation as well as the excellent long-time numerical stability.

An Arlequin poromechanics model is introduced to simulate the hydro-mechanical coupling effects of fluid-infiltrated porous media across different spatial scales within a concurrent computational framework. A two-field poromechanics problem is first recast as the two-fold saddle point of an incremental energy functional. We then introduce Lagrange multipliers and compatibility energy functionals to enforce the weak compatibility of hydro-mechanical responses in the overlapped domain. To examine the numerical stability of this hydro-mechanical Arlequin model, we derive a necessary condition for stability, the two-fold inf–sup condition for multi-field problems, and establish a modified inf–sup test formulated in the product space of the solution field. We verify the implementation of the Arlequin poromechanics model through benchmark problems covering the entire range of drainage conditions. Through these numerical examples, we demonstrate the performance, robustness, and numerical stability of the Arlequin poromechanics model. This article is protected by copyright. All rights reserved.

In the present paper strong form finite elements are employed for the free vibration study of laminated arbitrarily shaped plates. In particular, the stability and accuracy of three different Fourier expansion-based differential quadrature techniques are shown. These techniques are used to solve the partial differential system of equations inside each computational element. The three approaches are called Harmonic Differential Quadrature (HDQ), Fourier Differential Quadrature (FDQ) and Improved Fourier expansion-based Differential Quadrature (IFDQ) methods. IFDQ method implements auxiliary functions in order to approximate functional derivatives up to the fourth order, with respect to FDQ method that has a basis made of sines and cosines. All the present applications are related to literature comparisons and the presentation of new results for further investigation within the same topic. A study of such kind has never been proposed in the liteature and it could be useful as a reference for future investigation in this matter.

For the response analysis of engineering systems with uncertain-but-bounded parameters, three uncertain models have been considered according to the available probability distribution information. One is the bounded random model in which the uncertain parameters are well defined with sufficient probability distribution information and described as bounded random variables. The second one is the interval model in which the uncertain parameters are expressed as interval variables without giving any information of probability distribution. The last one is the bounded hybrid uncertain model which includes both bounded random variables and interval variables. On the basis theory of Gegenbauer polynomial approximation theory, a unified *Interval and Random Gegenbauer Series Expansion*(IRGSE) Method is proposed and extended for the response prediction of three uncertain models of acoustic fields with uncertain-but-bounded parameters. In IRGSE, the uncertain-but-bounded variables with different probability distribution information, including interval variables and bounded random variables with different *Probability Density Functions* (PDFs), are transformed into the function of unitary variables defined on associated with the corresponding polynomial parameter(

) of *Gegenbauer series expansion*(GSE). The coefficients of GSE is calculated by Gauss-Gegenbauer integration method. By using IRGSE, the responses of three uncertain acoustic models are approximated uniformly by GSE, through which the interval and random analysis can be easily implemented by many numerical solvers. Two numerical examples are applied to investigate the effectiveness of the proposed method.

Lobatto elements are high-order or spectral elements featured by non-equispaced Lobatto nodes and the Lobatto nodal quadrature. The tensor product nature of the element interpolation allows the derivatives at the node be computed thru the derivatives along the nodal line which involve only 1D interpolations. In this paper, generic schemes for formulating assumed natural strain (ANS) and stabilized Lobatto Lagrange C^{0} plate/shell elements of order ≥ 2 are presented. Two ANS schemes are devised. Both schemes sample the normal membrane and transverse shear natural strain components at the 1D reduced-order Gaussian quadrature points along the nodal lines. The difference is that the first and second schemes sample the membrane shear natural strain component at the Lobatto nodes and the 2D reduced-order Gaussian quadrature points, respectively. Meanwhile, all the other natural strain components in both schemes can be obtained by 1D interpolation along the nodal lines. In the stabilization scheme for reduced-integrated elements, only five stabilization vectors are required regardless of the element order. The new elements outperform the standard Lobatto element except when membrane-dominated curved thin shell problems with strong boundary layer effect are considered by coarse meshes.

Shape Memory Polymers (SMPs) belong to a class of smart materials which have shown promise for a wide range of applications. They are characterized by their ability to maintain a temporary deformed shape and return to an original parent permanent shape. In this paper, we consider the coupled photo-mechanical behavior of Light Activated Shape Memory Polymers (LASMPs), focusing on the numerical aspects for finite element simulations at the engineering scale. The photo-mechanical continuum framework is summarized, and some specific constitutive equations for LASMPs are described. Numerical implementation of the multiphysics governing partial differential equations takes the form of a user defined element subroutine within the commercial software package Abaqus. We verify our two- and three-dimensional finite element procedure for multiple analytically tractable cases. To show the robustness of the numerical implementation, simulations are performed under various geometries and complex photo-mechanical loading. This article is protected by copyright. All rights reserved.

The hyper-reduced-order model (HROM) is proposed for the thermal calculation with a constant moving thermal load. Firstly, the constant velocity transient process is simplified to a steady-state process in the moving frame. Secondly, the control volume is determined by the temperature rate, and the thermal equilibrium equation in the moving frame is derived by introducing an advective term containing the loading velocity. redThirdly, the HROM is performed on the control volume with a moving frame formulation. This HROM has been applied to the thermal loading on brick and ring disk specimens with a CPU gain of the order of 7 (10^{7}). red In addition, two strategies are proposed for the HROM to improve its precision. Moreover, the high efficiency and high accuracy are kept for the parametric studies on thermal conductivity and amplitude of heat flux based on the developed HROM. This article is protected by copyright. All rights reserved.

A mean-strain 10-node tetrahedral element (T10MS) is developed for the solution of geometrically nonlinear solid mechanics problems using the concept of energy-sampling stabilization. A uniform-strain tetrahedron for applications to linear elasticity was recently described. The formulation as extended here is able to solve large-strain hyperelasticity. The present 10-node tetrahedron is composed of several four-node linear tetrahedral elements, four tetrahedra in the corners and four tetrahedra that tile the central octahedron in three possible sets of four-node tetrahedra, corresponding to three different choices for the internal diagonal. We formulate a mean-strain element with stabilization energy evaluated on the four corner tetrahedra, which is shown to guarantee consistency and stability. The stabilization energy is expressed through a stored-energy function, and contact with input parameters in the small-strain regime is made. The neo-Hookean model is used to formulate the stabilization energy. As for small-strain elasticity, the stabilization parameters are determined by actual material properties and geometry of tetrahedra without any user input. The numerical tests demonstrate that the present element performs well for solid, shell, and nearly incompressible structures. Copyright © 2016 John Wiley & Sons, Ltd.

This paper details a semi-analytical procedure to efficiently integrate the product of a smooth function and a complex exponential over tetrahedral elements. These highly oscillatory integrals appear at the core of different numerical techniques. Here, the Partition of Unity Method (PUM) enriched with plane waves is used as motivation. The high computational cost or the lack of accuracy in computing these integrals is a bottleneck for their application to engineering problems of industrial interest. In this integration rule, the non-oscillatory function is expanded into a set of Lagrange polynomials. In addition, Lagrange polynomials are expressed as a linear combination of the appropriate set of monomials, whose product with the complex exponentials is analytically integrated, leading to 16 specific cases that are developed in detail. Finally, we present several numerical examples to assess the accuracy and the computational efficiency of the proposed method, compared to standard Gauss-Legendre quadratures. Copyright © 2016 John Wiley & Sons, Ltd.

In this work we propose a homogenization formulation to model transient heat conduction in heterogeneous media that takes into account thermal inertia contributions which arise from a finite description of the microscale. Rewriting the variational form of the transient heat conduction problem and making use of key assumptions, we arrive at a mathematical formulation that suggests an extension of the Hill-Mandel principle when considering non-null heat flux divergence in the RVE. Along the manuscript we highlight that the main results of the proposed formulation are in agreement with recent advances in the field of computational homogenization applied to transient mechanical and heat flow problems. The proposed extension of the Hill-Mandel principle contributes to the understanding of the microscale thermal inertia effects incorporation into the multiscale framework. We also present the calculations needed for implementing the model and numerical results which give support to the theoretical model developed. The numerical results highlight the importance of considering full transient aspects when dealing with multiscale heat conduction in heterogeneous media which are subjected to high thermal gradients. Copyright © 2016 John Wiley & Sons, Ltd.

The parametric analysis of electric grids requires carrying out a large number of Power Flow computations. The different parameters describe loading conditions and grid properties. In this framework, the Proper Generalized Decomposition (PGD) provides a numerical solution explicitly accounting for the parametric dependence. Once the PGD solution is available, exploring the multidimensional parametric space is computationally inexpensive. The aim of this paper is to provide tools to monitor the error associated with this significant computational gain and to guarantee the quality of the PGD solution. In this case, the PGD algorithm consists in three nested loops that correspond to 1) iterating algebraic solver, 2) number of terms in the separable greedy expansion and 3) the alternated directions for each term. In the proposed approach, the three loops are controlled by stopping criteria based on residual goal-oriented error estimates. This allows one for using only the computational resources necessary to achieve the accuracy prescribed by the end-user. The paper discusses how to compute the goal-oriented error estimates. This requires linearizing the error equation and the Quantity of Interest to derive an efficient error representation based on an adjoint problem. The efficiency of the proposed approach is demonstrated on benchmark problems. Copyright © 2016 John Wiley & Sons, Ltd.

The present work introduces an efficient technique for the deformation of block-structured grids occurring in simulations of fluid-structure interaction (FSI) problems relying on large-eddy simulations (LES). The proposed hybrid approach combines the advantages of the inverse distance weighting (IDW) interpolation with the simplicity and low computational effort of transfinite interpolations (TFI), while preserving the mesh quality in boundary layers. It is an improvement over the state-of-the-art currently in use.

To reach this objective, in a first step three elementary mesh deformation methods (TFI, IDW and radial basis functions) are investigated based on several test cases of different complexities analyzing their capabilities but also their computational costs. That allows to point out the advantages of each method but also demonstrates their drawbacks. Based on these specific properties of the different methods, a hybrid methodology is suggested which splits the entire grid deformation into two steps: First the movement of the block-boundaries of the block-structured grid and second the deformation of each block of the grid. Both steps rely on different methodologies, which allows to work out the most appropriate method for each step leading to a reasonable compromise between the grid quality achieved and the computational effort required. Finally, a hybrid IDW-TFI methodology is suggested which best fits to the specific requirements of coupled FSI-LES applications. This hybrid procedure is then applied to a real-life FSI-LES case. Copyright © 2016 John Wiley & Sons, Ltd.

This paper presents a discontinuous thick level set model for modeling delamination initiation and propagation in composites. Interface elements are widely applied in delamination models to define a discontinuity between layers of a laminate. In this paper the common damage definition in the constitutive relation of interface elements is replaced with a new definition developed using the thick level set approach. Following this approach a band of damage with predefined length is considered where the damage is defined as a function of distance to the damage front. The specifications of suitable damage functions for the developed method are investigated and an efficient damage function is introduced. A sensitivity analysis of numerical input parameters is performed which proves that the model is not sensitive to the length of the damage band. Since the required element size is linked to the length of the damage band, the insensitivity to this length provides freedom to use coarser meshes. Furthermore, the model provides a direct link between fracture mechanics and damage mechanics which enables further development of the model for fatigue analysis. Validation of this model is presented by conducting three-dimensional mode I, mode II and mixed-mode simulations and comparing the results with analytical solutions. Copyright © 2016 John Wiley & Sons, Ltd.

This paper investigates the question of the building of admissible stress field in a substructured context. More precisely we analyze the special role played by multiple points. This study leads to (1) an improved recovery of the stress field, (2) an opportunity to minimize the estimator in the case of heterogeneous structures (in the parallel and sequential case), (3) a procedure to build admissible fields for FETI-DP and BDDC methods leading to an error bound which separates the contributions of the solver and of the discretization. Copyright © 2016 John Wiley & Sons, Ltd.

We present an approach for controlling the undercut and the minimal overhang angle in density based topology optimization, which are useful for reducing support structures in additive manufacturing. We cast both the undercut control and the minimal overhang angle control that are inherently constraints on the boundary shape into a domain integral of Heaviside projected density gradient. Such a Heaviside projection based integral of density gradient leads to a single constraint for controlling the undercut or controlling the overhang angle in the optimization. It effectively corresponds to a constraint on the projected perimeter that has undercut or has slope smaller than the prescribed overhang angle. In order to prevent trivial solutions of intermediate density to satisfy the density gradient constraints, a constraint on density grayness is also incorporated into the formulations. Numerical results on MBB beams, cantilever beams, and 2D and 3D heat conduction problems demonstrate the proposed formulations are effective in controlling the undercut and the minimal overhang angle in the optimized designs. Copyright © 2016 John Wiley & Sons, Ltd.

Topology optimization formulations using multiple design variables per finite element have been proposed to improve the design resolution. This paper discusses the relation between the number of design variables per element and the order of the elements used for analysis. We derive that beyond a maximum number of design variables, certain sets of material distributions cannot be discriminated based on the corresponding analysis results. This makes the design description inefficient and the solution of the optimization problem non-unique. To prevent this, we establish bounds for the maximum number of design variables that can be used to describe the material distribution for any given finite element scheme without introducing non-uniqueness. Copyright © 2016 John Wiley & Sons, Ltd.

The representation of structural boundaries is of significant importance for the applications of DEM in industry. In recent decade, triangle meshes are extensively used for representing structural boundaries. Structural boundaries of industrial objects are composed by regular shapes and irregular shapes in many occasions. A method for representing structural boundaries with regular shapes and irregular shapes has been developed by combining mathematical equations and triangle meshes for computational efficiency. When structural boundaries are represented by mathematical equations and triangle meshes, gaps or protuberances may exist at the connection boundaries between regular shapes and irregular shapes. For the exactness of representation of structural boundaries, gaps or protuberances are identified, expressed and treated successively, so that the geometrical shapes composed by regular shapes described with mathematical equations and irregular shapes described with triangle meshes can replace the original structural shapes. Two series of numerical tests have been conducted to verify this method. The results showed that this method can effectively represent the complicated structural boundaries containing regular shapes and irregular shapes, and greatly reduce the computational cost in the contact detection between particles and structural boundaries.

Efficient algorithms are considered for the computation of a reduced-order model based on the proper orthogonal decomposition methodology for the solution of parameterized elliptic partial differential equations. The method relies on partitioning the parameter space into subdomains based on the properties of the solution space and then forming a reduced basis for each of the subdomains. This yields more efficient offline and online stages for the proper orthogonal decomposition method. We extend these ideas for inexpensive adjoint based a posteriori error estimation of both the expensive finite element method solutions and the reduced-order model solutions, for a single and multiple quantities of interest. Various numerical results indicate the efficacy of the approach. This article is protected by copyright. All rights reserved.

The paper presents the generalization of the modification of classical boundary integral equation (BIE) and obtaining parametric integral equation system (PIES) for 2D elastoplastic problems. The modification was made to obtain such equations for which numerical solving does not require application of finite or boundary elements. This was achieved through the use of curves and surfaces for modeling introduced at the stage of analytical modification of the classic BIE. For approximation of plastic strains the Lagrange polynomials with various number and arrangement of interpolation nodes were used. Reliability of the modification was verified on examples with analytical solutions. This article is protected by copyright. All rights reserved.

This paper aims at accounting for the uncertainties due to material structure and surface topology of micro-beams in a stochastic multiscale model. For micro-resonators made of anisotropic polycrystalline materials, micro-scale uncertainties are due to the grain size, grain orientation, and to the surface profile. First, micro-scale realizations of stochastic volume elements (SVEs) are obtained based on experimental measurements. To account for the surface roughness, the SVEs are defined as a volume element having the same thickness as the MEMS, with a view to the use of a plate model at the structural scale. The uncertainties are then propagated up to an intermediate scale, the meso-scale, through a second-order homogenization procedure. From the meso-scale plate resultant material property realizations, a spatially correlated random field of the in plane, out of plane, and cross resultant material tensors can be characterized. Owing to this characterized random field, realizations of MEMS-scale problems can be defined on a plate finite element model. Samples of the macro-scale quantity of interest can then be computed by relying on a Monte-Carlo simulation procedure. As a case study, the resonance frequency of MEMS micro-beams is investigated for different uncertainty cases, such as grain preferred orientations and surface roughness effects. Copyright © 2016 John Wiley & Sons, Ltd.

We propose a highly effective approach using a novel adaptive methodology to perform topology optimization with polygonal meshes, called polytree meshes. Polytree is a hierarchical data structure based on the principle of recursive spatial decomposition of each polygonal element with *n* nodes into (*n* + 1) arbitrary new polygonal elements; enabling more efficient utilization of unstructured meshes and arbitrary design domains in topology optimization. In order to treat hanging nodes after each optimization loop, we define the Wachspress coordinate on a reference element and then utilize an affine map to obtain shape functions and their gradients on arbitrary polygons with *n* vertices and *m* hanging nodes, called side-nodes, as the polygonal element with (*n* + *m*) vertices. The polytree meshes do not only improve the boundary description quality of the optimal result, but also reduce the computational cost of optimization process in comparison with the use of uniformly fine meshes. Several numerical examples are investigated to show the high effectiveness of the proposed method. This article is protected by copyright. All rights reserved.

A new algorithm called recursive absolute nodal coordinate formulation algorithm (REC-ANCF) is presented for dynamic analysis of multi-flexible-body system including nonlinear large deformation. This method utilizes the absolute nodal coordinate formulation (ANCF) to describe flexible bodies, and establishes a kinematic and dynamic recursive relationship for the whole system based on the articulated-body algorithm (ABA). In the ordinary differential equations (ODEs) obtained by the REC-ANCF, a simple form of the system generalized Jacobian matrix and generalized mass matrix is obtained. Thus, a recursive forward dynamic solution is proposed to solve the ODEs one element by one element through an appropriate matrix manipulation. Utilizing the parent array to describe the topological structure, the REC-ANCF is suitable for generalized tree multibody systems. Besides, the cutting joint method is used in simple closed-loop systems to make sure the O(n) algorithm complexity of the REC-ANCF. Compared with common ANCF algorithms, the REC-ANCF has several advantages: the optimal algorithm complexity (O(n)) under limited processors, simple derivational process, no location or speed constraint violation problem, higher algorithm accuracy. The validity and efficiency of this method are verified by several numerical tests. This article is protected by copyright. All rights reserved.

The quasicontinuum (QC) method is a concurrent scale-bridging technique that extends atomistic accuracy to significantly larger length scales by reducing the full atomic ensemble to a small set of representative atoms and using interpolation to recover the motion of all lattice sites where full atomistic resolution is not necessary. While traditional QC methods thereby create interfaces between fully-resolved and coarse-grained regions, the recently introduced fully-nonlocal QC framework does not fundamentally differentiate between atomistic and coarsened domains. Adding adaptive refinement enables us to tie atomistic resolution to evolving regions of interest such as moving defects. However, model adaptivity is challenging because large particle motion is described based on a reference mesh (even in the atomistic regions). Unlike in the context of, e.g., finite element meshes, adaptivity here requires that (i) all vertices lie on a discrete point set (the atomic lattice), (ii) model refinement is performed locally and provides sufficient mesh quality, and (iii) Verlet neighborhood updates in the atomistic domain are performed against a Lagrangian mesh. With the suite of adaptivity tools outlined here, the nonlocal QC method is shown to bridge across scales from atomistics to the continuum in a truly seamless fashion, as illustrated for nanoindentation and void growth. This article is protected by copyright. All rights reserved.

In nonlinear model order reduction, hyper reduction designates the process of approximating a projection-based reduced-order operator on a reduced mesh, using a numerical algorithm whose computational complexity scales with the small size of the Projection-based Reduced-Order Model (PROM). Usually, the reduced mesh is constructed by sampling the large-scale mesh associated with the high-dimensional model underlying the PROM. The sampling process itself is governed by the minimization of the size of the reduced mesh for which the hyper reduction method of interest delivers the desired accuracy for a chosen set of training reduced-order quantities. Because such a construction procedure is combinatorially hard, its key objective function is conveniently substituted with a convex approximation. Nevertheless, for large-scale meshes, the resulting mesh sampling procedure remains computationally intensive. In this paper, three different convex approximations that promote sparsity in the solution are considered for constructing reduced meshes that are suitable for hyper reduction, and paired with appropriate active set algorithms for solving the resulting minimization problems. These algorithms are equipped with carefully designed parallel computational kernels in order to accelerate the overall process of mesh sampling for hyper reduction, and therefore achieve practicality for realistic, large-scale, nonlinear structural dynamics problems. Conclusions are also offered as to what algorithm is most suitable for constructing a reduced mesh for the purpose of hyper reduction. This article is protected by copyright. All rights reserved.

A large amount of research in computational mechanics has biased toward atomistic simulations. This trend, on one hand, is due to the increased demand to perform computations in nanoscale and, on the other hand, is due to the rather simple applications of pairwise potentials in modeling the interactions between atoms of a given crystal. The Cauchy–Born (CB) hypothesis has been used effectively to model the behavior of crystals under different loading conditions, in which the comparison with molecular dynamics simulations presents desirable coincidence between the results. A number of research works have been devoted to the validity of CB hypothesis and its application in post-elastic limit. However, the range of application of CB hypothesis is limited, and it remains questionable whether it is still applicable beyond the validity limit. In this paper, a multi-scale technique is developed for modeling of plastic deformations in nanoscale materials. The deformation gradient is decomposed into the plastic and elastic parts, i.e., **F** = **F**^{p}**F**^{e}. This decomposition is in contrast to the conventional decomposition, **F** = **F**^{e}**F**^{p}, generally encountered in continuum and crystal plasticity. It is shown that the former decomposition is more appropriate for the problem dealt within this work. Inspired by crystal plasticity, the plastic part is determined from the slip on potential slip systems. Based on the assumption that the CB hypothesis remains valid in the homogeneous deformation, the elastic deformation gradient resulting from the aforementioned decomposition is employed in conjunction with the CB hypothesis to update the state variables for face-centered cubic crystals. The assumption of homogeneity of elastic deformation gradient is justified by the fact that elastic deformations are considerably smaller than the plastic deformations. The computational algorithms are derived in details, and numerical simulations are presented through several examples to demonstrate the capability of the proposed computational algorithm in the modeling of golden crystals under different loading conditions. Copyright © 2016 John Wiley & Sons, Ltd.

A wide variety of jump discontinuities, such as shock fronts, are abound in high-speed flows. An accurate approximation of these fronts may require higher order techniques either under mesh-based methods or mesh-free methods. In the latter class, the smoothed particle hydrodynamics (SPH) is becoming popular as a promising method. However, the standard approach in SPH (like any other discrete methods) can result in highly diffusive solutions because of the inevitable use of artificial viscosity to suppress numerical oscillations. On the other hand, the SPH formulation allows innovative ways to model complicated phenomena. In this paper, we introduce the novel idea of a skewed Gaussian kernel, to improve the shock capturing capability in high speed flows. Here, the standard Gaussian kernel function is modified, and its ‘shape’ is altered with a predesigned tunable skewness parameter, while the basic unity property of the kernel function is still preserved. The SPH with the proposed skewed Gaussian kernel is then applied on a number of benchmark problems in computational fluid dynamics, featuring shocks. The simulations have shown significantly better shock capture through the skewed kernel approach as against the standard techniques, with almost no increase in computational time. Copyright © 2016 John Wiley & Sons, Ltd.

Finite deformation contact problems with frictional effects and finite shape changes due to wear are investigated. To capture the finite shape changes, a third configuration besides the well-known reference and spatial configurations is introduced, which represents the time-dependent worn state. Consistent interconnections between these states are realized by an arbitrary Lagrangean–Eulerian formulation. The newly developed partitioned and fully implicit algorithm is based on a Lagrangean step and a shape evolution step. Within the Lagrangean step, contact constraints as well as the wear equations are weakly enforced following the well-established mortar framework. Additional unknowns due to the employed Lagrange multiplier method for contact constraint enforcement and due to wear itself are eliminated by condensation procedures based on the concept of biorthogonality and the so-called dual shape functions. Several numerical examples in both 2D and 3D are provided to demonstrate the performance and accuracy of the proposed numerical algorithm. Copyright © 2016 John Wiley & Sons, Ltd.

Dynamic fragmentation is a rapid and catastrophic failure of a material. During this process, the nucleation, propagation, branching, and coalescence of cracks results in the formation of fragments. The numerical modeling of this phenomenon is challenging because it requires high-performance computational capabilities. For this purpose, the finite-element method with dynamic insertion of cohesive elements was chosen. This paper describes the parallel implementation of its fundamental algorithms in the C++ open-source library Akantu. Moreover, a numerical instability that can cause the loss of energy conservation and possible solutions to it are illustrated. Finally, the method is applied to the dynamic fragmentation of a hollow sphere subjected to uniform radial expansion. Copyright © 2016 John Wiley & Sons, Ltd.

An adaptive mesh refinement (AMR) technique is proposed for level set simulations of incompressible multiphase flows. The present AMR technique is implemented for two-dimensional/three-dimensional unstructured meshes and extended to multi-level refinement. Smooth variation of the element size is guaranteed near the interface region with the use of multi-level refinement. A Courant–Friedrich–Lewy condition for zone adaption frequency is newly introduced to obtain a mass-conservative solution of incompressible multiphase flows. Finite elements around the interface are dynamically refined using the classical element subdivision method. Accordingly, finite element method is employed to solve the problems governed by the incompressible Navier–Stokes equations, using the level set method for dynamically updated meshes. The accuracy of the adaptive solutions is found to be comparable with that of non-adaptive solutions only if a similar mesh resolution near the interface is provided. Because of the substantial reduction in the total number of nodes, the adaptive simulations with two-level refinement used to solve the incompressible Navier–Stokes equations with a free surface are about four times faster than the non-adaptive ones. Further, the overhead of the present AMR procedure is found to be very small, as compared with the total CPU time for an adaptive simulation. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, a structure-preserving direct method for the optimal control of mechanical systems is developed. The new method accommodates a large class of one-step integrators for the underlying state equations. The state equations under consideration govern the motion of affine Hamiltonian control systems. If the optimal control problem has symmetry, associated generalized momentum maps are conserved along an optimal path. This is in accordance with an extension of Noether's theorem to the realm of optimal control problems. In the present work, we focus on optimal control problems with rotational symmetries. The newly proposed direct approach is capable of exactly conserving generalized momentum maps associated with rotational symmetries of the optimal control problem. This is true for a variety of one-step integrators used for the discretization of the state equations. Examples are the one-step theta method, a partitioned variant of the theta method, and energy-momentum (EM) consistent integrators. Numerical investigations confirm the theoretical findings. Copyright © 2016 John Wiley & Sons, Ltd.

This paper proposes a most probable point (MPP)-based dimension reduction method (DRM) using the Hessian matrix called HeDRM to improve accuracy of reliability analysis in existing MPP-based DRM methods. Conventional MPP-based DRMs contain two types of errors: (1) error due to eliminating cross-terms of a performance function by using the univariate DRM; (2) error because of dependency of an axis direction after a rotational transformation. The proposed method minimizes the aforementioned errors by utilizing the Hessian matrix of a performance function. By performing an orthogonal transformation using the eigenvectors of the Hessian matrix, the cross-term effect of the performance function is minimized and the axis direction that results in the most accurate calculation is obtained because the Gaussian quadrature points for numerical integration are arranged along the eigenvector directions. In this way, the error incurred by exiting MPP-based DRMs can be reduced that leads to more accurate probability of failure estimation. In addition, this paper proposes to allocate the Gaussian quadrature points using the magnitude of the eigenvalues of the Hessian matrix. This allocation makes it possible to predetermine the number of function evaluations required to estimate the probability of failure accurately and efficiently. Copyright © 2016 John Wiley & Sons, Ltd.

Numerical solution of nonlinear eigenvalue problems (NEPs) is frequently encountered in computational science and engineering. The applicability of most existing methods is limited by the matrix structures, properties of the eigen-solutions, sizes of the problems, etc. This paper aims to remove those limitations and develop robust and universal NEP solvers for large-scale engineering applications. The novelty lies in two aspects. First, a rational interpolation approach (RIA) is proposed based on the Keldysh theorem for holomorphic matrix functions. Comparing with the existing contour integral approach, the RIA provides the possibility to select sampling points in more general regions and has advantages in improving the accuracy and reducing the computational cost. Second, a resolvent sampling scheme using the RIA is proposed to construct reliable search spaces for the Rayleigh–Ritz procedure, based on which a robust eigen-solver, called resolvent sampling based Rayleigh–Ritz method (RSRR), is developed for solving general NEPs. The RSRR can be easily implemented and parallelized. The advantages of the RIA and the performance of the RSRR are demonstrated by a variety of benchmark and application examples. Copyright © 2016 John Wiley & Sons, Ltd.

We present a higher order discretization scheme for the compressible Euler and Navier–Stokes equations with immersed boundaries. Our approach makes use of a discontinuous Galerkin discretization in a domain that is implicitly defined by means of a level set function. The zero iso-contour of this level set function is considered as an additional domain boundary where we weakly enforce boundary conditions in the same manner as in boundary-fitted cells. In order to retain the full order of convergence of the scheme, it is crucial to perform volume and surface integrals in boundary cells with high accuracy. This is achieved using a linear moment-fitting strategy. Moreover, we apply a non-intrusive cell-agglomeration technique that averts problems with very small and ill-shaped cuts. The robustness, accuracy, and convergence properties of the scheme are assessed in several two-dimensional test cases for the steady compressible Euler and Navier–Stokes equations. Approximation orders range from 0 to 4, even though the approach directly generalizes to even higher orders. In all test cases with a sufficiently smooth solution, the experimental order of convergence matches the expected rate for discontinuous Galerkin schemes. Copyright © 2016 John Wiley & Sons, Ltd.

We present an efficient numerical method to solve for cyclic steady states of nonlinear electro-mechanical devices excited at resonance. Many electro-mechanical systems are designed to operate at resonance, where the ramp-up simulation to steady state is computationally very expensive – especially when low damping is present. The proposed method relies on a Newton–Krylov shooting scheme for the direct calculation of the cyclic steady state, as opposed to a naïve transient time-stepping from zero initial conditions. We use a recently developed high-order Eulerian–Lagrangian finite element method in combination with an energy-preserving dynamic contact algorithm in order to solve the coupled electro-mechanical boundary value problem. The nonlinear coupled equations are evolved by means of an operator split of the mechanical and electrical problem with an explicit as well as implicit approach. The presented benchmark examples include the first three fundamental modes of a vibrating nanotube, as well as a micro-electro-mechanical disk resonator in dynamic steady contact. For the examples discussed, we observe power law computational speed-ups of the form *S* = 0.6·*ξ*^{ − 0.8}, where *ξ* is the linear damping ratio of the corresponding resonance frequency. Copyright © 2016 John Wiley & Sons, Ltd.

The peridynamic theory reformulates the equations of continuum mechanics in terms of integro-differential equations instead of partial differential equations. It is not straightforward to apply the available artificial boundary conditions for continua to peridynamic modeling. We therefore develop peridynamic transmitting boundary conditions (PTBCs) for 1D wave propagation. Differently from the previous method where the matching boundary condition is constructed for only one boundary material point, the PTBCs are established by considering the interaction and exchange of information between a group of boundary material points and another group of inner material points. The motion of the boundary material points is recursively constructed in terms of their locations and is determined through matching the peridynamic dispersion relation. The effectiveness of the PTBCs is examined by reflection analyses, numerical tests, and numerical convergent conditions. Furthermore, two-way interfacial conditions are proposed. The PTBCs are then applied to simulations of wave propagation in a bar with a defect, a composite bar with interfaces, and a domain with a seismic source. All the analyses and applications demonstrate that the PTBCs can effectively remove undesired numerical reflections at artificial boundaries. The methodology may be applied to modeling of wave propagation by other nonlocal theories. Copyright © 2016 John Wiley & Sons, Ltd.

A fast multipole boundary element method (FMBEM) extended by an adaptive mesh refinement algorithm for solving acoustic problems in three-dimensional space is presented in this paper. The Collocation method is used, and the Burton–Miller formulation is employed to overcome the fictitious eigenfrequencies arising for exterior domain problems. Because of the application of the combined integral equation, the developed FMBEM is feasible for all positive wave numbers even up to high frequencies. In order to evaluate the hypersingular integral resulting from the Burton–Miller formulation of the boundary integral equation, an integration technique for arbitrary element order is applied. The fast multipole method combined with an arbitrary order h-p mesh refinement strategy enables accurate computation of large-scale systems. Numerical examples substantiate the high accuracy attainable by the developed FMBEM, while requiring only moderate computational effort at the same time. Copyright © 2016 John Wiley & Sons, Ltd.

We propose a method to couple smoothed particle hydrodynamics and finite elements methods for nonlinear transient fluid–structure interaction simulations by adopting different time-steps depending on the fluid or solid sub-domains. These developments were motivated by the need to simulate highly non-linear and sudden phenomena requiring the use of explicit time integrators on both sub-domains (explicit Newmark for the solid and Runge–Kutta 2 for the fluid). However, due to critical time-step required for the stability of the explicit time integrators in, it becomes important to be able to integrate each sub-domain with a different time-step while respecting the features that a previously developed mono time-step coupling algorithm offered. For this matter, a dual-Schur decomposition method originally proposed for structural dynamics was considered, allowing to couple time integrators of the Newmark family with different time-steps with the use of Lagrange multipliers. Copyright © 2016 John Wiley & Sons, Ltd.

In the present paper, we address the delicate balance between computational efficiency and level of detailing at the modelling of ductile fracture in thin-walled structures. To represent the fine-scale nature of the ductile process, we propose a new extended finite element method-based enrichment of the displacement field to allow for crack tips that end or kink within an element. The idea is to refine the crack tip element locally in a way such that the macroscale node connectivity is unaltered. This allows for a better representation of the discontinuous kinematics without affecting the macroscale solution procedure, which would be a direct consequence of a regular mesh refinement. The method is first presented in a general 3D setting, and thereafter, it is specialised to shell theory for the modelling of crack propagation in thin-walled structures. The paper is concluded by a number of representative examples showing the accuracy of the method. We conclude that the ideas proposed in the paper enhance the current methodology for the analysis of ductile fracture of thin-walled large-scale structures under high strain rates. Copyright © 2016 John Wiley & Sons, Ltd.

We present two accurate and efficient numerical schemes for a phase field dendritic crystal growth model, which is derived from the variation of a free-energy functional, consisting of a temperature dependent bulk potential and a conformational entropy with a gradient-dependent anisotropic coefficient. We introduce a novel Invariant Energy Quadratization approach to transform the free-energy functional into a quadratic form by introducing new variables to substitute the nonlinear transformations. Based on the reformulated equivalent governing system, we develop a first and a second order semi-discretized scheme in time for the system, in which all nonlinear terms are treated semi-explicitly. The resulting semi-discretized equations consist of a linear elliptic equation system at each time step, where the coefficient matrix operator is positive definite and thus, the semi-discrete system can be solved efficiently. We further prove that the proposed schemes are unconditionally energy stable. Convergence test together with 2D and 3D numerical simulations for dendritic crystal growth are presented after the semi-discrete schemes are fully discretized in space using the finite difference method to demonstrate the stability and the accuracy of the proposed schemes. Copyright © 2016 John Wiley & Sons, Ltd.

A new four-node quadrilateral membrane finite element with drilling rotational degree of freedom based on the enhanced assumed strain formulation is presented. A simple formulation is achieved by five incompatible modes that are added to the Allman-type interpolation. Furthermore, modified shape functions are used to improve the behaviour of distorted elements. Numerical results show that the proposed new element exhibits good numerical accuracy and improved performance, and in many cases, superior to existing elements. In particular, Poisson's locking in nearly incompressible elasticity fades and the element performs well when it becomes considerably distorted even when it takes almost triangular shape. Copyright © 2016 John Wiley & Sons, Ltd.

The automatic generation of meshes for the finite element (FE) method can be an expensive computational burden, especially in structural problems with localized stress peaks. The use of meshless methods can address such an issue, as these techniques do not require the existence of an underlying connection among the particles selected in a general domain. This study advances a numerical strategy that blends the FE method with the meshless local Petrov–Galerkin technique in structural mechanics, with the aim at exploiting the most attractive features of each procedure. The idea relies on the use of FEs to compute a background solution that is locally improved by enriching the approximation space with the basis functions associated to a few meshless points, thus taking advantage of the flexibility ensured by the use of particles disconnected from an underlying grid. Adding the meshless particles only where needed avoids the cost of mesh refining, or even of remeshing, without the prohibitive computational cost of a thoroughly meshfree approach. In the present implementation, an efficient integration strategy for the computation of the coefficients taking into account the mutual FE–meshless local Petrov–Galerkin interactions is introduced. Moreover, essential boundary conditions are enforced separately on both FEs and meshless particles, thus allowing for an overall accuracy improvement also when the enriched region is close to the domain boundary. Numerical examples in structural problems show that the proposed approach can significantly improve the solution accuracy at a local level, with no remeshing effort, and at a low computational cost. Copyright © 2016 John Wiley & Sons, Ltd.

The possibility of using free-slip conditions within the context of the particle finite element method (PFEM) is investigated. For high Reynolds number engineering applications in which tangential effects at the fluid–solid boundaries are not of primary interest, the use of free-slip conditions can alleviate the need for very fine boundary layer meshes. Two novel ways for the imposition of free-slip conditions in the framework of the PFEM are presented. The proposed approach emphasizes robustness and simplicity, while retaining a sufficient level of generality. These two methods are then tested in the case of dam break and sloshing problems, and their respective advantages and drawbacks are discussed. It is also shown how the use of free-slip conditions can indirectly improve mass conservation properties of the PFEM, even when coarse meshes are employed. Copyright © 2016 John Wiley & Sons, Ltd.

We propose the study of a posteriori error estimates for time-dependent generalized finite element simulations of heat transfer problems. A residual estimate is shown to provide reliable and practically useful upper bounds for the numerical errors, independent of the heuristically chosen enrichment functions. Two sets of numerical experiments are presented. First, the error estimate is shown to capture the decrease in the error as the number of enrichment functions is increased or the time discretization refined. Second, the estimate is used to predict the behaviour of the error where no exact solution is available. It also reflects the errors incurred in the poorly conditioned systems typically encountered in generalized finite element methods. Finally, we study local error indicators in individual time steps and elements of the mesh. This creates a basis towards the adaptive selection and refinement of the enrichment functions. Copyright © 2016 John Wiley & Sons, Ltd.

The need for remeshing when computing flow problems in domains suffering large deformations has motivated the implementation of a tool that allows the proper transmission of information between finite element meshes. Because the Lagrangian projection of results from one mesh to another is a dissipative method, a new conservative interpolation method has been developed. A series of constraints, such as the conservation of mass or energy, are applied to the interpolated arrays through Lagrange multipliers in an error minimization problem, so that the resulting array satisfies these physical properties while staying as close as possible to the original interpolated values in the *L*^{2} norm. Unlike other conservative interpolation methods that require a considerable effort in mesh generation and modification, the proposed formulation is mesh independent and is only based on the physical properties of the field being interpolated. Moreover, the performed corrections are neither coupled with the main calculation nor with the interpolation itself, for which reason the computational cost is very low. Copyright © 2016 John Wiley & Sons, Ltd.

We apply a combination of the transient scaled boundary finite element method (SBFEM) and quadtree-based discretization to model dynamic problems at high frequencies. We demonstrate that the current formulation of the SBFEM for dynamics tends to require more degrees of freedom than a corresponding spectral element discretization when dealing with smooth problems on regular domains. Thus, we improve the efficiency of the SBFEM by proposing a novel approach to reduce the number of auxiliary variables for transient analyses. Based on this improved SBFEM, we present a modified meshing procedure, which creates a quadtree mesh purely based on the geometry and allows arbitrary sizes and orders of elements, as well as an arbitrary number of different materials. The discretization of each subdomain is created automatically based on material parameters and the highest frequency of interest. The transition between regions of different properties is straightforward when using the SBFEM. The proposed approach is applied to image-based analysis with a particular focus on geological models. Copyright © 2016 John Wiley & Sons, Ltd.

The node-based or edge-based smoothed finite element method is extended to develop polyhedral elements that are allowed to have an arbitrary number of nodes or faces, and so retain a good geometric adaptability. The strain smoothing technique and implicit shape functions based on the linear point interpolation make the element formulation simple and straightforward. The resulting polyhedral elements are free from the excessive zero-energy modes and yield a robust solution very much insensitive to mesh distortion. Several numerical examples within the framework of linear elasticity demonstrate the accuracy and convergence behavior. The smoothed finite element method-based polyhedral elements in general yield solutions of better accuracy and faster convergence rate than those of the conventional finite element methods. Copyright © 2016 John Wiley & Sons, Ltd.

We present two multiscale approaches for fracture analysis of full scale femur. The two approaches are the reduced order homogenization (ROH) previously developed by the first author and his associates and a novel accelerated ROH (AROH). The AROH is based on utilizing ROH calibrated to limited experimental data as a training tool to calibrate a simpler, single-scale anisotropic damage model. For bone tissue orientation, we take advantage of so-called Wolff's law, which states that bone tissue orientation is well correlated with principal strain direction in a stance position. The meso-phase properties are identified by minimizing error between the overall cortical and trabecular bone properties resulting from the quantitative computed tomography (QCT) scans and those predicted by the two-scale homogenization. The overall elastic and inelastic properties of the cortical and trabecular bone microstructure are derived from bone density that can be estimated from the Hounsfield units, which represent the measured grey levels in the QCT scans. For model validation, we conduct ROH and AROH simulations of full scale finite element model of femur created from the QCT and compare the simulation results with available experimental data. Copyright © 2016 John Wiley & Sons, Ltd.

This paper explores advantages offered by the stochastic collocation method based on the Smolyak grids for the solution of differential equations with random inputs in the parameter space. We use sparse Smolyak grids and the Chebyshev polynomials to construct multidimensional basis and approximate decoupled stochastic differential equations via interpolation. Disjoint set of grid points and basis functions allow us to gain significant improvement to conventional Smolyak algorithm. Density function and statistical moments of the solution are obtained by means of quadrature rules if inputs are uncorrelated and uniformly distributed. Otherwise, the Monte Carlo analysis can run inexpensively using obtained sparse approximation. An adaptive technique to sample from a multivariate density function using sparse grid is proposed to reduce the number of required sampling points. Global sensitivity analysis is viewed as an extension of the sparse interpolant construction and is performed by means of the Sobol' variance-based or the Kullback–Leibler entropy methods identifying the degree of contribution from the individual inputs as well as the cross terms. Copyright © 2016 John Wiley & Sons, Ltd.

Uncertain static plane stress analysis of continuous structure involving interval fields is investigated in this study. Unlike traditional interval analysis of discrete structure, the interval field is adopted to model the uncertainty, as well as the dependency between the physical locations and degrees of variability, of all interval system parameters presented in the continuous structures. By implementing the flexibility properties of some common structural elements, a new computational scheme is proposed to reformulate the uncertain static plane stress analysis with interval fields into standard mathematical programming problems. Consequently, feasible upper and lower bounds of structural responses can be effectively yet efficiently determined. In addition, the proposed method is adequate to deal with situations involving one-dimensional and two-dimensional interval fields, which enhances the pertinence of the proposed approach by incorporating both discrete and continuous structures. In addition, the proposed computational scheme is able to establish the realizations of the uncertain parameters causing the extreme structural responses at zero computational cost. The applicability and credibility of the established computational framework are rigorously justified by various numerical investigations. Copyright © 2016 John Wiley & Sons, Ltd.

A general method is proposed to couple two subregions analyzed with finite element digital image correlation even when using a mechanical regularization (regularized digital image correlation). A Lagrange multiplier is introduced to stitch both displacements fields in order to recover continuity over the full region of interest. Another interface unknown is introduced to ensure, additionally, the equilibrium of the mechanical models used for regularization. As a first application, the method is used to perform a single measurement from images at two different resolutions. Secondly, the method is also extended to parallel computing in regularized digital image correlation. The problem is formulated at the interface and solved with a Krylov-type algorithm. A dedicated preconditioner is proposed to significantly accelerate convergence. The resulting method is a good candidate for the analysis of large data sets. Copyright © 2016 John Wiley & Sons, Ltd.

We introduce a method to mesh the boundary Γ of a smooth, open domain in
immersed in a mesh of tetrahedra. The mesh follows by mapping a specific collection of triangular faces in the mesh to Γ. Two types of surface meshes follow: (a) a mesh that *exactly* meshes Γ, and (b) meshes that approximate Γ to any order, by interpolating the map over the selected faces; that is, an isoparametric approximation to Γ. The map we use to deform the faces is the closest point projection to Γ. We formulate conditions for the closest point projection to define a homeomorphism between each face and its image. These are conditions on some of the tetrahedra intersected by the boundary, and they essentially state that each such tetrahedron should (a) have a small enough diameter, and (b) have two of its dihedral angles be acute. We provide explicit upper bounds on the mesh size, and these can be computed on the fly.

We showcase the quality of the resulting meshes with several numerical examples. More importantly, *all* surfaces in these examples were meshed with a *single* background mesh. This is an important feature for problems in which the geometry evolves or changes, because it could be possible for the background mesh to never change as the geometry does. In this case, the background mesh would be a universal mesh for all these geometries. We expect the method introduced here to be the basis for the construction of universal meshes for domains in three dimensions. Copyright © 2016 John Wiley & Sons, Ltd.

For thin elastic structures submerged in heavy fluid, e.g., water, a strong interaction between the structural domain and the fluid domain occurs and significantly alters the eigenfrequencies. Therefore, the eigenanalysis of the fluid–structure interaction system is necessary. In this paper, a coupled finite element and boundary element (FE–BE) method is developed for the numerical eigenanalysis of the fluid–structure interaction problems. The structure is modeled by the finite element method. The compressibility of the fluid is taken into consideration, and hence the Helmholtz equation is employed as the governing equation and solved by the boundary element method (BEM). The resulting nonlinear eigenvalue problem is converted into a small linear one by applying a contour integral method. Adequate modifications are suggested to improve the efficiency of the contour integral method and avoid missing the eigenvalues of interest. The Burton–Miller formulation is applied to tackle the fictitious eigenfrequency problem of the BEM, and the optimal choice of its coupling parameter is investigated for the coupled FE–BE method. Numerical examples are given and discussed to demonstrate the effectiveness and accuracy of the developed FE–BE method. Copyright © 2016 John Wiley & Sons, Ltd.

A computational method is developed for evaluating the plastic strain gradient hardening term within a crystal plasticity formulation. While such gradient terms reproduce the size effects exhibited in experiments, incorporating derivatives of the plastic strain yields a nonlocal constitutive model. Rather than applying mixed methods, we propose an alternative method whereby the plastic deformation gradient is variationally projected from the elemental integration points onto a smoothed nodal field. Crucially, the projection utilizes the mapping between Lie groups and algebras in order to preserve essential physical properties, such as orthogonality of the plastic rotation tensor. Following the projection, the plastic strain field is directly differentiated to yield the Nye tensor. Additionally, an augmentation scheme is introduced within the global Newton iteration loop such that the computed Nye tensor field is fed back into the stress update procedure. Effectively, this method results in a fully implicit evolution of the constitutive model within a traditional displacement-based formulation. An elemental projection method with explicit time integration of the plastic rotation tensor is compared as a reference. A series of numerical tests are performed for several element types in order to assess the robustness of the method, with emphasis placed upon polycrystalline domains and multi-axis loading. Copyright © 2016 John Wiley & Sons, Ltd.

Atomistic models, which are crucial for performing molecular dynamics simulations of carbon nanostructures, consist of virtual hexagonal meshes with defects properly distributed in the intersectional areas. Currently, atomistic models are created mostly by hand, which is a notably tedious and time-consuming process. In this paper, we develop a method that produces atomistic models automatically. Because a hexagonal mesh and triangulation represent dual graphs, our work focuses on the creation of proper triangulation. The edge lengths of the triangulation should be compatible with the lengths of the C–C bonds, and vertices with valences other than 6 (due to the defects in the hexagonal mesh) should be properly arranged around the boundaries of the different components of a carbon nanostructure. Two techniques play important roles in our method: (1) sphere packing is used to place the nodes for triangulation that produces nearly constant edge lengths of the triangles and (2) the movement and editing of defects is used to control the number and positions of the defects. We subsequently develop a computer program based on this method that can create models much easier and faster than the current handwork method, thereby reducing the operation time significantly. Copyright © 2016 John Wiley & Sons, Ltd.

This work focuses on providing accurate low-cost approximations of stochastic finite elements simulations in the framework of linear elasticity. In a previous work, an adaptive strategy was introduced as an improved Monte-Carlo method for multi-dimensional large stochastic problems. We provide here a complete analysis of the method including a new enhanced goal-oriented error estimator and estimates of CPU (computational processing unit) cost gain. Technical insights of these two topics are presented in details, and numerical examples show the interest of these new developments. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, an efficient local radial basis function collocation method (LRBFCM) is presented for computing the band structures of the two-dimensional (2D) solid/fluid and fluid/solid phononic crystals. Both systems of solid scatterers embedded in a fluid matrix (solid/fluid phononic crystals) and fluid scatterers embedded in a solid matrix (fluid/solid phononic crystals) are investigated. The solid–fluid interactions are taken into account by properly formulating and treating the continuity/equilibrium conditions on the solid–fluid interfaces, which require an accurate computation of the normal derivatives of the displacements and the pressure on the fluid–solid interfaces and the unit-cell boundaries. The developed LRBFCM for the mixed wave propagation problems in 2D solid/fluid and fluid/solid phononic crystals is validated by the corresponding results obtained by the finite element method (FEM). To the best knowledge of the authors, the present LRBFCM has yet not been applied to the band structure computations of 2D solid/fluid and fluid/solid phononic crystals. For different lattice forms, scatterers' shapes, acoustic impedance ratios, and material combinations (solid scatterers in fluid matrix or fluid scatterers in solid matrix), numerical results are presented and discussed to reveal the efficiency and the accuracy of the developed LRBFCM for calculating the band structures of 2D solid/fluid and fluid/solid phononic crystals. A comparison of the present numerical results with that of the FEM shows that the present LRBFCM is much more efficient than the FEM for the band structure computations of the considered 2D solid/fluid and fluid/solid phononic crystals. Copyright © 2016 John Wiley & Sons, Ltd.

This article is concerned with the identification of stiffness parameters for short fiber-reinforced plastic materials, produced by injection molding and modeled by an orientation averaged transversely isotropic material law. Experiments are carried out in the form of computer simulations for a representative volume element of the material. The goal is to robustly identify the parameters with the fewest number of representative volume element simulations and such that the macroscopic model response attains a minimal variance for arbitrary strains and fiber orientations. Our approach carries over to problems involving anisotropic viscosity, thermal expansion and conductivity, and structurally similar problems. Copyright © 2016 John Wiley & Sons, Ltd.

In isogeometric analysis, identical basis functions are used for geometrical representation and analysis. In this work, non-uniform rational basis splines basis functions are applied in an isoparametric approach. An isogeometric Reissner–Mindlin shell formulation for implicit dynamic calculations using the Galerkin method is presented. A consistent as well as a lumped matrix formulation is implemented. The suitability of the developed shell formulation for natural frequency analysis is demonstrated by a numerical example. In a second set of examples, transient problems of plane and curved geometries undergoing large deformations in combination with nonlinear material behavior are investigated. Via a zero-thickness stress algorithm for arbitrary material models, a *J*_{2}-plasticity constitutive law is implemented. In the numerical examples, the effectiveness, robustness, and superior accuracy of a continuous interpolation method of the shell director vector is compared with experimental results and alternative numerical approaches. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, we model crack discontinuities in two-dimensional linear elastic continua using the extended finite element method without the need to partition an enriched element into a collection of triangles or quadrilaterals. For crack modeling in the extended finite element, the standard finite element approximation is enriched with a discontinuous function and the near-tip crack functions. Each element that is fully cut by the crack is decomposed into two simple (convex or nonconvex) polygons, whereas the element that contains the crack tip is treated as a nonconvex polygon. On using Euler's homogeneous function theorem and Stokes's theorem to numerically integrate homogeneous functions on convex and nonconvex polygons, the exact contributions to the stiffness matrix from discontinuous enriched basis functions are computed. For contributions to the stiffness matrix from weakly singular integrals (because of enrichment with asymptotic crack-tip functions), we only require a one-dimensional quadrature rule along the edges of a polygon. Hence, neither element-partitioning on either side of the crack discontinuity nor use of any cubature rule within an enriched element are needed. Structured finite element meshes consisting of rectangular elements, as well as unstructured triangular meshes, are used. We demonstrate the flexibility of the approach and its excellent accuracy in stress intensity factor computations for two-dimensional crack problems. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, a new discrete elements generation method based on geometry is proposed to fill geometric domains with particles (disks or spheres). By generating particles each one with a random radius or with a radius calculated from the iteration to ensure no overlaps exist between particles and identifying unstable particles and changing them to stable ones, a dense and stable packing can be created. A partitioning particle radius interval method and a particle stability inspection and improvement method are introduced to guarantee the algorithm's success and the stability of the particles. Some packings were created to evaluate the performances of the new method. The results showed that the algorithm was very efficient and was able to create isotropic packings of low porosities and large coordinate numbers. The partitioning particle radius interval method improved the generation efficiency significantly and increased the packing densities. Through the comparisons with several existing methods proposed recently, the method proposed in this work is found to be more efficient and can fill geometric domains with the lowest porosities. In addition, the stability of the particles is guaranteed and no complex triangular or tetrahedral mesh is required in particle generation, thereby making the new method simpler. Copyright © 2016 John Wiley & Sons, Ltd.

This paper proposes an original method to simulate the electrical conduction in continuums with the Discrete Element Method (DEM). The proposed method is based on the graphs theory applied to electrical resistance network, where the resistance between two discrete elements is estimated through the notion of ‘transmission surface’ to assume the discrete domain as a continuous medium. In addition to the electrical conduction, the Joule heating of a DEM domain has also been developed to take full advantage of the electrical conduction. The proposed method has been implemented in the free DEM software named ‘GranOO’.

The numerical results were compared against analytical approaches when applicable, or against Finite Element Method if the geometries become more complex or in case of dynamic loadings. The results are found satisfactory with errors around 3*%* for the electrical conduction and Joule heating of reasonably complex domains and loading cases. When it comes to more complex domains, such as electrical constriction, whilst the results remain close to those obtained with reference solutions (around 6*%*), they highlight the importance of taking care about the domains discretization.

Finally, the proposed method is applied to detect cracks onset on a cylindrical rod torsion test to show how to take advantage of the proposed work. Copyright © 2016 John Wiley & Sons, Ltd.

This article presents a detailed study on the potential and limitations of performing higher-order multi-resolution topology optimization with the finite cell method. To circumvent stiffness overestimation in high-contrast topologies, a length-scale is applied on the solution using filter methods. The relations between stiffness overestimation, the analysis system, and the applied length-scale are examined, while a high-resolution topology is maintained. The computational cost associated with nested topology optimization is reduced significantly compared with the use of first-order finite elements. This reduction is caused by exploiting the decoupling of density and analysis mesh, and by condensing the higher-order modes out of the stiffness matrix. Copyright © 2016 John Wiley & Sons, Ltd.

This paper presents an alternative topology optimization method based on an efficient meshless smoothed particle hydrodynamics (SPH) algorithm. To currently calculate the objective compliance, the deficiencies in standard SPH method are eliminated by introducing corrective smoothed particle method and total Lagrangian formulation. The compliance is established relative to a designed density variable at each SPH particle which is updated by optimality criteria method. Topology optimization is realized by minimizing the compliance using a modified solid isotropic material with penalization approach. Some numerical examples of plane elastic structure are carried out and the results demonstrate the suitability and effectiveness of the proposed SPH method in the topology optimization problem. Copyright © 2016 John Wiley & Sons, Ltd.

An accurate and easy integration technique is desired for the meshless methods of weak form. As is well known, a sub-domain method is often used in computational mechanics. The conforming sub-domains, where the sub-domains are not separated nor overlapped each other, are often used, while the nonconforming sub-domains could be employed if needed. In the latter cases, the integrations of the sub-domains may be performed easily by choosing a simple configuration. And then, the meshless method with nonconforming sub-domains is considered one of the reasonable choices for computational mechanics without the troublesome integration. In this paper, we propose a new sub-domain meshless method. It is noted that, because the method can employ both the conforming and the nonconforming sub-domains, the integration for the weak form is necessarily accurate and easy by selecting the nonconforming sub-domains with simple configuration. The boundary value problems including the Poisson's equation and the Helmholtz's equation are analyzed by using the proposed method. The numerical solutions are compared with the exact solutions and the solutions of the collocation method, showing that the relative errors by using the proposed method are smaller than those by using the collocation method and that the proposed method possesses a good convergence. Copyright © 2016 John Wiley & Sons, Ltd.

Large-scale systems of nonlinear equations appear in many applications. In various applications, the solution of the nonlinear equations should also be in a certain interval. A typical application is a discretized system of reaction diffusion equations. It is well known that chemical species should be positive otherwise the solution is not physical and in general blow up occurs. Recently, a projected Newton method has been developed, which can be used to solve this type of problems. A drawback is that the projected Newton method is not globally convergent. This motivates us to develop a new feasible projected Newton–Krylov algorithm for solving a constrained system of nonlinear equations. Combined with a projected gradient direction, our feasible projected Newton–Krylov algorithm circumvents the non-descent drawback of search directions which appear in the classical projected Newton methods. Global and local superlinear convergence of our approach is established under some standard assumptions. Numerical experiments are used to illustrate that the new projected Newton method is globally convergent and is a significate complementarity for Newton–Krylov algorithms known in the literature. © 2016 The Authors. *International Journal for Numerical Methods in Engineering* Published by John Wiley & Sons Ltd.

In this paper, we develop a strategy that enables a consistent wave transfer when coupling mechanical models exhibiting various description scales under dynamic loading. Incompatibilities between concurrent models are addressed by a selective filtering based on the perfectly matched layer framework and involving a dedicated projection operator. This filtering, performed in the vicinity of the coupling interface, aims at avoiding spurious wave reflections by automatically killing part of the signal with wavelength range not compatible with the downstream (coarser scale) model while the relevant complementary part is transferred. Performances of the proposed approach are evaluated in details by conducting one-dimensional and two-dimensional numerical experiments. The definition of the projection operator, which is a key point of the method, is particularly tackled by analyzing and comparing several possible choices. Copyright © 2016 John Wiley & Sons, Ltd.

In this work, a novel time marching procedure is presented in which adaptive time integration parameters are locally computed, allowing different spatial and temporal distributions. As a consequence, a more versatile technique arises, enabling an enhanced performance. The main features of the proposed technique are: (1) it is simple; (2) it has guaranteed stability; (3) it is an efficient non-iterative adaptive single-step procedure; (4) it provides enhanced accuracy; (5) it enables advanced controllable algorithmic dissipation in the higher modes; (6) it is truly self-starting; and (7) it is entirely automatic, requiring no input parameter from the user. The proposed technique is very complete, providing the main positive attributes that are requested from an effective time marching procedure. Numerical results are presented along the paper, illustrating the excellent performance of the method. Copyright © 2016 John Wiley & Sons, Ltd.

A numerical procedure for analysis of general laminated plates under transverse load is developed utilizing the Mindlin plate theory, the finite volume discretization, and a segregated solution algorithm. The force and moment balance equations with the laminate constitutive relations are written in the form of a generic transport equation. In order to obtain discrete counterparts of the governing equations, the plate is subdivided into N control volumes by a Cartesian numerical mesh. As a result, five sets of N linear equations with N unknowns are obtained and solved using the conjugate gradient method with preconditioning.

For the method validation, a number of test cases are designed to cover thick and thin laminated plates with aspect ratio (width to thickness) from 4 to 100. Simply supported orthotropic, symmetric cross-ply, and angle-ply laminated plates under uniform and sinusoidal pressure loads are solved, and results are compared with available analytical solutions. The shear correction factor of 5/6 is utilized throughout the procedure, which is consistent with test cases used in the reviewed literature. Comparisons of the finite volume method results for maximum deflections at the center of the plate and the Navier solutions obtained for aspect ratios 10, 20, and 100 shows a very good agreement. Copyright © 2016 John Wiley & Sons, Ltd.

We consider a technique to estimate an approximate gradient using an ensemble of randomly chosen control vectors, known as Ensemble Optimization (EnOpt) in the oil and gas reservoir simulation community. In particular, we address how to obtain accurate approximate gradients when the underlying numerical models contain uncertain parameters because of geological uncertainties. In that case, ‘robust optimization’ is performed by optimizing the expected value of the objective function over an ensemble of geological models. In earlier publications, based on the pioneering work of Chen *et al*. (2009), it has been suggested that a straightforward one-to-one combination of random control vectors and random geological models is capable of generating sufficiently accurate approximate gradients. However, this form of EnOpt does not always yield satisfactory results. In a recent article, Fonseca *et al*. (2015) formulate a modified EnOpt algorithm, referred to here as a Stochastic Simplex Approximate Gradient (StoSAG; in earlier publications referred to as ‘modified robust EnOpt’) and show, via computational experiments, that StoSAG generally yields significantly better gradient approximations than the standard EnOpt algorithm. Here, we provide theoretical arguments to show why StoSAG is superior to EnOpt. © 2016 The Authors. *International Journal for Numerical Methods in Engineering* Published by John Wiley & Sons, Ltd.

Represented by the element free Galerkin method, the meshless methods based on the Galerkin variational procedure have made great progress in both research and application. Nevertheless, their shape functions free of the Kronecker delta property present great troubles in enforcing the essential boundary condition and the material continuity condition. The procedures based on the relaxed variational formulations, such as the Lagrange multiplier-based methods and the penalty method, strongly depend on the problem in study, the interpolation scheme, or the artificial parameters. Some techniques for this issue developed for a particular method are hard to extend to other meshless methods. Under the framework of partition of unity and strict Galerkin variational formulation, this study, taking Poisson's boundary value problem for instance, proposes a unified way to treat exactly both the material interface and the nonhomogeneous essential boundary as in the finite element analysis, which is fit for any partition of unity-based meshless methods. The solution of several typical examples suggests that compared with the Lagrange multiplier method and the penalty method, the proposed method can be always used safely to yield satisfactory results. Copyright © 2016 John Wiley & Sons, Ltd.

This paper is concerned with the development of mesh-free models for the static analysis of smart laminated composite beams. The overall smart composite beam is composed of a laminated substrate composite beam and a piezoelectric layer attached partially or fully at the top surface of the substrate beam. The piezoelectric layer acts as the distributed actuator layer of the smart beam. A layer-wise displacement theory and an equivalent single-layer theory have been used to derive the models. Several cross-ply substrate beams are considered for presenting the numerical results. The responses of the smart composite beams computed by the present new mesh-free model based on the layer-wise displacement theory excellently match with those obtained by the exact solutions. The mesh-free model based on the equivalent single-layer theory cannot accurately compute the responses due to transverse actuation by the piezoelectric actuator. The models derived here suggest that the mesh-free method can be efficiently used for the numerical analysis of smart structures. Copyright © 2016 John Wiley & Sons, Ltd.

This work deals with the efficient time integration of mechanical systems with elastohydrodynamic (EHD) lubricated joints. Two novel approaches are presented. First, a projection function is used to formulate the well-known Swift–Stieber cavitation condition and the mass-conservative cavitation condition of Elrod as an unconstrained problem. Based on this formulation, the pressure variable from the EHD problem is added to the dynamic equations of a multi-body system in a monolithic manner so that cavitation is solved within a global iteration. Compared with a partitioned state-of-the-art formulation, where the pressure is solved locally in a nonlinear force element, this global search reduces simulation time.

Second, a Quasi-Newton method of DeGroote is applied during time integration to solve the nonlinear relation between pressure and deformation. Compared with a simplified Newton method, the calculation of the time-consuming parts of the Jacobian are avoided, and therefore, simulation time is reduced significantly, when the Jacobian is calculated numerically. Solution strategies with the Quasi-Newton method are presented for the partitioned formulation as well as for the new DAE formulations with projection function. Results are given for a simulation example of a rigid shaft in a flexible bearing. Copyright © 2016 John Wiley & Sons, Ltd.

In this work, we use Nitsche's formulation to weakly enforce kinematic constraints at an embedded interface in Helmholtz problems. Allowing embedded interfaces in a mesh provides significant ease for discretization, especially when material interfaces have complex geometries. We provide analytical results that establish the well-posedness of Helmholtz variational problems and convergence of the corresponding finite element discretizations when Nitsche's method is used to enforce kinematic constraints. As in the analysis of conventional Helmholtz problems, we show that the inf-sup constant remains positive provided that the Nitsche's stabilization parameter is judiciously chosen. We then apply our formulation to several 2D plane-wave examples that confirm our analytical findings. Doing so, we demonstrate the asymptotic convergence of the proposed method and show that numerical results are in accordance with the theoretical analysis. Copyright © 2016 John Wiley & Sons, Ltd.

We present a framework to efficiently solve large deformation contact problems with nearly incompressible materials by implementing adaptive remeshing. Specifically, nodally integrated elements are employed to avoid mesh locking when linear triangular or tetrahedral elements are used to facilitate mesh re-generation. Solution variables in the bulk and on contact surfaces are transferred between meshes such that accuracy is maintained and re-equilibration on the new mesh is possible. In particular, the displacement transfer in the bulk is accomplished through a constrained least squares problem based on nodal integration, while the transfer of contact tractions relies on parallel transport. Finally, a residual-based error indicator is chosen to drive adaptive mesh refinement. The proposed strategies are applicable to both two-dimensional or three-dimensional analysis and are demonstrated to be robust by a number of numerical examples. Copyright © 2016 John Wiley & Sons, Ltd.

The computation of elastic continua is widely used in today's engineering practice, and finite element models yield a reasonable approximation of the physically observed behaviour. In contrast, the failure of materials due to overloading can be seen as a sequence of discontinuous effects finally leading to a system failure. Until now, it has not been possible to sufficiently predict this process with numerical simulations. It has recently been shown that discrete models like lattice spring models are a promising alternative to finite element models for computing the breakdown of materials because of static overstress and fatigue. In this paper, we will address one of the downsides of current lattice spring models, the need for a periodic mesh leading to a mesh-induced anisotropy of material failure in simulations. We will show how to derive irregular cells that still behave as part of a homogeneous continuum irrespectively of their shape and which should be able to eliminate mesh-induced anisotropy. In addition, no restraints concerning the material stiffness tensor are introduced, which allows the simulation of non-isotropic materials. Finally, we compare the elastic response of the presented model with present lattice spring models using mechanical systems with a known analytical stress field. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, we propose a fully smoothed extended finite element method for axisymmetric problems with weak discontinuities. The salient feature of the proposed approach is that all the terms in the stiffness and mass matrices can be computed by smoothing technique. This is accomplished by combining the Gaussian divergence theorem with the evaluation of indefinite integral based on smoothing technique, which is used to transform the domain integral into boundary integral. The proposed technique completely eliminates the need for isoparametric mapping and the computing of Jacobian matrix even for the mass matrix. When employed over the enriched elements, the proposed technique does not require sub-triangulation for the purpose of numerical integration. The accuracy and convergence properties of the proposed technique are demonstrated with a few problems in elastostatics and elastodynamics with weak discontinuities. It can be seen that the proposed technique yields stable and accurate solutions and is less sensitive to mesh distortion. Copyright © 2016 John Wiley & Sons, Ltd.

The treatments of heterogeneities and periodic boundary conditions are explored to properly perform isogeometric analysis (IGA) based on NURBS basis functions in solving homogenization problems for heterogeneous media with omni-directional periodicity and composite plates with in-plane periodicity. Because the treatment of the combination of different materials in IGA models is not trivial especially for periodicity constraints, the first priority is to clearly specify points at issue in the numerical modeling, or equivalently mesh generation, for IG homogenization analysis (IGHA). The most awkward, but important issue is how to generate patches for NURBS representation of the geometry of a rectangular parallelepiped unit cell to realize appropriate deformations in consideration of the convex-hull property of IGA. The issue arises from the introduction of overlapped control points located at angular points in the heterogeneous unit cell, which must satisfy multiple point constraint (MPC) conditions associated with periodic boundary conditions (PBCs). Although two measures may be conceivable, we suggest the use of multiple patches along with double MPC that imposes PBCs and the continuity conditions between different patches simultaneously. Several numerical examples of numerical material and plate tests are presented to demonstrate the validity of the proposed strategy of IG modeling for IGHA. Copyright © 2016 John Wiley & Sons, Ltd.

The discrete element method, developed by Cundall and Strack, typically uses some variations of the central difference numerical integration scheme. However, like all explicit schemes, the scheme is only conditionally stable, with the stability determined by the size of the time-step. The current methods for estimating appropriate discrete element method time-steps are based on many assumptions; therefore, large factors of safety are usually applied to the time-step to ensure stability, which substantially increases the computational cost of a simulation. This work introduces a general framework for estimating critical time-steps for any planar rigid body subject to linear damping and forcing. A numerical investigation of how system damping, coupled with non-collinear impact, affects the critical time-step is also presented. It is shown that the critical time-step is proportional to
if a linear contact model is adopted, where *m* and *k* represent mass and stiffness, respectively. The term which multiplies this factor is a function of known physical parameters of the system. The stability of a system is independent of the initial conditions. Copyright © 2016 John Wiley & Sons, Ltd.

We present a novel methodology to effectively localize radial basis function approximation methods in three dimensions. The local scheme requires shape parameter-dependent functions that can be used to approximate gradients of scattered data and to solve partial differential operators. The optimum shape parameter is obtained from the highest gradient of interest, where a known analytical function, when boundary conditions are not present, or a shape parameter-free global approximation are used to educate the localized scheme. The later option is applicable to problems where the operator needs to be solved multiple times, like in time evolution or stochastic integration. Past shape parameter's optimizations, for two-dimensional domains, based on the condition number of the interpolant matrix, were unable to provide satisfactory approximations. The applicability of our method is illustrated in the context of an analytical expression interpolation and during a Ginzburg–Landau relaxation of a free energy functional. In general, the optimum shape parameter depends on geometry, node distribution, and density, whereas the approximation errors decrease as the node density and the local stencil size increase. The effective localization of radial basis functions motivates its use in moving boundary problems and accelerates solutions through sparse matrix solvers. Copyright © 2016 John Wiley & Sons, Ltd.

The research work extends the *scaled boundary finite element method* to non-deterministic framework defined on random domain wherein random behaviour is exhibited in the presence of the random-field uncertainties. The aim is to blend the scaled boundary finite element method into the Galerkin spectral stochastic methods, which leads to a proficient procedure for handling the stress singularity problems and crack analysis. The Young's modulus of structures is considered to have random-field uncertainty resulting in the stochastic behaviour of responses. Mathematical expressions and the solution procedure are derived to evaluate the statistical characteristics of responses (displacement, strain, and stress) and stress intensity factors of cracked structures. The feasibility and effectiveness of the presented method are demonstrated by particularly chosen numerical examples. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, we present a solution framework for high-order discretizations of conjugate heat transfer problems on non-body-conforming meshes. The framework consists of and leverages recent developments in discontinuous Galerkin discretization, simplex cut-cell techniques, and anisotropic output-based adaptation. With the cut-cell technique, the mesh generation process is completely decoupled from the interface definitions. In addition, the adaptive scheme combined with the discontinuous Galerkin discretization automatically adjusts the mesh in each sub-domain and achieves high-order accuracy in outputs of interest. We demonstrate the solution framework through several multi-domained conjugate heat transfer problems consisting of laminar and turbulent flows, curved geometry, and highly coupled heat transfer regions. The combination of these attributes yield nonintuitive coupled interactions between fluid and solid domains, which can be difficult to capture with user-generated meshes. Copyright © 2016 John Wiley & Sons, Ltd.

This paper proposes a new methodology to guarantee the accuracy of the homogenisation schemes that are traditionally employed to approximate the solution of PDEs with random, fast evolving diffusion coefficients. More precisely, in the context of linear elliptic diffusion problems in randomly packed particulate composites, we develop an approach to strictly bound the error in the expectation and second moment of quantities of interest, without ever solving the fine-scale, intractable stochastic problem. The most attractive feature of our approach is that the error bounds are computed without any integration of the fine-scale features. Our computations are purely macroscopic, deterministic and remain tractable even for small scale ratios. The second contribution of the paper is an alternative derivation of modelling error bounds through the Prager–Synge hypercircle theorem. We show that this approach allows us to fully characterise and optimally tighten the interval in which predicted quantities of interest are guaranteed to lie. We interpret our optimum result as an extension of Reuss–Voigt approaches, which are classically used to estimate the homogenised diffusion coefficients of composites, to the estimation of macroscopic engineering quantities of interest. Finally, we make use of these derivations to obtain an efficient procedure for multiscale model verification and adaptation. Copyright © 2016 John Wiley & Sons, Ltd.

We study the simultaneous analysis and design (SAND) formulation of the ‘classical’ topology optimization problem subject to linear constraints on material density variables. Based on a dual method in theory, and a primal-dual method in practice, we propose a separable and strictly convex quadratic Lagrange–Newton subproblem for use in sequential approximate optimization of the SAND-formulated classical topology design problem.

The SAND problem is characterized by a large number of nonlinear equality constraints (the equations of equilibrium) that are linearized in the approximate convex subproblems. The availability of cheap second-order information is exploited in a Lagrange–Newton sequential quadratic programming-like framework. In the spirit of efficient structural optimization methods, the quadratic terms are restricted to the diagonal of the Hessian matrix; the subproblems have minimal storage requirements, are easy to solve, and positive definiteness of the diagonal Hessian matrix is trivially enforced.

Theoretical considerations reveal that the dual statement of the proposed subproblem for SAND minimum compliance design agrees with the ever-popular optimality criterion method – which is a nested analysis and design formulation. This relates, in turn, to the known equivalence between rudimentary dual sequential approximate optimization algorithms based on reciprocal (and exponential) intervening variables and the optimality criterion method. Copyright © 2016 John Wiley & Sons, Ltd.

Nonlocal integral and/or gradient enhancements are widely used to resolve the mesh dependency issue with standard continuum damage models. However, it is reported that whereas the structural response is mesh independent, a spurious damage growth is observed. Accordingly, a class of modified nonlocal enhancements is developed in literature, where the interaction domain increases with damage. In this contribution, we adopt a contrary view that the interaction domain *decreases* with damage. This is motivated by the fact that the fracture of quasi-brittle materials typically starts as a diffuse network of microcracks, before localizing into a macroscopic crack. To ensure thermodynamics consistency, the micromorphic theory is adopted in the model development. The ensuing microforce balance resembles closely the Helmholtz expression in a conventional gradient damage model. The superior performance of the localizing gradient damage model is demonstrated through a one-dimensional problem, as well as mode I and II failure in plane deformation. For all three cases, a localized deformation band at material failure is obtained. Copyright © 2016 John Wiley & Sons, Ltd.

No abstract is available for this article.

]]>No abstract is available for this article.

]]>Modal derivative is an approach to compute a reduced basis for model order reduction of large-scale nonlinear systems that typically stem from the discretization of partial differential equations. In this way, a complex nonlinear simulation model can be integrated into an optimization problem or the design of a controller, based on the resulting small-scale state-space model. We investigate the approximation properties of modal derivatives analytically and thus lay a theoretical foundation of their use in model order reduction, which has been missing so far. Concentrating on the application field of structural mechanics and structural dynamics, we show that the concept of modal derivatives can also be applied as nonlinear extension of the Craig–Bampton family of methods for substructuring. We furthermore generalize the approach from a pure projection scheme to a novel reduced-order modeling method that replaces all nonlinear terms by quadratic expressions in the reduced state variables. This complexity reduction leads to a frequency-preserving nonlinear quadratic state-space model. Numerical examples with carefully chosen nonlinear model problems and three-dimensional nonlinear elasticity confirm the analytical properties of the modal derivative reduction and show the potential of the proposed novel complexity reduction methods, along with the current limitations. Copyright © 2016 John Wiley & Sons, Ltd.

A computational framework for assisting in the development of novel textiles is presented. Electronic textiles are key in the rapidly growing field of wearable electronics for both consumer and military uses. There are two main challenges to the modeling of electronic textiles: the discretization of the textile microstructure and the interaction between electromagnetic and mechanical fields. A director-based beam formulation with an assumed electrical current is used to discretize the fabric at the level of individual fibrils. The open-source package FEniCS was used to implement the finite element model. Contact integrals were added into the FEniCS framework so that multiphysics contact laws can be incorporated in the same framework, leveraging the code generation and automated differentiation capabilities of FEniCS to produce the tangents needed by the implicit solution method. The computational model is used to construct and determine the mechanical, thermal, and electrical properties of a representative volume elements of a plain woven textile. Dynamic relaxation to solve the mechanical fields and the electrical and thermal fields is solved statically for a given mechanical state. The simulated electrical responses are fit to a simplified Kirchhoff network model to determine effective resistances of the textile. Copyright © 2016 John Wiley & Sons, Ltd.

Computational modelling of fracture has been attempted in the past with a range of numerical approaches including finite element, extended finite element and meshless methods. The cracking particle method (CPM) of Rabczuk is a pragmatic alternative to explicit modelling of crack surfaces in which a crack is represented by a set of cracking particles that can be easily updated when the crack propagates. The change of cracking angle is recorded in discrete segments of broken lines, which makes this methodology suitable to model discontinuous cracks. In this paper, a new CPM is presented that improves on two counts: firstly, crack path curvature modelling is improved by the use of bilinear segments centred at each particle and secondly, efficiency for larger problems is improved via an adaptive process of both refinement and recovery. The system stiffness is calculated and stored in local matrices, so only a small influenced domain should be recalculated for each step while the remainder can be read directly from storage, which greatly reduces the computational expense. The methodology is applied to several 2D crack problems, and good agreement to analytical solutions and previous work is obtained. Copyright © 2016 John Wiley & Sons, Ltd.

A computational framework for scale-bridging in multi-scale simulations is presented. The framework enables seamless combination of at-scale models into highly dynamic hierarchies to build a multi-scale model. Its centerpiece is formulated as a standalone module capable of fully asynchronous operation. We assess its feasibility and performance for a two-scale model applied to two challenging test problems from impact physics. We find that the computational cost associated with using the framework may, as expected, become substantial. However, the framework has the ability of effortlessly combining at-scale models to render complex multi-scale models. The main source of the computational inefficiency of the framework is related to poor load balancing of the lower-scale model evaluation We demonstrate that the load balancing can be efficiently addressed by recourse to conventional load-balancing strategies. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, we present a generalized prismatic hybrid meshing method for viscous flow simulations. One major difficulty in implementing a robust prismatic hybrid meshing tool is to handle boundary layer mesh collisions, and normally an extra data structure (e.g. quadtree in two-dimensional and octree in three-dimensional) is required. The proposed method overcomes this difficulty via an heuristic approach, and it only relies on constrained delaunay triangulation/tetrahedralization(CDT). No extra data structures are required. Geometrical reasoning is used to approximate the maximum marching distance of each point by walking through the CDT. This is combined with post-processing of marching vectors and distance and prohibition of multilevel differences to form an automatic and robust mechanism to remove boundary layer mesh collisions. Benefiting from the matureness of CDT techniques, the proposed method is robust, efficient and simple to implement. Its capability is demonstrated by generating quality prismatic hybrid meshes for industrial models with complex geometries. The proposed method is believed to be able considerably reduce the effort to implement a robust hybrid prismatic mesh generator for viscous flow simulations. © 2016 The Authors. International Journal for Numerical Methods in Engineering. Published by John Wiley & Sons Ltd.

This paper presents a 3D formulation for quasi-kinematic limit analysis, which is based on a radial point interpolation meshless method and numerical optimization. The velocity field is interpolated using radial point interpolation shape functions, and the resulting optimization problem is cast as a standard second-order cone programming problem. Because the essential boundary conditions can be only guaranteed at the position of the nodes when using radial point interpolation, the results obtained with the proposed approach are not rigorous upper bound solutions. This paper aims to improve the computing efficiency of 3D upper bound limit analysis and large problems, with tens of thousands of nodes, can be solved efficiently. Five numerical examples are given to confirm the effectiveness of the proposed approach with the von Mises yield criterion: an internally pressurized cylinder; a cantilever beam; a double-notched tensile specimen; and strip, square and rectangular footings. Copyright © 2016 John Wiley & Sons, Ltd.

In this article, a black-box higher-order fast multipole method for solving boundary integral equations on parametric surfaces in three spatial dimensions is proposed. Such piecewise smooth surfaces are the topic of recent studies in isogeometric analysis. Due to the exact surface representation, the rate of convergence of higher-order methods is not limited by approximation errors of the surface. An element-wise clustering strategy yields a balanced cluster tree and an efficient numerical integration scheme for the underlying Galerkin method. By performing the interpolation for the fast multipole method directly on the reference domain, the cost complexity in the polynomial degree is reduced by one order. This gain is independent of the application of either - or -matrices. In fact, several simplifications in the construction of -matrices are pointed out, which are a by-product of the surface representation. Extensive numerical examples are provided in order to quantify and qualify the proposed method. Copyright © 2016 John Wiley & Sons, Ltd.