Three‐dimensional superconvergent gradient recovery on tetrahedral meshes

In this paper, finite element superconvergence phenomenon based on centroidal Voronoi Delaunay tessellations (CVDT) in three‐dimensional space is investigated. The Laplacian operator with the Dirichlet boundary condition is considered. A modified superconvergence patch recovery (MSPR) method is established to overcome the influence of slivers on CVDT meshes. With these two key preconditions, a CVDT mesh and the MSPR, the gradients recovered from the linear finite element solutions have superconvergence in the l₂ norm at nodes of a CVDT mesh for an arbitrary three‐dimensional bounded domain. Numerous numerical examples are presented to demonstrate this superconvergence property and good performance of the MSPR method. Copyright © 2016 John Wiley & Sons, Ltd.     

Robust and efficient monolithic fluid‐structure‐interaction solvers

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

An interpolation‐based fast multipole method for higher‐order boundary elements on parametric surfaces

In this article, a black‐box higher‐order fast multipole method for solving boundary integral equations on parametric surfaces in three spatial dimensions is proposed. Such piecewise smooth surfaces are the topic of recent studies in isogeometric analysis. Due to the exact surface representation, the rate of convergence of higher‐order methods is not limited by approximation errors of the surface. An element‐wise clustering strategy yields a balanced cluster tree and an efficient numerical integration scheme for the underlying Galerkin method. By performing the interpolation for the fast multipole method directly on the reference domain, the cost complexity in the polynomial degree is reduced by one order. This gain is independent of the application of either ‐ or ‐matrices. In fact, several simplifications in the construction of ‐matrices are pointed out, which are a by‐product of the surface representation. Extensive numerical examples are provided in order to quantify and qualify the proposed method. Copyright © 2016 John Wiley & Sons, Ltd.     

A stable nonconforming finite element on hexahedra

A new nonconforming brick element is introduced, which only has 13 DOFs locally and takes as its shape functions space. The vector‐valued version generates, together with a discontinuous approximation, an inf‐sup stable finite element pair of order 2 for the Stokes problem in the energy norm. Copyright © 2016 John Wiley & Sons, Ltd.     

Dispersion and isogeometric analyses of second‐order and fourth‐order implicit gradient‐enhanced plasticity models

Implicit gradient plasticity models incorporate higher‐order spatial gradients via an additional Helmholtz type equation for the plastic multiplier. So far, the enrichment has been limited to second‐order spatial gradients, resulting in a formulation that can be discretised using ‐continuous finite elements. Herein, an implicit gradient plasticity model is formulated that includes a fourth‐order gradient term as well. A comparison between the localisation properties of both the implicit gradient plasticity formulations and the explicit second‐order gradient plasticity model is made using a dispersion analysis. The higher‐order continuity requirement for the fourth‐order implicit gradient plasticity model has been met by exploiting the higher‐order continuity property of isogeometric analysis, which uses nonuniform rational B‐splines as shape functions instead of Lagrange polynomials. The discretised variables, displacements, and plastic multiplier may require different orders of interpolation, an issue that is also addressed. Numerical results show that both formulations can be used as a localisation limiter, but that quantitative differences occur, and a different evolution of the localisation band is obtained for 2‐dimensional problems.     

Powell–Sabin B‐splines and unstructured standard T‐splines for the solution of the Kirchhoff–Love plate theory exploiting Bézier extraction

The equations that govern Kirchhoff–Love plate theory are solved using quadratic Powell–Sabin B‐splines and unstructured standard T‐splines. Bézier extraction is exploited to make the formulation computationally efficient. Because quadratic Powell–Sabin B‐splines result in ‐continuous shape functions, they are of sufficiently high continuity to capture Kirchhoff–Love plate theory when cast in a weak form. Unlike non‐uniform rational B‐splines (NURBS), which are commonly used in isogeometric analysis, Powell–Sabin B‐splines do not necessarily capture the geometry exactly. However, the fact that they are defined on triangles instead of on quadrilaterals increases their flexibility in meshing and can make them competitive with respect to NURBS, as no bending strip method for joined NURBS patches is needed. This paper further illustrates how unstructured T‐splines can be modified such that they are ‐continuous around extraordinary points, and that the blending functions fulfil the partition of unity property. The performance of quadratic NURBS, unstructured T‐splines, Powell–Sabin B‐splines and NURBS‐to‐NURPS (non‐uniform rational Powell–Sabin B‐splines, which are obtained by a transformation from a NURBS patch) is compared in a study of a circular plate. Copyright © 2015 John Wiley & Sons, Ltd.     

An isogeometric analysis approach to gradient‐dependent plasticity

Gradient‐dependent plasticity can be used to achieve mesh‐objective results upon loss of well‐posedness of the initial/boundary value problem because of the introduction of strain softening, non‐associated flow, and geometric nonlinearity. A prominent class of gradient plasticity models considers a dependence of the yield strength on the Laplacian of the hardening parameter, usually an invariant of the plastic strain tensor. This inclusion causes the consistency condition to become a partial differential equation, in addition to the momentum balance. At the internal moving boundary, one has to impose appropriate boundary conditions on the hardening parameter or, equivalently, on the plastic multiplier. This internal boundary condition can be enforced without tracking the elastic‐plastic boundary by requiring ‐continuity with respect to the plastic multiplier. In this contribution, this continuity has been achieved by using nonuniform rational B‐splines as shape functions both for the plastic multiplier and for the displacements. One advantage of this isogeometric analysis approach is that the displacements can be interpolated one order higher, making it consistent with the interpolation of the plastic multiplier. This is different from previous approaches, which have been exploited. The regularising effect of gradient plasticity is shown for 1‐ and 2‐dimensional boundary value problems.     

Establishing stable time‐steps for DEM simulations of non‐collinear planar collisions with linear contact laws

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

On the convergence analysis,stability, and implementation of meshless local radial point interpolation on a class of three‐dimensional wave equations

In this article, the meshless local radial point interpolation method is applied to analyze three space dimensional wave equations of the form subject to given initial and Dirichlet boundary conditions. The main difficulty of the great number of methods in full 3‐D problems is the large computational costs. In meshless local radial point interpolation method, it does not require any background integration cells, so that all integrations are carried out locally over small quadrature domains of regular shapes such as circles or squares in two dimensions and spheres or cubes in three dimensions. The point interpolation method with the help of radial basis functions is proposed to construct shape functions that have Kronecker delta function property. A weak formulation with the Heaviside step function converts the set of governing equations into local integral equations on local subdomains. A two‐step time discretization method is employed to evaluate the time derivatives. This suggests Crank‐Nicolson technique to be applied on the right hand side of the equation. The convergence analysis and stability of the method are fully discussed. Three illustrative examples are presented, and satisfactory agreements are achieved. It is shown theoretically that the proposed method is unconditionally stable for the second example whereas it is not for the first and third ones. Copyright © 2015 John Wiley & Sons, Ltd.     

Virtual and smoothed finite elements: A connection and its application to polygonal/polyhedral finite element methods

We show both theoretically and numerically a connection between the smoothed finite element method (SFEM) and the virtual element method and use this approach to derive stable, cheap and optimally convergent polyhedral FEM. We show that the stiffness matrix computed with one subcell SFEM is identical to the consistency term of the virtual element method, irrespective of the topology of the element, as long as the shape functions vary linearly on the boundary. Using this connection, we propose a new stable approach to strain smoothing for polygonal/polyhedral elements where, instead of using sub‐triangulations, we are able to use one single polygonal/polyhedral subcell for each element while maintaining stability. For a similar number of degrees of freedom, the proposed approach is more accurate than the conventional SFEM with triangular subcells. The time to compute the stiffness matrix scales with the in case of the conventional polygonal FEM, while it scales as in the proposed approach. The accuracy and the convergence properties of the SFEM are studied with a few benchmark problems in 2D and 3D linear elasticity. Copyright © 2015 John Wiley & Sons, Ltd.     

High‐order space‐time finite element schemes for acoustic and viscodynamic wave equations with temporal decoupling

We revisit a method originally introduced by Werder et al. (in Comput. Methods Appl. Mech. Engrg., 190:6685–6708, 2001) for temporally discontinuous Galerkin FEMs applied to a parabolic partial differential equation. In that approach, block systems arise because of the coupling of the spatial systems through inner products of the temporal basis functions. If the spatial finite element space is of dimension D and polynomials of degree r are used in time, the block system has dimension (r + 1)D and is usually regarded as being too large when r > 1. Werder et al. found that the space‐time coupling matrices are diagonalizable over for r ?100, and this means that the time‐coupled computations within a time step can actually be decoupled. By using either continuous Galerkin or spectral element methods in space, we apply this DG‐in‐time methodology, for the first time, to second‐order wave equations including elastodynamics with and without Kelvin–Voigt and Maxwell–Zener viscoelasticity. An example set of numerical results is given to demonstrate the favourable effect on error and computational work of the moderately high‐order (up to degree 7) temporal and spatio‐temporal approximations, and we also touch on an application of this method to an ambitious problem related to the diagnosis of coronary artery disease. Copyright © 2014 The Authors. International Journal for Numerical Methods in Engineering published by John Wiley & Sons Ltd.     

Efficient adaptive integration of functions with sharp gradients and cusps in n‐dimensional parallelepipeds

In this paper, we study the efficient numerical integration of functions with sharp gradients and cusps. An adaptive integration algorithm is presented that systematically improves the accuracy of the integration of a set of functions. The algorithm is based on a divide and conquer strategy and is independent of the location of the sharp gradient or cusp. The error analysis reveals that for a C⁰ function (derivative discontinuity at a point), a rate of convergence of n + 1 is obtained in . Two applications of the adaptive integration scheme are studied. First, we use the adaptive quadratures for the integration of the regularized Heaviside function—a strongly localized function that is used for modeling sharp gradients. Then the adaptive quadratures are employed in the enriched finite element solution of the all‐electron Coulomb problem in crystalline diamond. The source term and enrichment functions of this problem have sharp gradients and cusps at the nuclei. We show that the optimal rate of convergence is obtained with only a marginal increase in the number of integration points with respect to the pure finite element solution with the same number of elements. The adaptive integration scheme is simple, robust, and directly applicable to any generalized finite element method employing enrichments with sharp local variations or cusps in n‐dimensional parallelepiped elements. Copyright © 2012 John Wiley & Sons, Ltd.     

Nonconforming rotated Q1 tetrahedral element with explicit time stepping for elastodynamics

In this paper, we apply a rotated bilinear tetrahedral element to elastodynamics in . This element performs superior to the constant strain element in bending and, unlike the conforming linear strain tetrahedron, allows for row‐sum lumping of the mass matrix. We study the effect of different choices of approximation (pointwise continuity versus edge average continuity) as well as lumping versus consistent mass in the setting of eigenvibrations. We also use the element in combination with the leapfrog method for time domain computations and make numerical comparisons with the constant strain and linear strain tetrahedra. Copyright © 2012 John Wiley & Sons, Ltd.     

Extended finite element method for plastic limit load computation of cracked structures

The extended finite element method is extended to allow computation of the limit load of cracked structures. In the paper, it is demonstrated that the linear elastic tip enrichment basis with and without radial term may be used in the framework of limit analysis, but the six‐function enrichment basis based on the well‐known Hutchinson–Rice–Rosengren asymptotic fields appears to be the best. The discrete kinematic formulation is cast in the form of a second‐order cone problem, which can be solved using highly efficient interior‐point solvers. Finally, the proposed numerical procedure is applied to various benchmark problems, showing that the present results are in good agreement with those in the literature. Copyright © 2015 John Wiley & Sons, Ltd.     

A second‐order decoupled implicit/explicit method of the 3D primitive equations of ocean II: finite element spatial discretization

A fully discrete second‐order decoupled implicit/explicit method is proposed for solving 3D primitive equations of ocean in the case of Dirichlet boundary conditions on the side, where a second‐order decoupled implicit/explicit scheme is used for time discretization, and a finite element method based on the P₁(P₁) ? P₁?P₁(P₁) elements for velocity, pressure and density is used for spatial discretization of these primitive equations. Optimal H¹?L²?H¹ error estimates for numerical solution and an optimal L² error estimate for are established under the convergence condition of 0 < h≤β₁,0 < τ≤β₂, and τ≤β₃h for some positive constants β₁,β₂, and β₃. Furthermore, numerical computations show that the H¹?L²?H¹ convergence rate for numerical solution is of O(h + τ²) and an L² convergence rate for is O(h²+τ²) with the assumed convergence condition, where h is a mesh size and τ is a time step size. More practical calculations are performed as a further validation of the numerical method. Copyright © 2016 John Wiley & Sons, Ltd.     

Two‐sided Grassmann manifold algorithm for optimal model reduction

We consider an optimal model reduction problem for large‐scale dynamical systems. The problem is formulated as a minimization problem over Grassmann manifold with two variables. This formulation allows us to develop a two‐sided Grassmann manifold algorithm, which is numerically efficient and suitable for the reduction of large‐scale systems. The resulting reduced system preserves the stability of the original system. Numerical examples are presented to show that the proposed algorithm is computationally efficient and robust with respect to the selection of initial projection matrices. Copyright © 2015 John Wiley & Sons, Ltd.     

A new multi‐scale dispersive gradient elasticity model with micro‐inertia: Formulation and ‐finite element implementation

Motivated by nano‐scale experimental evidence on the dispersion characteristics of materials with a lattice structure, a new multi‐scale gradient elasticity model is developed. In the framework of gradient elasticity, the simultaneous presence of acceleration and strain gradients has been denoted as dynamic consistency. This model represents an extension of an earlier dynamically consistent model with an additional micro‐inertia contribution to improve the dispersion behaviour. The model can therefore be seen as an enhanced dynamic extension of the Aifantis' 1992 strain‐gradient theory for statics obtained by including two acceleration gradients in addition to the strain gradient. Compared with the previous dynamically consistent model, the additional micro‐inertia term is found to improve the prediction of wave dispersion significantly and, more importantly, requires no extra computational cost. The fourth‐order equations are rewritten in two sets of symmetric second‐order equations so that ‐continuity is sufficient in the finite element implementation. Two sets of unknowns are identified as the microstructural and macrostructural displacements, thus highlighting the multi‐scale nature of the present formulation. The associated energy functionals and variationally consistent boundary conditions are presented, after which the finite element equations are derived. Considerable improvements over previous gradient models are observed as confirmed by two numerical examples. Copyright © 2016 John Wiley & Sons, Ltd.     

Energy conservation during remeshing in the analysis of dynamic fracture

The analysis of (dynamic) fracture often requires multiple changes to the discretisation during crack propagation. The state vector from the previous time step must then be transferred to provide the initial values of the next time step. A novel methodology based on a least-squares fit is proposed for this mapping. The energy balance is taken as a constraint in the mapping, which results in a complete energy preservation. Apart from capturing the physics better, this also has advantages for numerical stability. To further improve the accuracy, Powell-Sabin B-splines, which are based on triangles, have been used for the discretisation. Since continuity of the displacement field holds at crack tips for Powell-Sabin B-splines, the stresses at and around crack tips are captured much more accurately than when using elements with a standard Lagrangian interpolation, or with NURBS and T-splines. The versatility and accuracy of the approach to simulate dynamic crack propagation are assessed in two case studies, featuring mode-I and mixed-mode crack propagation.     

Considerations on a phase-field model for adhesive fracture

A recently proposed phase-field model for cohesive fracture is examined. Previous investigations have shown stress oscillations to occur when using unstructured meshes. It is now shown that the use of nonuniform rational B-splines (NURBS) as basis functions rather than traditional Lagrange polynomials significantly reduces this oscillatory behavior. Moreover, there is no effect on the global structural behavior, as evidenced through load-displacement curves. The phase-field model imposes restrictions on the interpolation order of the NURBS used for the three different fields: displacement, phase field, and crack opening. This holds within the Bézier element, but also at the boundaries, where a reduction to ${𝒞}^0$-continuity yields optimal results. Application to a range of cases, including debonding of a hard fiber embedded in a soft matrix, illustrates the potential of the cohesive phase-field model.