Formulations in rates and increments for elastic-plastic analysis

Minimum principles in velocities, stress rates and plastic strain rates are extended in order to derive formulations for finite increments of displacement, stress and plastic strain fields defining complete numerical methods. Kinematical, statical and mixed principles are developed from a new variational formulation of the elastic-plastic work-hardening constitutive relation. The consequences of this time discretization are discussed independently of any discretization of the continuum. In particular, the incremental formulations derived from extended rate principles account for local elastic unloading and produce stress field approximations complying with equilibrium and plastic admissibility without any additional procedure, at least for piecewise linear yield functions. These properties are not fulfilled when the incremental analysis is based on direct discrete versions of classical rate principles. Finally, FEM approximations are formally introduced and the solution of the resulting finite dimensional quadratic optimization problem is considered.     

Adaptive multivlevel methods in three space dimensions

We consider the approximate solution of self-adjoint elliptic problems in three space dimensions by piecewise linear finite elements with respect to a highly non-uniform tetrahedral mesh which is generated adaptively. The arising linear systems are solved iteratively by the conjugate gradient method provided with a multilevel preconditioner. Here, the accuracy of the iterative solution is coupled with the discretization error. As the performance of hierarchical bases preconditioners deteriorates in three space dimensions, the BPX preconditioner is used, taking special care of an efficient implementation. Reliable a posteriori estimates for the discretization error are derived from a local comparison with the approximation resulting from piecewise quadratic elements. To illustrate the theoretical results, we consider a familiar model problem involving reentrant corners and a real-life problem arising from hyperthermia, a recent clinical method for cancer therapy.     

Transient dynamic crack analysis in piecewise homogeneous, anisotropic and linear elastic composites by a spatial symmetric Galerkin-BEM

Transient elastodynamic analysis of two-dimensional, piecewise homogeneous, anisotropic and linear elastic solids containing interior and interface cracks is presented in this paper. To solve the initial boundary value problem, a spatial symmetric time-domain boundary element method is developed. Stationary cracks subjected to impact loading conditions are considered. Elastodynamic fundamental solutions for homogenous, anisotropic and linear elastic solids are implemented. The piecewise homogeneous, anisotropic and linear elastic solids are modeled by the multi-domain technique. The spatial discretization is performed by a symmetric Galerkin-method, while a collocation method is utilized for the temporal discretization. An explicit time-stepping scheme is obtained for computing the unknown boundary data. Numerical examples are presented and discussed to show the effects of the interface cracks, the material anisotropy, the material combination and the dynamic loading on the dynamic stress intensity factors.     

Film flow of a suspension down an inclined plane

A method is developed for simulating the film flow of a suspension of rigid particles with arbitrary shapes down an inclined plane in the limit of vanishing Reynolds number. The problem is formulated in terms of a system of integral equations of the first and second kind for the free-surface velocity and the traction distribution along the particle surfaces involving the a priori unknown particle linear velocity of translation and angular velocity of rotation about designated centres. The problem statement is completed by introducing scalar constraints that specify the force and torque exerted on the individual particles. A boundary-element method is implemented for solving the governing equations for the case of a two-dimensional periodic suspension. The system of linear equations arising from numerical discretization is solved using a preconditioner based on a particle-cluster iterative method recently developed by Pozrikidis (2000 Engng Analysis Bound. Elem. 25, 19-30). Numerical investigations show that the generalized minimal residual (GMRES) method with this preconditioner is significantly more efficient than the plain GMRES method used routinely in boundary-element implementations. Extensive numerical simulations for solitary particles and random suspensions illustrate the effect of the particle shape, size and aspect ratio in semi-finite shear flow, and the effect of free-surface deformability in film flow.     

A general non‐linear optimization algorithm for lower bound limit analysis

The non‐linear programming problem associated with the discrete lower bound limit analysis problem is treated by means of an algorithm where the need to linearize the yield criteria is avoided. The algorithm is an interior point method and is completely general in the sense that no particular finite element discretization or yield criterion is required. As with interior point methods for linear programming the number of iterations is affected only little by the problem size. Some practical implementation issues are discussed with reference to the special structure of the common lower bound load optimization problem, and finally the efficiency and accuracy of the method is demonstrated by means of examples of plate and slab structures obeying different non‐linear yield criteria. Copyright © 2002 John Wiley & Sons, Ltd.     

Non-linear vibration of frames by the incremental dynamic stiffness method

The dynamic stiffness method is extended to large amplitude free and forced vibrations of frames. When the steady state vibration is concerned, the time variable is replaced by the frequency parameter in the Fourier series sense and the governing partial differential equations are replaced by a set of ordinary differential equations in the spatial variables alone. The frequency-dependent shape functons are generated approximately for the spatial discretization. These shape functions are the exact solutions of a beam element subjected to mono-frequency excitation and constant axial force to minimize the spatial discretization errors. The system of ordinary differential equations is replaced by a system of non-linear algebraic equations with the Fourier coefficients of the nodal displacements as unknowns. The Fourier nodal coefficients are solved by the Newtonian algorithm in an incremental manner. When an approximate solution is available, an improved solution is obtained by solving a system of linear equations with the Fourier nodal increments as unknowns. The method is very suitable for parametric studies. When the excitation frequency is taken as a parameter, the free vibration response of various resonances can be obtained without actually computing the linear natural modes. For regular points along the response curves, the accuracy of the gradient matrix (Jacobian or tangential stiffness matrix) is secondary (cf. the modified Newtonian method). However, at the critical positions such as the turning points at resonances and the branching points at bifurcations, the gradient matrix becomes important. The minimum number of harmonic terms required is governed by the conditions of completeness and balanceability for predicting physically realistic response curves. The evaluations of the newly introduced mixed geometric matrices and their derivatives are given explicitly for the computation of the gradient matrix.     

Finite Deflection of Slender Cantilever with Predefined Load Application Locus using an Incremental Formulation

In this paper, a class of problems involving space constrained loading on thin beams with large deflections is considered. The loading is such that, the locus of the force application point moves along an arbitrarily predefined path, fixed in space. Both linear elastic as well as elastic-perfectly plastic materials are considered. A simplification is realized using the moment-curvature relationship directly. The governing equation obtained is highly non-linear owing to inclusion of both material and geometric non-linearity. A general algorithm is described to solve the governing equation using an incremental formulation coupled with Runge Kutta 4^th order initial value explicit solver. Additionally, the presented method is capable of handling unloading and reverse loading conditions. An example problem where the load application point locus is an inclined straight line is solved to demonstrate the performance of the method. It is found that, the force response due to the inclined locus is stiffer than the vertical locus. This response is akin to dry friction condition on a vertical locus case.     

A comparative study of the behaviour of a direct solver and a preconditioned iterative solver for the equations arising from the discretization by Chebyshev collocation of a second-order partial differential equation on a square

The behaviour is compared of two solvers for the discrete equations arising from the discretization using Chebyshev collocation of a second-order linear partial differential equation on a square. The alternative solvers considered are a direct solver and an iterative solver based on preconditioning with the matrix arising from finite-difference discretization of the governing equation. The total error of the collocation derivatives and the separate contributions from round-off and discretization error are examined. The efficiency of the two solvers is compared. The iterative solver is more efficient than the direct solver on fine grids for equations similar to the Poisson equation, provided that there are Dirichlet boundary conditions on at least three of the sides of the square.     

Normalized explicit approximate inverse preconditioning for solving 3D boundary value problems on uniprocessor and distributed systems

Normalized explicit approximate inverse matrix techniques, based on normalized approximate factorization procedures, for solving sparse linear systems resulting from the finite difference discretization of partial differential equations in three space variables are introduced. Normalized explicit preconditioned conjugate gradient schemes in conjunction with normalized approximate inverse matrix techniques are presented for solving sparse linear systems. The convergence analysis with theoretical estimates on the rate of convergence and computational complexity of the normalized explicit preconditioned conjugate gradient method are also derived. A Parallel Normalized Explicit Preconditioned Conjugate Gradient method for distributed memory systems, using message passing interface (MPI) communication library, is also given along with theoretical estimates on speedups, efficiency and computational complexity. Application of the proposed method on a three‐dimensional boundary value problem is discussed and numerical results are given for uniprocessor and multicomputer systems. Copyright © 2005 John Wiley & Sons, Ltd.     

A co-rotational weak-form quadrature planar beam element for geometric nonlinear static and dynamic analysis

The co-rotational formulation of quadrature planar beam element undergoing large displacement and large rotation is presented. A local frame co-rotates with the differential element and decomposes the motion into a rigid body movement and a strain-producing deformation. General explicit formulations of elemental vectors and matrices, including internal force vector, external force vector, tangent stiffness matrix, and mass matrix, are derived via the numerical integration together with the differential quadrature law. Thus, the element nodes and numerical integration method can be chosen arbitrarily based on the accuracy requirement and problem type. A number of case studies on the static, postbuckling, and dynamic response of beams and frame structures are conducted. The convergence study shows that the co-rotational quadrature element has an exponential rate of convergence and the reduced Gauss integration yield the highest accuracy. It is seen that the proposed co-rotational quadrature beam element is simple in formulations, computationally efficient, and capable of capturing the complex nonlinear behavior of beam and frame structures with high precision.     

Sparse hierarchical solvers with guaranteed convergence

Solving sparse linear systems from discretized partial differential equations is challenging. Direct solvers have, in many cases, quadratic complexity (depending on geometry), while iterative solvers require problem dependent preconditioners to be robust and efficient. Approximate factorization preconditioners such as incomplete LU factorization provide cheap approximations to the system matrix. However, even a highly accurate preconditioner may have deteriorating performance when the condition number of the system matrix increases. By increasing the accuracy on low-frequency errors, we propose a novel hierarchical solver with improved robustness with respect to the condition number of the linear system. This solver retains the linear computational cost and memory footprint of the original algorithm.     

Dynamic response of a slab of elastic-viscoplastic material that exhibits induced plastic anisotropy

Various 2-dimensional problems of the dynamic loading of a slab are solved for a material characterization that is elastic-viscoplastic and exhibits anisotropic work-hardening. The governing constitutive equations are based on a unified formulation which requires neither a yield criterion nor loading or unloading conditions. They include multi-dimensional anisotropic effects induced by the plastic deformation history. The theory also considers plastic compressibility which depends on the extent of the anisotropy. A numerical procedure for solving the equations is developed which incorporates the history dependent anisotropic hardening effects. Cases considered are the dynamic penetration of a slab by a rigid cylindrical indenter, and a distributed force rapidly applied over part of the slab surface. Both conditions of fixed and free rear surfaces of the slab are examined. A uniaxial problem is also considered in which different bases for the anisotropic hardening law are examined.     

Non-linear steady state vibration of frames by finite element method

A finite element method based on the virtual work principle to determine the steady state response of frams in free or forced periodic vibration is introduced. The axial and flexural deformations are coupled by mean of the induced axial force along the element. The spatial discretization of the deformations is achieved by the usual finite element method and the time discretization by Fourier coefficients of the nodal displacements. No unconventional element matrices are needed. After applying the harmonic balance method, a set of non-linear algebraic equations of the Fourier coefficients is obtained. These equations are solved by the Newtonian iteration method in terms of the Fourier coefficient increments. Nodal damping can easily be included by a diagonal damping matrix. The direct numerical determination of the Fourier coefficient increments is difficult owing to the presence of peaks, loops and discontinuities of slope along the amplitude-frequency response curves. Parametric construction of the response curves using the phase difference between the response and excitation is recommended to provide more points during the rapid change of the phase (i.e. at resonance). For undamped natural vibration, the method of selective coefficients adopted. Numerical examples on the Duffing equation, a hinged–hinged beam, a clamped–hinged beam, a ring and a frame are given. For reasonably accurate results, it is shown that the number of finite elements must be sufficient to predict at least the linear mode at the frequency of interest and the number of harmones considered must satisfy the conditions of completeness and balanceability, which are discussed in detail.     

Analysis of a row of elliptical inclusions in a plate using singular integral equations

This paper deals with the interaction problem of a row of elliptical inclusions under uniaxial tension. The body force method is used to formulate the problem as a system of singular integral equations with Cauchy--type and logarithmic singularities, where the unknowns are densities of body forces distributed in infinite plates that have the same elastic constants as those of the matrix and inclusion. In order to satisfy the boundary conditions along the elliptical boundaries, eight kinds of fundamental density functions, proposed in a previous paper, are applied. In the analysis, the number, shape, and position of inclusions are varied systematically; after which the magnitude and position of the maximum stress are examined. For any fixed shape and size of inclusions, the maximum stress is shown to be linear with the reciprocal of the number of inclusions. The present method is found to yield rapidly converging numerical results for various geometrical conditions of inclusions. This revised version was published online in August 2006 with corrections to the Cover Date.     

Numerical shakedown analysis of axisymmetric sandwich shells: An upper bound formulation

Shakedown analysis of axisymmetric elastic–perfectly plastic sandwich shells is performed here using a new upper bound formulation based on a special form of Koiter's theorem concerning piecewise linearized yield surfaces. Starting from finite element techniques and the Tresca sandwich yield condition, shakedown analysis is reduced to a linear programming problem which is solved by a powerful simplex algorithm. Numerical results are given for a number of examples and a comparison is made with a previously computed lower bound formulation.     

High‐order finite volume schemes for the advection–diffusion equation

A high‐order finite volume method based on piecewise interpolant polynomials is proposed to discretize spatially the one‐dimensional and two‐dimensional advection–diffusion equation. Evolution equations for the mean values of each control volume are integrated in time by a classical fourth‐order Runge–Kutta. Since our work focuses on the behaviour of the spatial discretization, the time step is chosen small enough to neglect the time integration error. Two‐dimensional interpolants are built by means of one‐dimensional interpolants. It is shown that when the degree of the one‐dimensional interpolant q is odd, the proper selection of a fixed stencil gives rise to centred schemes of order q+1. In order not to lose precision due to the change of stencil near boundaries, the degree of the interpolants close to boundaries is raised to q+1. Four test cases with small values of diffusion are integrated with high‐order methods. It is shown that the spatial discretization of the advection–diffusion equation with periodic boundary conditions leads to normal discretization matrices, and asymptotic stability must be assured to bound the spatial discretization error. Once the asymptotic stability is assured by means of the spectra of the discretization matrix, the spatial error is of the order of the truncation error. However, it is shown that the discretization of the advection–diffusion equation with arbitrary boundary conditions gives rise to non‐normal matrices. If asymptotic stability is assured, the spatial order of steady solutions is of the order of the truncation error. But, for transient processes, the order of the spatial error is determined by both the truncation error and the norm of the exponential matrix of the spatial discretization. The use of the pseudospectra of the discretization matrix is proposed as a valuable tool to analyse the transient error of the numerical solution. Copyright © 2001 John Wiley & Sons, Ltd.     

On computational procedures for the force method

It is known that the matrix force method has certain advantages over the displacement method for a class of structural problems. It is also known that the force method, when carried out by the conventional Gauss-Jordan procedure, tends to fill in the problem data, making the method unattractive for large size, sparse problems. This poor fill-in property, however, is not necessarily inherent to the method, and the sparsity may be maintained if one uses what we call the Turn-Back LU Procedure. The purpose of this paper is two-fold. First, it is shown that there exist some close relationships between the force method and the least squares problem, and that many existing algebraic procedures to perform the force method can be regarded as applications/extensions of certain well-known matrix factorization schemes for the least squares problem. Secondly, it is demonstrated that these algebraic procedures for the force method can be unified form the matrix factorization viewpoint. Included in this unification is the Turn-Back LU Procedure, which was originally proposed by Topçu in his thesis.⁸ It is explained why this procedure tends to produce sparse and banded ‘self-stress’ and flexibility matrices with small band width. Some computational results are presented to demonstrate the superiority of the Turn-Back LU Procedure over the other schemes considered in this paper.     

Topology optimization of compliant circular path mechanisms based on an aggregated linear system and singular value decomposition

This paper proposes a topology optimization method for the design of compliant circular path mechanisms, or compliant mechanisms having a set of output displacement vectors with a constant norm, which is induced by a given set of input forces. To perform the optimization, a simple linear system composed of an input force vector, an output displacement vector and a matrix connecting them is constructed in the context of a discretized linear elasticity problem using FEM. By adding two constraints: 1, the dimensions of the input and the output vectors are equal; 2, the Euclidean norms of all local input force vectors are constant; from the singular value decomposition of the matrix connecting the input force vector and the output displacement vector, the optimization problem, which specifies and equalizes the norms of all output vectors, is formulated. It is a minimization problem of the weighted summation of the condition number of the matrix and the least square error of the second singular value and the specified value. This methodology is implemented as a topology optimization problem using the solid isotropic material with penalization method, sensitivity analysis and method of moving asymptotes. The numerical examples illustrate mechanically reasonable compliant circular path mechanisms and other mechanisms having multiple outputs with a constant norm. Copyright © 2011 John Wiley & Sons, Ltd.     

A new soft-particle DEM model of micro-particle impact integrated adhesive,elastoplastic and microslip behaviors

Micro-particle impact is a problem of solid mechanics that is common in many applications. To address this problem, a new soft-particle DEM model of micro-particle impact is proposed, which incorporates adhesive, elastoplastic and microslip behaviors. The normal force model is developed as two contiguous loading stages: the elastic stage and the elastoplastic stage in which the transition is from the elastic deformation to fully plastic deformation. Most innovative in unloading, the normal force model is also evolved into two contiguous stages: unloading under elastic loading and unloading under elastoplastic loading in which it combines Hertz elastic model and Mesarovic-Johnson plastic model. The normal force model is further assumed as the one-way coupling with pressure-based Maw tangential model with the micro-slip behavior. Further model validations are performed by employing the experimental results in literatures. The validation results indicate that model predictions agree with the experimental data, and are demonstrated to be incredibly accurate than other models, particularly for restitution coefficients and critical sticking velocity. Furthermore we can find that the smaller size particle has a longer period of nonlinear loading, while the larger size particle has a longer period of linear loading. For tangential restitution coefficient at the small incident angle, a down trend may be due to the oscillation of the tangential force.     

The effect of ordering on preconditioned GMRES algorithm,for solving the compressible Navier-Stokes equations

We investigate the effect of the ordering of the blocks of unknowns on the rate of convergence of a preconditioned non-linear GMRES algorithm, for solving the Navier-Stokes equations for compressible flows, using finite element methods on unstructured grids. The GMRES algorithm is preconditioned by an incomplete LDU block factorization of the Jacobian matrix associated with the non-linear problem to solve. We examine a wide range of ordering methods including minimum degree, (reverse) Cuthill-McKee and snake, and consider preconditionings without fill-in. We show empirically that there can be a significant difference in the number of iterations required by the preconditioned non-linear GMRES method and suggest a criterion for choosing a good ordering algorithm, according to the problem to solve. We also consider the effect of orderings when an incomplete factorization which allows some fill-in is performed. We consider the effect of automatically controlling the sparsity of the incomplete factorization through the level of fill-in. Finally, following the principal ideas of non-linear GMRES algorithm, we suggest other inexact Newton methods.