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.

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.

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.

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.

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.

A direct solver for symmetric sparse matrices from finite element problems is presented. The solver is supposed to work as a local solver of domain decomposition methods for hybrid parallelization on cluster systems of multi-core CPUs, and then it is required to run on shared memory computers and to have an ability of kernel detection. Symmetric pivoting with a given threshold factorizes a matrix with a decomposition introduced by a nested bisection and selects suspicious null pivots from the threshold. The Schur complement constructed from the suspicious null pivots is examined by a factorization with 1 × 1 and 2 × 2 pivoting and by a robust kernel detection algorithm based on measurement of residuals with orthogonal projections onto supposed image spaces. A static data structure from the nested bisection and a block sub-structure for Schur complements at all bisection levels can use level 3 BLAS routines efficiently. Asynchronous task execution for each block can reduce idle time of processors drastically, and as a result, the solver has high parallel efficiency. Competitive performance of the developed solver to Intel Pardiso on shared memory computers is shown by numerical experiments. 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.

This paper is devoted to the formulation of a plane scaled boundary finite element with initially constant thickness for physically and geometrically nonlinear material behavior. Special two-dimensional element shape functions are derived by using the analytical displacement solution of the standard scaled boundary finite element method, which is originally based on linear material behavior and small strains. These 2D shape functions can be constructed for an arbitrary number of element nodes and allow to capture singularities (e.g., at a plane crack tip) analytically, without extensive mesh refinement. Mapping these proposed 2D shape functions to the 3D case, a formulation that is compatible with standard finite elements is obtained. The resulting physically and geometrically nonlinear scaled boundary finite element formulation is implemented into the framework of the finite element method for bounded plane domains with and without geometrical singularities. The numerical realization is shown in detail. To represent the physically and geometrically nonlinear material and structural behavior of elastomer specimens, the extended tube model and the Yeoh model are used. Numerical studies on the convergence behavior and comparisons with standard Q1P0 finite elements demonstrate the correct implementation and the advantages of the developed scaled boundary finite element. Copyright © 2014 John Wiley & Sons, Ltd.

The use of cohesive zone models is an efficient way to treat the damage especially when the crack path is known *a priori*. It is the case in the modeling of delamination in composite laminates. However, the simulations using cohesive zone models are expensive in a computational point of view. When using implicit time integration or when solving static problems, the non-linearity related to the cohesive model requires many iteration before reaching convergence. In explicit approaches, an important number of iterations are also needed because of the time step stability condition. In this article, a new approach based on a separated representation of the solution is proposed. The proper generalized decomposition is used to build the solution. This technique coupled with a cohesive zone model allows a significant reduction of the computational cost. The results approximated with the proper generalized decomposition are very close the ones obtained using the classical finite element approach.Copyright © 2014 John Wiley & Sons, Ltd.

In this work, a linear hexahedral element based on an assumed strain finite element technique is presented for the solution of plasticity problems. The element stems from the Nodally Integrated Continuum Element (NICE) formulation and its extensions. Assumed gradient operators are derived via nodal integration from the kinematic-weighted residual; the degrees of freedom are only the displacements at the nodes. The adopted constitutive model is the classical associative von Mises plasticity model with isotropic and kinematic hardening; in particular, a double-step midpoint integration algorithm is adopted for the integration and solution of the relevant nonlinear evolution equations. Efficiency of the proposed method is assessed through simple benchmark problems and comparison with reference solutions. Copyright © 2014 John Wiley & Sons, Ltd.

Phase-field approaches to fracture offer new perspectives toward the numerical solution of crack propagation. In this paper, a phase-field method for finite deformations and general nonlinear material models is introduced using a novel multiplicative split of the principal stretches to account for the different behavior of fracture in tension and compression. An energy-momentum consistent integrator is developed and applied to the arising nonlinear coupled phase-field model. This approach is thermodynamically consistent in the sense that the first law of thermodynamics if fulfilled with respect to the dissipation function. The capabilities and the performance of the proposed approach is demonstrated in several representative examples. 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 two-level domain decomposition method is introduced for general shape optimization problems constrained by the incompressible Navier–Stokes equations. The optimization problem is first discretized with a finite element method on an unstructured moving mesh that is implicitly defined without assuming that the computational domain is known and then solved by some one-shot Lagrange–Newton–Krylov–Schwarz algorithms. In this approach, the shape of the domain, its corresponding finite element mesh, the flow fields and their corresponding Lagrange multipliers are all obtained computationally in a single solve of a nonlinear system of equations. Highly scalable parallel algorithms are absolutely necessary to solve such an expensive system. The one-level domain decomposition method works reasonably well when the number of processors is not large. Aiming for machines with a large number of processors and robust nonlinear convergence, we introduce a two-level inexact Newton method with a hybrid two-level overlapping Schwarz preconditioner. As applications, we consider the shape optimization of a cannula problem and an artery bypass problem in 2D. Numerical experiments show that our algorithm performs well on a supercomputer with over 1000 processors for problems with millions of unknowns. Copyright © 2014 John Wiley & Sons, Ltd.

A new method is developed here for the real-time integration of the equations of solid dynamics based on the use of proper orthogonal decomposition (POD)–proper generalized decomposition (PGD) approaches and direct time integration. The method is based upon the formulation of solid dynamics equations as a parametric problem, depending on their initial conditions. A sort of black-box integrator that takes the resulting displacement field of the current time step as input and (via POD) provides the result for the subsequent time step at feedback rates on the order of 1 kHz is obtained. To avoid the so-called *curse of dimensionality* produced by the large amount of parameters in the formulation (one per degree of freedom of the full model), a combined POD–PGD strategy is implemented. Examples that show the promising results of this technique are included. Copyright © 2014 John Wiley & Sons, Ltd.

A simple and compact representation framework and the corresponding efficient numerical integration algorithm are developed for constitutive equations of isotropic elastoplasticity. Central to this work is the utilization of a set of mutually orthogonal unit tensor bases and the corresponding invariants. The set of bases can be regarded equivalently as a local cylindrical coordinate system in the three-dimensional coaxial tensor subspace, namely, the principal space. The base tensors are given in the global coordinate system. Similar to the principal space approach, the proposed method reduces the problem dimension from six to three. In contrast to the conventional approach, the transformation procedure between the principal space and the general space is avoided and explicit computation of the principal axes is bypassed. With the proposed technique, the matrices, which need to be inverted during iteration, take a simple form for the great majority of constitutive equations in use. The tangent operator consistent with the proposed algorithm can be decomposed into the direct sum of two linear maps over the coaxial tensor subspace and the subspace orthogonal to it. Consequently, its closed form is derived in an extremely simple manner. Finally, numerical examples demonstrate the high quality performances of the proposed scheme. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents an extension to the existing dynamic relaxation method to include equality constraint conditions in the process. The existing dynamic relaxation method is presented as a general, gradient-based, minimization technique. This representation allows for the introduction of the projected gradient, discrete parallel transportation and pull back operators that enable the formulation of the geodesic dynamic relaxation method, a method that accounts for equality constraint conditions. The characteristics of both the existing and geodesic dynamic relaxation methods are discussed in terms of the system's conservation of energy, damping (viscous, kinetic, and drift), and geometry generation. Particular attention is drawn to the introduction of a novel damping approach named drift damping. This technique is essentially a combination of viscous and kinetic damping. It allows for a smooth and fast convergence rate in both the existing and geodesic dynamic relaxation processes. The case study was performed on the form-finding of an iconic, ridge-and-valley, pre-stressed membrane system, which is supported by masts. The study shows the potential of the proposed method to account for specified (total) length requirements. The geodesic dynamic relaxation technique is widely applicable to the form-finding of force-modeled systems (including mechanically and pressurized pre-stressed membranes) where equality constraint control is desired. Copyright © 2014 John Wiley & Sons, Ltd.