Quadrilateral and triangular elements with curved edges are developed in the framework of spectral, discontinuous, hybrid control-volume/finite-element method for elliptic problems. In order to accommodate hybrid meshes, encompassing both triangular and quadrilateral elements, one single mapping is used. The scheme is applied to two-dimensional problems with discontinuous, anisotropic diffusion coefficients, and the exponential convergence of the method is verified in the presence of curved geometries. Copyright © 2014 John Wiley & Sons, Ltd.

We describe a new mesh smoothing method that consists of minimizing the sum of squared element volumes over the free vertex positions. To the extent permitted by the fixed vertices and mesh topology, the resulting mesh elements have uniformly distributed volumes. In the case of a triangulation, uniform volume implies well-shaped triangles. If a graded mesh is required, the element volumes may be weighted by centroidal values of a sizing function, resulting in a mesh that conforms to the required vertex density. The method has both a local and a global implementation. In addition to smoothing, the method serves as a simple parameter-free means of untangling a mesh with inverted elements. It applies to all types of meshes, but we present test results here only for planar triangle meshes. Our test results show the new method to be fast, superior in uniformity or conformity to a sizing function, and among the best methods in terms of triangle shape quality. We also present a new angle-based method that is simpler and more effective than alternatives. This method is directly aimed at producing well-shaped triangles and is particularly effective when combined with the volume-based method. It also generalizes to anisotropic mesh smoothing. Copyright © 2014 John Wiley & Sons, Ltd.

The purpose of this paper is to study the effect of the bulk modulus on the iterative solution of free surface quasi-incompressible fluids using a mixed partitioned scheme. A practical rule to set up the value of a pseudo-bulk modulus a priori in the tangent bulk stiffness matrix for improving the conditioning of the linear system of algebraic equations is also given. The efficiency of the proposed strategy is tested in several problems analyzing the advantage of the modified bulk tangent matrix with regard to the stability of the pressure field, the convergence rate and the computational speed of the analyses. The technique has been tested on a finite calculus/particle finite element method Lagrangian formulation, but it can be easily extended to other quasi-incompressible stabilized finite element formulations. Copyright © 2014 John Wiley & Sons, Ltd.

Dynamic crack-branching instabilities in a brittle material are studied numerically by using a non-local damage model. PMMA is taken as our model brittle material. The simulated crack patterns, crack velocities, and dissipated energies compare favorably with experimental data gathered from the literature, as long as the critical strain for damage initiation as well as the parameters for a rate-dependent damage law are carefully selected. Nonetheless, the transition from a straight crack propagation to the emergence of crack branches is very sensitive to the damage initiation threshold. The transition regime is thus a particularly interesting challenge for numerical approaches. We advocate using the present numerical study as a benchmark to test the robustness of alternative non-local numerical approaches. Copyright © 2014 John Wiley & Sons, Ltd.

We adopt a numerical method to solve Poisson's equation on a fixed grid with embedded boundary conditions, where we put a special focus on the accurate representation of the normal gradient on the boundary. The lack of accuracy in the gradient evaluation on the boundary is a common issue with low-order embedded boundary methods. Whereas a direct evaluation of the gradient is preferable, one typically uses post-processing techniques to improve the quality of the gradient. Here, we adopt a new method based on the discontinuous-Galerkin (DG) finite element method, inspired by the recent work of [A.J. Lew and G.C. Buscaglia. A discontinuous-Galerkin-based immersed boundary method. *International Journal for Numerical Methods in Engineering*, 76:427-454, 2008]. The method has been enhanced in two aspects: firstly, we approximate the boundary shape locally by higher-order geometric primitives. Secondly, we employ higher-order shape functions within intersected elements. These are derived for the various geometric features of the boundary based on analytical solutions of the underlying partial differential equation. The development includes three basic geometric features in two dimensions for the solution of Poisson's equation: a straight boundary, a circular boundary, and a boundary with a discontinuity. We demonstrate the performance of the method via analytical benchmark examples with a smooth circular boundary as well as in the presence of a singularity due to a re-entrant corner. Results are compared to a low-order extended finite element method as well as the DG method of [1]. We report improved accuracy of the gradient on the boundary by one order of magnitude, as well as improved convergence rates in the presence of a singular source. In principle, the method can be extended to three dimensions, more complicated boundary shapes, and other partial differential equations. Copyright © 2014 John Wiley & Sons, Ltd.

In this work, we develop an isogeometric non-uniform rational B-spline (NURBS)-based solid-shell element for the geometrically nonlinear static analysis of elastic shell structures. A single layer of continuous 3D elements through the thickness of the shell is considered, and the order of approximation in that direction is chosen to be equal to two. A complete 3D constitutive relation is assumed. The objective is to develop a highly accurate low-order element for coarse meshes. We propose an extension of the mixed method of Bouclier *et al.* [11] to deal with locking in the context of large rotations and large displacements. The main idea is to modify the interpolation of the average through the thickness of the stress components. It is also necessary to stabilize the element in order to avoid the occurrence of spurious zero-energy modes. This was achieved, for the quadratic version, through the adjunction of artificial elementary stabilization stiffnesses. The result is an element of order 2, which is at least as accurate as standard NURBS shell elements of order 4. Linear and nonlinear test calculations have been carried out along with comparisons with other published NURBS and classical techniques in order to assess the performance of the element. Copyright © 2014 John Wiley & Sons, Ltd.

A novel nonlinear parametric model order reduction technique for the solution of contact problems in flexible multibody dynamics is presented. These problems are characterized by significant variations in the location and size of the contact area and typically require high-dimensional finite element models having multiple inputs and outputs to be solved. The presented technique draws from the fields of nonlinear and parametric model reduction to construct a reduced-order model whose dimensions are insensitive to the dimensions of the full-order model. The solution of interest is approximated in a lower-dimensional subspace spanned by a constant set of eigenvectors augmented with a parameter-dependent set of *global contact shapes*. The latter represent deformation patterns of the interacting bodies obtained from a series of static contact analyses. The set of global contact shapes is parameterized with respect to the system configuration and therefore continuously varies in time. An energy-consistent formulation is assured by explicitly taking into account the dynamic parameter variability in the derivation of the equations of motion. The performance of the novel technique is demonstrated by simulating a dynamic gear contact problem and comparing results against traditional model reduction techniques as well as commercial nonlinear finite element software. Copyright © 2014 John Wiley & Sons, Ltd.

Effective simulation of the solid-liquid-gas coupling effect in unsaturated porous media is of great significance in many diverse areas. Because of the strongly nonlinear characteristics of the fully coupled formulations for the three-phase porous media, an effective numerical solution scheme, such as the finite element method with an efficient iterative algorithm, has to be employed. In this paper, an efficient finite element procedure based on the adaptive relaxed Picard method is developed for analyzing the coupled solid-liquid-gas interactions in porous media. The coupled model and the finite element analysis procedure are implemented into a computer code PorousH2M, and the proposed procedure is validated through comparing the numerical simulations with the experimental benchmarks. It is shown that the adaptive relaxed Picard method has salient advantage over the traditional one with respect to both the efficiency and the robustness, especially for the case of relatively large time step sizes. Compared with the Newton-Raphson scheme, the Picard method successfully avoids the unphysical ‘spurious unloading’ phenomenon under the plastic deformation condition, although the latter shows a better convergence rate. The proposed procedure provides an important reference for analyzing the fully coupled problems related to the multi-phase, multi-field coupling in porous media. Copyright © 2014 John Wiley & Sons, Ltd.

An adaptive approach to using reduced-order models (ROMs) as surrogates in partial differential equations (PDE)-constrained optimization is introduced that breaks the traditional offline-online framework of model order reduction. A sequence of optimization problems constrained by a given ROM is defined with the goal of converging to the solution of a given PDE-constrained optimization problem. For each reduced optimization problem, the constraining ROM is trained from sampling the high-dimensional model (HDM) at the solution of some of the previous problems in the sequence. The reduced optimization problems are equipped with a nonlinear trust-region based on a residual error indicator to keep the optimization trajectory in a region of the parameter space where the ROM is accurate. A technique for incorporating sensitivities into a reduced-order basis is also presented, along with a methodology for computing sensitivities of the ROM that minimizes the distance to the corresponding HDM sensitivity, in a suitable norm.

The proposed reduced optimization framework is applied to subsonic aerodynamic shape optimization and shown to reduce the number of queries to the HDM by a factor of 4-5, compared with the optimization problem solved using only the HDM, with errors in the optimal solution far less than 0.1%. Copyright © 2014 John Wiley & Sons, Ltd.

This paper primarily deals with the computational aspects of chemical dissolution-front instability problems in two-dimensional fluid-saturated porous media under non-isothermal conditions. After the dimensionless governing partial differential equations of the non-isothermal chemical dissolution-front instability problem are briefly described, the formulation of a computational procedure, which contains a combination of using the finite difference and finite element method, is derived for simulating the morphological evolution of chemical dissolution fronts in the non-isothermal chemical dissolution system within two-dimensional fluid-saturated porous media. To ensure the correctness and accuracy of the numerical solutions, the proposed computational procedure is verified through comparing the numerical solutions with the analytical solutions for a benchmark problem. As an application example, the verified computational procedure is then used to simulate the morphological evolution of chemical dissolution fronts in the supercritical non-isothermal chemical dissolution system. The related numerical results have demonstrated the following: (1) the proposed computational procedure can produce accurate numerical solutions for the planar chemical dissolution-front propagation problem in the non-isothermal chemical dissolution system consisting of a fluid-saturated porous medium; (2) the Zhao number has a significant effect not only on the dimensionless propagation speed of the chemical dissolution front but also on the distribution patterns of the dimensionless temperature, dimensionless pore-fluid pressure, and dimensionless chemical-species concentration in a non-isothermal chemical dissolution system; (3) once the finger penetrates the whole computational domain, the dimensionless pore-fluid pressure decreases drastically in the non-isothermal chemical dissolution system.

A displacement and rotation-based dynamic finite element formulation for Cosserat plates is discussed in detail in this paper. Special attention is given to the validation of the element through adequate benchmarks and patch tests. The choice of the interpolation functions and the order of integration of the stiffness and the mass matrices are extensively argued. The possibility of local system deficiencies is explored, and a shear locking investigation specifically elaborated for Cosserat plates is carried out. It is shown how the present formulation has interesting computational properties as compared to a classical continuum-based formulation and how it can provide suitable results despite the use of reduced integration. An example of application of the finite element is given, in which the natural frequencies of a masonry panel modelled by means of discrete elements are computed and compared with the finite elements solution. Copyright © 2014 John Wiley & Sons, Ltd.

In this paper, we are interested in the dynamical response of a material body subjected to impact loadings. Such loadings are brutal and intense and may damage the material, leading to strain localization and rupture. Before strain localization occurs, computation of such problems is often accurate enough and very efficient when an explicit time integration scheme is applied. However, after strain localization occurs, the mathematical relevance of a model is preserved only if non-locality is introduced. This is often resulting in a dramatic increase of computational costs. We propose in this work to introduce non-locality through the Thick Level Set approach (TLS). It is the first time the TLS approach has been presented in a dynamical context. In this approach, additional computational efforts are limited in space to a domain slightly bigger than the strain localization region and the time discretization is explicit. The non-local computation is based on a new technique where basis functions are built on the damaged band. The resulting function space has needed properties to compute non-local fields. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents a coupling method between a discrete element code CeaMka3D and a finite element code Sem. The coupling is based on a least-squares method, which adds terms of forces to finite element code and imposes the velocity at coupling particles. For each coupling face, a small linear system with a constant matrix is solved. This method remains conservative in energy and shows good results in applications. Copyright © 2014 John Wiley & Sons, Ltd.

Computational structural dynamics plays an essential role in the simulation of linear and nonlinear systems. Indeed, the characteristics of the time integration procedure have a critical impact on the feasibility of the calculation. In order to go beyond the classical approach (a unique time integrator and a unique timescale), the pioneer approach of Belytschko and co-workers consisted in developing mixed implicit–explicit time integrators for structural dynamics. In a first step, the implementation and stability analyses of partitioned integrators with one time step have been achieved for a large class of time integrators. In a second step, the implementation and stability analyses of partitioned integrators with different time steps were studied in detail for particular cases. However, stability results involving different time steps and different time integrators in different parts of the mesh is still an open question in the general case for structural dynamics. The aim of this paper is to propose a state-of-the art of heterogeneous (different time schemes) asynchronous (different time steps) time integrators (HATI) for computational structural dynamics. Finally, an alternative approach based on energy considerations (with velocity continuity at the interface) is proposed in order to develop a general class of HATI for structural dynamics. Copyright © 2014 John Wiley & Sons, Ltd.

A general topology optimization method, which is capable of simultaneous design of density and orientation of anisotropic material, is proposed by introducing orientation design variables in addition to the density design variable. In this work, the Cartesian components of the orientation vector are utilized as the orientation design variables. The proposed method supports continuous orientation design, which is out of the scope of discrete material optimization approaches, as well as design using discrete angle sets. The advantage of this approach is that vector element representation is less likely to fail into local optima because it depends less on designs of former steps, especially compared with using the angle as a design variable (Continuous Fiber Angle Optimization) by providing a flexible path from one angle to another with relaxation of orientation design space. An additional advantage is that it is compatible with various projection or filtering methods such as sensitivity filters and density filters because it is free from unphysical bound or discontinuity such as the one at *θ* = 2*π* and *θ* = 0 seen with direct angle representation. One complication of Cartesian component representation is the point-wise quadratic bound of the design variables; that is, each pair of element values has to reside in a given circular bound. To overcome this issue, we propose an isoparametric projection method, which transforms box bounds into circular bounds by a coordinate transformation with isoparametric shape functions without having the singular point that is seen at the origin with polar coordinate representation. A new topology optimization method is built by taking advantage of the aforementioned features and modern topology optimization techniques. Several numerical examples are provided to demonstrate its capability. Copyright © 2014 John Wiley & Sons, Ltd.

This contribution investigates the performance of a least-squares finite element method based on non-uniform rational B-splines (NURBS) basis functions. The least-squares functional is formulated directly in terms of the strong form of the governing equations and boundary conditions. Thus, the introduction of auxiliary variables is avoided, but the order of the basis functions must be higher or equal to the order of the highest spatial derivatives. The methodology is applied to the incompressible Navier–Stokes equations and to linear as well as nonlinear elastic solid mechanics. The numerical examples presented feature convective effects and incompressible or nearly incompressible material. The numerical results, which are obtained with equal-order interpolation and without any stabilisation techniques, are smooth and accurate. It is shown that for *p* and *h* refinement, the theoretical rates of convergence are achieved. Copyright © 2014 John Wiley & Sons, Ltd.

A new continuous-discontinuous strategy for the simulation of failure is presented. The continuous bulk is regularised by means of a gradient-enhanced damage model, where non-locality is introduced at the level of displacements. As soon as the damage parameter is close or equal to 1, a traction-free crack is introduced. To determine the direction of crack growth, a new criterion is proposed. In contrast to traditional techniques, where mechanical criteria are used to define the crack path, here, a geometrical approach is used. More specifically, given a regularised damage field *D*(** x**), we propose to propagate the discontinuity following the direction dictated by the medial axis of the isoline (or isosurface in 3D)

In this work, we present a stable proper orthogonal decomposition–Galerkin approximation for parametrized steady incompressible Navier–Stokes equations with low Reynolds number. Copyright © 2014 John Wiley & Sons, Ltd.

This paper establishes the basic framework for the traction-based equilibrium finite element method (traction-based EFEM). Stable linear traction-based equilibrium elements are formulated using the macro-element technique. An arbitrary internal macro-point renders a linear triangular element stable, while a stable linear quadrilateral element requires the macro-point to locate at the intersection of diagonals. Then, a Lagrangian formulation is utilized to minimize the complementary energy under equilibrium constraints, and consequently, tractions as well as additional Lagrange multipliers are obtained. Linear statically admissible (SA) stresses are thereafter acquired from tractions. As for Lagrange multipliers, they turn out to coincide well with rigid-body displacements in each element after simple modifications. With rigid-body displacements and linear tractions known, quadratic displacements and the kinematically admissible (KA) counterpart thereof by recovery are obtainable. The knowledge of both SA stresses and KA displacements renders dual analysis directly applicable. That is to say, the traction-based EFEM is featured with direct access to strict upper and lower bounds of strain energy and other quantities of interest. Copyright © 2014 John Wiley & Sons, Ltd.

This paper proposes a level-set based topology optimization method incorporating a boundary tracking mesh generating method and nonlinear programming. Because the boundary tracking mesh is always conformed to the structural boundary, good approximation to the boundary is maintained during optimization; therefore, structural design problems are solved completely without grayscale material. Previously, we introduced the boundary tracking mesh generating method into level-set based topology optimization and updated the design variables by solving the level-set equation. In order to adapt our previous method to general structural optimization frameworks, the incorporation of the method with nonlinear programming is investigated in this paper. To successfully incorporate nonlinear programming, the optimization problem is regularized using a double-well potential. Furthermore, the sensitivities with respect to the design variables are strictly derived to maintain consistency in mathematical programming. We expect the investigation to open up a new class of grayscale-free topology optimization. The usefulness of the proposed method is demonstrated using several numerical examples targeting two-dimensional compliant mechanism and metallic waveguide design problems. Copyright © 2014 John Wiley & Sons, Ltd.

This paper deals with data transfer between two meshes as it happens in a finite element context when a remeshing has to be performed. We propose a finite-volume-based data transfer method for an efficient remeshing of three-dimensional solid mechanics problems. The originality of this transfer method stems from a linear reconstruction of the fields to be transferred on an auxiliary finite volume mesh, a fast computation of the transfer operator and the application to the complete remeshing of 3D problems. This procedure is applicable to both nodal values and discrete fields defined at quadrature points. In addition, a data transfer method using mortar elements is presented. The main improvement made to this second method comes from a fast computation of mortar elements. These two data transfer methods are compared with the simplest transfer method, which consists of a classical interpolation. After some academic examples, we present 2D forging and 3D friction stir welding applications. Copyright © 2014 John Wiley & Sons, Ltd.

In this paper, the effects of element shape on the critical time step are investigated. The common rule-of-thumb, used in practice, is that the critical time step is set by the shortest distance within an element divided by the dilatational (compressive) wave speed, with a modest safety factor. For regularly shaped elements, many analytical solutions for the critical time step are available, but this paper focusses on distorted element shapes. The main purpose is to verify whether element distortion adversely affects the critical time step or not. Two types of element distortion will be considered, namely aspect ratio distortion and angular distortion, and two particular elements will be studied: four-noded bilinear quadrilaterals and three-noded linear triangles. The maximum eigenfrequencies of the distorted elements are determined and compared to those of the corresponding undistorted elements. The critical time steps obtained from single element calculations are also compared to those from calculations based on finite element patches with multiple elements. Copyright © 2014 John Wiley & Sons, Ltd.

An efficient procedure for embedding kinematic boundary conditions in the biharmonic equation, for problems such as the pure streamfunction formulation of the Navier–Stokes equations and thin plate bending, is based on a stabilized variational formulation, obtained by Nitsche's approach for enforcing boundary constraints. The absence of kinematic admissibility constraints allows the use of non-conforming meshes with non-interpolatory approximations, thereby providing added flexibility in addressing the higher continuity requirements typical of these problems. Variationally conjugate pairs weakly enforce kinematic boundary conditions. The use of a scaling factor leads to a formulation with a single stabilization parameter. For plates, the enforcement of tangential derivatives of deflections obviates the need for pointwise enforcement of corner values in the presence of corners. The single stabilization parameter is determined from a local generalized eigenvalue problem, guaranteeing coercivity of the discrete bilinear form. The accuracy of the approach is verified by representative computations with bicubic B-splines, providing guidance to the determination of the scaling and exhibiting optimal rates of convergence and robust performance with respect to values of the stabilization parameter. Copyright © 2014 John Wiley & Sons, Ltd.

A parallel localization scheme is presented to enable solution transfers between unstructured grids. The scheme relies on neighbor walks to reduce the number of candidate elements that are visited to find the enclosing element. An advancing front method efficiently allows a subset of nodes to efficiently sweep through the grid, progressively reducing search spaces. The algorithm is parallelized permitting solution transfers over arbitrary grid decompositions. A hierarchical localization process helps prevent the neighbor walk algorithm from failing when encountering the boundaries of a concave domain by localizing the boundaries before the interior of the domain is localized. Random selections of the next step interrupt cyclic loops that may occur during a neighbor walk. The complexity of the search algorithm is verified over parallel decompositions and is effectively independent of the number of partitions. Copyright © 2014 John Wiley & Sons, Ltd.

We present an original algorithm and accompanying mathematical formulation for topology optimization of structures that can sustain material damage and are subject to multiple load cases with varying configurations. Damage accumulation is simulated using a coupled, non-linear brittle damage model. The structures are optimized for minimum mass subject to stiffness constraints defined as the compliance evaluated at the end of each loading sequence. To achieve robustness of the optimized structures, the respective damage fields caused by each individual load case are computed and combined using superposition to simulate a *worst-case* damage field. All load cases are then run a second time using the worst-case damage distribution as a starting point. In this way, one effectively accounts for the spectrum of possible load sequences to which the structure may be subjected. Results from this method are compared with an exhaustive, brute-force approach in which all non-repeating load sequences are analyzed individually. For each method, the corresponding sensitivities are derived and implemented analytically using a path-dependent adjoint method. The two approaches are implemented on a series of numerical examples, which demonstrate that the superposition method produces structures that are as robust as those obtained using the exhaustive method but require significantly less computational effort. Copyright © 2014 John Wiley & Sons, Ltd.

Ignoring crack tip effects, the stability of the X-FEM discretizations is trivial for open cracks but remains a challenge if we constrain the crack to be closed (i.e., the bi-material problem). Here, we develop a formulation for general cohesive interactions between crack faces within the X-FEM framework. The stability of the new formulation is demonstrated for any cohesive crack stiffness (including the closed crack) and illustrated for a nonlinear cohesive softening law. A benchmark of the new model is carried out with simpler approaches for a closed crack (i.e., Lagrange multipliers) and for a cohesive crack (i.e., penalty approach). Due to the analogies between stable cohesive X-FEM and Nitsche's methods, the new method simplifies the implementation and is attractive in dynamic explicit codes. Copyright © 2014 John Wiley & Sons, Ltd.

Nonlinear elastic materials are of great engineering interest, but challenging to model with standard finite elements. The challenges arise because nonlinear elastic materials are characterized by non-convex stored-energy functions as a result of their ability to undergo large reversible deformations, are incompressible or nearly incompressible, and often times possess complex microstructures. In this work, we propose and explore an alternative approach to model finite elasticity problems in two dimensions by using polygonal discretizations. We present both lower order displacement-based and mixed polygonal finite element approximations, the latter of which consist of a piecewise constant pressure field and a linearly-complete displacement field at the element level. Through numerical studies, the mixed polygonal finite elements are shown to be stable and convergent. For demonstration purposes, we deploy the proposed polygonal discretization to study the nonlinear elastic response of rubber filled with random and periodic distributions of rigid particles, as well as the development of cavitation instabilities in elastomers containing vacuous defects. These physically-based examples illustrate the potential of polygonal finite elements in studying and modeling nonlinear elastic materials with complex microstructures under finite deformations. Copyright © 2014 John Wiley & Sons, Ltd.

This work presents a method to adaptively refine reduced-order models *a posteriori* without requiring additional full-order-model solves. The technique is analogous to mesh-adaptive *h*-refinement: it enriches the reduced-basis space online by ‘splitting’ a given basis vector into several vectors with disjoint support. The splitting scheme is defined by a tree structure constructed offline via recursive *k*-means clustering of the state variables using snapshot data. The method identifies the vectors to split online using a dual-weighted-residual approach that aims to reduce error in an output quantity of interest. The resulting method generates a hierarchy of subspaces online without requiring large-scale operations or full-order-model solves. Further, it enables the reduced-order model to satisfy *any prescribed error tolerance* regardless of its original fidelity, as a completely refined reduced-order model is mathematically equivalent to the original full-order model. Experiments on a parameterized inviscid Burgers equation highlight the ability of the method to capture phenomena (e.g., moving shocks) not contained in the span of the original reduced basis. Copyright © 2014 John Wiley & Sons, Ltd.

In this work, we propose a method to prescribe essential boundary conditions in the finite element approximation of elliptic problems when the boundary of the computational domain does not match with the element boundaries. The problems considered are the Poisson problem, the Stokes problem, and the Darcy problem, the latter both in the primal and in the dual formulation. The formulation proposed is of variational type. The key idea is to start with the variational form that defines the problem and treat the boundary condition as a constraint. The particular feature is that the Lagrange multiplier is not defined on the boundary where the essential condition needs to be prescribed but is taken as a certain trace of a field defined in the computational domain, either in all of it or just in a region surrounding the boundary. When approximated numerically, this may allow one to condense the DOFs of the new field and end up with a problem posed only in terms of the original unknowns. The nature of the field used to weakly impose boundary conditions depends on the problem being treated. For the Poisson problem, it is a flux; for the Stokes problem, a stress; for the Darcy problem in primal form, a velocity field; and for the Darcy problem in dual form, it is a potential. Copyright © 2014 John Wiley & Sons, Ltd.

This paper consists of an extension of simulation with direct estimation of stress intensity factors to the three-dimensional case. Here, it combines X-FEM with localized multigrids and direct estimation of quantities of interest along the crack front (SIF, *T*-stress, etc.) based on crack tip asymptotic series expansion. In practice, a three-dimensional patch is introduced locally with a truncated basis of Williams series expansion and is linked in a weak sense with the X-FEM localized multigrids. Some examples (with available analytical solutions) illustrate the efficiency and the robustness of the method. These examples consider planar cracks with curved front, but the proposed method aims to apply to any continuously curved crack. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents a continuum formulation based on the theory of porous media for the mechanics of liquid unsaturated porous media. The hysteresis of the liquid retention model is carefully modelled, including the derivation of the corresponding consistent tangent moduli. The quadratic convergence of Newton's method for solving the highly nonlinear system with an implicit finite element code is demonstrated. A u-p formulation is proposed where the time discretisation is carried out prior to the space discretisation. In this way, the derivation of all consistent moduli is fairly straightforward. Time integration is approximated with the Theta and Newmark's methods, and hence the fully coupled nonlinear dynamics of porous media is considered. It is shown that the liquid retention model requires also the consistent second-order derivative for quadratic convergence. Some predictive simulations are presented illustrating the capabilities of the formulation, in particular to the modelling of complex porous media behaviour. Copyright © 2014 John Wiley & Sons, Ltd.

The method of weighted residuals can efficiently enforce time-periodic solutions of flexible structures experiencing unilateral contact. The harmonic balance method (HBM) based on Fourier expansion of the sought solution is a common formulation, although wavelet bases that can sparsely define nonsmooth solutions may be superior. This hypothesis is investigated using a full three-dimensional blade with unilateral contact conditions on a set of *N*_{c} discrete contact points located at its tip. The unilateral contact conditions are first regularized, and a distributional formulation in time is introduced, allowing trial functions to properly approximate in the time domain the solution to the governing equations. The mixed wavelet Petrov–Galerkin solutions are found to yield consistent or better results than HBM, with higher convergence rates and seemingly more accurate contact force prediction. Copyright © 2014 John Wiley & Sons, Ltd.

A general technique is proposed to improve the stability property of a structure-dependent integration method, which is very computationally efficient in solving inertia-type problems when compared with conventional integration methods, by introducing a stability amplification factor into the coefficient matrices of difference equations. As a result, an improved structure-dependent integration method with a better stability property can be achieved. The concept, derivation, and validation of this technique are intensively studied and are presented in this work. It is evident that the technique can be applied to any structure-dependent integration method to enhance its stability property. Copyright © 2014 John Wiley & Sons, Ltd.

Classical explicit finite element formulations rely on lumped mass matrices. A diagonalized mass matrix enables a trivial computation of the acceleration vector from the force vector. Recently, non-diagonal mass matrices for explicit finite element analysis (FEA) have received attention due to the selective mass scaling (SMS) technique. SMS allows larger time step sizes without substantial loss of accuracy. However, an expensive solution for accelerations is required at each time step. In the present study, this problem is solved by directly constructing the inverse mass matrix. First, a consistent and sparse inverse mass matrix is built from the modified Hamiltons principle with independent displacement and momentum variables. Usage of biorthogonal bases for momentum allows elimination of momentum unknowns without matrix inversions and directly yields the inverse mass matrix denoted here as reciprocal mass matrix (RMM). Secondly, a variational mass scaling technique is applied to the RMM. It is based on the penalized Hamiltons principle with an additional velocity variable and a free parameter. Using element-wise bases for velocity and a local elimination yields variationally scaled RMM. Thirdly, examples illustrating the efficiency of the proposed method for simplex elements are presented and discussed. Copyright © 2014 John Wiley & Sons, Ltd.

A novel extended variational multiscale method for incompressible two-phase flow is proposed. In this approach, the level-set method, which allows for accurately representing complex interface evolutions, is combined with an extended finite element method for the fluid field. Sharp representation of the discontinuities at the interface related to surface-tension effects and large material–parameter ratios are enabled by this approach. To capture the discontinuities, jump enrichments are applied for both velocity and pressure field. Nitsche's method is then used to weakly impose the continuity of the velocity field. For a stable formulation on the entire domain, residual-based variational multiscale terms are supported by appropriate face-oriented ghost-penalty and fluid stabilization terms in the region of enriched elements. Both face-oriented terms and interfacial terms related to Nitsche's method are introduced such that it is accounted for viscous-dominated and convection-dominated transient flows. As a result, stability and well-conditioned systems are guaranteed independent of the interface location. The proposed method is applied to four numerical examples of increasing complexity: two-dimensional Rayleigh–Taylor instabilities, a two-dimensional collapsing water column, three-dimensional rising bubbles, and a three-dimensional bubble coalescence. Excellent agreement with either analytical solutions or numerical and experimental reference data as well as robustness for all flow regimes is demonstrated for all examples. Copyright © 2014 John Wiley & Sons, Ltd.

For the delamination and matrix crack prediction of composite laminated structures, the methods based on the damage mechanics and fracture mechanics are most commonly used. However, there are very few methods that can accurately simulate the delaminations together with matrix cracks, although the in-plane matrix cracks always exist alongside the delaminations under impact loading. In this work, an extended layerwise method is developed to model the composite laminated beam with multiple delaminations and matrix cracks. In the displacement field, the nodes in the thickness direction are located at the middle surface of each single layer, the top surface and the bottom surface of the composite beams. The displacement field contains the linear Lagrange interpolation functions, the one-dimensional weak discontinuous function and strong discontinuous function. The strong and weak discontinuous function are applied to model the displacement discontinuity induced by delaminations and the strain discontinuity induced by the interface between the layers, respectively. Because the nodes in the thickness direction are located at the middle surface of each single layer, the extended layerwise method can be conveniently employed to deal with the in-plane matrix cracks combined with the extend FEM. Copyright © 2014 John Wiley & Sons, Ltd.

An embedded mesh method using piecewise constant multipliers originally proposed by Puso *et al.* (CMAME, 2012) is analyzed here to determine effects of the pressure stabilization term and small cut cells. The approach is implemented for transient dynamics using the central difference scheme for the time discretization. It is shown that the resulting equations of motion are a stable linear system with a condition number independent of mesh size. Next, it is shown that the constraints and the stabilization terms can be recast as non-proportional damping such that the time integration of the scheme is provably stable with a critical time step computed from the undamped equations of motion. Effects of small cuts are discussed throughout the presentation. A mesh study is conducted to evaluate the effects of the stabilization on the discretization error and conditioning and is used to recommend an optimal value for stabilization scaling parameter. Several nonlinear problems are also analyzed and compared with comparable conforming mesh results. Finally, several demanding problems highlighting the robustness of the proposed approach are shown. Copyright © 2014 John Wiley & Sons, Ltd.

A new computational contact formulation is presented and analyzed for large deformation frictional contact. The new formulation uses an unbiased treatment of the two neighboring contact surfaces considering the two-half-pass contact algorithm, originally derived for frictionless contact. The presented work thus introduces several novelties to unbiased friction algorithms. The new algorithm does not enforce traction continuity at the contact interface explicitly but rather satisfies it intrinsically to high accuracy, as is shown. A new 3D friction formulation is also proposed, which is a direct extension of the 1D setup, expressing the friction variables in the parameter space used for the curvilinear surface description. The new formulation resorts to classical expressions in the continuum limit. The current approach uses *C*^{1}-smooth contact surface representations based on either Hermite or non-uniform rational B-spline interpolation. A penalty regularization is considered for the impenetrability and tangential sticking constraints. The new, unbiased friction formulation is illustrated by several 2D and 3D examples, which include an extensive analysis of the model parameters, a convergence study, and the comparison with a classical biased master/slave contact algorithm. Copyright © 2014 John Wiley & Sons, Ltd.

In the recent paper, Fish and Kuznetsov introduced the so-called computational continua (C^{2}) approach, which is a variant of the higher order computational homogenization that does not require higher order continuity, introduces no new degrees of freedom, and is free of higher order boundary conditions. In a follow-up paper on reduced order computational continua, the C^{2} formulation has been enhanced with a model reduction scheme based on construction of residual-free fields to yield a computationally efficient framework coined as RC^{2}. The original C^{2} formulations were limited to rectangular and box elements. The primary objectives of the present manuscript is to revisit the original formulation in three respects: (i) consistent formulation of boundary conditions for unit cells subjected to higher order coarse scale fields, (ii) effective solution of the unit cell problem for lower order approximation of eigenstrains, and (iii) development of nonlocal quadrature schemes for general two-dimensional (quad and triangle) and three-dimensional (hexahedral and tetrahedral) elements. Copyright © 2014 John Wiley & Sons, Ltd.

This work presents a simple technique for real-time monitoring of thermal processes. Real-time simulation-based control of thermal processes is a big challenge because high-fidelity numerical simulations are costly and cannot be used, in general, for real-time decision making. Very often, processes are monitored or controlled with a few measurements at some specific points. Thus, the strategy presented here is centered on fast evaluation of the response only where it is needed. To accomplish this, classical harmonic analysis is combined with recent model reduction techniques. This leads to an advanced harmonic methodology, which solves in real time the transient heat equation at the monitored point.

In order to apply the *reciprocity principle*, harmonic analysis is used in the space-frequency domain. Then, Proper Generalized Decomposition, a reduced order approach, pre-computes a transfer function able to produce the output response for a given excitation. This transfer function is computed offline and only once. The response at the monitoring point can be recovered performing a computationally inexpensive post-processing step. This last step can be performed online for real-time monitoring of the thermal process. Examples show the applicability of this approach for a wide range of problems ranging from fast temperature evaluation to inverse problems. Copyright © 2014 John Wiley & Sons, Ltd.

We demonstrate the potential of collocation methods for efficient higher-order analysis on standard nodal finite element meshes. We focus on a collocation method that is variationally consistent and geometrically flexible, converges optimally, embraces concepts of reduced quadrature, and leads to symmetric stiffness and diagonal consistent mass matrices. At the same time, it minimizes the evaluation cost per quadrature point, thus reducing formation and assembly effort significantly with respect to standard Galerkin finite element methods. We provide a detailed review of all components of the technology in the context of elastodynamics, that is, weighted residual formulation, nodal basis functions on Gauss–Lobatto quadrature points, and symmetrization by averaging with the ultra-weak formulation. We quantify potential gains by comparing the computational efficiency of collocated and standard finite elements in terms of basic operation counts and timings. Our results show that collocation is significantly less expensive for problems dominated by the formation and assembly effort, such as higher-order elastostatic analysis. Furthermore, we illustrate the potential of collocation for efficient higher-order explicit dynamics. Throughout this work, we advocate a straightforward implementation based on simple modifications of standard finite element codes. We also point out the close connection to spectral element methods, where many of the key ideas are already established. Copyright © 2014 John Wiley & Sons, Ltd.

We present the first a priori error analysis for the first hybridizable discontinuous Galerkin method for linear elasticity proposed in Internat. J. Numer. Methods Engrg. **80** (2009), no. 8, 1058–1092. We consider meshes made of polyhedral, shape-regular elements of arbitrary shape and show that, whenever piecewise-polynomial approximations of degree are used and the exact solution is smooth enough, the antisymmetric part of the gradient of the displacement converges with order *k*, the stress and the symmetric part of the gradient of the displacement converge with order *k* + 1/2, and the displacement converges with order *k* + 1. We also provide numerical results showing that the orders of convergence are actually sharp. Copyright © 2014 John Wiley & Sons, Ltd.

The optimization of a vibro-acoustic problem is of main importance for passengers' comfort in transportation vehicles in terms of interior noise. Engineers use numerical tools to predict the response of this coupled problem, but it may lead to a prohibitive computational time. Based on FEM, this work aims at reducing the computational time. The first idea is to keep the same mesh of the acoustic cavity for all the structure configurations and to enrich the pressure approximation by using the extended FEM (XFEM). The enrichment is based on a Heaviside function completed at the structure tip by a continuous ramp function. The second idea is to build reduced basis. The structure basis is composed of its eigenmodes, whereas a modal synthesis method with a fixed interface is used to build the fluid basis. The interface DOFs are the enriched DOF of the XFEM, whereas the internal domain corresponds to the acoustic cavity with no structure. These two combined ideas enable to minimize the computational time in the study of the influence of the structure positions in an acoustic cavity. The method is implemented for shell structures embedded in a 3D acoustic domain. Copyright © 2014 John Wiley & Sons, Ltd.

Topology optimization of large scale structures is computationally expensive, notably because of the cost of solving the equilibrium equations at each iteration. Reduced order models by projection, also known as reduced basis models, have been proposed in the past for alleviating this cost. We propose here a new method for coupling reduced basis models with topology optimization to improve the efficiency of topology optimization of large scale structures. The novel approach is based on constructing the reduced basis on the fly, using previously calculated solutions of the equilibrium equations. The reduced basis is thus adaptively constructed and enriched, based on the convergence behavior of the topology optimization. A direct approach and an approach with adjusted sensitivities are described, and their algorithms provided. The approaches are tested and compared on various 2D and 3D minimum compliance topology optimization benchmark problems. Computational cost savings by up to a factor of 12 are demonstrated using the proposed methods. Copyright © 2014 John Wiley & Sons, Ltd.

In this work, we propose a novel model order reduction approach for two-phase flow in porous media by introducing a global pressure formulation in which the total mobility, which realizes the coupling between saturation and pressure, is regarded as a parameter to the pressure equation. In our approach, the mobility itself is an optimal fit of mobility profiles that are precomputed using the time-of-flight for the initial saturation. Using this formulation and the localized reduced basis multiscale method, we obtain a low-dimensional surrogate of the high-dimensional pressure equation. By applying ideas from model order reduction for parametrized partial differential equations, we are able to split the computational effort for solving the pressure equation into a costly *offline* step that is performed only once and an inexpensive *online* step that is carried out in every time step of the two-phase flow simulation, which is thereby largely accelerated. Usage of elements from numerical multiscale methods allows us to displace the computational intensity between the offline and online steps to reach an ideal runtime at acceptable error increase for the two-phase flow simulation. This is achieved by constructing reduced-dimensional local spaces that lead to a non-conforming global approximation space. As one example for a coupling on the global space, we introduce a discontinuous Galerkin scheme. Copyright © 2014 John Wiley & Sons, Ltd.

An adaptive and efficient approach for constructing reduced-order models (ROMs) that are robust to changes in parameters is developed. The approach is based on a greedy sampling of the underlying high-dimensional model (HDM) together with an efficient procedure for exploring the configuration space and identifying parameters for which the error is likely to be high. Because this exploration is based on a surrogate model for an error indicator, it is amenable to a fast training phase. Furthermore, a model for the exact error based on the error indicator is constructed and used to determine when the greedy procedure reaches a desired error tolerance. An efficient procedure for updating the reduced-order basis constructed by proper orthogonal decomposition is also introduced in this paper, avoiding the cost associated with computing large-scale singular value decompositions. The proposed procedure is then shown to require less evaluations of *a posteriori* error estimators than the classical procedure in order to identify locations of the parameter space to be sampled. It is illustrated on the training of parametric ROMs for three linear and nonlinear mechanical systems, including the realistic prediction of the response of a V-hull vehicle to underbody blasts. 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.

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.

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.

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.

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

We are concerned with the numerical simulation of wave motion in arbitrarily heterogeneous, elastic, perfectly-matched-layer-(PML)-truncated media. We extend in three dimensions a recently developed two-dimensional formulation, by treating the PML via an unsplit-field, but mixed-field, displacement-stress formulation, which is then coupled to a standard displacement-only formulation for the interior domain, thus leading to a computationally cost-efficient hybrid scheme. The hybrid treatment leads to, at most, third-order in time semi-discrete forms. The formulation is flexible enough to accommodate the standard PML, as well as the multi-axial PML.

We discuss several time-marching schemes, which can be used à la carte, depending on the application: (a) an extended Newmark scheme for third-order in time, either unsymmetric or fully symmetric semi-discrete forms; (b) a standard implicit Newmark for the second-order, unsymmetric semi-discrete forms; and (c) an explicit Runge–Kutta scheme for a first-order in time unsymmetric system. The latter is well-suited for large-scale problems on parallel architectures, while the second-order treatment is particularly attractive for ready incorporation in existing codes written originally for finite domains.

We compare the schemes and report numerical results demonstrating stability and efficacy of the proposed formulations. Copyright © 2014 John Wiley & Sons, Ltd.

We introduce a new class of explicit coupling schemes for the numerical solution of fluid-structure interaction problems involving a viscous incompressible fluid and an elastic structure. These methods generalize the arguments reported in [*Comput. Methods Appl. Mech. Engrg., 267:566–593, 2013*, *Numer. Math., 123(1):21–65, 2013*] to the case of the coupling with thick-walled structures. The basic idea lies in the derivation of an intrinsic interface Robin consistency at the space semi-discrete level, using a lumped-mass approximation in the structure. The fluid–solid splitting is then performed through appropriate extrapolations of the solid velocity and stress on the interface. Based on these methods, a new, parameter-free, Robin–Neumann iterative procedure is also proposed for the partitioned solution of implicit coupling. A priori energy estimates, guaranteeing the stability of the schemes and the convergence of the iterative procedure, are established within a representative linear setting. The accuracy and performance of the methods are illustrated in several numerical examples. Copyright © 2014 John Wiley & Sons, Ltd.

In this paper, a novel constitutive model combining continuum damage with embedded discontinuity is developed for explicit dynamic analyses of quasi-brittle failure phenomena. The model is capable of describing the rate-dependent behavior in dynamics and the three phases in failure of quasi-brittle materials. The first phase is always linear elastic, followed by the second phase corresponding to fracture-process zone creation, represented with rate-dependent continuum damage with isotropic hardening formulated by utilizing consistency approach. The third and final phase, involving nonlinear softening, is formulated by using an embedded displacement discontinuity model with constant displacement jumps both in normal and tangential directions. The proposed model is capable of describing the rate-dependent ductile to brittle transition typical of cohesive materials (e.g., rocks and ice). The model is implemented in the finite element setting by using the CST elements. The displacement jump vector is solved for implicitly at the local (finite element) level along with a viscoplastic return mapping algorithm, whereas the global equations of motion are solved with explicit time-stepping scheme. The model performance is illustrated by several numerical simulations, including both material point and structural tests. The final validation example concerns the dynamic Brazilian disc test on rock material under plane stress assumption. Copyright © 2014 John Wiley & Sons, Ltd.