The LATIN (acronym of LArge Time INcrement) method was originally devised as a non-incremental procedure for the solution of quasi-static problems in continuum mechanics with material nonlinearity. In contrast to standard incremental methods like Newton and modified Newton, LATIN is an iterative procedure applied to the entire loading path. In each LATIN iteration two problems are solved: a local problem, which is nonlinear but algebraic and miniature, and a global problem, which involves the entire loading process but is linear. The convergence of these iterations, which has been shown to occur for a large class of nonlinear problems, provides an approximate solution to the original problem. In this paper the LATIN method is presented from a different viewpoint, taking advantage of the causality principle. In this new view, LATIN is an *incremental method*, and the LATIN iterations are performed within each load step, similarly to the way that Newton iterations are performed. The advantages of the new approach are discussed. In addition, LATIN is extended for the solution of time-dependent wave problems. As a relatively simple model for illustrating the new formulation, lateral wave propagation in a flat membrane made of a nonlinear material is considered. Numerical examples demonstrate the performance of the scheme, in conjunction with finite element discretization in space and the Newmark trapezoidal algorithm in time. This article is protected by copyright. All rights reserved.

We introduce the use of hybridizable discontinuous Galerkin (HDG) finite element methods on overlapping (overset) meshes. Overset mesh methods are advantageous for solving problems on complex geometrical domains. We combine geometric flexibility of overset methods with the advantages of HDG methods: arbitrarily high-order accuracy, reduced size of the global discrete problem, and the ability to solve elliptic, parabolic, and/or hyperbolic problems with a unified form of discretization. Our approach to developing the ‘overset HDG’ method is to couple the global solution from one mesh to the local solution on the overset mesh. We present numerical examples for steady convection-diffusion and static elasticity problems. The examples demonstrate optimal order convergence in all primal fields for an arbitrary amount of overlap of the underlying meshes. This article is protected by copyright. All rights reserved.

Lack of conservation has been the biggest drawback in meshfree generalized finite difference methods (GFDMs). In this paper, we present a novel modification of classical meshfree GFDMs to include local balances which produce an approximate conservation of numerical fluxes. This numerical flux conservation is done within the usual moving least squares framework. Unlike Finite Volume Methods, it is based on locally defined control cells, rather than a globally defined mesh. We present the application of this method to an advection diffusion equation and the incompressible Navier–Stokes equations. Our simulations show that the introduction of flux conservation significantly reduces the errors in conservation in meshfree GFDMs. This article is protected by copyright. All rights reserved.

This paper presents a numerical investigation of microscopic lubricant flows from the cavities to the plateaus of the surface roughness of metal sheets during forming processes. This phenomenon, called micro-plasto-hydrodynamic (MPH) lubrication, was observed experimentally in various situations such as compression sliding tests, strip drawing and cold rolling. It leads to local friction drop and wear reduction. It is therefore critical to achieve a good understanding of this phenomenon.

To move towards that goal, a multiscale fluid-structure interaction (FSI) model is developed to model lubricant flows at the microscopic scale. These simulations are made possible through the use of the Arbitrary Lagrangian Eulerian (ALE) formalism.

In this paper, this methodology is used to study plane strip drawing. The numerical model is able to predict the onset of lubricant escape and the amount of lubricant flowing on the plateaus. Numerical results exhibit good agreement with experimental measurements.

This paper presents a novel formulation based on Hellinger-Reissner variational principle in the framework of quasi-conforming method for static and free vibration analysis of Reissner-Mindlin plates. The formulation starts from polynomial approximation of stresses, which satisfy the equilibrium equations of Reissner-Mindlin plate theory. Then the stress matrix is treated as the weighted function to weaken the strain-displacement equations after the strains are derived by using the constitutive equations. Finally the string-net functions are introduced to calculate strain integration. As examples, two new plate bending elements, a 4-node quadrilateral element QC-P4-11*β* and a 3-node triangular element QC-P3-7*β*, are proposed. Several benchmark examples are demonstrated to show the performance of the elements, and the results obtained are compared with other available ones. Numerical results have proved that both elements possess excellent precision. In particular, the quadrilateral element performs well even when the element shape degenerates into a triangle or concave quadrangle. This article is protected by copyright. All rights reserved.

A true integration-free meshless method based on Taylor series named Taylor Meshless Method (TMM) has been proposed to solve two dimensional Partial Differential Equations (PDE). In this framework, the shape functions are approximated solutions of the PDE and the discretization concerns only the boundary. In this paper, the applicability of TMM to solve large scale problems is discussed under two aspects. First, as in some other meshless methods, ill-conditioned matrices and round-off error propagation could lead to a loss of accuracy when the number of unknowns increases. This point will be investigated in the case of large scale problems. Second the computation time and its distribution are analyzed from numerical experiments for PDE's in a 3D domain. It is established that the TMM method is efficient and robust, even in the case of large scale problems while the finite element numerical model involves more than 3 millions degrees of freedom. This article is protected by copyright. All rights reserved.

We consider model reduction for magneto-quasistatic field equations in the vector potential formulation. A finite element discretization of such equations leads to large-scale differential-algebraic equations of special structure. For model reduction of linear systems, we employ a balanced truncation approach, whereas nonlinear systems are reduced using a proper orthogonal decomposition method combined with a discrete empirical interpolation technique. We will exploit the special block structure of the underlying problem to improve the performance of the model reduction algorithms. Furthermore, we discuss an efficient evaluation of the Jacobi matrix required in nonlinear time integration of the reduced models. This article is protected by copyright. All rights reserved.

Multi-scale problems are often solved by decomposing the problem domain into multiple subdomains, solving them independently using different levels of spatial and temporal refinement, and coupling the subdomain solutions back to obtain the global solution. Most commonly finite elements are used for spatial discretization and finite difference time-stepping is used for time integration. Given a finite element mesh for the global problem domain, the number of possible decompositions into subdomains and the possible choices for associated time-steps is exponentially large and the computational costs associated with different decompositions can vary by orders of magnitude. The problem of finding an optimal decomposition and the associated time discretization that minimizes computational costs while maintaining accuracy is non-trivial. Existing mesh partitioning tools, such as METIS, overlook the constraints posed by multi-scale methods, and lead to sub-optimal partitions with a high performance penalty. We present a multi-level mesh partitioning approach that exploits domain-specific knowledge of multi-scale methods to produce nearly optimal mesh partitions and associated time-steps automatically. Results show that for multi-scale problems, our approach produces decompositions that outperform those produced by state-of-the-art partitioners like METIS and even those that are manually constructed by domain experts. This article is protected by copyright. All rights reserved.

Because of the complexity of fluid flow solvers, non-intrusive uncertainty quantification techniques have been developed in aerodynamic simulations in order to compute the quantities of interest required in an optimization process, for example. The objective function is commonly expressed in terms of moments of these quantities, such as the mean, standard deviation, or even higher-order moments. Polynomial surrogate models based on polynomial chaos expansions have often been implemented in this respect. The original approach of uncertainty quantification using polynomial chaos is however intrusive. It is based on a Galerkin-type formulation of the model equations to derive the governing equations for the polynomial expansion coefficients. Third-order, indeed fourth-order moments of the polynomials are needed in this analysis. Besides, both intrusive and non-intrusive approaches call for their computation provided that higher-order moments of the quantities of interest need be post-processed. In most applications they are evaluated by Gauss quadratures, and eventually stored for use throughout the computations. In this paper analytical formulas are rather considered for the moments of the continuous polynomials of the Askey scheme, so that they can be evaluated by quadrature-free procedures instead. Matlab^{©} codes have been developed for this purpose and tested by comparisons with Gauss quadratures. This article is protected by copyright. All rights reserved.

The optimal design of mechanical structures subject to periodic excitations within a large frequency interval is quite challenging. In order to avoid bad performances for non-discretized frequencies, it is necessary to finely discretize the frequency interval, leading to a very large number of state equations. Then, if a standard adjoint-based approach is used for optimization, the computational cost (both in terms of CPU and memory storage) may be prohibitive for large problems, especially in three space dimensions. The goal of the present work is to introduce two new non-adjoint approaches for dealing with frequency response problems in shape and topology optimization. In both cases, we rely on a classical modal basis approach to compute the states, solutions of the direct problems. In the first method, we do not use any adjoint but rather directly compute the shape derivatives of the eigenmodes in the modal basis. In the second method, we compute the adjoints of the standard approach by using again the modal basis. The numerical cost of these two new strategies are much smaller than the usual ones if the number of modes in the modal basis is much smaller than the number of discretized excitation frequencies. We present numerical examples for the minimization of the dynamic compliance in two and three space dimensions. This article is protected by copyright. All rights reserved.

This paper presents a level set-based shape and topology optimization method for conceptual design of cast parts. In order to be successfully manufactured by the casting process, the geometry of cast parts should satisfy certain moldability conditions, which poses additional constraints in the shape and topology optimization of cast parts. Instead of using the originally point-wise constraint statement, we propose a casting constraint in the form of domain integration over a narrowband near the material boundaries. This constraint is expressed in terms of the gradient of the level set function defining the structural shape and topology. Its explicit and analytical form facilitates the sensitivity analysis and numerical implementation. As compared with the standard implementation of the level set method based on the steepest descent algorithm, the proposed method uses velocity field design variables and combines the level set method with the gradient-based mathematical programming algorithm on the basis of the derived sensitivity scheme of the objective function and the constraints. This approach is able to simultaneously account for the casting constraint and the conventional material volume constraint in a convenient way. In this method, the optimization process can be started from an arbitrary initial design, without the need for an initial design satisfying the cast constraint. Numerical examples in both 2D and 3D design domain are given to demonstrate the validity and effectiveness of the proposed method.

This paper develops an enhanced response sensitivity approach for structural damage identification. The whole work is mainly two-fold. Firstly, the general response sensitivity approach has shown to perform well for small damage, but may not work for moderate or large damage. Therefore, to overcome this drawback, the trust-region restriction is additionally enforced to enhance the general response sensitivity approach. In doing so, the Tikhonov regularization is invoked which is simple for practical manipulation and is shown equivalent to trust-region considerations. Secondly, concrete convergence analysis is presented to guarantee the performance of the enhanced response sensitivity approach. Numerical examples are studied to verify the proposed approach. Copyright © 2016 John Wiley & Sons, Ltd.

An accurate spectral-sampling surface method for the vibration analysis of 2-D curved beams with variable curvatures and general boundary conditions is presented. The method combines the advantages of the sampling surface method and spectral method. The formulation is based on the 2-D elasticity theory, which provides complete accuracy and efficiency for curved beams with arbitrary thicknesses and variable curvatures because no other assumptions on the deformations and stresses along the thickness direction are introduced. Specifically, a set of non-equally spaced sampling surfaces parallel to the beam's middle surface are primarily collocated along the thickness direction, and the displacements of these surfaces are chosen as fundamental beam unknowns. This fact provides an opportunity to derive elasticity solutions for thick beams with a prescribed accuracy by selecting sufficient sampling surfaces. Each of the fundamental beam unknowns is then invariantly expanded as Chebyshev polynomials of the first kind, and the problems are stated in variational form with the aid of the penalty technique and Lagrange multipliers which provides complete flexibility to describe any arbitrary boundary conditions. Finally, the desired solutions are obtained by the variational operation.

This paper focuses on a new framework for obtaining a non-intrusive (i.e. not requiring projecting of the governing equations onto the reduced basis modes) reduced order model for 2*D* fluid problems. To overcome the shortcomings of intrusive model order reduction usually derived by combining the POD and the Galerkin projection methods, we developed a novel technique based on Randomized Dynamic Mode Decomposition as a fast and accurate option in model order reduction. Our approach utilizes an adaptive randomized dynamic mode decomposition (ARDMD) to obtain a reduced basis in the offline stage, then the temporal values of the ROM are obtained in the online stage through an interpolation using Radial Basis Functions (RBF). The rank of the reduced DMD model is given as the unique solution of a constrained optimization problem. The Saint-Venant (shallow water) equations in a channel on the rotating earth are employed to provide the numerical data. We emphasize the excellent behavior of the non-intrusive reduced order model by performing a qualitative analysis. In addition, we gain a significantly reduction of CPU time in computation of the reduced order models (ROMs) compared to the classical DMD method. This article is protected by copyright. All rights reserved.

This work deals with the formulation and implementation of a mixed finite element formulation for nonlinear thermoelasticity. In the literature, a ‘consistent’ mixed formulation to reduce spurious oscillations in thermal stresses has been developed which involves the use of a temperature interpolation that is one order lower than the displacement interpolation functions. However, such a strategy would require a very fine mesh in problems where thermal gradients are high. In order to reduce the computational cost, we propose a new hybrid formulation for both the mechanical and thermal parts of the stress tensor, that not only overcomes membrane, shear and volumetric locking but also eliminates thermal stress oscillations, even with the use of a coarse mesh. For transient problems, a new energy-momentum conserving time stepping scheme is also proposed so that linear and angular momenta and energy are conserved exactly in the fully discrete hybrid framework in the absence of loading and dissipation. Several examples for the St. Venant–Kirchhoff and Ogden material models where the solutions are compared against either analytical or other numerical strategies show the efficacy of the developed procedure. This article is protected by copyright. All rights reserved.

In this paper, we present an effective 6-node triangular solid-shell element (MITC-S6), with particular attention on shear locking and thickness locking. To alleviate shear locking, the assumed transverse strain field of the MITC3+ shell element is used while modifying the bending enhancement mechanism. Thickness locking is treated using the assumed and enhanced strain methods for thickness strain. Two independent enhancements of strains are applied: The in-plane and transverse shear strain fields are enhanced using the strain fields obtained from a bubble interpolation function for in-plane translations and the thickness strain field is enhanced for linear variation in the thickness direction. The general three dimensional material law is employed. The proposed element passes all the basic tests including zero-energy mode, patch, and isotropy tests. Excellent performance is observed in various linear and nonlinear benchmark tests, wherein its performance is compared to that of existing 6-node triangular and 8-node quadrilateral solid-shell elements.

The response of a random dynamical system is totally characterized by its probability density function (pdf). However determining a pdf by a direct approach requires a high numerical cost; similarly, surrogate models such as direct polynomial chaos expansions, are not generally efficient, especially around the eigenfrequencies of the dynamical system. In the present study a new approach based on Padé approximants to obtain moments and pdf of the dynamic response in the frequency domain is proposed. A key difference between the direct polynomial chaos representation and Padé representation is that the Padé approach has polynomials in both numerator and denominator. For frequency response functions, the denominator plays a vital role as it contains the information related to resonance frequencies, which are uncertain. A Galerkin approach in conjunction with polynomial chaos is proposed for the Padé approximation. Another physics based approach, utilizing polynomial chaos expansions of the random eigenmodes is proposed and compared with the proposed Padé approach. It is shown that both methods give accurate results even if a very low degree of the polynomial expansion is used. The methods are demonstrated for two degree of freedom system with one and two uncertain parameters. This article is protected by copyright. All rights reserved.

For a Mindlin-Reissner plate subjected to transverse loadings, the distributions of the rotations and some resultant forces may vary very sharply within a narrow district near certain boundaries. This edge effect is indeed a great challenge for conventional finite element analysis. Recently, an effective hybrid displacement function (HDF) finite element method was successfully developed for solving such difficulty [1, 2]. Although good performances can be obtained in most cases, the distribution continuity of some resulting resultants is destroyed when coarse meshes are employed. Moreover, an additional local coordinate system must be used for avoiding singular problem in matrix inversion, which makes the derivations more complicated. Based on a modified complementary energy functional containing Lagrangian multipliers, an improved HDF (IHDF) element scheme is proposed in this work. And two new special IHDF elements, named by IHDF-P4-Free and IHDF-P4-SS1, are constructed for modeling plate behaviors near free and soft simply-supported (SS1) boundaries, respectively. The present modeling scheme not only greatly improves the precision of the numerical results, but also avoids usage of the additional local Coordinate system. The numerical tests demonstrate that the new IHDF element scheme is an effective way for solving the challenging edge effect problem in Mindlin-Reissner plates.

An upscaling approach for multi-porosity media is developed and demonstrated on a three-scale porous medium. The three-scale porous medium problem is analyzed using the so-called MultiSPH approach, which is a sequential SPH solver at multiple scales. By this approach, upscaling introduces drag forces on the SPH particles at coarser scales. The resulting permeabilities obtained by the MultiSPH approach are validated against available analytical solutions for simple microstructures. For complex microstructures, the MultiSPH approach is validated against detailed single-scale SPH simulations and experimental data. Sensitivities of fluid-structure interaction on the permeability is investigated.

In this paper we propose a numerical method based on an active structural control to dampen the discrete critical frequencies from a fe/be coupled system, which models a time harmonic fluid-structure interaction problem. The active structural acoustic control (ASAC) used consists of applying external forces directly on the structure, so that, in the presence of critical frequencies the linear system remains stable, and in turn, in the presence of non-critical frequencies, the active control does not alter the system. We consider the problem in the framework of the theory of optimal control and present bi-dimensional numerical simulations to show the behavior of the scheme in some vibroacoustic structures. This article is protected by copyright. All rights reserved.

Linear elasticity problems posed on cracked domains, or domains with re-entrant corners, yield singular solutions that deteriorate the optimality of convergence of finite element methods. In this work, we propose an optimally convergent finite element method for this class of problems. The method is based on approximating a much smoother function obtained by locally reparameterizing the solution around the singularities. This reparameterized solution can be approximated using standard finite element procedures yielding optimal convergence rates for any order of interpolating polynomials, without additional degrees of freedom or special shape functions. Hence the method provides optimally convergent solutions for the same computational complexity of standard finite element methods. Furthermore, the sparsity and the conditioning of the resulting system are preserved. The method handles body forces and crack-face tractions, as well as multiple crack tips and re-entrant corners. The advantages of the method are showcased for four different problems: a straight crack with loaded faces, a circular arc crack, an L-shaped domain undergoing anti-plane deformation, and lastly a crack along a bimaterial interface. Optimality in convergence is observed for all the examples. A proof of optimal convergence is accomplished mainly by proving the regularity of the reparameterized solution. This article is protected by copyright. All rights reserved.

A new approach based on the use of the Newton and level set methods allows to follow the motion of interfaces with surface tension immersed in an incompressible Newtonian fluid. Our method features the use of a high order fully implicit time integration scheme that circumvents the stability issues related to the explicit discretization of the capillary force when capillary effects dominate. A strategy based on a consistent Newton-Raphson linearization is introduced, and performances are enhanced by using an exact Newton variant which guarantees a third-order convergence behavior without requiring second order derivatives. The problem is approximated by mixed finite elements, while the anisotropic adaptive mesh refinements enable us to increase the computational accuracy. Numerical investigations of the convergence properties and comparisons with benchmark results provide evidence regarding the efficacy of the methodology. The robustness of the method is tested with respect to the standard explicit method, and stability is maintained for significantly larger time steps compared to those allowed by the stability condition. This article is protected by copyright. All rights reserved.

This work investigates the use of hierarchical mesh decomposition strategies for topology optimization using bi-directional evolutionary structural optimization (BESO) algorithm. The proposed method uses a dual mesh system which decouples the design variables from the finite element analysis mesh. The investigation focuses on previously unexplored areas of these techniques to investigate the effect of five meshing parameters on the analysis solving time (i.e. computational effort) and the analysis quality (i.e. solution optimality). The foreground mesh parameters, including adjacency ratio, and minimum and maximum element size, were varied independently across solid and void domain regions. Within the topology optimization, strategies for controlling the mesh parameters were investigated. The differing effects of these parameters on the efficiency and efficacy of the analysis and optimization stages are discussed and recommendations made for parameter combinations. Some of the key findings were that increasing the adjacency ratio increased the efficiency only modestly; the largest effect was for the minimum and maximum element size parameters, and that the most dramatic reduction in solve time can be achieved by not setting the minimum element size too low, assuming mapping onto a background mesh with a minimum element size of 1. This article is protected by copyright. All rights reserved.

The most used method for calculation of mechanical properties of the fibre reinforced composites on the scale of unit cell of the reinforcement consists in continuous meshing of the reinforcement and matrix volumes. Often this leads to bad quality or large size of the generated mesh. Mesh superposition techniques (MSP) allow overcoming these issues by independent meshing of both components with coupling them together based on constraint equations. In present work two different implementations of MSP techniques are investigated in detail and results are compared with common method and experiment for a simple unidirectional reinforcement geometry and a complex 3D weave composite. An MSP formulation for large deformations case is described. Results show good agreement in local strain distributions and in predicting of the homogenized mechanical properties.

The study of molecular flows at low Knudsen numbers (~0.1-0.5), over nano-scaled objects of 20-100 nm size is becoming an important area of research. The simulation of fluid-structure interaction at nano-scale is important for understanding the adsorption and drag resistance characteristics of nano devices in the fields of drug delivery, surface cleaning and protein movement.

A novel formulation has been proposed that calculates localised values for both the kinetic and configurational parts of the Irwin-Kirkwood stress tensor at given fixed positions within the computational domain.

Macroscopic properties, such as streaming velocity, pressure and drag coefficients are predicted by modelling the fluid-structure interaction using a moving least-squares method. The gravitation driven molecular flow is examined over three different cross sectional shapes, i.e. diamond, circular and square shaped cylinders; confined within parallel walls, has been simulated for rough and smooth surfaces.

The molecular dynamics formulation has allowed, for the first time, the calculation of localised drag forces over nano-cylinders. The computational simulation has shown that existing methods, including continuum based approaches, significantly underestimate drag coefficients over nano-cylinders. The proposed molecular dynamics formulation has been verified on simulation based tests, as experimental and analytical results are unavailable at this scale.

For the two-dimensional three-temperature (2-D 3-T) radiative heat conduction problem appearing in the inertial confinement numerical stimulations, we choose the *Freezing coefficient method* (FCM) to linearize the nonlinear equations, and initially apply the well-known mixed finite element scheme with the lowest order Raviart-Thomas element associated with the triangulation to the linearized equations, and obtain the convergence with one order with respect to the space direction for the temperature and flux function approximations, and design a simple but efficient algorithm for the discrete system. Three numerical examples are displayed. The former two verify theoretical results and show the super-convergence for temperature and flux functions at the barycenter of the element, which is helpful for solving the radiative heat conduction problems. The third validates the robustness of this scheme with small energy conservative error and one order convergence for the time discretization. This article is protected by copyright. All rights reserved.

This work investigates a model reduction method applied to coupled multi-physics systems. The case in which a system of interest interacts with an external system is considered. An approximation of the PoincaréSteklov operator is computed by simulating, in an offline phase, the external problem when the inputs are the Laplace-Beltrami eigenfunctions defined at the interface. In the online phase, only the reduced representation of the operator is needed to account for the influence of the external problem on the main system. An online basis enrichment is proposed in order to guarantee a precise reduced-order computation. Several test-cases are proposed on different fluid-structure couplings. This article is protected by copyright. All rights reserved.

We propose several algorithms to recover the location and intensity of a radiation source located in a simulated 250 m ×180 m block of an urban center based on synthetic measurements. Radioactive decay and detection are Poisson random processes, so we employ likelihood functions based on this distribution. Due to the domain geometry and the proposed response model, the negative logarithm of the likelihood is only piecewise continuous differentiable, and it has multiple local minima. To address these difficulties, we investigate three hybrid algorithms comprised of mixed optimization techniques. For global optimization, we consider Simulated Annealing (SA), Particle Swarm (PS) and Genetic Algorithm (GA), which rely solely on objective function evaluations; i.e., they do not evaluate the gradient in the objective function. By employing early stopping criteria for the global optimization methods, a pseudo-optimum point is obtained. This is subsequently utilized as the initial value by the deterministic Implicit Filtering method (IF), which is able to find local extrema in non-smooth functions, to finish the search in a narrow domain. These new hybrid techniques, combining global optimization and Implicit Filtering address, difficulties associated with the non-smooth response, and their performances are shown to significantly decrease the computational time over the global optimization methods. To quantify uncertainties associated with the source location and intensity, we employ the Delayed Rejection Adaptive Metropolis (DRAM) and DiffeRential Evolution Adaptive Metropolis (DREAM) algorithms. Marginal densities of the source properties are obtained, and the means of the chains' compare accurately with the estimates produced by the hybrid algorithms. This article is protected by copyright. All rights reserved.

This work provides a robust variational-based numerical implementation of a phase field model of ductile fracture in elastic-plastic solids undergoing large strains. This covers a computationally efficient micromorphic regularization of the coupled gradient plasticity-damage formulation recently proposed in Miehe et al. [1, 2, 3, 4]. The phase field approach regularizes sharp crack surfaces within a pure continuum setting by a specific gradient damage modeling with geometric features rooted in fracture mechanics. It has proven immensely successful with regard to the analysis of complex crack topologies without the need for fracture-specific computational structures such as finite element design of crack discontinuities or intricate crack-tracking algorithms. The proposed gradient-extended plasticity-damage formulation includes two independent length scales which regularize both the plastic response as well as the crack discontinuities. This ensures that the damage zones of ductile fracture are inside of plastic zones or vice versa, and guarantees on the computational side a mesh objectivity in post-critical ranges. The proposed setting is rooted in a canonical variational principle. The coupling of gradient plasticity to gradient damage is realized by a constitutive work density function that includes the stored elastic energy and the dissipated work due to plasticity and fracture. The latter represents a coupled resistance to plasticity and damage, depending on the gradient-extended internal variables which enter plastic yield functions and fracture threshold functions. With this viewpoint on the generalized internal variables at hand, the thermodynamic formulation is outlined for gradient-extended dissipative solids with generalized internal variables which are passive in nature. It is specified for a conceptual model of von Mises-type elasto-plasticity at finite strains coupled with fracture. The canonical theory proposed is shown to be governed by a rate-type minimization principle, which fully determines the coupled multi-field evolution problem. This is exploited on the numerical side by a fully symmetric monolithic finite element implementation. An important aspect of this work is the regularization towards a micromorphic gradient plasticity-damage setting by taking into account additional internal variable fields linked to the original ones by penalty terms. This enhances the robustness of the finite element implementation, in particular on the side of gradient plasticity. The performance of the formulation is demonstrated by means of some representative examples. This article is protected by copyright. All rights reserved.

In this contribution a three-dimensional frictional contact formulation in the form of an interface solid finite element is discussed. The interface element incorporates normal as well as tangential contact conditions and aims at the modeling of phenomena occurring in engineering problems on the interface between different materials. The focus of this work lies on the description of the tangential interaction of bodies on their contact interface. Instead of the widely used Coulomb friction law, a simple elastoplastic material law is applied in order to model the interaction mechanism on the tangent plane. Appropriate local coordinate systems are introduced on the surface and point of contact and all contact operations are performed in these coordinate systems by adopting a fully covariant description. In this way we derive a symmetric element stiffness matrix not only for the case of sticking but also for sliding. Exploiting the advantages of the covariant description we show that the fact that the derived stiffness matrix discretely consists of constitutive and geometrical parts enables us to introduce rather straightforwardly new appropriate interface laws describing the tangential interaction of the contacting bodies. The presented numerical results show very good correlation of the developed approach with “pull-out” experiments. This article is protected by copyright. All rights reserved.

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.

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.

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.

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 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.

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.

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. Because 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.

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 simulation (LES). The proposed hybrid approach combines the advantages of the inverse distance weighting (IDW) interpolation with the simplicity and low computational effort of transfinite interpolation (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 not only their capabilities but also their computational costs. That not only 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 that 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 that 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.

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 Messerschmitt–Bolkow–Blohm 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.

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. Copyright © 2016 John Wiley & Sons, Ltd.

The representation of structural boundaries is of significant importance for the applications of the discrete element method in an industry. In recent decades, 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 complicated structural boundaries containing regular shapes and irregular shapes and greatly reduce computational cost in the contact detection between particles and structural boundaries. Copyright © 2016 John Wiley & Sons, Ltd.

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, for example, 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. 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. Copyright © 2016 John Wiley & Sons, Ltd.

Shape-memory polymers (SMPs) belong to a class of smart materials that 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 photomechanical behavior of light activated shape-memory polymers (LASMPs), focusing on the numerical aspects for finite element simulations at the engineering scale. The photomechanical 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-dimensional 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 photomechanical loading. Copyright © 2016 John Wiley & Sons, Ltd.

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 through the derivatives along the nodal line which involve only 1D interpolations. In this paper, generic schemes for formulating assumed natural strain and stabilized Lobatto Lagrange C^{0} plate/shell elements of order
2 are presented. Two assumed natural strain 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. Copyright © 2016 John Wiley & Sons, Ltd.

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. 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.

No abstract is available for this article.

]]>We propose a numerical method for a fluid–structure interaction problem. The material of the structure is homogeneous, isotropic, and it can be described by the compressible neo-Hookean constitutive equation, while the fluid is governed by the Navier–Stokes equations. Our study does not use turbulence model. Updated Lagrangian method is used for the structure and fluid equations are written in Arbitrary Lagrangian–Eulerian coordinates. One global moving mesh is employed for the fluid–structure domain, where the fluid–structure interface is an ‘interior boundary’ of the global mesh. At each time step, we solve a monolithic system of unknown velocity and pressure defined on the global mesh. The continuity of velocity at the interface is automatically satisfied, while the continuity of stress does not appear explicitly in the monolithic fluid–structure system. This method is very fast because at each time step, we solve only one linear system. This linear system was obtained by the linearization of the structure around the previous position in the updated Lagrangian formulation and by the employment of a linear convection term for the fluid. Numerical results are presented. Copyright © 2016 John Wiley & Sons, Ltd.

The identification of the geological structure from seismic data is formulated as an inverse problem. The properties and the shape of the rock formations in the subsoil are described by material and geometric parameters, which are taken as input data for a predictive model. Here, the model is based on the Helmholtz equation, describing the acoustic response of the system for a given wave length. Thus, the inverse problem consists in identifying the values of these parameters such that the output of the model agrees the best with observations. This optimization algorithm requires multiple queries to the model with different values of the parameters. Reduced order models are especially well suited to significantly reduce the computational overhead of the multiple evaluations of the model.

In particular, the proper generalized decomposition produces a solution explicitly stating the parametric dependence, where the parameters play the same role as the physical coordinates. A proper generalized decomposition solver is devised to inexpensively explore the parametric space along the iterative process. This exploration of the parametric space is in fact seen as a post-process of the generalized solution. The approach adopted demonstrates its viability when tested in two illustrative examples. Copyright © 2016 John Wiley & Sons, Ltd.

This paper presents a new efficient monolithic finite element solution scheme to treat the set of PDEs governing a 2D, biphasic, saturated theory of porous media model with intrinsically coupled and incompressible solid and fluid constituents for infinitesimal and large elastic deformation. Our approach, which inherits some of its techniques from CFD, is characterized by the following aspects: (1) it only performs operator evaluation with no additional Gateaux derivatives. In particular, the computations of the time-consuming material tangent matrix are not involved here; (2) it solves the non-linear dynamic problem with no restriction on the strength of coupling; (3) it is more efficient than the linear **u****v***p* solver discussed in previous works; (4) it requires weaker derivatives, and hence, lower-order FE can be tested; and (5) the boundary conditions are reduced, solution independent and more convenient to apply than in the old **u****v***p* formulation. For the purpose of validation and comparison, prototypical simulations including analytical solutions are carried out, and at the end, an adaptive time stepping procedure is introduced to handle the rapid change in the numbers of nonlinear iterations that may occur. Copyright © 2016 John Wiley & Sons, Ltd.

Among all 3D 8-node hexahedral solid elements in current finite element library, the ‘best’ one can produce good results for bending problems using coarse regular meshes. However, once the mesh is distorted, the accuracy will drop dramatically. And how to solve this problem is still a challenge that remains outstanding. This paper develops an 8-node, 24-DOF (three conventional DOFs per node) hexahedral element based on the virtual work principle, in which two different sets of displacement fields are employed simultaneously to formulate an unsymmetric element stiffness matrix. The first set simply utilizes the formulations of the traditional 8-node trilinear isoparametric element, while the second set mainly employs the analytical trial functions in terms of 3D oblique coordinates (*R*, *S*, *T*). The resulting element, denoted by US-ATFH8, contains no adjustable factor and can be used for both isotropic and anisotropic cases. Numerical examples show it can strictly pass both the first-order (constant stress/strain) patch test and the second-order patch test for pure bending, remove the volume locking, and provide the invariance for coordinate rotation. Especially, it is insensitive to various severe mesh distortions. Copyright © 2016 John Wiley & Sons, Ltd.

This paper presents a novel class of preconditioners for the iterative solution of the sequence of symmetric positive-definite linear systems arising from the numerical discretization of transient parabolic and self-adjoint partial differential equations. The preconditioners are obtained by nesting appropriate projections of reduced-order models into the classical iteration of the preconditioned conjugate gradient (PCG). The main idea is to employ the reduced-order solver to project the residual associated with the conjugate gradient iterations onto the space spanned by the reduced bases. This approach is particularly appealing for transient systems where the full-model solution has to be computed at each time step. In these cases, the natural reduced space is the one generated by full-model solutions at previous time steps. When increasing the size of the projection space, the proposed methodology highly reduces the system conditioning number and the number of PCG iterations at every time step. The cost of the application of the preconditioner linearly increases with the size of the projection basis, and a trade-off must be found to effectively reduce the PCG computational cost. The quality and efficiency of the proposed approach is finally tested in the solution of groundwater flow models. © 2016 The Authors. *International Journal for Numerical Methods in Engineering* Published by John Wiley & Sons Ltd.

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.