Localized Jacobian ILU preconditioners for hydraulic fractures

We discuss the properties of a class of sparse localized approximations to the Jacobian operator that arises in modelling the evolution of a hydraulically driven fracture in a multi‐layered elastic medium. The governing equations involve a highly non‐linear coupled system of integro‐partial differential equations along with the fracture front free boundary problem. We demonstrate that an incomplete LU factorization of these localized Jacobians yields an efficient preconditioner for the fully populated, stiff, non‐symmetric system of algebraic equations that need to be solved multiple times for every growth increment of the fracture. The performance characteristics of this class of preconditioners is explored via spectral analysis and numerical experiment. Copyright © 2005 John Wiley & Sons, Ltd.     

Analysis of two‐grid methods for reaction‐diffusion equations by expanded mixed finite element methods

We present two efficient methods of two‐grid scheme for the approximation of two‐dimensional semi‐linear reaction‐diffusion equations using an expanded mixed finite element method. To linearize the discretized equations, we use two Newton iterations on the fine grid in our methods. Firstly, we solve an original non‐linear problem on the coarse grid. Then we use twice Newton iterations on the fine grid in our first method, and while in second method we make a correction on the coarse grid between two Newton iterations on the fine grid. These two‐grid ideas are from Xu's work (SIAM J. Sci. Comput. 1994; 15 :231–237; SIAM J. Numer. Anal. 1996; 33 :1759–1777) on standard finite element method. We extend the ideas to the mixed finite element method. Moreover, we obtain the error estimates for two algorithms of two‐grid method. It is showed that coarse space can be extremely coarse and we achieve asymptotically optimal approximation as long as the mesh sizes satisfy H =??(h^¼) in the first algorithm and H =??(h^?) in second algorithm. Copyright © 2006 John Wiley & Sons, Ltd.     

Robust and efficient domain decomposition preconditioners for adaptive hp finite element approximations of linear elasticity with and without discontinuous coefficients

Adaptive finite element methods (FEM) generate linear equation systems that require dynamic and irregular patterns of storage, access, and computation, making their parallelization difficult. Additional difficulties are generated for problems in which the coefficients of the governing partial differential equations have large discontinuities. We describe in this paper the development of a set of iterative substructuring based solvers and domain decomposition preconditioners with an algebraic coarse‐grid component that address these difficulties for adaptive hp approximations of linear elasticity with both homogeneous and inhomogeneous material properties. Our solvers are robust and efficient and place no restrictions on the mesh or partitioning. Copyright © 2003 John Wiley & Sons, Ltd.     

A two‐grid method for expanded mixed finite‐element solution of semilinear reaction–diffusion equations

We present a scheme for solving two‐dimensional semilinear reaction–diffusion equations using an expanded mixed finite element method. To linearize the mixed‐method equations, we use a two‐grid algorithm based on the Newton iteration method. The solution of a non‐linear system on the fine space is reduced to the solution of two small (one linear and one non‐linear) systems on the coarse space and a linear system on the fine space. It is shown that the coarse grid can be much coarser than the fine grid and achieve asymptotically optimal approximation as long as the mesh sizes satisfy H=O(h^1/3). As a result, solving such a large class of non‐linear equation will not be much more difficult than solving one single linearized equation. Copyright © 2003 John Wiley & Sons, Ltd.     

Decoupling joint and elastic accelerations in deformable mechanical systems using non‐linear recursive formulations

The aim of this paper is to develop non‐linear recursive formulations for decoupling joint and elastic accelerations, while maintaining the non‐linear inertia coupling between rigid body motion and elastic deformation in deformable mechanical systems. The inertia projection schemes used in most existing recursive formulations for the dynamic analysis of deformable mechanisms lead to dense coefficient matrices in the equations of motion. Consequently, there are strong dynamic couplings between the joint and elastic coordinates. When the number of elastic degrees of freedom increases, the size of the coefficient matrix in the equations of motion becomes large. Consequently, the use of these recursive formulations for solving the joint and elastic accelerations becomes less efficient. In this paper, the non‐linear recursive formulations have been used to decouple the elastic and joint accelerations in deformable mechanical systems. The relationships between the absolute, elastic and joint variables and generalized Newton–Euler equations are used to develop systems of loosely coupled equations that have sparse matrix structure. By using the inertia matrix structure of deformable mechanical systems and the fact that joint reaction forces associated with elastic coordinates do represent independent variables, a reduced system of equations whose dimension is dependent of the number of elastic degrees of freedom is obtained. This system can be solved for the joint accelerations as well as for the joint reaction forces. The use of the approaches developed in this investigation is illustrated using deformable open‐loop serial robot and closed‐loop four‐bar mechanical systems. Copyright © 2007 John Wiley & Sons, Ltd.     

An efficient method for solving implicit and explicit stiff differential equations

When the s‐stage fully implicit Runge–Kutta (RK) method is used to solve a system of n ordinary differential equations (ODE) the resulting algebraic system has a dimension ns. Its solution by Gauss elimination is expensive and requires 2s³n³/3 operations. In this paper we present an efficient algorithm, which differs from the traditional RK method. The formal procedure for uncoupling the algebraic system into a block‐diagonal matrix with s blocks of size n is derived for any s. Its solution is s²/2 times faster than the original, nondiagonalized system, for s even, and s³/(s−1) for s odd in terms of number of multiplications, as well as s² times in terms of number of additions/multiplications. In particular, for s=3 the method has the same precision and stability properties as the well‐known RK‐based RadauIIA quadrature of Ehle, implemented by Hairer and Wanner in RADAU5 algorithm. Unlike RADAU5, however, the method is applicable with any s and not only to the explicit ODEs My′=f(x, y), where M=const., but also to the general implicit ODEs of the form f(x, y, y′)=0. The block‐diagonal form of the algebraic system allows parallel processing. The algorithm formally differs from the implicit RK methods in that the solution for y is assumed to have a form of the algebraic polynomial whose coefficients are found by enforcing y to satisfy the differential equation at the collocation points. Locations of those points are found from the derived stability function such as to guarantee either A‐ or L‐stability properties as well as a superior precision of the algorithm. If constructed such as to be L‐stable the method is a good candidate for solving differential‐algebraic equations (DAEs). Although not limited to any specific field, the application of the method is illustrated by its implementation in the multibody dynamics described by both ODEs and DAEs. Copyright © 2000 John Wiley & Sons, Ltd.     

Hybrid parallel multigrid preconditioner based on automatic mesh coarsening for 3D metal forming simulations

A parallel multigrid (MG) method is developed to reduce the large computational costs involved by the finite element simulation of highly viscous fluid flows, especially those resulting from metal forming applications, which are characterized by using a mixed velocity/pressure implicit formulation, unstructured meshes of tetrahedra, and frequent remeshings. The developed MG method follows a hybrid approach where the different levels of nonnested meshes are geometrically constructed by mesh coarsening, while the linear systems of the intermediate levels result from the Galerkin algebraic approach. A linear O(N) convergence rate is expected (with N being the number of unknowns), while keeping software parallel efficiency. These objectives lead to selecting unusual MG smoothers (iterative solvers) for the upper grid levels and to developing parallel mesh coarsening algorithms along with parallel transfer operators between the different levels of partitioned meshes. Within the utilized PETSc library, the developed MG method is employed as a preconditioner for the usual conjugate residual algorithm because of the symmetric undefinite matrix of the system to solve. It shows a convergence rate close to optimal, an excellent parallel efficiency, and the ability to handle the complex forming problems encountered in 3‐dimensional hot forging, which involve large material deformations and frequent remeshings.     

Algebraic grid generation on trimmed parametric surface using non‐self‐overlapping planar Coons patch

Using a Coons patch mapping to generate a structured grid in the parametric region of a trimmed surface can avoid the singularity of elliptic PDE methods when only C¹ continuous boundary is given; the error of converting generic parametric C¹ boundary curves into a specified representation form is also avoided. However, overlap may happen on some portions of the algebraically generated grid when a linear or naïve cubic blending function is used in the mapping; this severely limits its usage in most of engineering and scientific applications where a grid system of non‐self‐overlapping is strictly required. To solve the problem, non‐trivial blending functions in a Coons patch mapping should be determined adaptively by the given boundary so that self‐overlapping can be averted. We address the adaptive determination problem by a functional optimization method. The governing equation of the optimization is derived by adding a virtual dimension in the parametric space of the given trimmed surface. Both one‐ and two‐parameter blending functions are studied. To resolve the difficulty of guessing good initial blending functions for the conjugate gradient method used, a progressive optimization algorithm is then proposed which has been shown to be very effective in a variety of practical examples. Also, an extension is added to the objective function to control the element shape. Finally, experiment results are shown to illustrate the usefulness and effectiveness of the presented method. Copyright © 2004 John Wiley & Sons, Ltd.     

Fast multipole method as an efficient solver for 2D elastic wave surface integral equations

The fast multipole method (FMM) is very efficient in solving integral equations. This paper applies the method to solve large solid-solid boundary integral equations for elastic waves in two dimensions. The scattering problem is first formulated with the boundary element method. FMM is then introduced to expedite the solution process. By using the FMM technique, the number of floating-point operations of the matrix-vector multiplication in a standard conjugate gradient algorithm is reduced from O(N ²) to O(N ^1.5), where N is the number of unknowns. The matrix-filling time and the memory requirement are also of the order N ^1.5. The computational complexity of the algorithm is further reduced to O(N ^4/3) by using a ray propagation technique. Numerical results are given to show the accuracy and efficiency of FMM compared to the boundary element method with dense matrix.     

A two‐dimensional ordinary,state‐based peridynamic model for linearly elastic solids

Peridynamics is a non‐local mechanics theory that uses integral equations to include discontinuities directly in the constitutive equations. A three‐dimensional, state‐based peridynamics model has been developed previously for linearly elastic solids with a customizable Poisson's ratio. For plane stress and plane strain conditions, however, a two‐dimensional model is more efficient computationally. Here, such a two‐dimensional state‐based peridynamics model is presented. For verification, a 2D rectangular plate with a round hole in the middle is simulated under constant tensile stress. Dynamic relaxation and energy minimization methods are used to find the steady‐state solution. The model shows m‐convergence and δ‐convergence behaviors when m increases and δ decreases. Simulation results show a close quantitative matching of the displacement and stress obtained from the 2D peridynamics and a finite element model used for comparison. Copyright © 2014 John Wiley & Sons, Ltd.     

A numerical study of a semi‐algebraic multilevel preconditioner for the local discontinuous Galerkin method

In this work, we consider the local discontinuous Galerkin (LDG) method applied to second‐order elliptic problems arising in the modeling of single‐phase flows in porous media. It has been recently proven that the spectral condition number of the stiffness matrix exhibits an asymptotic behavior of ??(h^?2) on structured and unstructured meshes, where h is the mesh size. Thus, efficient preconditioners are mandatory. We present a semi‐algebraic multilevel preconditioner for the LDG method using local Lagrange‐type interpolatory basis functions. We show, numerically, that its performance does not degrade, or at least the number of iterations increases very slowly, as the number of unknowns augments. The preconditioner is tested on problems with high jumps in the coefficients, which is the typical scenario of problems arising in porous media. Copyright © 2007 John Wiley & Sons, Ltd.     

Solution of low p'eclet number heat transfer with viscous dissipation in the infinite axial region between parallel plates via functional analytic methods

Abstract

Functional analytic methods have been applied to the analysis of the extended Graetz problem with prescribed wall flux and viscous dissipation between parallel plates. First, the non‐self‐adjoint elliptic energy equation is decomposed into a set of first order partial differential equations to obtain a self‐adjoint formulism. Next, the induced eigenvalue problems are solved by applying an approximation method in a Hilbert space, and an algebraic characteristic equation is obtained. In addition, the expansion coefficients of the solutions on upstream and downstream regions can be explicitly obtained and unnecessary to match at the entrance.     

Dynamic modelling of the flat 2‐D crack by a semi‐analytic BIEM scheme

We present an efficient numerical method for solving indirect boundary integral equations that describe the dynamics of a flat two‐dimensional (2‐D) crack in all modes of fracture. The method is based on a piecewise‐constant interpolation, both in space and time, of the slip‐rate function, by which the original equation is reduced to a discrete convolution, in space and time, of the slip‐rate and a set of analytically obtained coefficients. If the time‐step interval is set sufficiently small with respect to the spatial grid size, the discrete equations decouple and can be solved explicitly. This semi‐analytic scheme can be extended to the calculation of the wave field off the crack plane. A necessary condition for the numerical stability of this scheme is investigated by way of an exhaustive set of trial runs for a kinematic problem. For the case investigated, our scheme is very stable for a fairly wide range of control parameters in modes III and I, whereas, in mode II, it is unstable except for some limited ranges of the parameters. The use of Peirce and Siebrits' ε‐scheme in time collocation is found helpful in stabilizing the numerical calculation. Our scheme also allows for variable time steps. Copyright © 2001 John Wiley & Sons, Ltd.     

A comparison between two solution techniques to solve the equations of glacially induced deformation of an elastic Earth

In this paper we compare two modeling frameworks within which we compute the Earth's elastic response to an external load. The main difference between the two frameworks lies in their capability to solve the partial differential equations (PDE) for a pre‐stressed elastic body, and especially how the term describing the pre‐stress advection is represented. Furthermore, the two frameworks differ in their applicability in the incompressible limit, the choice of the finite elements (FEs) used for the discretization, and the solution strategy for the arising algebraic problem. The two frameworks are denoted as Framework I and II. Framework I is confined to the features of a commercial FE package (ABAQUS), whereas Framework II is based on a code developed by the authors. The efficiency and accuracy of the two frameworks are compared via theoretical analysis and numerical experiments for two model problems in 2D and 3D. We find that the preconditioned iterative solution method used in Framework II is more efficient than the direct solver used in Framework I. Furthermore, we find that the solutions obtained from the two frameworks differ significantly for compressible solids. This is due to how the pre‐stress advection term is handled. Copyright © 2008 John Wiley & Sons, Ltd.     

A mixed finite element formulation of triphasic mechano‐electrochemical theory for charged,hydrated biological soft tissues

An equivalent new expression of the triphasic mechano‐electrochemical theory [9] is presented and a mixed finite element formulation is developed using the standard Galerkin weighted residual method. Solid displacement u ^s, modified electrochemical/chemical potentials ϵ^w, ϵ⁺and ϵ⁻ (with dimensions of concentration) for water, cation and anion are chosen as the four primary degrees of freedom (DOFs) and are independently interpolated. The modified Newton–Raphson iterative procedure is employed to handle the non‐linear terms. The resulting first‐order Ordinary Differential Equations (ODEs) with respect to time are solved using the implicit Euler backward scheme which is unconditionally stable. One‐dimensional (1‐D) linear isoparametric element is developed. The final algebraic equations form a non‐symmetric but sparse matrix system. With the current choice of primary DOFs, the formulation has the advantage of small amount of storage, and the jump conditions between elements and across the interface boundary are satisfied automatically. The finite element formulation has been used to investigate a 1‐D triphasic stress relaxation problem in the confined compression configuration and a 1‐D triphasic free swelling problem. The formulation accuracy and convergence for 1‐D cases are examined with independent finite difference methods. The FEM results are in excellent agreement with those obtained from the other methods. Copyright © 1999 John Wiley & Sons, Ltd.     

Elasto‐plastic finite element analysis of shells with damage due to microvoids

This paper presents a non‐linear finite element analysis for the elasto‐plastic behaviour of thick/thin shells and plates with large rotations and damage effects. The refined shell theory given by Voyiadjis and Woelke (Int. J. Solids Struct. 2004; 41 :3747–3769) provides a set of shell constitutive equations. Numerical implementation of the shell theory leading to the development of the C⁰ quadrilateral shell element (Woelke and Voyiadjis, Shell element based on the refined theory for thick spherical shells. 2006, submitted) is used here as an effective tool for a linear elastic analysis of shells. The large rotation elasto‐plastic model for shells presented by Voyiadjis and Woelke (General non‐linear finite element analysis of thick plates and shells. 2006, submitted) is enhanced here to account for the damage effects due to microvoids, formulated within the framework of a micromechanical damage model. The evolution equation of the scalar porosity parameter as given by Duszek‐Perzyna and Perzyna (Material Instabilities: Theory and Applications, ASME Congress, Chicago, AMD‐Vol. 183/MD‐50, 9–11 November 1994; 59–85) is reduced here to describe the most relevant damage effects for isotropic plates and shells, i.e. the growth of voids as a function of the plastic flow. The anisotropic damage effects, the influence of the microcracks and elastic damage are not considered in this paper. The damage modelled through the evolution of porosity is incorporated directly into the yield function, giving a generalized and convenient loading surface expressed in terms of stress resultants and stress couples. A plastic node method (Comput. Methods Appl. Mech. Eng. 1982; 34 :1089–1104) is used to derive the large rotation, elasto‐plastic‐damage tangent stiffness matrix. Some of the important features of this paper are that the elastic stiffness matrix is derived explicitly, with all the integrals calculated analytically (Woelke and Voyiadjis, Shell element based on the refined theory for thick spherical shells. 2006, submitted). In addition, a non‐layered model is adopted in which integration through the thickness is not necessary. Consequently, the elasto‐plastic‐damage stiffness matrix is also given explicitly and numerical integration is not performed. This makes this model consistent mathematically, accurate for a variety of applications and very inexpensive from the point of view of computer power and time. Copyright © 2006 John Wiley & Sons, Ltd.     

A stream function implicit finite difference scheme for 2D incompressible flows of Newtonian fluids

The development of a new algorithm to solve the Navier–Stokes equations by an implicit formulation for the finite difference method is presented, that can be used to solve two‐dimensional incompressible flows by formulating the problem in terms of only one variable, the stream function. Two algebraic equations with 11 unknowns are obtained from the discretized mathematical model through the ADI method. An original algorithm is developed which allows a reduction from the original 11 unknowns to five and the use of the Pentadiagonal Matrix Algorithm (PDMA) in each one of the equations. An iterative cycle of calculations is implemented to assess the accuracy and speed of convergence of the algorithm. The relaxation parameter required is analytically obtained in terms of the size of the grid and the value of the Reynolds number by imposing the diagonal dominancy condition in the resulting pentadiagonal matrixes. The algorithm developed is tested by solving two classical steady fluid mechanics problems: cavity‐driven flow with Re=100, 400 and 1000 and flow in a sudden expansion with expansion ratio H/h=2 and Re=50, 100 and 200. The results obtained for the stream function are compared with values obtained by different available numerical methods, to evaluate the accuracy and the CPU time required by the proposed algorithm. Copyright © 2002 John Wiley & Sons, Ltd.     

A finite‐volume method for deformation analysis of woven fabrics

Efficient numerical methods for simulating cloth deformations have been identified as the key to the development of successful Computer‐Aided Design systems for clothing products. This paper presents the formulation of a new finite‐volume method for the simulation of complex deformations of initially flat woven fabric sheets under self‐weight or externally applied loading. The fabric sheet is assumed to undergo very large displacements and rotations but small strains during the process of deformation. The fabric material is assumed to be linear elastic and orthotropic. The fabric sheet is discretized into many small structured patches called finite volumes (or control volumes), each containing one grid node and several face nodes. The bending and membrane deformations of a typical volume can be defined using the global co‐ordinates of its grid node and surrounding face nodes. The equilibrium equations governing the complex deformations are derived employing the principle of stationary total potential energy. These equations are solved using a single‐step full Newton–Raphson method which is found to be capable of predicting the final deformed shape, the only result of interest in a fabric drape analysis. To speed up convergence, the line search technique is adopted with good effect. This single‐step approach is more efficient than the step‐by‐step incremental approach employed in conventional non‐linear finite element analysis of load‐bearing structures. A number of example simulations of fabric drape/buckling deformations are included in the paper, which demonstrate the efficiency and validity of the proposed method. Copyright © 1999 John Wiley & Sons, Ltd.     

Solution of low péclet number heat transfer with viscous dissipation in the semi‐infinite axial region of a tube via functional analytic methods

Abstract

A solution of the extended Graetz problem with prescribed wall heat flux and viscous dissipation in a semi‐infinite axial region of a tube is obtained by functional analytic methods. The energy equation is split into a set of partial differential equations to obtain a self‐adjoint formulism. Then, an algebraic characteristic equation of the eigenvalue problem for an arbitrary velocity profile is obtained by an approximation method in L ₂[0, 1]. In addition, a backward recursive formula for calculating the expansion coefficients of the solution is developed.