Strain-softening Damage Modeling UsingBoundary Element Method

A dual reciprocity boundary element method based on global shapes functions belonging to a complete set to transform domain integrals into boundary ones is presented. Governing equations of Poisson and Helmholtz type are tested. It is demonstrated that the complete set shape functions are more efficient and accurate in approximating domain integrals than those based on incomplete set. INTRODUCTION The application of Boundary Element Method to linear partial differential operators, or to nonlinear and time-dependent ones, often leads to domain integrals that are cumbersome to handle. The dual reciprocity method pioneered by Nardini and Brebbia* is a powerful numerical technique which converts domain integrals into boundary ones using a set of global shape functions to approximate the integrand. It has been successfully applied in many problems involving dynamic, diffusion and nonlinear type operators. See Partridge, et al. for a comprehensive bibliographical review. The success of the dual reciprocity hinges on the ability to express the global shape functions into the Laplacian of a set of known functions such that the divergence theorem can be applied to resolve the domain integral into a boundary one. Such a constraint has so far limited the selection of shape functions to a set of radial functions of the form 1 + ar -f br^ -f • • •• Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X 344 Boundary Elements These radial functions do not form a complete set hence may lack convergence properties. Although practical evidences have indicated that those shape functions work quite satisfactorily, it is still desirable to use shape functions from a set that is complete in the mathematical sense. This paper present a dual reciprocity BEM based on global shape functions derived from a polynomial series. Perhaps due to its mathematical completeness, we are able to demonstrate that the numerical results are more efficient and accurate than the traditional radial set. DUAL RECIPROCITY For many non-standard, transient, or nonlinear partial differential equations it is possible to extract a standard, time-independent, linear partial differential operator C and lump the remainder to the right-hand-side: jC% = 6(z, ?/, Z, f , %, . . . , %,., %„ , %,, %,, %= . • •) (1) For the purpose of BEM implementation, these non-standard portions can be formally regarded as body forces. By the application of generalized Green's theorem = / Jr (2) where v is a conjugate variable, £* is the adjoint operator of £, B» and B* are operators resulting from integration by parts, which may be viewed as generalized normal derivatives, F is the solution boundary, and fi the solution domain, a boundary integral equation can be formed as: cu = / (u*q uq*) j bu* dSl (3) where quantities denoted by asterisk are free-space Green's functions satisfying £*u* = —6, with 6 a Dirac delta function, c is a constant dependent on the Cauchy principal value integration of the singularity at the base point, q and q* are generalized normal derivatives of u and u*. To circumvent the domain integration, the function b can be approximated by an expansion based on a set of global shape functions /, such that &*fX£ (4) j=l in which a, are coefficients determined by collocation, and n^ is the number of collocation nodes. The shape functions are chosen such that particular solutions of the set of inhomogeneous equations Ciij = /,• (5) Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X Boundary Elements 345 are explicitly known. The domain integral in (3) can be transformed as [ / j^bu da = i ŷ /ju (6) j=l L ^1 J The last expression was a consequence of the application of (2). Substitution of the above into (3) yields / (u*q uq} dY = Y a, \cu, I (%*& titf) dT (7) •/r fr(. L 7r J cu We hence in principal have achieved the elimination of the domain integral. However, the detail of the implementation is somewhat dependent on the type of body force terms. We shall illustrate below the implementation by examples. EXAMPLES Operator of V^u = b(x, y) type For the governing equation of Poisson type, q in (7) is du/dn, and u* is — lnr/2%for 2-D or l/4%r for 3-D geometries. For simplicity, only two-dimensional problems will be studied herein. If the right-hand-side 6(x, y) is a mathematical expression whose particular solution is known as Up, the solution of (8) can be expressed as % = %& + %? (9) where u^ is the homogeneous solution and V2%& = 0 (10) It is apparent that Uh = u — Up satisfies the boundary integral equation c(u «„) = / [u(q q,} -(uu,)q'] (11) If we approximate the right-hand-side b(x,y) as (4), and find the particular solution term by term as V = /,(12) Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X 346 Boundary Elements it is clear that j=l By the above interpretation, we have provided an alternative way for deriving (7), which is rewritten below with explicit reference to coordinates as c(x)u(x)f(u*(x',x)q(x')-u(x')q-(x',x)] dT(x') j r

A boundary element method for analysing non-linear phenomena on gas-liquid two-phase free surface at middle and higher Reynolds numbers has been developed. Employing this method we analyze two kinds of non-linear phenomena on gas-liquid two-phase free surface (l) Solitary wave on free falling thin liquid film. (2) Surface tension drives convective flows in open-boat type crystal-growth techniques in microgravity environment. INTRODUCTION Although some new schemes of discretizing convective terms are developed as QUICK and QUICK ER in finite difference methods[FDM] and finite element methods[FEM], but it is very difficult to avoid numerical diffusion introduced by convective discretizations, especially, at high Reynolds number. In [1], they developed a boundary element iterative scheme(BEIS) to calculate convective terms and compared the advantages of using upwind or central difference schemes and BEIS. As shown as numerical results, BEIS of the convective terms was the most accurate of the two approaches. However, computation was complex. The main objective of this research is to develop a boundary element method for analysing some non-linear phenomena on gas-liquid two-phase free surface at middle and higher Reynolds numbers. On the basis of the boundary element formulation with the pressure penalty function, the divergence theorem is applied to the non-linear convective volume integral of the boundary element formulation. Consequently, velocity gradients are eliminated, the complete formulation is written in terms of velocity. This avoids the difficulty of convective discretizations, and provides considerable reducing storage and computational requirement while improving accuracy. In order to calculate free surface, we construct shooting iterative method of free surface on the basis of normal stress balance condition Project supported by the National Natural Science Foundation of China Transactions on Modelling and Simulation vol 6, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X 168 Free and Moving Boundary Problems similar to [2] [3]. The method is applied to analyze two kinds of free surface problems. (1) Solitary wave on free falling thin liquid film. This problem was researched by Wasden and Kheshgi etc. with finite difference method [4] and finite element method [5]. However, in finite difference method[4], one applied an adjustment procedure for wave velocity in order to obtain reasonable flow fields, in finite element methods[5] numerical solution was limited to low flow rates. Moreover, the upwind schemes of discretizing convective momentum terms common to FDM and FEM introduce numerical diffusion or oscillatory behavior of the wave peak region, more importantly, the upwind scheme lacks sensitivity to cross stream diffusion and source terms which are especially important in the case of a thin film. Therefore it is valuable to research a non-traditional numerical method such as BEM. (2) Surface tension drives convective flows in open-boat type crystal growth techniques in microgravity environment. The analytical or numerical study of free surface problem in crystal-growth techniques is very difficult. We shall restrict ourselves to model a problem of open-boat type crystal growth technology, and we only study the melt-gas free boundary, assume the melt-crystal interface shape is flat. The free surface problem of open-boat type is simplified as the free surface problem in open rectangular cavity. The left vertical boundary of the cavity can be interpreted as the crucible wall. The right vertical boundary of the cavity can be interpreted as the melt-crystal interface. The temperatures equal #/ and 0,., respectively. Some important numerical results are presented in paper. BASIC EQUATIONS We consider a system of steady incompressible viscous liquid in a open rectangular cavity with differentially heated lateral walls (the temperature difference is A0*). There is a gradient of surface tension on the free surface since the distribution of temperature is non-homogeneous. The gradient of surface tension drives the convection of liquid. On the basis of the basic law of fluid mechanics and Boussinesq approximation, the governing equations of these problems can be described by a tensor notation and a non-dimensional form as follow: continuity equation «,-„• = o (i) momentum equations energy equation V,-^,=0 (3) where, #,„% is the Kronecker delta symbol, m = x for solitary wave, m = y for thermocapillary convection. The lengths, velocities are scaled with respect to H, Transactions on Modelling and Simulation vol 6, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X

Kussmaul's (1969) formulation of acoustic scattering problems involves the superposition of a double-layer distribution and a simple-layer distribution on the boundary.This ensured the existence of a unique solution for all wave numbers of the incoming wave. However,the boundary integral equation involves a highly singular kernel which proves difficult to treat numerically. An alternative approach is to place the double-layer distribution on an interior boundary similar (similarly situated) to the given boundary.This avoids the singularity problems and enables excellent numerical results to be achieved.Some numerical results are presented. INTRODUCTION Problems of hard acoustic scattering may be formulated in terms of boundary integral equations.Direct formulations are based upon the Helmholtz formula;indirect formulations utilise layer potentials.This is a Neumann boundary-value problem which has a unique solution subject to regular behaviour at infinity.However the relevant classical integral equations,whether direct or indirect,fail at wave numbers corresponding with the eigenfrequencies of the interior Dirichlet problem.In this paper we show that a relatively simple, nonsingular,supplement to the classical indirect integral equation ensures a a unique solution at all wave numbers. An immediate consequence is that meaningful numerical solutions can be achieved,for a wide range of wavenumbers,by computationally effective boundary-element methods (BEMs). Some results for scattering by a hard sphere are provided.For completeness we mention some recent papers presenting alternative resolutions of the difficulty. Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X 282 Boundary Elements Rego Silva et al [1] Numerical Implementation of a Hypersingular BEM formulation for Acoustic Radiation Problems,following the approach of Burton and Miller [2].Here the hypersingular integral is evaluated numerically by an algorithm of Guiggiani et al [3] building essentially upon Hadamard's finite-part integration. Liu and Rizzo [4] A weakly-singular form of the hypersingular boundary integral equation applied to 3-D acoustic wave problems,again following the approach of Burton and Miller. However here the hypersingular integral is weakened by a regularisation technique,which enables it to be evaluated numerically by straightforward methods. Wu et al [8] Vectorization and parallelization of the acoustic boundary element code BEMAP on the IBM ES/3090 VF,exploits the supercomputing facilities to improve the CHIEF method originated by Schenck [15]. Kirkup and Henwood [5] Computational Solution of Acoustic Radiation Problems by Kussmaul's Boundary Element Method,adopting Kussmaul's representation [6] of the wave function as the superposition of a simple-layer and a double-layer potential generated from sources on the boundary.This formulation also involves a hypersingular integral,which is evaluated numerically by the method of Terai [7],an early exploitation of Hadamard's finite part integration. Kussmaul's approach may be modified by locating the double-layer sources on an interior boundary,similar and similarly situated to the given boundary.This eliminates the hypersingular integral and at the same time provides a unique solution at all wave numbers. Brief details are given in the following section. KUSSMAUL FORMULATION Any solution, ,of the Helmholtz equation existing on the infinite domain B exterior to a closed (simply connected) surface dB (fig. 1) of order 0(r~') as /•co can generally be represented by a simple-layer potential dq ; q € dB , p € B+dB (1) dB where o(q) is the source density at q and dq is the area element at q ; also Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X Boundary Elements 283 /-',-'» ; r=|f-g| (2) where k is the wave number associated with , defined by

The Boundary Element Method (BEM) is not local and generates a full and nonsymmetric matrix, therefore the matrix solution time could easily grow beyond acceptable limits. Furthermore, in many cases a decoupling between geometry description and solution approximationwould be advantageous and desirable. This work aims at offering the BEM a new approximation capability, which could reduce the number of degrees of freedom (dofs) needed, and at providing the desired decoupling. The method proposed here takes advantage of the less restrictive requirements in BEM compared with FEM for the field approximation. The base of the new approximation becomes what has been called “Spectral Element” (SE): a super-entity gathering many elements, which are used only for geometry description and numerical integration. Spectral shape functions are generated using smooth positive definite compact supported Radial Basis Functions (RBF). The method was applied to linear elasticity in 3D and numerical tests were performed for different typical cases in order to highlight advantages, disadvantages and limits of the proposed methodology. 1 Approximation field in the BEM 1.1 Some general properties of the “finite element approximation” Both BEM and FEM, in order to reduce the infinite dimensional space of the solution to a finite one, use the so called finite transform. If x represents the vector of space coordinates and u the approximation field, a set of “generalised degrees of freedom” in the form of a vector α must be chosen so that: u(x) = A(x)α (1) Boundary Elements XXVII 165 © 2005 WIT Press WIT Transactions on Modelling and Simulation, Vol 39, www.witpress.com, ISSN 1743-355X (on-line) In equation (1) A(x) represents the coordinate function matrix, which contains functions fixed in space. The finite elements appear when the continuum under analysis is ideally cut into “small” pieces, the elements of the mesh. In order to ensure that the values of the field in few specific points (nodes) belonging to an element govern the whole behaviour of u inside it, equation (1) must be imposed simultaneously in each of the n nodes of an element:   u1(x) u2(x) .. un(x)   = U = Cα =   A(x1) A(x2) .. A(xn)  α (2) Now if C is square, i.e. if the number of nodes n is equal to the number of α dofs, and is not singular, the following equations hold: u(x) = A(x)α = A(x)C−1U = Φ(x)U where Φ = A(C)−1 (3) Generally the Φ matrix is called shape matrix. It interpolates the dof of the element, represented by the nodal values, in between the nodes. In the FEM four general properties are needed for the functions in Φ, called shape functions, each of which belongs to a different node: 1. Completeness condition: in order to be able to represent a constant field, (first term in a Taylor series development of any field) the sum of all functions must be one in every point of the element:

In the boundary element method (BEM) for stress analysis, it is well known that thermal loads give rise to an additional volume integral in the primary form of the boundary integral equation (BIE). This volume integral needs to be further transformed to surface ones in order to retain the characteristic of the BEM as a boundary solution technique. In this study of the BEM for 3D thermoelasticity in general anisotropy, the fundamental solutions are expressed as Fourier series with coefficients calculated using an explicit-form Green’s function. In the exact volume-to-surface integral transformation associated with the term for the thermal effects in the BIE, a new kernel function is constructed. All formulations are implemented in an existing BEM code for 3D elastostatic analysis. Some numerical examples are presented to demonstrate the veracity of the formulations and the implementation, where the numerical results are compared with those obtained using the finite element method (FEM).

In this paper, formulations for homogenization method using the boundary element method (BEM) for unit cell analysis are presented. Unlike the finite element method (FEM), the application of BEM for the homogenization analyses is not fully explored. Possible formulations, a full knowledge of numerical implementation and appropriate expressions for the effective elastic moduli, when BEM is employed in unit cell analysis, are not fully explored yet. The purpose of present paper is to discuss possible equation formulations, numerical implementations and formulations for the effective elastic moduli, which are suited for BEM unit cell analysis. Two different BEM formulations for unit cell analysis (multi-region and single-region with volume integral) are presented. Periodic boundary conditions and numerical procedures are discussed, and then formulas for effective elastic moduli are presented. Finally some numerical examples are presented.

The present work deals the development of a nonlinear numerical model for structural analysis of solids composed by multi-domains considering cohesive discontinuities along its interfaces. The numerical method adopted is the boundary element method (BEM), through its singular and hyper–singular integral equations. Due to the mesh dimensionality reduction provided by BEM, this numerical method is robust and accurate for analyzing the fracture process in solids, as well as physical nonlinearities that occurs along the body’s boundaries. Multi-domain structures are modelled considering the sub-region technique, in which both equilibrium of forces and compatibility of displacements are enforced along all interfaces. The crack propagation process is simulated by the fictitious crack model, in which the residual resistance of the region ahead the crack tip is represented by cohesive tractions. It leads to a nonlinear problem relating the tractions at cohesive interface cracks to its crack opening displacements. The implemented formulation is applied to analysis of three examples. The numerical responses achieved are compared to numerical and experimental solutions available in literature in order to show the robustness and accuracy of the formulation.

The direct boundary integral equations method has been applied to analyze stresses in a fracture process zone (a crack bridged zone) and to calculate stress intensity factors module for structures with bridged interfacial cracks under mechanical loading. Bridged zones at interfacial cracks are considered as parts of these cracks with assumption that surfaces of interfacial cracks are connected by distributed spring-like bonds with given bond deformation law. For numerical analysis of piecewise structures with bridged interfacial cracks the multi-domain formulation of the boundary elements method is used. The stress intensity factors module evaluation is performed on the basis of displacements and stresses computed at nodal points of special quadratic boundary elements adjoined to a crack tip. The comparative study between the results obtained by the boundary elements method and the results obtained previously by the singular integral---differential equations method is performed and the validity of the presented numerical formulation is demonstrated. The new problem for a bridged circumferential crack between a cylindrical inclusion and a matrix in plate of finite size is also solved. Stresses distributions along the bridged zone and the stress intensity factors modulus dependencies versus the bridged zone length and bonds stiffness are presented and discussed for this problem.

A direct traction BIE approach for three-dimensional crack problems

Optimal design of machine components using notch correction and plasticity models

An approach in solving defect shape identification problems is proposed in this paper. This approach is developed with the boundary element method in conjunction with the Kalman filtering technique, especially the extended Kalman filtering technique. The boundary element approach is utilized by reason of easyness on lemeshing of the defect shape in identification process. Applicability and effectiveness of the proposed method are examined with numerical simulation of several identification examples on a circular or an elliptic defect in 2-D elastostatics.

A numerical solution of general frictional contact problems by the direct boundary element and mathematical programming approach

The influence of the magnetic configuration on the behavior of magneto-electro-elastic bimorph beams is analyzed by using a boundary element approach. The problem is formulated by using the generalized displacements and generalized tractions. The boundary integral equation formulation is obtained by extending the reciprocity theorem to magneto-electro-elastic problems; it is numerically implemented by using the boundary element method multidomain technique to address problems involving nonhomogeneous configurations. Results under different magnetic configurations are compared highlighting the characteristic features of magnetopiezoelectric behavior particularly focusing on the link between interlaminar stress and magnetic induction.

Numerical analysis of acoustic radiation responses of shear deformable laminated composite shell panel in hygrothermal environment

Solution of some 2D transport problems by a high order A N–SP 2 N−1 method