In this work, a hybrid numerical approach, which combines elements of SPH and coarse-grained molecular dynamics, is used to investigate the effect of various flow conditions to deformable and breakable shell-structures such as capsules, vesicles or biological cells. Copyright © 2014 John Wiley & Sons, Ltd.

A computational contact homogenization framework is established for the modeling and simulation of soft matter friction. The main challenges toward the realization of the framework are (1) the establishment of a frictional contact algorithm that displays an optimal combination of accuracy, efficiency, and robustness and plays a central role in (2) the construction of a micromechanical contact test within which samples of arbitrary size may be embedded and which is not restricted to a single deformable body. The former challenge is addressed through the extension of mixed variational formulations of contact mechanics to a mortar-based isogeometric setting where the augmented Lagrangian approach serves as the constraint enforcement method. The latter challenge is addressed through the concept of periodic embedding, with which a periodically replicated -continuous interface topography is realized across which not only pending but also ensuing contact among simulation cells will be automatically captured. Two-dimensional and three-dimensional investigations with unilateral/bilateral periodic/random roughness on two elastic micromechanical samples demonstrate the overall framework and the nature of the macroscopic frictional response. Copyright © 2014 John Wiley & Sons, Ltd.

Most engineering applications involving solutions by numerical methods are dependent on several parameters, whose impact on the solution may significantly vary from one to the other. At times an evaluation of these multivariate solutions may be required at the expense of a prohibitively high computational cost. In the present paper, an adaptive approach is proposed as a way to estimate the solution of such multivariate finite element problems. It is based upon the integration of so-called nested Padé approximants within the finite element procedure. This procedure includes an effective control of the approximation error, which enables adaptive refinements of the converged intervals upon reconstruction of the solution. The main advantages lie in a potential reduction of the computational effort and the fact that the level of *a priori* knowledge required about the solution in order to have an accurate, efficient, and well-sampled estimate of the solution is low. The approach is introduced for bivariate problems, for which it is validated on elasto-poro-acoustic problems of both academic and more industrial scale. It is argued that the methodology in general holds for more than two variables, and a discussion is opened about the truncation refinements required in order to generalize the results accordingly. Copyright © 2014 John Wiley & Sons, Ltd.

The present work addresses the lower bound limit analysis (or yield design) of thick plates under shear-bending interaction. Equilibrium finite elements are used to discretize the bending moment and the shear force fields. Different strength criteria, formulated in the five-dimensional space of bending moment and shear force, are considered, one of them taking into account the interaction between bending and shear resistances. The criteria are chosen to be sufficiently simple so that the resulting optimization problem can be formulated as a second-order cone programming problem (SOCP), which is solved by the dedicated solver MOSEK. The efficiency of the proposed finite element is illustrated by means of numerical examples on different plate geometries, for which the thin plate solutions as well as the pure shear solutions are accurately obtained as two different limit cases of the plate slenderness ratio. In particular, the proposed element exhibits a good behavior in the thin plate limit. Copyright © 2014 John Wiley & Sons, Ltd.

First-order reliability method (FORM) has been mostly utilized for solving reliability-based design optimization (RBDO) problems efficiently. However, second-order reliability method (SORM) is required in order to estimate a probability of failure accurately in highly nonlinear performance functions. Despite accuracy of SORM, its application to RBDO is quite challenging due to unaffordable numerical burden incurred by a Hessian calculation. For reducing the numerical efforts, a quasi-Newton approach to approximate the Hessian is introduced in this study instead of calculating the true Hessian. The proposed SORM with the approximated Hessian requires computations only used in FORM, leading to very efficient and accurate reliability analysis. The proposed SORM also utilizes a generalized chi-squared distribution in order to achieve better accuracy. Furthermore, SORM-based inverse reliability method is proposed in this study. An accurate reliability index corresponding to a target probability of failure is updated using the proposed SORM. Two approaches in terms of finding an accurate most probable point using the updated reliability index are proposed. The proposed SORM-based inverse analysis is then extended to RBDO in order to obtain a reliability-based optimum design satisfying probabilistic constraints more accurately even for a highly nonlinear system. The numerical study results show that the proposed reliability analysis and RBDO achieve efficiency of FORM and accuracy of SORM at the same time. Copyright © 2014 John Wiley & Sons, Ltd.

This contribution presents a mesh adaptive crack propagation scheme for the evaluation of the viscoelastic fracture response of elastomers at large strains and up to high loading rates. The approach accounts for micromechanical based features of both elastic and viscoelastic bulk responses of idealized polymer networks. To this end, the Bergstörm–Boyce model is considered to introduce hyperelastic and nonlinear finite viscoelastic responses. Moreover, the crack driving force and the crack driving direction are predicted by the material force approach. A consistent thermodynamic framework for the combined configurational motion in viscoelastic continua at finite strain regime is discussed.

The fracture toughness of non-strain-crystallizing elastomers shows strong rate dependency and the energy release rate versus the rate of tearing to be a fundamental material property. Therefore, in this contribution, a dynamic fracture criterion, which is a function of the rate of crack growth, is shown to be adequate in numerical simulations. The use of the presented method enables to study fracture behaviour of any material nonlinearity within the implicit time integration. Main feature of the proposed algorithm is restructuring the overall discrete system by duplication of crack front DOFs based on minimization of the overall energy via the Griffith criterion. Copyright © 2014 John Wiley & Sons, Ltd.

We compare the computational efficiency of isogeometric Galerkin and collocation methods for partial differential equations in the asymptotic regime. We define a metric to identify when numerical experiments have reached this regime. We then apply these ideas to analyze the performance of different isogeometric discretizations, which encompass *C*^{0} finite element spaces and higher-continuous spaces. We derive convergence and cost estimates in terms of the total number of degrees of freedom and then perform an asymptotic numerical comparison of the efficiency of these methods applied to an elliptic problem. These estimates are derived assuming that the underlying solution is smooth, the full Gauss quadrature is used in each non-zero knot span and the numerical solution of the discrete system is found using a direct multi-frontal solver. We conclude that under the assumptions detailed in this paper, higher-continuous basis functions provide marginal benefits. Copyright © 2014 John Wiley & Sons, Ltd.

It is important to design robust and reliable systems by accounting for uncertainty and variability in the design process. However, performing optimization in this setting can be computationally expensive, requiring many evaluations of the numerical model to compute statistics of the system performance at every optimization iteration. This paper proposes a multifidelity approach to optimization under uncertainty that makes use of inexpensive, low-fidelity models to provide approximate information about the expensive, high-fidelity model. The multifidelity estimator is developed based on the control variate method to reduce the computational cost of achieving a specified mean square error in the statistic estimate. The method optimally allocates the computational load between the two models based on their relative evaluation cost and the strength of the correlation between them. This paper also develops an information reuse estimator that exploits the autocorrelation structure of the high-fidelity model in the design space to reduce the cost of repeatedly estimating statistics during the course of optimization. Finally, a combined estimator incorporates the features of both the multifidelity estimator and the information reuse estimator. The methods demonstrate 90% computational savings in an acoustic horn robust optimization example and practical design turnaround time in a robust wing optimization problem. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents a novel mathematical programming approach for the static stability analysis of structures with uncertainties within the framework of FEM. The considered uncertain parameters are material properties, geometry of element cross section, and loading conditions, all of which are described by an interval model. The proposed method formulates the two cases of interest, namely, worst and best buckling load calculation, into a pair of mathematical programming problems. Two straightforward advantages are exhibited by such formulations. The first advantage is that the proposed formulation can overcome the interference on the sharpness of bounds of the buckling load due to the interval dependence issue. The second benefit is that the information of uncertain parameters causing the extremities of buckling load can always be retrieved as by-products of the uncertain stability analysis. Some numerical examples are presented to illustrate the capability of the proposed method on various structures and the sharpness of the bounds of the buckling load factors. The efficiency and effectiveness of the proposed method are also demonstrated through comparison with the classical Monte Carlo simulation method. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents a numerical scheme based on a fictitious domain framework for the numerical modeling of earthquake-induced ground motion in the presence of realistic surface topography of the Earth's crust. We show that by adopting a non-conforming octree-based meshing approach associated with a virtual representation of the surficial irregularity, we can obtain accurate representations of ground motion. From the computational point of view, our methodology proves to be also efficient, and more importantly, it allows us to preserve the salient features of multi-resolution cubic-shaped finite elements for wave propagation applications. We implemented the non-conforming meshing scheme for the treatment of realistic topographies into Hercules, the octree-based finite-element earthquake simulator developed by the Quake Group at Carnegie Mellon University. We tested the benefits of the strategy by benchmarking its results against reference examples, and by means of numerical error estimate analyses. Our qualitative and quantitative comparisons showed a close agreement between our numerical results and the reference results, and also, that the order of convergence of the displacement field is preserved in the presence of surface topography. Moreover, this performance was obtained by using the same mesh refinement techniques with cubic elements as in traditional flat free-surface simulations. Copyright © 2014 John Wiley & Sons, Ltd.

High Order Numerical Manifold Method (HONMM) is a powerful method to solve static problems. A development of HONMM to achieve a dynamic solution with high accuracy and less computational cost is addressed in the current paper. In the developed method, the global approximation is obtained through increasing the order of local approximation functions without any Linear Dependence (LD) of the unknowns. The weighted residual formulations are modified to be used in dynamic high order simulation. Moreover, a modified Newmark method formulation is adjusted for time integration of high order equations. The superiority of the proposed method over the conventional NMM is demonstrated through a special beam example. The dynamic free fall block example is used to exhibit the removal of mass matrix singularity. As cases of dynamic analysis, beam free and forced vibrations are illustrated which include a moving load. Finally, a non-uniform cross-section beam under dynamic variable loads with accelerated motion is solved while demonstrating the capability of the new method such as simplicity, accuracy and time efficiency for simulation of complex dynamic problems. Copyright © 2014 John Wiley & Sons, Ltd.

Micro-mechanical and macro-mechanical behavior of face-centered cubic (FCC) crystals is investigated by using different forms of strain energy functions in hyperelastic material models in crystal plasticity finite element framework. A quadratic strain energy function with anisotropic elastic constants, a polyconvex strain energy function with invariants associated with the cubic symmetry, and a strain energy function from an inter-atomic potential are considered in hyperelastic material models to describe the elastic deformation of FCC crystals. In our numerical experiments, the trajectories of {111} poles in the pole figure and the accumulated plastic slips of FCC coppers under uniaxial tension and simple shear depend on the choice of strain energy functions when the slip resistance of the slip systems is high. The ability of strain energy functions in this study to represent elastic lattice distortions in crystals varies with the amount of elastic deformation and the shape of deformed lattice. However, numerical results show that the change of macroscopic mechanical behavior of FCC coppers is not significant for the choice of strain energy functions, compared with the change of crystallographic texture evolution. Copyright © 2014 John Wiley & Sons, Ltd.

A numerical technique that is based on the integration of the asymptotic solution in the numerical framework for computing the local singular behavior of Stokes flow near a sharp corner is presented. Moffat's asymptotic solution is used, and special enriched shape functions are developed and integrated in the extended finite element method (X-FEM) framework to solve the Navier–Stokes equations. The no-slip boundary condition on the walls of the corner is enforced via the use of Lagrange multipliers. Flows around corners with different angles are simulated, and the results are compared with both those of the known analytic solution and the X-FEM with no special enrichment near the corner. The results of the present technique are shown to greatly reduce the error made in computing the pressure and velocity fields near a corner tip without the need for mesh refinement near the corner. The method is then applied to the estimation of the permeability of a network of fibers, where it is shown that the local small-scale pressure singularities have a large impact on the large-scale network permeability. Copyright © 2014 John Wiley & Sons, Ltd.

A new Molecular Dynamics Finite Element Method (MDFEM) with a coupled mechanical-charge/dipole formulation is proposed. The equilibrium equations of Molecular Dynamics (MD) are embedded exactly within the computationally more favourable Finite Element Method (FEM). This MDFEM can readily implement any force field because the constitutive relations are explicitly uncoupled from the corresponding geometric element topologies. This formal uncoupling allows to differentiate between chemical-constitutive, geometric and mixed-mode instabilities. Different force fields, including bond-order reactive and polarisable fluctuating charge–dipole potentials, are implemented exactly in both explicit and implicit dynamic commercial finite element code. The implicit formulation allows for larger length and time scales and more varied eigenvalue-based solution strategies.

The proposed multi-physics and multi-scale compatible MDFEM is shown to be equivalent to MD, as demonstrated by examples of fracture in carbon nanotubes (CNT), and electric charge distribution in graphene, but at a considerably reduced computational cost. The proposed MDFEM is shown to scale linearly, with concurrent continuum FEM multi-scale couplings allowing for further computational savings. Moreover, novel conformational analyses of pillared graphene structures (PGS) are produced. The proposed model finds potential applications in the parametric topology and numerical design studies of nano-structures for desired electro-mechanical properties (e.g. stiffness, toughness and electric field induced vibrational/electron-emission properties). Copyright © 2014 John Wiley & Sons, Ltd.

We formulate a class of delicately controlled problems to model the kink-free evolution of quasistatic cracks in brittle, isotropic, linearly elastic materials in two dimensions. The evolving crack satisfies familiar principles—Griffith's criterion, local symmetry, and irreversibility. A novel feature of the formulation is that in addition to the crack path, the loading is also treated as an unknown. Specifically, a scaling factor for prescribed Dirichlet and Neumann boundary conditions is computed as part of the solution to yield an always-propagating and ostensibly kink-free crack and a continuous loading history beyond the initial step. A dimensionless statement of the problem depends only on the Poisson's ratio of a homogeneous material, and is in particular, independent of its Young's modulus and fracture toughness.

Numerical resolution of the formulated problem relies on two new ideas. The first is an algorithm to compute triangulations conforming to cracked domains by locally deforming a *given* background mesh in the vicinity of evolving cracks. The algorithm is robust under mild assumptions on the sizes and angles of triangles near the crack and its smoothness. Hence, a subset, if not the entire family, of cracked domains realized during the course of a simulation can be discretized with the same background mesh; we term the latter a *universal mesh* for such a family of domains. Universal meshes facilitate adopting a discrete representation for the crack (as splines in our examples), preclude the need for local splitting/retriangulation operations, liberate the crack from following directions prescribed by the mesh, and enable the adoption of standard FEMs to compute the elastic fields in the cracked solid. Second, we employ a method specifically designed to approximate the stress intensity factors for curvilinear cracks. We examine the performance of the resulting numerical method with detailed examples, including comparisons with an exact solution of a crack propagating along a circular arc (which we construct) and comparisons with experimental fracture paths. In all cases, we observe convergence of computed paths, their derivatives, and loading histories with refinement of the universal mesh. Copyright © 2014 John Wiley & Sons, Ltd.

We develop a three-dimensional, hierarchically parallel, finite strain multiscale solver capable of computing nonlinear multiscale solutions with over 1 billion finite elements and over 574 million nonlinear equations on 1552 computing cores. In the vein of FE^{2}, we use the nested iterative procedure and devote the solver to multiscale cohesive modeling of heterogeneous hyperelastic layers. The hierarchically parallel multiscale solver takes advantage of a client-server non-blocking communication matrix that limits latency, starvation, and overhead by overlaying computations at different scales. We perform simulations of real-scale engineering devices and bridge in length-scales, spanning from mm to nm in spatial resolution. Verification of the hierarchically parallel solver is provided together with a mesh convergence study. Moreover, we report on the scaling performance. Copyright © 2014 John Wiley & Sons, Ltd.

We present a parameterized-background data-weak (PBDW) formulation of the variational data assimilation (state estimation) problem for systems modeled by partial differential equations. The main contributions are a constrained optimization weak framework informed by the notion of experimentally observable spaces; *a priori* and *a posteriori* error estimates for the field and associated linear-functional outputs; weak greedy construction of prior (background) spaces associated with an underlying potentially high-dimensional parametric manifold; stability-informed choice of observation functionals and related sensor locations; and finally, output prediction from the optimality saddle in operations, where *M* is the number of experimental observations. We present results for a synthetic Helmholtz acoustics model problem to illustrate the elements of the methodology and confirm the numerical properties suggested by the theory. To conclude, we consider a physical raised-box acoustic resonator chamber: we integrate the PBDW methodology and a Robotic Observation Platform to achieve real-time *in situ* state estimation of the time-harmonic pressure field; we demonstrate the considerable improvement in prediction provided by the integration of a best-knowledge model and experimental observations; we extract, even from these results with real data, the numerical trends indicated by the theoretical convergence and stability analyses. Copyright © 2014 John Wiley & Sons, Ltd.

Reduced-order models for linear time-invariant dynamical systems are considered, and the error between the full-order model and the reduced-order model solutions is characterized. Based on the analytical representation of the error, an *a posteriori* error indicator is proposed that combines a Krylov-based exponential integrator and an *a posteriori* residual-based estimate. Numerical experiments illustrate the quality of the error estimator. Copyright © 2014 John Wiley & Sons, Ltd.

One of the major challenges in the Bayesian solution of inverse problems governed by partial differential equations (PDEs) is the computational cost of repeatedly evaluating numerical PDE models, as required by Markov chain Monte Carlo (MCMC) methods for posterior sampling. This paper proposes a data-driven projection-based model reduction technique to reduce this computational cost. The proposed technique has two distinctive features. First, the model reduction strategy is tailored to inverse problems: the snapshots used to construct the reduced-order model are computed adaptively from the posterior distribution. Posterior exploration and model reduction are thus pursued simultaneously. Second, to avoid repeated evaluations of the full-scale numerical model as in a standard MCMC method, we couple the full-scale model and the reduced-order model together in the MCMC algorithm. This maintains accurate inference while reducing its overall computational cost. In numerical experiments considering steady-state flow in a porous medium, the data-driven reduced-order model achieves better accuracy than a reduced-order model constructed using the classical approach. It also improves posterior sampling efficiency by several orders of magnitude compared with a standard MCMC method. Copyright © 2014 John Wiley & Sons, Ltd.

The volumetric locking issue is critical in finite element analysis of nearly incompressible problems. In this paper, a new five-node quadrilateral element (Q5) is proposed by enriching the four-node quadrilateral element (Q4) with a centroid node to solve the volumetric locking problem in FEM. The cell-based smoothed FEM is employed with Q5 element (referred as Q5-CS-SC4) to soften the stiffness in order to obtain a better solution. To eliminate pressure oscillation in near-incompressible problems, an edge-based area-weighted smoothing scheme incorporated with the cell-wise divergence-free Q5 element is carried out (referred as Q5-EAW), and an adjustable area-weighted strain smoothing scheme using a parameter *p* is proposed to improve the performance of the Q5 element in dealing with incompressible media (referred as Q5-*p*EAW). The formulation of Q5-*p*EAW is a simple combination of Q5-CS-SC4 and Q5-EAW by an adjustable area weight. It can search the exact strain energy of the problem in near-incompressible cases. We also introduce another node-based strain smoothing technique (Q5-NAW) into the domain-based selective scheme to obtain a Q5-*p*EAW/NAW model to solve the pressure oscillation, which gives a much smoother pressure solution than Q5-EAW. Finally, the Q5 element is extended into hexahedral element to develop a nine-node hexahedral (H9) element shape function by enriching the eight-node hexahedral element (Q8) with a centroid node. An H9-Gi/NAW model similar to Q5-*p*EAW/NAW is proposed by using H9 element to solve the 3D volumetric locking. A series of benchmark problems are provided to demonstrate that the proposed Q5-*p*EAW/NAW for 2D plane strain problems, and H9-Gi/NAW model for 3D cases are locking-free for nearly incompressible problems. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents a bubble-enhanced smoothed finite element formulation for the analysis of volume-constrained problems in two-dimensional linear elasticity. The new formulation is derived based on the variational multi-scale approach in which unequal order displacement-pressure pairs are used for the mixed finite element approximation and hierarchical bubble function is selected for the fine-scale displacement approximation. An area-weighted averaging scheme is employed for the two-scale smoothed strain calculation under the framework of edge-based smoothed FEM. The smoothed fine-scale solution is shown to naturally contain the stress field jump of the smoothed coarse-scale solution across the boundary of edge-based smoothing domain and thus provides the possibility to stabilize the global solution for volume-constrained problems. A global monolithic solution strategy is employed, and the fine-scale solution is solved without the consideration of approximating the strong form of the fine-scale equation. Several numerical examples are analyzed to demonstrate the accuracy of the present formulation. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents new developments on a weakly intrusive approach for the simplified implementation of space and time multiscale methods within an explicit dynamics software. The ‘substitution’ method proposed in previous works allows to take advantage of a global coarse model, typically used in an industrial context, running separate, refined in space and in time, local analyses only where needed. The proposed technique is iterative, but the explicit character of the method allows to perform the global computation only once per global time step, while a repeated solution is required for the small local problems only. Nevertheless, a desirable goal is to reach convergence with a reduced number of iterations. To this purpose, we propose here a new iterative algorithm based on an improved interface inertia operator. The new operator exploits a combined property of velocity Hermite time interpolation on the interface and of the central difference integration scheme, allowing the consistent upscaling of interface inertia contributions from the lower scale. This property is exploited to construct an improved mass matrix operator for the interface coupling, allowing to significantly enhance the convergence rate. The efficiency and robustness of the procedure are demonstrated through several examples of growing complexity. Copyright © 2014 John Wiley & Sons, Ltd.

A global format is developed for momentum and energy consistent time integration of second-order dynamic systems with general nonlinear stiffness. The algorithm is formulated by integrating the state-space equations of motion over the time increment. The internal force is first represented in fourth-order form consisting of the end-point mean value plus a term containing the stiffness matrix increment. This form gives energy conservation for systems with internal energy as a quartic function of the displacement components. This representation is then extended to general energy conservation via a discrete gradient representation. The present procedure works directly with the internal force and the stiffness matrix at the time integration interval end-points, and in contrast to previous energy-conserving algorithms, it does not require any special form of the energy function nor use of mean value products at the element level or explicit use of a geometric stiffness matrix. An optional monotonic algorithmic damping, increasing with response frequency, is developed in terms of a single damping parameter. In the solution procedure, the velocity is eliminated and the nonlinear iterations are based on the displacement components alone. The procedure represents an energy consistent alternative to available collocation methods, with an equally simple implementation. Copyright © 2014 John Wiley & Sons, Ltd.

We present a method to reduce mesh bias in dynamic fracture simulations using the finite element method with adaptive insertion of extrinsic cohesive zone elements along element boundaries. The geometry of the domain discretization is important in this setting because cracks are only allowed to propagate along element facets and can potentially bias the crack paths. To reduce mesh bias, we consider unstructured polygonal finite elements in this work. The meshes are generated with centroidal Voronoi tessellations to ensure element quality. However, the possible crack directions at each node are limited, making this discretization a poor candidate for dynamic fracture simulation. To overcome this problem, and significantly improve crack patterns, we propose adaptive element splitting, whereby the number of potential crack directions is increased at each crack tip. *Thus, the crack is allowed to propagate through the polygonal element.* Geometric studies illustrate the benefits of polygonal element discretizations employed with element splitting over other structured and unstructured discretizations for crack propagation applications. Numerical examples are performed and demonstrate good agreement with previous experimental and numerical results in the literature. Copyright © 2014 John Wiley & Sons, Ltd.

The development of cell-centered finite volume discretizations for deformation is motivated by the desire for a compatible approach with the discretization of fluid flow in deformable porous media. We express the conservation of momentum in the finite volume sense, and introduce three approximations methods for the cell-face stresses. The discretization method is developed for general grids in one to three spatial dimensions, and leads to a global discrete system of equations for the displacement vector in each cell, after which the stresses are calculated based on a local expression. The method allows for anisotropic, heterogeneous and discontinuous coefficients.

The novel finite volume discretization is justified through numerical validation tests, designed to investigate classical challenges in discretization of mechanical equations. In particular our examples explore the stability with respect to the Poisson ratio and spatial discontinuities in the material parameters. For applications, logically Cartesian grids are prevailing, and we also explore the performance on perturbations on such grids, as well as on unstructured grids. For reference, comparison is made in all cases with the lowest-order Lagrangian finite elements, and the finite volume methods proposed herein is comparable for approximating displacement, and is superior for approximating stresses. Copyright © 2014 © 2014 The Authors. *International Journal for Numerical**Methods in Engineering* published by John Wiley & Sons Ltd.

Model order reduction (MOR) techniques for enriched reproducing kernel meshfree methods are proposed for analysis of Poisson problems with mild and strong singularities. The employment of an integrated singular basis function method (ISBFM), in conjunction with the selection of harmonic near-tip asymptotic basis functions, leads to a Galerkin formulation in which the non-smooth near-tip basis functions appear only on the boundaries away from the singularity point. This approach avoids the need of integrating the derivatives of non-smooth functions near the singularity point and yields a discrete system that allows effective MOR procedures. Under this framework, a decomposed reduction method equipped with two distinct projections for smooth and non-smooth parts of the finite-dimensional space is proposed. Compared with the uniform reduction approach using a single projection operator, the decomposed projection on ISBFM discrete system preserves the singularity behavior of the fine-scale solution in its lower-dimensional approximation. Analytical error estimation and stability analysis show that ISBFM with the decomposed projection can achieve better accuracy with only a slight increase of condition number compared with the uniform reduction approach in the reduced-order solution of singularity problems. The numerical tests validate the effectiveness of the proposed methods. Copyright © 2014 John Wiley & Sons, Ltd.

The goal of this work is to examine the efficacy of interpolatory model order reduction on frequency sweep problems with many forcing vectors. The adaptive method proposed relies on the implicit interpolatory properties of subspace projection with basis vectors spanning the forced response of the system and its derivatives. The algorithm is similar to a recently proposed adaptive scheme in that it determines both interpolation location and order within the frequency domain of interest. The bounds of efficiency of the approach as the number of forcing vectors grows are explored through the use of rough floating operation counts. This, in turn, guides choices made in the adaptive algorithm, including a new technique that prevents excessive subspace growth by capitalizing on subspace direction coupling. In order to demonstrate the algorithm's utility, a series of frequency sweep problems drawn from the field of structural acoustics is analyzed. It is demonstrated that increased order of interpolation can actually degrade the efficiency of the algorithm as the number of forcing vectors grows but that this can be limited by the subspace size throttling procedure developed here. Copyright © 2014 John Wiley & Sons, Ltd.

Mass transport processes are known to play an important role in many fields of biomechanics such as respiratory, cardiovascular, and biofilm mechanics. In this paper, we present a novel computational model considering the effect of local solid deformation and fluid flow on mass transport. As the transport processes are assumed to influence neither structure deformation nor fluid flow, a sequential one-way coupling of a fluid–structure interaction (FSI) and a multi-field scalar transport model is realized. In each time step, first the non-linear monolithic FSI problem is solved to determine current local deformations and velocities. Using this information, the mass transport equations can then be formulated on the deformed fluid and solid domains. At the interface, concentrations are related depending on the interfacial permeability. First numerical examples demonstrate that the proposed approach is suitable for simulating convective and diffusive scalar transport on coupled, deformable fluid and solid domains. Copyright © 2014 John Wiley & Sons, Ltd.

A natural kernel contact (NKC) algorithm under the framework of the semi-Lagrangian reproducing kernel particle method (semi-Lagrangian RKPM) is proposed to model multi-body contact with specific consideration for impact and penetration modeling. The NKC algorithm utilizes the interaction of the semi-Lagrangian kernel functions associated with contacting bodies to serve as the non-penetration condition. The effects of friction are represented by introducing a layer of the friction-like elasto-plasticity material between contacting bodies. This approach allows the frictional contact conditions and the associated kinematics to be naturally embedded in the semi-Lagrangian RKPM inter-particle force calculation. The equivalence in the Karush–Kuhn–Tucker conditions between the proposed NKC algorithm and the conventional contact kinematic constraints as well as the associated state variable relationships are identified. A level set method is further introduced in the NKC algorithm to represent the contact surfaces without pre-defined potential contact surfaces. The stability analysis performed in this work shows that temporal stability in the semi-Lagrangian RKPM with NKC algorithms is related to the velocity gradient between contacting bodies. The proposed methods have been verified by several benchmark problems and applied to the simulation of impact and penetration processes. Copyright © 2014 John Wiley & Sons, Ltd.

A two-phase flow eXtended Finite Element Method (XFEM) model is presented to analyse the injection and sequestration of carbon dioxide (CO_{2}) in deep saline aquifers. XFEM is introduced to accurately approximate near-injection well pressure behaviour with elements significantly larger than the injection well diameter. We present a vertically averaged multiphase flow model that combines XFEM to approximate the pressure field, with a Streamline Upwind/Finite Element Method/Finite Difference Method (SU-FEM-FDM) to approximate the distribution of CO_{2} in the aquifer. Near-well enrichment functions are presented along with the solution procedure for the coupled problem. Two examples are presented: in the first, CO_{2} injection into a perfectly horizontal aquifer is modelled with both XFEM and FEM-based methods. The results suggest that XFEM is able to provide low relative errors in the pressure near the well at a reduced computational cost compared with FEM. The impact and selection of the stabilization coefficient of the SU-FEM-FDM is also discussed. In the second example, the XFEM and SU-FEM-FDM model is applied to a more realistic problem of an inclined aquifer to demonstrate the ability of the model to capture the buoyancy-driven migration of CO_{2} in a deep saline aquifer. Copyright © 2014 John Wiley & Sons, Ltd.

We present a novel multiscale algorithm for nondestructive detection of multiple flaws in structures, within an inverse problem type setting. The key idea is to apply a two-step optimization scheme, where first rough flaw locations are quickly determined, and then, fine tuning is applied in these localized subdomains to obtain global convergence to the true flaws.

The two-step framework combines the strengths of heuristic and gradient-based optimization methods. The first phase employs a discrete-type optimization in which the optimizer is limited to specific flaw locations and shapes, thus converting a continuous optimization problem in the entire domain into a coarse discrete optimization problem with limited number of choices. To this end, we develop a special algorithm called discrete artificial bee colony. The second phase employs a gradient-based optimization of the Broyden–Fletcher–Goldfarb–Shanno type on local well-defined and bounded subdomains determined in the previous phase. A semi-analytical approach is developed to compute the stiffness derivative associated with the evaluation of objective function gradients.

The eXtended FEM (XFEM), with both circular and elliptical void enrichment functions, is used to solve the forward problem and alleviate the costly remeshing of every candidate flaw, in both optimization steps.

The multiscale algorithm is tested on several benchmark examples to identify various numbers and types of flaws with arbitrary shapes and sizes (e.g., cracks, voids, and their combination), without knowing the number of flaws beforehand. We study the size effect of the pseudo grids in the first optimization step and consider the effect of modeling error and measurement noise. The results are compared with the previous work that employed a single continuous optimization scheme (XFEM–genetic algorithm and XFEM–artificial bee colony methods). We illustrate that the proposed methodology is robust, yields accurate flaw detection results, and in particular leads to significant improvements in convergence rates compared with the previous work. Copyright © 2014 John Wiley & Sons, Ltd.

Many practical applications require the analysis of elastic wave propagation in a homogeneous isotropic media in an unbounded domain. One widely used approach for truncating the infinite domain is the so-called method of perfectly matched layers (PMLs). Most existing PML formulations are developed for finite difference methods based on the first-order velocity-stress form of the elasticity equations, and they are not straight-forward to implement using standard finite element methods (FEMs) on unstructured meshes. Some of the problems with these formulations include the application of boundary conditions in half-space problems and in the treatment of edges and/or corners for time-domain problems. Several PML formulations, which do work with FEMs have been proposed, although most of them still have some of these problems and/or they require a large number of auxiliary nodal history/memory variables. In this work, we develop a new PML formulation for time-domain elastodynamics on a spherical domain, which reduces to a two-dimensional formulation under the assumption of axisymmetry. Our formulation is well-suited for implementation using FEMs, where it requires lower memory than existing formulations, and it allows for natural application of boundary conditions. We solve example problems on two-dimensional and three-dimensional domains using a high-order discontinuous Galerkin (DG) discretization on unstructured meshes and explicit time-stepping. We also study an approach for stabilization of the discrete equations, and we show several practical applications for quality factor predictions of micromechanical resonators along with verifying the accuracy and versatility of our formulation. Copyright © 2014 John Wiley & Sons, Ltd.

We develop an efficient semi-local method for speeding up the solution of linear systems arising in spectral/hp element discretization of the linear elasticity equations. The main idea is to approximate the element-wise residual distribution with a localization operator we introduce in this paper, and subsequently solve the local linear system. Additionally, we decouple the three directions of displacement in the localization operator, hence enabling the use of an efficient low energy preconditioner for the conjugate gradient solver. This approach is effective for both nodal and modal bases in the spectral/hp element method, but here, we focus on the modal hierarchical basis. In numerical tests, we verify that there is no loss of accuracy in the semi-local method, and we obtain good parallel scalability and substantial speed-up compared to the original formulation. In particular, our tests include both structure-only and fluid-structure interaction problems, with the latter modeling a 3D patient-specific brain aneurysm. Copyright © 2014 John Wiley & Sons, Ltd.

Of fundamental interest are multidisciplinary interactions encompassing: (1) first order systems such as those encountered in parabolic heat conduction, first order hyperbolic systems such as fluid flow, and so on, and (2) second order systems such as those encountered in hyperbolic heat conduction, hyperbolic second order systems such as elastodynamics and wave propagation, and so on. After space discretization using methods such as finite differences, finite volumes, finite elements, and the like, the consequent proper integration of the time continuous ordinary differential equations is extremely important. In particular, the physical quantities of interest may need to be mostly preserved and/or the equations should be optimally integrated so that there is minimal numerical dissipation, dispersion, algorithm overshoot, capture shocks without too much dissipation, solve stiff problems and enable the completion of the analysis, and so on. To-date, practical methods in most commercial and research software include the trapezoidal family (Euler forward/backward, Galerkin, and Crank Nicholson) for first order systems and the other counterpart trapezoidal family (Newmark family and variants with controllable numerical dissipation) for second order systems. For the respective first/second order systems, they are totally separate families of algorithms and are derived from altogether totally different numerical approximation techniques. Focusing on the class of the linear multistep (LMS) methods, algorithms by design was first utilized to develop GS4-2 framework for time integration of second order systems. We have also recently developed the GS4-1 framework for integrating first order systems. In contrast to all past efforts over the past 50 years or so, we present the formalism of a generalized unified framework, termed GS4 (generalized single step single solve), that unifies GS4-1 (first-order systems) and GS4-2 (second-order systems) frameworks for simultaneous use in both first and second order systems with optimal algorithms, numerical and order preserving attributes (in particular, second-order time accuracy) as well. The principal contribution emanating from such an integrated framework is the practicality and convenience of using the same computational framework and implementation when solving first and/or second order systems without having to resort to the individual framework. All that is needed is a single novel GS4-2 framework for either second- and/or first-order systems, and we show how to switch from one to the other for illustrative applications to thermo-mechanical problems influenced by first/second order systems, respectively. Copyright © 2014 John Wiley & Sons, Ltd.

In this paper, we target more advanced fluid–structure interaction (FSI) simulations of wind turbines than reported previously. For this, we illustrate how the recent advances in isogeometric analysis of thin structures may be used for efficient structural mechanics modeling of full wind turbine structures, including tower, nacelle, and blades. We consider both horizontal axis and vertical axis wind turbine designs. We enhance the sliding–interface formulation of aerodynamics, previously developed to handle flows about mechanical components in relative motion such as rotor–tower interaction to allow nonstationary sliding interfaces. To accommodate the nonstationary sliding interfaces, we propose a new mesh moving technique and present its mathematical formulation. The numerical examples include structural mechanics verification for the new offshore wind turbine blade design, FSI simulation of a horizontal axis wind turbine undergoing yawing motion as it turns into the wind and FSI simulation of a vertical axis wind turbine. The FSI simulations are performed at full scale and using realistic wind conditions and rotor speeds. Copyright © 2014 John Wiley & Sons, Ltd.

Crack propagation in brittle materials with anisotropic surface energy is important in applications involving single crystals, extruded polymers, or geological and organic materials. Furthermore, when this anisotropy is strong, the phenomenology of crack propagation becomes very rich, with forbidden crack propagation directions or complex sawtooth crack patterns. This problem interrogates fundamental issues in fracture mechanics, including the principles behind the selection of crack direction. Here, we propose a variational phase-field model for strongly anisotropic fracture, which resorts to the extended Cahn-Hilliard framework proposed in the context of crystal growth. Previous phase-field models for anisotropic fracture were formulated in a framework only allowing for weak anisotropy. We implement numerically our higher-order phase-field model with smooth local maximum entropy approximants in a direct Galerkin method. The numerical results exhibit all the features of strongly anisotropic fracture and reproduce strikingly well recent experimental observations. Copyright © 2014 John Wiley & Sons, Ltd.

We present an adaptive B-spline FEM that employs hierarchical refinement and coarsening techniques. Hierarchical B-spline refinement augments the approximation space through the use of finer, more compactly supported B-splines. With this approach compatibility is automatic because problematic T-vertices do not arise. The process of refinement does not change the regularity of the basis functions, and can be easily generalized to arbitrary dimensions. Data transfer during both refinement and coarsening is addressed. The techniques are then applied to several time-dependent problems for which adaptivity is advantageous, and the performance of the method is evaluated. Copyright © 2014 John Wiley & Sons, Ltd.

A posteriori estimation of *defeaturing error*, that is, engineering analysis error caused by defeaturing, remains a key challenge in computer-aided design/computer-aided engineering integration; indeed, it is essential if analysis based on automatic simplification of the geometry of a complex design is to be reliable. Previous work has mainly focused on removing a *single, negative, internal* design feature (i.e. a void within the model's interior); effects of removal of multiple features are found by simple summation. In this paper, we give a *second-order* defeaturing error estimator that can handle *multiple boundary* features, either *positive* or *negative* (cutouts or additions on the model's boundary). This is the first such method to take into account *interactions* between different features. (Removing internal holes in 2D or through holes in 3D will change the model's topology and is not addressed here.) The proposed estimator is also the first to handle positive features without heuristics. Our approach uses *second-order shape sensitivity*, on the basis of reformulating defeaturing as a shape transformation process. The reference model used for sensitivity computation is carefully chosen, as is the associated shape variation velocity, to ensure efficient computation and simple sensitivity expressions. Mixed partial derivatives in a second-order Taylor expansion are used to account for interaction between features. The advantages of our approach are demonstrated using various numerical examples. Copyright © 2014 John Wiley & Sons, Ltd.

A continuum shell element based on the isogeometric analysis concept is extended to model propagating delaminations that can occur in composite materials and structures. The interpolation in the thickness direction is carried out using a quadratic B-spline, and delamination is modelled by a double-knot insertion to reduce the inter-layer continuity. Within the discontinuity, the traction is derived from the relative displacement between the layers by a cohesive relation. A range of examples, including delamination propagation in straight and curved planes and buckling-delamination, illustrate the versatility and the potential of the approach. Copyright © 2014 John Wiley & Sons, Ltd.

Metal stamping is one of the most cost effective manufacturing methods for producing precision parts. The precision of the parts is a function of the uniformity of the material used for the blank and the design of the dies. During the metal stamping process, the blank undergoes large, inelastic deformations, and once it is removed from the dies, springs back to a shape different from the die. This *spring back* must be accounted for in the design of the die, and its prediction is a major challenge in computational mechanics. Isogeometric analysis, which uses the same basis functions as the computer-aided design programs used to design the shape of the part, is an attractive alternative to traditional finite element analysis for metal stamping. Mass scaling, and the underlying stable time step estimates, that are commonly used in metal stamping simulations are presented for isogeometric analysis. Example calculations are presented to demonstrate the effectiveness of isogeometric analysis in metal stamping. Copyright © 2014 John Wiley & Sons, Ltd.

A method for stabilizing the mean-strain hexahedron is described that differs from the currently known approaches. For simplicity, the developments are limited to linear elasticity but with an arbitrarily anisotropic elasticity matrix. The technique relies on a sampling of the stabilization energy using two quadrature rules, the mean-strain quadrature and the full Gaussian integration rule. The use of two quadrature rules is shown to guarantee consistency and stability. The stabilization energy is assumed to be generated by a modified constitutive matrix based on the spectral decomposition. The spectral decomposition of the constitutive matrix identifies the stiff and flexible modes of deformation. The stiff modes of deformation are only sampled by the mean-strain integration, which eliminates volumetric locking for isotropic materials as well as locking due to strongly anisotropic material properties. The accuracy and convergence characteristics of the present formulations compare favorably with the capabilities of mean-strain and other high-performance hexahedral elements as implemented in ABAQUS. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents a generalized FEM based on the solution of interdependent coarse-scale (global) and fine-scale (local) problems in order to resolve multiscale effects due to fine-scale heterogeneities. Overall structural behavior is captured by the global problem, while local problems focus on the resolution of fine-scale solution features in regions where material heterogeneities may govern the structural response. Fine-scale problems are accurately solved in parallel, and, to address the intrinsic coupling of scales, these solutions are embedded into the global solution space using a partition of unity approach. This method is demonstrated on representative heat transfer examples in order to examine its accuracy, efficiency, and flexibility. Copyright © 2014 John Wiley & Sons, Ltd.

Stress concentrations near grain boundaries, precipitates, and similar micro-heterogeneities nucleate instabilities leading to macroscale fracture. As it is not practical to model each flaw explicitly, their ensemble effect is modeled statistically. Accounting for this aleatory uncertainty requires smaller specimens (e.g., small finite elements) to have generally higher and more variable strengths, which is necessary for the initial failure probability of a finite domain to be unaffected by its discretization into elements. Localization itself, which might be attributed to constitutive instability, requires realistic numerical perturbations to predict bifurcations such as radial cracking in axisymmetric problems. These perturbations, stemming from microscale heterogeneity, are incorporated in simulations by imposing statistical spatial variability in the parameters of an otherwise conventional (deterministic and scale-independent) damage model. This approach is attractive for its algorithmic simplicity and straightforward calibration from standard strength tests. In addition, it results in virtually no loss of efficiency or robustness relative to deterministic models and accommodates general three-dimensional loading. Despite these advantages, some significant challenges remain and are discussed. However, it is demonstrated that including aleatory uncertainty with associated scale effects significantly improves predictiveness on large-scale computational domains, where it is impractical to resolve each crack or localization zone.Copyright © 2014 John Wiley & Sons, Ltd.

This paper is dedicated to Professor Ted Belytscho's 70th birthday. His seminal contributions to Fracture Mechanics and the extended finite element method have significantly inspired us in this work. The current paper extends our recent work on the extraction of stress intensity factors (SIFs) from Irwin's integral, using a high-order extended FEM (XFEM) formulation. By matching leading *r* terms (*r* being the distance from the crack tip to any other material point) in XFEM with analytical expansion of Irwin's Integral, a closed-form solution for SIFs is obtained. Hence, SIFs can directly be obtained upon solution of the XFEM discrete system, which were shown to be accurate on structured rectangular meshes. However, extension to unstructured quadrilateral elements may not be trivial. To this end, we extend the formulation to triangular elements, where significant simplification of the closed-form expression is obtained. Moreover, the formulation is directly applicable to unstructured meshes without extra work. Hence, this paper is considered a significant step toward automating and extending this method to more general applications. Numerical results show that accurate and consistent SIFs can be obtained on some pure mode and inclined crack benchmark problems, on structured and unstructured meshes. Examples of a crack approaching a hole and two cracks approaching each other are also investigated. The latter illustrates the advantage of this method over a J-integral class of methods, as SIFs can still be calculated when cracks are in proximity and no remeshing is required. Hence, potentially, this method can address crack coalescence, branching and proximity to boundaries more rigorously. Copyright © 2014 John Wiley & Sons, Ltd.

One of the advantages of partition-of-unity FEMs, like the extended FEM, is the ability of modeling discontinuities independent of the mesh structure. The enrichment of the element functional space with discontinuous or non-differentiable functions requires, when the element stiffness is computed, partitioning into subdomains for quadrature. However, the arbitrary intersection between the base mesh and the discontinuity plane generates quadrature subdomains of complex shape. This is particularly true in three-dimensional problems, where quite sophisticate methodologies have been presented in the literature for the element stiffness evaluation.

The present work addresses the problem of Heaviside function enrichments and is based on the replacement of the discontinuous enrichment function with the limit of an equivalent polynomial defined on the entire element domain. This allows for the use of standard Gaussian quadrature in the elements crossed by the discontinuity. The work redefines conceptually the first version of the equivalent polynomial methodology introduced in 2006, allowing a much broader applicability. As a consequence, equivalent polynomials can be computed for all continuum element families in one, two, and three dimensions. Copyright © 2014 John Wiley & Sons, Ltd.

In the present study, a general dynamic data-driven application system (DDDAS) is developed for real-time monitoring of damage in composite materials using methods and models that account for uncertainty in experimental data, model parameters, and in the selection of the model itself. The methodology involves (i) data data from uniaxial tensile experiments conducted on a composite material; (ii) continuum damage mechanics based material constitutive models; (iii) a Bayesian framework for uncertainty quantification, calibration, validation, and selection of models; and (iv) general Bayesian filtering, as well as Kalman and extended Kalman filters. A software infrastructure is developed and implemented in order to integrate the various parts of the DDDAS. The outcomes of computational analyses using the experimental data prove the feasibility of the Bayesian-based methods for model calibration, validation, and selection. Moreover, using such DDDAS infrastructure for real-time monitoring of the damage and degradation in materials results in results in an improved prediction of failure in the system. Copyright © 2014 John Wiley & Sons, Ltd.

In this work, an explicit family of time marching procedures with adaptive dissipation control is introduced. The proposed technique is conditionally stable, second-order accurate, and has controllable algorithm dissipation, which adapts according to the properties of the governing system of equations. Thus, spurious modes can be more effectively dissipated and accuracy is improved. Because this is an explicit time integration technique, the new family is quite efficient, requiring no system of equations to be dealt with at each time step. Moreover, the technique is simple and very easy to implement. Numerical results are presented along the paper, illustrating the good performance of the proposed method, as well as its potentialities. Copyright © 2014 John Wiley & Sons, Ltd.

The problem of representing random fields describing the material and boundary properties of the physical system at discrete points of the spatial domain is studied in the context of linear stochastic finite element method. A randomly parameterized diffusion system with a set of independent identically distributed stochastic variables is considered. The discretized parametric fields are interpolated within each element with multidimensional Lagrange polynomials and integrated into the weak formulation. The proposed discretized random-field representation has been utilized to express the random fluctuations of the domain boundary with nodal position coordinates and a set of random variables. The description of the boundary perturbation has been incorporated into the weak stochastic finite element formulation using a stochastic isoparametric mapping of the random domain to a deterministic master domain. A method for obtaining the linear system of equations under the proposed mapping using generic finite element weak formulation and the stochastic spectral Galerkin framework is studied in detail. The treatment presents a unified way of handling the parametric uncertainty and random boundary fluctuations for dynamic systems. The convergence behavior of the proposed methodologies has been demonstrated with numerical examples to establish the validity of the numerical scheme. Copyright © 2014 John Wiley & Sons, Ltd.

The present work aims to look into the contribution of the extended finite element method for large deformation of cracked bodies in plane strain approximation. The unavailability of sufficient mathematical tools and proofs for such problem makes the study exploratory. First, the asymptotic solution is presented. Then, a numerical analysis is realized to verify the pertinence of solution given by the asymptotic procedure, because it serves as an eXtended finite element method enrichment basis. Finally, a convergence study is carried out to show the contribution of the exploitation of such method. Copyright © 2014 John Wiley & Sons, Ltd.