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. 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, since 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. This article is protected by copyright. All rights reserved.

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, while 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. This article is protected by copyright. All rights reserved.

In this paper, a new discrete elements generation method based on geometry is proposed to fill geometric domains with particles (discs 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 (PPRIM) and a particle stability inspection and improvement method (PSIIM) 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 PPRIM 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. This article is protected by copyright. All rights reserved.

In this paper, we model crack discontinuities in two-dimensional linear elastic continua using the extended finite element mssethod without the need to partition an enriched element into a collection of triangles or quadrilaterals. For crack modeling in the X-FEM, 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 (due to 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. This article is protected by copyright. All rights reserved.

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. This article is protected by copyright. All rights reserved.

In isogeometric analysis identical basis functions are used for geometrical representation and analysis. In this work non-uniform rational B-splines (NURBS) 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 are 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 to experimental results and alternative numerical approaches. This article is protected by copyright. All rights reserved.

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. This article is protected by copyright. All rights reserved.

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 (CSPM) and Total Lagrangian Formulation (TLF). The compliance is established relative to a designed density variable at each SPH particle which is updated by Optimality Criteria (OC) method. Topology optimization is realized by minimizing the compliance using a modified Solid Isotropic Material with Penalization (SIMP) 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. This article is protected by copyright. All rights reserved.

This article is concerned with the identification of stiffness parameters for short fiber-reinforced plastic (SFRP) 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 (RVE) of the material. The goal is to robustly identify the parameters with the fewest number of RVE 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 & conductivity and structurally similar problems. This article is protected by copyright. All rights reserved.

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. In [1] 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. This article is protected by copyright. All rights reserved.

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.

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 one-dimensional 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. This article is protected by copyright. All rights reserved.

Represented by the element free Galerkin method (EFGM), 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.

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 XFEM based enrichment of the displacement field to allow for cracks 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. This article is protected by copyright. All rights reserved.

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 to 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 to 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. This article is protected by copyright. All rights reserved.

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. This article is protected by copyright. All rights reserved.

This work focuses on providing accurate low-cost approximations of stochastic finite elements simulations in the framework of linear elasticity. In [E. Florentin, P. Diez, Adaptive reduced basis strategy based on goal oriented error assessment for stochastic problems, Comput. Methods Appl. Mech. Engrg. 225-228 (2012) 116-127], an adaptive strategy has been 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 cost gain. Technical insight of these two topics are presented in details and numerical examples show the interest of these new developments. This article is protected by copyright. All rights reserved.

In this paper, an effective 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. This article is protected by copyright. All rights reserved.

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.

Atomistic models, which are crucial for performing molecular dynamics (MD) simulations of carbon nanostructures (CNS), 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 defects in the hexagonal mesh) should be properly arranged around the boundaries of the different components of a CNS. Two techniques play important roles in our method: (1) sphere packing is used to place the nodes for triangulation, which 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 hand-work method, thereby reducing the operation time significantly.

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. This article is protected by copyright. All rights reserved.

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 (FEM). 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. This article is protected by copyright. All rights reserved.

We introduce a method to mesh the boundary *Γ* of a smooth, open domain in
*R*^{3} 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; i.e., 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, since 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 (see [1] ) 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. This article is protected by copyright. All rights reserved.

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 (DG) 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 zero to four, 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 DG schemes. This article is protected by copyright. All rights reserved.

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. This article is protected by copyright. All rights reserved.

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, 5 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 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. This article is protected by copyright. All rights reserved.

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. This article is protected by copyright. All rights reserved.

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

We propose a method to couple Smoothed Particle Hydrodynamics (SPH) and Finite Elements (FEM) methods for non-linear 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. This article is protected by copyright. All rights reserved.

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: (i) it is simple; (ii) it has guaranteed stability; (iii) it is an efficient non-iterative adaptive single-step procedure; (iv) it provides enhanced accuracy; (v) it enables advanced controllable algorithmic dissipation in the higher modes; (vi) it is truly self-starting; (vii) 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.

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 nano-scale 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 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 nano-scale 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 FCC 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 modeling of Golden crystals under different loading conditions. This article is protected by copyright. All rights reserved.

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

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

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

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

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

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

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

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

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

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

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

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

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

A reduced order model (ROM) is presented for the long-term calculation of subsurface oil/water flows. As in several previous ROMs in the field, the Newton iterations in the full model (FM) equations, which are implicit in time, are projected onto a set of modes obtained by applying proper orthogonal decomposition (POD) to a set of snapshots computed by the FM itself. The novelty of the present ROM is that the POD modes are (i) first calculated from snapshots computed by the FM in a short initial stage, and then (ii) updated on the fly along the simulation itself, using new sets of snapshots computed by the FM in even shorter additional runs. Thus, the POD modes adapt themselves to the local dynamics along the simulation, instead of being completely calculated at the outset, which requires a computationally expensive preprocess. This strategy is robust and computationally efficient, which is tested in 10- and 30-year simulations for a realistic reservoir model taken from the SAIGUP project. Copyright © 2016 John Wiley & Sons, Ltd.

This article presents a family of variational integrators from a continuous time point of view. A general procedure for deriving symplectic integration schemes preserving an energy-like quantity is shown, which is based on the principle of virtual work. The framework is extended to incorporate holonomic constraints without using additional regularization. In addition, it is related to well-known partitioned Runge–Kutta methods and to other variational integration schemes. As an example, a concrete integration scheme is derived for the planar pendulum using both polar and Cartesian coordinates. Copyright © 2016 John Wiley & Sons, Ltd.

We study the alternative ‘simultaneous analysis and design’ (SAND) formulation of the local stress-constrained and slope-constrained topology design problem. It is demonstrated that a standard trust-region Lagrange–Newton sequential quadratic programming-type algorithm—based, in this case, on strictly convex and separable approximate subproblems—may converge to singular optima of the local stress-constrained problem without having to resort to relaxation or perturbation techniques. Moreover, because of the negation of the sensitivity analyses—in SAND, the density and displacement variables are independent—and the immense sparsity of the SAND problem, solutions to large-scale problem instances may be obtained in a reasonable amount of computation time. Copyright © 2016 John Wiley & Sons, Ltd.

Finite element method (FEM) is a well-developed method to solve real-world problems that can be modeled with differential equations. As the available computational power increases, complex and large-size problems can be solved using FEM, which typically involves multiple degrees of freedom (DOF) per node, high order of elements, and an iterative solver requiring several sparse matrix-vector multiplication operations. In this work, a new storage scheme is proposed for sparse matrices arising from FEM simulations with multiple DOF per node. A sparse matrix-vector multiplication kernel and its variants using the proposed scheme are also given for CUDA-enabled GPUs. The proposed scheme and the kernels rely on the mesh connectivity data from FEM discretization and the number of DOF per node. The proposed kernel performance was evaluated on seven test matrices for double-precision floating point operations. The performance analysis showed that the proposed GPU kernel outperforms the ELLPACK (ELL) and CUSPARSE Hybrid (HYB) format GPU kernels by an average of 42% and 32%, respectively, on a Tesla K20c card. Copyright © 2016 John Wiley & Sons, Ltd.

The current article presents a Lagrangian cell-centred finite volume solution methodology for simulation of metal forming processes. Details are given of the mathematical model in updated Lagrangian form, where a hyperelastoplastic *J*_{2} constitutive relation has been employed. The cell-centred finite volume discretisation is described, where a modified discretised is proposed to alleviate erroneous hydrostatic pressure oscillations; an outline of the memory efficient segregated solution procedure is given. The accuracy and order of accuracy of the method are examined on a number of 2-D and 3-D elastoplastic benchmark test cases, where good agreement with available analytical and finite element solutions is achieved. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, we develop a robust reduced-order modeling method, named algebraic dynamic condensation, which is based on the improved reduced system method. Using algebraic substructuring, the global mass and stiffness matrices are divided into many small submatrices without considering the physical domain, and substructures and interface boundary are defined in the algebraic perspective. The reduced model is then constructed using three additional procedures: substructural stiffness condensation, interface boundary reduction, and substructural inertial effect condensation. The formulation of the reduced model is simply expressed at a submatrix level without using a transformation matrix that induces huge computational cost. Through various numerical examples, the performance of the proposed method is demonstrated in terms of accuracy and computational cost. Copyright © 2016 John Wiley & Sons, Ltd.

This work proposes novel stability-preserving model order reduction approaches for vibro-acoustic finite element models. As most research in the past for these systems has focused on noise attenuation in the frequency-domain, stability-preserving properties were of low priority. However, as the interest for time-domain auralization and (model based) active noise control increases, stability-preserving model order reduction techniques are becoming indispensable. The original finite element models for vibro-acoustic simulation are already well established but require too much computational load for these applications. This work therefore proposes two new global approaches for the generation of stable reduced-order models. Based on proven conditions for stability preservation under one-sided projection, a reformulation of the displacement-fluid velocity potential (*u* − *ϕ*) formulation is proposed. In contrast to the regular formulation, the proposed approach leads to a new asymmetric structure for the system matrices which is proven to preserve stability under one-sided projection. The second approach starts from a displacement-pressure (*u* − *p*) description where the system level projection space is decoupled for the two domains, for which we also prove the preservation of stability. Two numerical validation cases are presented which demonstrate the inadequacy of straightforward model order reduction on typical vibro-acoustic models for time-domain simulation and compare the performance of the proposed approaches. Both proposed approaches effectively preserve the stability of the original system. Copyright © 2016 John Wiley & Sons, Ltd.

The strain smoothing technique over higher order elements and arbitrary polytopes yields less accurate solutions than other techniques such as the conventional polygonal finite element method. In this work, we propose a linear strain smoothing scheme that improves the accuracy of linear and quadratic approximations over convex polytopes. The main idea is to subdivide the polytope into simplicial subcells and use a linear smoothing function in each subcell to compute the strain. This new strain is then used in the computation of the stiffness matrix. The convergence properties and accuracy of the proposed scheme are discussed by solving a few benchmark problems. Numerical results show that the proposed linear strain smoothing scheme makes the approximation based on polytopes able to deliver the same optimal convergence rate as traditional quadrilateral and hexahedral approximations. The accuracy is also improved, and all the methods tested pass the patch test to machine precision. Copyright © 2016 John Wiley & Sons, Ltd.

When geometric uncertainties arising from manufacturing errors are comparable with the characteristic length or the product responses are sensitive to such uncertainties, the products of deterministic design cannot perform robustly. This paper presents a new level set-based framework for robust shape and topology optimization against geometric uncertainties. We first propose a stochastic level set perturbation model of uncertain topology/shape to characterize manufacturing errors in conjunction with Karhunen–Loève (K–L) expansion. We then utilize polynomial chaos expansion to implement the stochastic response analysis. In this context, the mathematical formulation of the considered robust shape and topology optimization problem is developed, and the adjoint-variable shape sensitivity scheme is derived. An advantage of this method is that relatively large shape variations and even topological changes can be accounted for with desired accuracy and efficiency. Numerical examples are given to demonstrate the validity of the present formulation and numerical techniques. In particular, this method is justified by the observations in minimum compliance problems, where slender bars vanish when the manufacturing errors become comparable with the characteristic length of the structures. Copyright © 2016 John Wiley & Sons, Ltd.

Discretization-induced oscillations in the load–displacement curve are a well-known problem for simulations of cohesive crack growth with finite elements. The problem results from an insufficient resolution of the complex stress state within the cohesive zone ahead of the crack tip. This work demonstrates that the *hp*-version of the finite element method is ideally suited to resolve this complex and localized solution characteristic with high accuracy and low computational effort. To this end, we formulate a local and hierarchic mesh refinement scheme that follows *dynamically* the propagating crack tip. In this way, the usually applied *static* a priori mesh refinement along the complete potential crack path is avoided, which significantly reduces the size of the numerical problem. Studying systematically the influence of *h*-refinement, *p*-refinement, and *h**p*-refinement, we demonstrate why the suggested *h**p*-formulation allows to capture accurately the complex stress state at the crack front preventing artificial snap-through and snap-back effects. This allows to decrease significantly the number of degrees of freedom and the simulation runtime. Furthermore, we show that by combining this idea with the finite cell method, the crack propagation within complex domains can be simulated efficiently without resolving the geometry by the mesh. Copyright © 2016 John Wiley & Sons, Ltd.

The FFT-based homogenization method of Moulinec–Suquet has recently emerged as a powerful tool for computing the macroscopic response of complex microstructures for elastic and inelastic problems. In this work, we generalize the method to problems discretized by trilinear hexahedral elements on Cartesian grids and physically nonlinear elasticity problems. We present an implementation of the basic scheme that reduces the memory requirements by a factor of four and of the conjugate gradient scheme that reduces the storage necessary by a factor of nine compared with a naive implementation.

For benchmark problems in linear elasticity, the solver exhibits mesh- and contrast-independent convergence behavior and enables the computational homogenization of complex structures, for instance, arising from computed tomography computed tomography (CT) imaging techniques. There exist 3D microstructures involving pores and defects, for which the original FFT-based homogenization scheme does not converge. In contrast, for the proposed scheme, convergence is ensured. Also, the solution fields are devoid of the spurious oscillations and checkerboarding artifacts associated to conventional schemes. We demonstrate the power of the approach by computing the elasto-plastic response of a long-fiber reinforced thermoplastic material with 172 × 10^{6} (displacement) degrees of freedom. Copyright © 2016 John Wiley & Sons, Ltd.

This work introduces a semi-analytical formulation for the simulation and modeling of curved structures based on the scaled boundary finite element method (SBFEM). This approach adapts the fundamental idea of the SBFEM concept to scale a boundary to describe a geometry. Until now, scaling in SBFEM has exclusively been performed along a straight coordinate that enlarges, shrinks, or shifts a given boundary. In this novel approach, scaling is based on a polar or cylindrical coordinate system such that a boundary is shifted along a curved scaling direction. The derived formulations are used to compute the static and dynamic stiffness matrices of homogeneous curved structures. The resulting elements can be coupled to general SBFEM or FEM domains. For elastodynamic problems, computations are performed in the frequency domain. Results of this work are validated using the global matrix method and standard finite element analysis. Copyright © 2016 John Wiley & Sons, Ltd.

We propose a new method to obtain contact forces under a non-smoothed contact problem between arbitrarily-shaped bodies which are discretized by finite element method. Contact forces are calculated by the specific contact algorithm between two particles of smoothed particle hydrodynamics, which is a meshfree method, and that are applied to each colliding body. This approach has advantages that accurate contact forces can be obtained within an accelerated collision without a jump problem in a discrete time increment. Also, this can be simply applied into any contact problems like a point-to-point, a point-to-line, and a point-to-surface contact for complex shaped and deformable bodies. In order to describe this method, an impulse based method, a unilateral contact method and smoothed particle hydrodynamics method are firstly introduced in this paper. Then, a procedure about the proposed method is handled in great detail. Finally, accuracy of the proposed method is verified by a conservation of momentum through three contact examples. Copyright © 2016 John Wiley & Sons, Ltd.

This paper considers stochastic hybrid stress quadrilateral finite element analysis of plane elasticity equations with stochastic Young's modulus and stochastic loads. Firstly, we apply Karhunen–Loève expansion to stochastic Young's modulus and stochastic loads so as to turn the original problem into a system containing a finite number of deterministic parameters. Then we deal with the stochastic field and the space field by *k*−version/*p*−version finite element methods and a hybrid stress quadrilateral finite element method, respectively. We derive a priori error estimates, which are uniform with respect to the Lamè constant *λ*∈(0,+*∞*). Finally, we provide some numerical results. Copyright © 2016 John Wiley & Sons, Ltd.

This paper builds on recent work developed by the authors for the numerical analysis of large strain solid dynamics, by introducing an upwind cell centred hexahedral finite volume framework implemented within the open source code OpenFOAM [http://www.openfoam.com/]. In Lee, Gil and Bonet (2013), a first-order hyperbolic system of conservation laws was introduced in terms of the linear momentum and the deformation gradient tensor of the system, leading to excellent behaviour in two-dimensional bending dominated nearly incompressible scenarios. The main aim of this paper is the extension of this algorithm into three dimensions, its tailor-made implementation into OpenFOAM and the enhancement of the formulation with three key novelties. First, the introduction of two different strategies in order to ensure the satisfaction of the underlying involutions of the system, that is, that the deformation gradient tensor must be curl-free throughout the deformation process. Second, the use of a discrete angular momentum projection algorithm and a monolithic Total Variation Diminishing Runge–Kutta time integrator combined in order to guarantee the conservation of angular momentum. Third, and for comparison purposes, an adapted Total Lagrangian version of the hyperelastic-GLACE nodal scheme of Kluth and Després (2010) is presented. A series of challenging numerical examples are examined in order to assess the robustness and accuracy of the proposed algorithm, benchmarking it against an ample spectrum of alternative numerical strategies developed by the authors in recent publications. Copyright © 2016 John Wiley & Sons, Ltd.

This paper describes a novel methodology that combines smoothed discrete particle hydrodynamics (SDPH) and finite volume method (FVM) to enhance the effective performance in solving the problems of gas-particle multiphase flow. To describe the collision and fluctuation of particles, this method also increases a new parameter, namely, granular temperature, according to the kinetic theory of granular flow. The coupled framework of SDPH–FVM has been established, in which the drag force and pressure gradient act on the SDPH particles and the momentum sources of drag force are added back onto the FVM mesh. The proposed technique is a coupled discrete-continuum method based on the two-fluid model. To compute for the discrete phase, its SDPH is developed from smoothed particle hydrodynamics (SPH), in which the properties of SPH are redefined with some new physical quantities added into the traditional SPH parameters, so that it is more beneficial for SDPH in representing the particle characteristics. For the continuum phase, FVM is employed to discretize the continuum flow field on a stationary grid by capturing fluid characteristics. The coupled method exhibits strong efficiency and accuracy in several two-dimensional numerical simulations. Copyright © 2016 John Wiley & Sons, Ltd.

We present an adaptive variant of the measure-theoretic approach for stochastic characterization of micromechanical properties based on the observations of quantities of interest at the coarse (macro) scale. The salient features of the proposed nonintrusive stochastic inverse solver are identification of a nearly optimal sampling domain using enhanced ant colony optimization algorithm for multiscale problems, incremental Latin-hypercube sampling method, adaptive discretization of the parameter and observation spaces, and adaptive selection of number of samples. A complete test data of the TORAY T700GC-12K-31E and epoxy #2510 material system from the National Institute for Aviation Research report is employed to characterize and validate the proposed adaptive nonintrusive stochastic inverse algorithm for various unnotched and open-hole laminates. Copyright © 2016 John Wiley & Sons, Ltd.

A new smoothed finite element method (S-FEM) with tetrahedral elements for finite strain analysis of nearly incompressible solids is proposed. The proposed method is basically a combination of the F-bar method and edge-based S-FEM with tetrahedral elements (ES-FEM-T4) and is named ‘F-barES-FEM-T4’. F-barES-FEM-T4 inherits the accuracy and shear locking-free property of ES-FEM-T4. At the same time, it also inherits the volumetric locking-free property of the F-bar method. The isovolumetric part of the deformation gradient (*F*^{iso}) is derived from the ** F** of ES-FEM-T4, whereas the volumetric part (

In this study, a new mean-strain 10-node tetrahedral element is developed using energy-sampling stabilization. The proposed 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 linear tetrahedra, corresponding to three different choices for the internal diagonal. The assumed strains are calculated from mean ‘basis function gradients.’ The energy-sampling technique introduced previously for removing zero-energy modes in the mean-strain hexahedron is adapted for the present element: the stabilization energy is evaluated on the four-corner tetrahedra. The proposed element naturally leads to a lumped-mass matrix and does not have unphysical low-energy vibration modes. For simplicity, we limit our developments to linear elasticity with compressible and nearly incompressible material. The numerical tests demonstrate that the present element performs well compared with the classical 10-node tetrahedral elements for shell and plate structures, and nearly incompressible materials. Copyright © 2016 John Wiley & Sons, Ltd.

Couple stress formulations have been given much attention lately because of the possibility to explain cases, when the classical theory of elasticity fails to describe adequately the mechanical behavior. Such cases may include size-dependent stiffness, high stress gradients, and the response of materials with microstructure. Here, a new mixed Lagrangian formulation is developed for elastodynamic response within consistent size-dependent skew-symmetric couple stress theory. With a specific choice of mixed variables, the formulation can be written with only C^{0} continuity requirements, without the need to introduce a Lagrange multiplier or penalty method. Furthermore, this formulation permits, for the first time, the determination of natural frequencies within the consistent couple stress theory. Details for the strong form of the equilibrium equations, constitutive model relations, boundary conditions, and the corresponding weak form are provided. In addition, the discrete forms are also discussed with two approaches for reducing variables. Several simple two-dimensional computational example problems are then examined, along with a brief investigation of the effect of couple stress on natural frequencies, which exhibit size-dependence for most, but not all, modes. Copyright © 2016 John Wiley & Sons, Ltd.

Contact and fracture in the material point method require grid-scale enrichment or partitioning of material into distinct velocity fields to allow for displacement or velocity discontinuities at a material interface. A new method is presented in which a kernel-based damage field is constructed from the particle data. The gradient of this field is used to dynamically repartition the material into contact pairs at each node. This approach avoids the need to construct and evolve explicit cracks or contact surfaces and is therefore well suited to problems involving complex 3-D fracture with crack branching and coalescence. A straightforward extension of this approach permits frictional ‘self-contact’ between surfaces that are initially part of a single velocity field, enabling more accurate simulation of granular flow, porous compaction, fragmentation, and comminution of brittle materials. Numerical simulations of self contact and dynamic crack propagation are presented to demonstrate the accuracy of the approach. Copyright © 2016 John Wiley & Sons, Ltd.

We present three new sets of *C*^{1} hierarchical high-order tensor-product bases for conforming finite elements. The first basis is a high-order extension of the Bogner–Fox–Schmit basis. The edge and face functions are constructed using a combination of cubic Hermite and Jacobi polynomials with *C*^{1} global continuity on the common edges of elements. The second basis uses the tensor product of fifth-order Hermite polynomials and high-order functions and achieves global *C*^{1} continuity for meshes of quadrilaterals and *C*^{2} continuity on the element vertices. The third basis for triangles is also constructed using the tensor product of one-dimensional functions defined in barycentric coordinates. It also has global *C*^{1} continuity on edges and *C*^{2} continuity on vertices. A patch test is applied to the three considered elements. Projection and plate problems with smooth fabricated solutions are solved, and the performance of the *h*- and *p*-refinements are evaluated by comparing the approximation errors in the *L*_{2}- and energy norms. A plate with singularity is then studied, and *h*- and *p*-refinements are analysed. Finally, a transient problem with implicit time integration is considered. The results show exponential convergence rates with increasing polynomial order for the triangular and quadrilateral meshes of non-distorted and distorted elements. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, a simple locking-free triangular plate element, labeled here as Mindlin-type triangular plate element with nine degrees of freedom (MTP9), is presented. The element employs an incompatible approximation with nine degrees of freedom (DOFs) independent of the nodes and the shape of the triangle to define the displacements *u*/*v*/*w*(which is similar to a general solid element) along the *x*/*y*/*z* axes. It is free of shear locking, has a proper rank, and provides stable solutions for thick and thin plates. Moreover, the paper provides a new way to develop simple and efficient locking-free thick–thin-plate/shell elements. A variety of numerical examples demonstrate the convergence, accuracy, and robustness of the present element MTP9. Copyright © 2016 John Wiley & Sons, Ltd.

The paper is concerned with the modeling of the planar motion of a horizontal sheet of metal in a rolling mill. Inhomogeneous velocity profiles, with which the material is generated at one roll stand and enters the next one, lead to the time evolution of the deformation of the metal strip. We propose a novel kinematic description in which the axial coordinate is an Eulerian one, while the transverse motion of the sheet is modeled in a Lagrangian framework. The material volume travels across a finite element mesh, whose boundaries are in contact with the roll stands.

The concise mathematical formulation of the method is different from the classical form of the arbitrary Lagrangian–Eulerian approach with a rate form of constitutive relations.

The undeformed state of the strip is incompatible owing to the varying time rate of the generation of material. We treat this phenomenon by introducing the notion of intrinsic strains, which are mathematically described using the multiplicative decomposition of the deformation gradient.

We currently present the approach for quasistatic simulations of in-plane elastoplastic deformations of the strip. A practically relevant problem with two strip segments and three roll stands is studied in a numerical example. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, we develop a dual-horizon peridynamics (DH-PD) formulation that naturally includes varying horizon sizes and completely solves the ‘ghost force’ issue. Therefore, the concept of dual horizon is introduced to consider the unbalanced interactions between the particles with different horizon sizes. The present formulation fulfills both the balances of linear momentum and angular momentum exactly. Neither the ‘partial stress tensor’ nor the ‘slice’ technique is needed to ameliorate the ghost force issue. We will show that the traditional peridynamics can be derived as a special case of the present DH-PD. All three peridynamic formulations, namely, bond-based, ordinary state-based, and non-ordinary state-based peridynamics, can be implemented within the DH-PD framework. Our DH-PD formulation allows for *h*-adaptivity and can be implemented in any existing peridynamics code with minimal changes. A simple adaptive refinement procedure is proposed, reducing the computational cost. Both two-dimensional and three-dimensional examples including the Kalthoff–Winkler experiment and plate with branching cracks are tested to demonstrate the capability of the method. Copyright © 2016 John Wiley & Sons, Ltd.

The variational approach to fracture is effective for simulating the nucleation and propagation of complex crack patterns but is computationally demanding. The model is a strongly nonlinear non-convex variational inequality that demands the resolution of small length scales. The current standard algorithm for its solution, alternate minimization, is robust but converges slowly and demands the solution of large, ill-conditioned linear subproblems. In this paper, we propose several advances in the numerical solution of this model that improve its computational efficiency. We reformulate alternate minimization as a nonlinear Gauss–Seidel iteration and employ over-relaxation to accelerate its convergence; we compose this accelerated alternate minimization with Newton's method, to further reduce the time to solution, and we formulate efficient preconditioners for the solution of the linear subproblems arising in both alternate minimization and in Newton's method. We investigate the improvements in efficiency on several examples from the literature; the new solver is five to six times faster on a majority of the test cases considered. © 2016 The Authors International Journal for Numerical Methods in Engineering Published by 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 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.

Accurate sizing functions are crucial for efficient generation of high-quality meshes, but to define the sizing function is often the bottleneck in complicated mesh generation tasks because of the tedious user interaction involved. We present a novel algorithm to automatically create high-quality sizing functions for surface mesh generation. First, the tessellation of a Computer Aided Design (CAD) model is taken as the background mesh, in which an initial sizing function is defined by considering geometrical factors and user-specified parameters. Then, a convex nonlinear programming problem is formulated and solved efficiently to obtain a smoothed sizing function that corresponds to a mesh satisfying necessary gradient constraint conditions and containing a significantly reduced element number. Finally, this sizing function is applied in an advancing front mesher. With the aid of a walk-through algorithm, an efficient sizing-value query scheme is developed. Meshing experiments of some very complicated geometry models are presented to demonstrate that the proposed sizing-function approach enables accurate and fully automatic surface mesh generation. Copyright © 2016 John Wiley & Sons, Ltd.

This contribution presents a novel approach to structural shape optimization that relies on an embedding domain discretization technique. The evolving shape design is embedded within a uniform finite element background mesh which is then used for the solution of the physical state problem throughout the course of the optimization. We consider a boundary tracking procedure based on adaptive mesh refinement to separate between interior elements, exterior elements, and elements intersected by the physical domain boundary. A selective domain integration procedure is employed to account for the geometric mismatch between the uniform embedding domain discretization and the evolving structural component. Thereby, we avoid the need to provide a finite element mesh that conforms to the structural component for every design iteration, as it is the case for a standard Lagrangian approach to structural shape optimization. Still, we adopt an explicit shape parametrization that allows for a direct manipulation of boundary vertices for the design evolution process. In order to avoid irregular and impracticable design updates, we consider a geometric regularization technique to render feasible descent directions for the course of the optimization. Copyright © 2016 John Wiley & Sons, Ltd.

The paper deals with two main advantages in the analysis of slender elastic structures both achieved through the mixed (stress and displacement) format with respect to the more commonly used displacement one: (i) the smaller error in the extrapolations usually employed in the solution strategies of nonlinear problems and (ii) the lower polynomial dependence of the problem equations on the finite element degrees of freedom when solid finite elements are used. The smaller extrapolation error produces a lower number of iterations and larger step length in path-following analysis and a greater accuracy in Koiter asymptotic method. To focus on the origin of the phenomenon, the two formats are derived for the same finite element interpolation. The reduced polynomial dependence improves the Koiter asymptotic strategy in terms of both computational efficiency, accuracy and simplicity. 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 shear-flexible isogeometric Reissner–Mindlin shell formulation using non-uniform rational B-splines basis functions is introduced, which is used for the demonstration of a coupling approach for multiple non-conforming patches. The six degrees of freedom formulation uses the exact surface normal vectors and curvature. The shell formulation is implemented in an isogeometric analysis framework for computation of structures composed of multiple geometric entities. To enable local model refinement as well as non-matching domains coupling, a conservative multi-patch approach using Lagrange multipliers for structured non-uniform rational B-splines patches is presented. Here, an additional border frame mesh is used to couple geometries during structural analyses. This frame interface approach avoids the problem of excessive constraints when multiple patches are coupled at one point. First, the shell formulation is verified with several reference cases. Then the influence of the frame interface discretization and frame penalty stiffness on the smoothness of the results is investigated. The effects of the perturbed Lagrangian method in combination with the frame interface approach is shown. In addition, results of models with T-joint interface connections and perpendicular stiffener patches are presented. Copyright © 2016 John Wiley & Sons, Ltd.

A non-gradient-based approach for topology optimization using a genetic algorithm is proposed in this paper. The genetic algorithm used in this paper is assisted by the Kriging surrogate model to reduce computational cost required for function evaluation. To validate the non-gradient-based topology optimization method in flow problems, this research focuses on two single-objective optimization problems, where the objective functions are to minimize pressure loss and to maximize heat transfer of flow channels, and one multi-objective optimization problem, which combines earlier two single-objective optimization problems. The shape of flow channels is represented by the level set function. The pressure loss and the heat transfer performance of the channels are evaluated by the Building-Cube Method code, which is a Cartesian-mesh CFD solver. The proposed method resulted in an agreement with previous study in the single-objective problems in its topology and achieved global exploration of non-dominated solutions in the multi-objective problems. © 2016 The Authors International Journal for Numerical Methods in Engineering Published by John Wiley & Sons Ltd

Reproducing kernel particle method (RKPM) has been applied to many large deformation problems. RKPM relies on polynomial reproducing conditions to yield desired accuracy and convergence properties but requires appropriate kernel support coverage of neighboring nodes to ensure kernel stability. This kernel stability condition is difficult to achieve for problems with large particle motion such as the fragment-impact processes that commonly exist in extreme events. A new reproducing kernel formulation with ‘quasi-linear’ reproducing conditions is introduced. In this approach, the first-order polynomial reproducing conditions are approximately enforced to yield a nonsingular moment matrix. With proper error control of the first-order completeness, nearly second-order convergence rate in L2 norm can be achieved while maintaining kernel stability. The effectiveness of this quasi-linear RKPM formulation is demonstrated by modeling several extremely large deformation and fragment-impact problems. Copyright © 2016 John Wiley & Sons, Ltd.

We propose a formulation for linear elastic fracture mechanics in which the stress intensity factors are found directly from the solution vector of an extended boundary element method formulation. The enrichment is embedded in the boundary element method formulation, rather than adding new degrees of freedom for each enriched node. Therefore, a very limited number of new degrees of freedom is added to the problem, which contributes to preserving the conditioning of the linear system of equations. The Stroh formalism is used to provide boundary element method fundamental solutions for any degree of anisotropy, and these are used for both conventional and enriched degrees of freedom. Several numerical examples are shown with benchmark solutions to validate the proposed method. Copyright © 2016 John Wiley & Sons, Ltd.

The fractional step method (FSM) is an efficient solution technique for the particle finite element method, a Lagrangian-based approach to simulate fluid–structure interaction (FSI). Despite various refinements, the applicability of the FSM has been limited to low viscosity flow and FSI simulations with a small number of equations along the fluid–structure interface. To overcome these limitations, while incorporating nonlinear response in the structural domain, an FSM that unifies structural and fluid response in the discrete governing equations is developed using the quasi-incompressible formulation. With this approach, fluid and structural particles do not need to be treated separately, and both domains are unified in the same system of equations. Thus, the equations along the fluid–structure interface do not need to be segregated from the fluid and structural domains. Numerical examples compare the unified FSM with the non-unified FSM and show that the computational cost of the proposed method overcomes the slow convergence of the non-unified FSM for high values of viscosity. As opposed to the non-unified FSM, the number of iterations required for convergence with the unified FSM becomes independent of viscosity and time step, and the simulation run time does not depend on the size of the FSI interface. Copyright © 2016 John Wiley & Sons, Ltd.

Reduced order models are useful for accelerating simulations in many-query contexts, such as optimization, uncertainty quantification, and sensitivity analysis. However, offline training of reduced order models (ROMs) can have prohibitively expensive memory and floating-point operation costs in high-performance computing applications, where memory per core is limited. To overcome this limitation for proper orthogonal decomposition, we propose a novel adaptive selection method for snapshots in time that limits offline training costs by selecting snapshots according an error control mechanism similar to that found in adaptive time-stepping ordinary differential equation solvers. The error estimator used in this work is related to theory bounding the approximation error in time of proper orthogonal decomposition-based ROMs, and memory usage is minimized by computing the singular value decomposition using a single-pass incremental algorithm. Results for a viscous Burgers' test problem demonstrate convergence in the limit as the algorithm error tolerances go to zero; in this limit, the full-order model is recovered to within discretization error. A parallel version of the resulting method can be used on supercomputers to generate proper orthogonal decomposition-based ROMs, or as a subroutine within hyperreduction algorithms that require taking snapshots in time, or within greedy algorithms for sampling parameter space. 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 proposes an efficient metamodeling approach for uncertainty quantification of complex system based on Gaussian process model (GPM). The proposed GPM-based method is able to efficiently and accurately calculate the mean and variance of model outputs with uncertain parameters specified by arbitrary probability distributions. Because of the use of GPM, the closed form expressions of mean and variance can be derived by decomposing high-dimensional integrals into one-dimensional integrals. This paper details on how to efficiently compute the one-dimensional integrals. When the parameters are either uniformly or normally distributed, the one-dimensional integrals can be analytically evaluated, while when parameters do not follow normal or uniform distributions, this paper adopts the effective Gaussian quadrature technique for the fast computation of the one-dimensional integrals. As a result, the developed GPM method is able to calculate mean and variance of model outputs in an efficient manner independent of parameter distributions. The proposed GPM method is applied to a collection of examples. And its accuracy and efficiency is compared with Monte Carlo simulation, which is used as benchmark solution. Results show that the proposed GPM method is feasible and reliable for efficient uncertainty quantification of complex systems in terms of the computational accuracy and efficiency. Copyright © 2016 John Wiley & Sons, Ltd.

A nonparametric probabilistic approach for modeling uncertainties in projection-based, nonlinear, reduced-order models is presented. When experimental data are available, this approach can also quantify uncertainties in the associated high-dimensional models. The main underlying idea is twofold. First, to substitute the deterministic reduced-order basis (ROB) with a stochastic counterpart. Second, to construct the probability measure of the stochastic reduced-order basis (SROB) on a subset of a compact Stiefel manifold in order to preserve some important properties of a ROB. The stochastic modeling is performed so that the probability distribution of the constructed SROB depends on a small number of hyperparameters. These are determined by solving a reduced-order statistical inverse problem. The mathematical properties of this novel approach for quantifying model uncertainties are analyzed through theoretical developments and numerical simulations. Its potential is demonstrated through several example problems from computational structural dynamics. Copyright © 2016 John Wiley & Sons, Ltd.

By separation of scales and the homogenization of a flow through porous media, a two-scale problem arises where a Darcy-type flow is present on the macroscale and a Stokes flow on the subscale. In this paper, the problem is given as the minimization of a potential. Additional constraints imposing periodicity in a weak sense are added using Lagrange multipliers. In particular, the upper and lower energy bounds for the corresponding strongly periodic problem are produced, quantifying the accuracy of the weakly periodic boundary conditions. A numerical example demonstrates the evaluation of energy bounds and the performance of weakly periodic boundary conditions on a representative volume element. © 2016 The Authors. International Journal for Numerical Methods in Engineering Published by John Wiley & Sons Ltd

In this paper, a four-node quadrilateral flat shell element is proposed for geometrically nonlinear analysis based on updated Lagrangian formulation with the co-rotational kinematics concept. The flat shell element combines the membrane element with drilling degrees of freedom and the plate element with shear deformation. By means of these linearized elements, a simplified nonlinear analysis procedure allowing for warping of the flat shell element and large rotation is proposed. The tangent stiffness matrix and the internal force recovery are formulated in this paper. Several classic benchmark examples are presented to validate the accuracy and efficiency of the proposed new and more proficient element for practical engineering analysis of shell structures. Copyright © 2016 John Wiley & Sons, Ltd.

It is common in solving topology optimization problems to replace an integer-valued characteristic function design field with the material volume fraction field, a real-valued approximation of the design field that permits ‘fictitious’ mixtures of materials during intermediate iterations in the optimization process. This is reasonable so long as one can interpolate properties for such materials and so long as the final design is integer valued. For this purpose, we present a method for smoothly thresholding the volume fractions of an arbitrary number of material phases which specify the design. This method is trivial for two-material design problems, for example, the canonical topology design problem of specifying the presence or absence of a single material within a domain, but it becomes more complex when three or more materials are used, as often occurs in material design problems. We take advantage of the similarity in properties between the volume fractions and the barycentric coordinates on a simplex to derive a thresholding, method which is applicable to an arbitrary number of materials. As we show in a sensitivity analysis, this method has smooth derivatives, allowing it to be used in gradient-based optimization algorithms. We present results, which show synergistic effects when used with Solid Isotropic Material with Penalty and Rational Approximation of Material Properties material interpolation functions, popular methods of ensuring integerness of solutions. Copyright © 2016 John Wiley & Sons, Ltd.

This work discusses a discontinuous Galerkin (DG) discretization for two-phase flows. The fluid interface is represented by a level set, and the DG approximation space is adapted such that jumps and kinks in pressure and velocity fields can be approximated sharply. This adaption of the original DG space, which can be performed ‘on-the-fly’ for arbitrary interface shapes, is referred to as extended discontinuous Galerkin. By combining this ansatz with a special quadrature technique, one can regain spectral convergence properties for low-regularity solutions, which is demonstrated by numerical examples. This work focuses on the aspects of spatial discretization, and special emphasis is devoted on how to overcome problems related to quadrature, small cut cells, and condition number of linear systems. Temporal discretization will be discussed in future works. Copyright © 2016 John Wiley & Sons, Ltd.

Digital imaging technologies such as X-ray scans and ultrasound provide a convenient and non-invasive way to capture high-resolution images. The colour intensity of digital images provides information on the geometrical features and material distribution which can be utilised for stress analysis. The proposed approach employs an automatic and robust algorithm to generate quadtree (2D) or octree (3D) meshes from digital images. The use of polygonal elements (2D) or polyhedral elements (3D) constructed by the scaled boundary finite element method avoids the issue of hanging nodes (mesh incompatibility) commonly encountered by finite elements on quadtree or octree meshes. The computational effort is reduced by considering the small number of cell patterns occurring in a quadtree or an octree mesh. Examples with analytical solutions in 2D and 3D are provided to show the validity of the approach. Other examples including the analysis of 2D and 3D microstructures of concrete specimens as well as of a domain containing multiple spherical holes are presented to demonstrate the versatility and the simplicity of the proposed technique. Copyright © 2016 John Wiley & Sons, Ltd.

A new isogeometric density-based approach for the topology optimization of multi-material structures is presented. In this method, the density fields of multiple material phases are represented using the isogeometric non-uniform rational B-spline-based parameterization leading to exact modeling of the geometry, removing numerical artifacts and full analytical computation of sensitivities in a cost-effective manner. An extension of the perimeter control technique is introduced where restrictions are imposed on the perimeters of density fields of all phases. Consequently, not only can one control the complexity of the optimal design but also the minimal lengths scales of all material phases. This leads to optimal designs with significantly enhanced manufacturability and comparable performance. Unlike the common element-wise or nodal-based density representations, owing to higher order continuity of density fields in this method, their gradients required for perimeter control restrictions are calculated exactly without additional computational cost. The problem is formulated with constraints on either (1) volume fractions of different material phases or (2) the total mass of the structure. The proposed method is applied for the minimal compliance design of two-dimensional structures consisting of multiple distinct materials as well as functionally graded ones. Copyright © 2016 John Wiley & Sons, Ltd.

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.

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

No abstract is available for this article.

]]>No abstract is available for this article.

]]>Inhomogeneous flows involving dense particulate media display clear size effects, in which the particle length scale has an important effect on flow fields. Hence, nonlocal constitutive relations must be used in order to predict these flows. Recently, a class of nonlocal fluidity models has been developed for emulsions and subsequently adapted to granular materials. These models have successfully provided a quantitative description of experimental flows in many different flow configurations. In this work, we present a finite element-based numerical approach for solving the nonlocal constitutive equations for granular materials, which involve an additional, non-standard nodal degree-of-freedom – the granular fluidity, which is a scalar state parameter describing the susceptibility of a granular element to flow. Our implementation is applied to three canonical inhomogeneous flow configurations: (1) linear shear with gravity, (2) annular shear flow without gravity, and (3) annular shear flow with gravity. We verify our implementation, demonstrate convergence, and show that our results are mesh independent. Copyright © 2016 John Wiley & Sons, Ltd.

We construct new robust and efficient preconditioned generalized minimal residual solvers for the monolithic linear systems of algebraic equations arising from the finite element discretization and Newton's linearization of the fully coupled fluid–structure interaction system of partial differential equations in the arbitrary Lagrangian–Eulerian formulation. We admit both linear elastic and nonlinear hyperelastic materials in the solid model and cover a large range of flows, for example, water, blood, and air, with highly varying density. The preconditioner is constructed in form of
, where
,
, and
are proper approximations to the matrices *L*, *D*, and *U* in the *LDU* block factorization of the fully coupled system matrix, respectively. The inverse of the corresponding Schur complement is approximated by applying a few cycles of a special class of algebraic multigrid methods to the perturbed fluid sub-problem, which is obtained by modifying corresponding entries in the original fluid matrix with an explicitly constructed approximation to the exact perturbation coming from the sparse matrix–matrix multiplications. The numerical studies presented impressively demonstrate the robustness and the efficiency of the preconditioner proposed in the paper. Copyright © 2016 John Wiley & Sons, Ltd.

Stress-related problems have not been given the same attention as the minimum compliance topological optimization problem in the literature. Continuum structural topological optimization with stress constraints is of wide engineering application prospect, in which there still are many problems to solve, such as the stress concentration, an equivalent approximate optimization model and etc. A new and effective topological optimization method of continuum structures with the stress constraints and the objective function being the structural volume has been presented in this paper. To solve the stress concentration issue, an approximate stress gradient evaluation for any element is introduced, and a total aggregation normalized stress gradient constraint is constructed for the optimized structure under the *r*−th load case. To obtain stable convergent series solutions and enhance the control on the stress level, two p-norm global stress constraint functions with different indexes are adopted, and some weighting p-norm global stress constraint functions are introduced for any load case. And an equivalent topological optimization model with reduced stress constraints is constructed,being incorporated with the rational approximation for material properties, an active constraint technique, a trust region scheme, and an effective local stress approach like the qp approach to resolve the stress singularity phenomenon. Hence, a set of stress quadratic explicit approximations are constructed, based on stress sensitivities and the method of moving asymptotes. A set of algorithm for the one level optimization problem with artificial variables and many possible non-active design variables is proposed by adopting an inequality constrained nonlinear programming method with simple trust regions, based on the primal-dual theory, in which the non-smooth expressions of the design variable solutions are reformulated as smoothing functions of the Lagrange multipliers by using a novel smoothing function. Finally, a two-level optimization design scheme with active constraint technique, i.e. varied constraint limits, is proposed to deal with the aggregation constraints that always are of loose constraint (non active constraint) features in the conventional structural optimization method. A novel structural topological optimization method with stress constraints and its algorithm are formed, and examples are provided to demonstrate that the proposed method is feasible and very effective. © 2016 The Authors. International Journal for Numerical Methods in Engineering published by John Wiley & Sons Ltd.