USE OF TREFFTZ FUNCTIONS IN NON-LINEAR BEM/FEM

Vladimír Kompiš, Milan Toma and Marián Handrik

Faculty of Mechanical Engineering, University of Žilina, Veľký Diel, 010 26 Žilina, Slovakia

Abstract - Simulation of many practical problems requires to use non-linear formulations with large displacements, large strains and large rotations. It is well known that the use of Trefftz (T-) functions (i. e. the functions satisfying the governing equations inside the domain) as weighting, or interpolation functions leads to more efficient formulations than those obtained by classical methods. In this paper we will show the use of T-functions and especially T-polynomials, Kupradze’s and Boussinesq functions (Green functions with singularity points defined outside of the domain) and their combination in connection with the total Lagrangian formulation for multi-domain BEM (reciprocity based FEM) analysis of displacements and for the post-processing phase in the analysis (evaluation of both gradient of displacements and stress fields). The formulation results in non-singular boundary integrals which has numerical advantages over other formulations using singular boundary integral equations.


Key words - Large strain-large rotation, multi-domain BEM, non-singular BIE, stress evaluation in the post-processing.


1. Introduction

   Trefftz (T-) functions [1] are functions satisfying all governing equations usually in the form of differential equations in the whole domain. The T-functions can be in the form of polynomials, harmonic, exponential, Bessel, Hankel, etc. functions, or they can be also so called Kupradze’s functions, i. e. the functions of the fundamental solution with source points outside the domain. They can be found only for linear problems. The T-functions are very efficient tool for the solution of linear problems. Two main formulations based on T-functions are used for the solution of the problems of solid and fluid mechanics: 1) hybrid FEM, in which the T-functions are used as interpolation functions of internal fields [2] and 2) reciprocity based BEM or FEM (Multi-Domain BEM), where the T-functions are used for reciprocity relations, which can be shown to be identical with definition of weighting functions for the governing equations satisfied in the weak (integral) sense.
   The simplest way to satisfy the boundary conditions is the collocation method, where the boundary conditions are satisfied in the discrete points of the boundary only. Such approximation does not guarantee the convergence of the solution in the multi-domain (MD) formulations. Many other approaches can be found in which the inter-domain continuity of the primary variables (e. g. displacements) is satisfied in the strong sense and the continuity of secondary variables (e. g. tractions) is enforced in a weak sense (e. g. weighted residuals, variational form, integral least squares, etc.).
    In the hybrid FEM, the internal and boundary fields are chosen independently, and the boundary field enforces both the inter-domain continuity and the inter-domain equilibrium are satisfied in a weak sense. Equivalent formulation cannot be found for non-linear problems, as corresponding T-functions are not known for internal variables.
   The most frequently used reciprocity formulations are well known as BEM formulations, in which the fundamental solution satisfies the governing equations in all points of the domain except one point of the boundary and lead to the solution by singular or hyper-singular boundary equations, which complicates the computational algorithms. On the other side, the formulations using the non-singular T-functions can cause the resulting system to be bad conditioned, which results in large numerical errors. However, using the formulation for simpler sub-domains makes the resulting system of equations sparse and good conditioned, as well.
   In this paper, we will show the use of reciprocity formulation for high non-linear problems arising by large strain, large rotation elasticity solved by the MD BEM, which can be named the Reciprocity Based (RB) FEM, too. The total Lagrangian (TL) formulation is used in the model and T-functions can be used as weighting functions for the weak equilibrium formulation similarly, as it is used by similar BEM formulation [3]. T-functions in polynomial form, Kupradze’s functions and a combination of both are used in our formulations. The formulations obtained in this way are compatible with classical FEM formulations [4] and the element matrices of both formulations can be arbitrarily combined into the total stiffness matrix.We demonstrated in our previous works [5] that for linear problems the stresses interpolated by T-polynomials from the patch of discrete nodal displacements and from the boundary conditions (the last ones are included for the points on or near boundaries only) are obtained with higher accuracy and converge faster than the stresses computed for the elements and smoothed out by some other procedures. When the relation between displacement and stress fields is non-linear, T-polynomials do not exist, and polynomials satisfying the equilibrium in the point of interest only (i. E. in the point, where the stress are to be computed) are used for the interpolation.
   Simple examples, where analytic solution can be found for comparison, demonstrate the formulation.


2. Total Lagrangian formulation for finite displacement problems

   In the total Lagrangian formulation the basic (equilibrium) equations refer to undeformed configuration of the body.
   Let Xi (i = 1, 2, 3 for 3D, i = 1, 2 for 2D problems) denotes the co-ordinate of a material particle X in the undeformed body. After the deformation, the co-ordinate of this particle will be xi. The cartesian components fij of the deformation gradient are defined by

(1)

   Using the displacement ui of a material particle X

(2)

leads to the alternative expression for the deformation gradient

(3)

where δij is Kronecker delta. Since the used formulation is presented in the undeformed configuration, partial derivatives denoted by (.),i are taken with respect to the undeformed co-ordinates Xi. The deformation gradient can be used to define the Green strain tensor

(4)

and the symmetric 2nd Piola-Kirchhoff stress tensor

   with   
(5)

both referring to the undeformed configuration. δij is the Cuachy stress tensor.
    The equilibrium equation

(6)

being expressed in the deformed configuration, can be transformed [6] to the initial (undeformed) configuration

(7)

where bi0 and bi denote the body force with respect to the initial and deformed cofiguration, respectively. The relation between the tractions, ti0, which express the force per unit undeformed area dΓ0, and the tractions, ti, in the deformed configuration dΓ is given by

(8)

where ni0 and ni is the outer surface normal in the initial and deformed configuration, respectively.
   Now, the equilibrium equation of the deformed body relative to initial configuration (7) will be used in the weak (integral) sense and T-displacements, i. e. the displacements satisfying the governing (equilibrium) equation in its homogeneous form (without any loads) in the whole domain, will be taken as the weight functions

(9)

The T-functions are taken in co-ordinates of the initial configuration for infinitesimal displacements, Ui, i.e. they will satisfy the Lame-Navier equations in their linear, homogeneous form (see Appendix A).
   Due to the total Lagrangian approach, Ω is the domain of the undeformed body, we omit its upper index 0 (similarly we do for the surface Γ) and all derivatives are taken with respect to this configuration.
    Applying integration by parts and the Gauss' theorem to equation (9) we obtain

(10)

Substituting the displacement gradients for the deformation gradient from (3) into (10) results in

(11)

Using the reciprocity relation in the form [3]

(12)

where Eij and Σij denote T-strain and -stress components corresponding to the T-displacements. Using this relation, the integration by parts and the Gauss' theorem once again, the equation (11) can be written in the form

(13)

where Ti are T-tractions derived from the T-stress ( Ti = Σij ni0 ). This equation expresses the reciprocity of works done by two systems of forces: the one, denoted by small letters, which is looked for and concerns the finite displacements, and the other, reference state (for which all, displacements, strain, stresses and tractions are known inside and on the boundaries), denoted by capital letters.
   In the classical BEM, the reference state defined by the Kelvin fundamental solution leads to the singular integral equation problem. When T-polynomials (see Appendix A), or Kupradze's functions (fundamental solution with singularity point located outside of the domain) are used for the reference state, all integrals are regular and simpler and more efficient algorithms can be used for their evaluation. Another possibility is to use Boussinesq functions (a force acting on a surface of a half space for 3D or a half plane for 2D) with the singularity outside the domain, or a combination of the Kupradze's or Boussinesq (both given shown in Appendix B) functions with the T-polynomials.
   With complexity of the shape and/or load conditions, higher order polynomials with large number of terms have to be used and the efficiency of the formulation decreases. Splitting the whole domain into the sub-domains, a MD BEM formulation is obtained with the sub-domain, or element matrices defined by the boundary integral equation (13). Such procedure enables to obtain a sparsely populated stiffness matrix for the whole structure similarly, as it is known from other FEM formulations [4], [6]. This improves also the numerical stability of the solution.
   For this purpose, the displacements between the sub-domains (elements) are chosen to be C0 compatible, however, the tractions are incompatible between such elements and the inter-element equilibrium and the natural boundary conditions will be satisfied in the weak, integral sense only:

(14)

where Γi, Γt and Γe are the inter-element boundaries, boundaries with prescribed tractions and element boundaries, respectively. The upper indices A and B denote the neighbouring elements.
   For the discretization, the boundary displacements and tractions are to be expressed by their values in nodal points (denoted by dash) and by shape functions, N, as

   or   
(15)
   or   
(16)

where ξ is the local co-ordinate of an element boundary point and capital letter index denotes the nodal point. The upper index e denotes the correspondence to the element.
   Similarly, displacements, their gradients and stresses in the volume integrals are interpolated from the nodal values of corresponding elements.
   The tractions are discontinuous in the edge and corner points, and thus a double and triple nodes for tractions is defined in those points.
   The displacement gradients in (13) are obtained from the derivatives of shape functions and corresponding nodal displacements as

(17)

   The last two integrals in equation (13) are non-linear functions of nodal displacements and thus, the dependence between displacements and tractions have to be solved iteratively. In order to obtain a tangential matrix, the non-linear terms have to be linearized. The integrand of the first of them will be

(18)

or, if the integrand is written in the form ( for the Hook's material)

(19)

and the material is isotropic, with

(20)

with λ and μ being the Lame material constants, then the linearized form is

(21)

The second integral will be linearized for Hook's isotropic material, for which

(22)

as follows

(23)

The field variables of the previous, (N - 1)-th, iteration step are computed from the nodal displacements of that step, whereas the field variables of the current, N-th, step are defined by the shape functions and corresponding nodal displacements of this step.
   The relation between the nodal displacements and nodal tractions (13) can be formally written as

(24)

or in the matrix form

(25)

where ue and te are vectors of element nodal displacements and tractions, respectively and

(26)

ξ(j) and ω(j) are co-ordinates and weights of the Gauss quadrature points and J is the Jacobian. The lower case index I denotes the I-th Trefftz function and the lower case index J corresponds to the nodal value. In order to distinguish the integration over the volume from that over the element boundaries, the volume Gauss points are denoted by α.
   We do not give the detail expressions for the non-linear term TiIJNL here, but its obtaining is straightforward from the linerized form of the integrands (18), (21) and (23).
   The discretized form of equation (14) is

(27)

or in the matrix form

(28)

with summation over all elements. The lower case indices K and L correspond to displacements and tractions, respectively.
   Setting from (28) for tractions into (25), the resulting system of non-linear equations is obtained in discretized form

(29)

or in the short form

(30)

   The matrix TL in eq. (29) is computed only once at the beginning of iteration process and contains the integration over the element boundaries, whereas the matrix TNL is computed in each following iteration step (but not in the first step) and contains the integration over the element volume only.


3. Displacement gradient and stress smoothing procedure in large displacement problems

   It is known that the stress smoothing using interpolation functions which satisfy the governing equations improves the accuracy and the rate of convergence of the stress field [10], [11]. The stress smoothing is usually applied in the post-processing stage of the computation. We have shown [11] that the rate of convergence of RB FEM in linear problems improves the rate of convergence of computed stresses from 1.5 in the worst (corner) points of the element up to 4, when T-polynomials interpolators are used together with boundary conditions.
   In non-linear finite displacement problems, the T-functions are non-linear in the stress-displacement relation and cannot be find exactly. The smoothing property of the inter-domain (inter-element) field inside the domain is well known and so, the moving least square (MLS) technique using classical interpolation polynomials can be applied to derive the smooth stress field inside the domain. Much larger errors are detected on the domain boundaries, especially those with prescribed tractions and in the regions closed to them. Inclusion of the local static boundary conditions into the smoothing technique enables to improve the accuracy of the stress fields. This technique leads to non-linear MLS procedure for obtaining the stress field on and near the domain boundaries.
   We assume here that the displacement formulation were used, i. e. we obtained the displacements as primary variables in the nodal points from the numerical solution. We will suppose that the discrete displacements obtained in this way are accurate enough comparing to the accuracy of their derivatives. It is well known that the stress and strain fields (secondary variables) computed from these displacements in classical way (i.e. element-wise) are of lower accuracy, discontinuous between the elements or sub-regions, for which the displacements were supposed to be continuous and that they converge slower than the primary variables.
   Following stress smoothing procedure will be supposed here: We will denote by the point of interest a discrete point where the displacement gradients and stress will be computed by the smoothing procedure. Usually, it will be a nodal point of an element (Fig.1). By the patch domain we will denote the domain containing all nodal points which will be considered to define the local polynomial displacement field in the least square (LS) sense. For linear problems, moreover the T-polynomials are used for this purpose. In large displacement problems full polynomials are used instead of the Trefftz polynomials in the smoothing procedures for displacements and the gradient fields derived from them. We suppose the vector of displacements approximated by

(31)

where [P(x)] is the full polynomial of co-ordinates {X} = (X1, X2, X3) in the original, undeformed configuration of the body and {c} is a vector of unknown coefficients.
    The numerical experiments show that the approach using full polynomials gives satisfactory results inside the domain. The components of the deformation gradient in each point of interest are obtained from the derivatives of the displacements (represented by the coefficients of the linear terms of the approximation polynomials (31), if the local co-ordinate system with the origin located in the point of interest is chosen) as

(32)

The components of the second Piola Kirchhoff stress tensor for isotropic material are given by

(33)

or

(34)

   If the point of interest is a point of the domain boundary, or a point near the domain boundaries, such approximation does not give satisfactory accuracy and the boundary conditions have to be included into the approximation. Only the traction boundary conditions, which are most important in practical problems are considered here. The tractions related to the undeformed configuration obtained from the displacement gradient are


(35)

or, those obtained from deformation gradient in shorter form

(36)

   The tractions are non-linear because of large strains and large rotations [6]. The MLS approximation of displacements by polynomial form defined over the discrete nodal displacements in the patch of nodes (domain of influence) and simultaneously satisfying the traction boundary conditions in the point of interest is non-linear and has to be solved in an iteration process as follows.
   The tractions can be expanded into the Taylor series

(37)

   The upper index 0 denoted the relation of tractions to the undeformed configuration was omitted for simplicity and the upper index (n) is now related to iteration step. The expansion (34) will be used in the Newton-Raphson (N-R) iteration process of numerical realisation of non-linear MLS smoothing of displacements when the boundary tractions are taken into account. Only first order terms of the Taylor expansion (34) will be considered.
   The necessary derivatives in (34) can be derived from (33) as follows:


(38)

with

(39)

   Having chosen the local origin of co-ordinates in the point of interest and realising from the equation (3) that the displacement gradients are

(40)

the following N-R iteration procedure for finding the constants {c} in (31) can be defined as follows (2D problem will be described for the sake of abbreviation of the record):
   In the n-th iteration step minimise:

(41)

with

(42)
(43)
(44)

where i = 1,2 and {ci}, {di} and {t} are the sub-vector of the vector of coefficients {c} corresponding to the i-th component of displacements, the vector of the i-th component of the nodal displacement and the vector of prescribed tractions, respectively. [Bi] are matrices corresponding to the definition (34) and (35).
   Solution of the problem (38) is given by:

(45)

Figure 1: Definitions

   Note that in the case of prescribed tractions in the boundary points, the boundary conditions define normal and shear components of the Cauchy stress tensor in the surface plane. We did not include weighting factors into the formulation (38) and (42) for simplicity. However, taking the prescribed boundary traction conditions with larger weighting function will define the lower importance of approximately computed displacements for evaluation of corresponding stress components and of the rotations and displacement gradients in such points. The displacements are important for evaluation of stresses and displacement gradients components in the tangential plane. In the corner points, the only quantities for which the displacements will be decisive are rotations and corresponding displacement gradients and not the stresses, if the tractions are known for the point.

4. Examples (click to continue)