

USE OF TREFFTZ FUNCTIONS IN NONLINEAR BEM/FEMAbstract  Simulation of many practical problems requires to use nonlinear 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 Tfunctions and especially Tpolynomials, 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 multidomain BEM (reciprocity based FEM) analysis of displacements and for the postprocessing phase in the analysis (evaluation of both gradient of displacements and stress fields). The formulation results in nonsingular boundary integrals which has numerical advantages over other formulations using singular boundary integral equations. Key words  Large strainlarge rotation, multidomain BEM, nonsingular BIE, stress evaluation in the postprocessing. 1. Introduction Trefftz (T) functions [1] are functions satisfying all governing equations usually in the form of differential equations in the whole domain. The Tfunctions 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 Tfunctions are very efficient tool for the solution of linear problems. Two main formulations based on Tfunctions are used for the solution of the problems of solid and fluid mechanics: 1) hybrid FEM, in which the Tfunctions are used as interpolation functions of internal fields [2] and 2) reciprocity based BEM or FEM (MultiDomain BEM), where the Tfunctions 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 multidomain (MD) formulations. Many other approaches can be found in which the interdomain 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 interdomain continuity and the interdomain equilibrium are satisfied in a weak sense. Equivalent formulation cannot be found for nonlinear problems, as corresponding Tfunctions 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 hypersingular boundary equations, which complicates the computational algorithms. On the other side, the formulations using the nonsingular Tfunctions can cause the resulting system to be bad conditioned, which results in large numerical errors. However, using the formulation for simpler subdomains 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 nonlinear 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 Tfunctions can be used as weighting functions for the weak equilibrium formulation similarly, as it is used by similar BEM formulation [3]. Tfunctions 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 Tpolynomials 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 nonlinear, Tpolynomials 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 X_{i} (i = 1, 2, 3 for 3D, i = 1, 2 for 2D problems) denotes the coordinate of a material particle X in the undeformed body. After the deformation, the coordinate of this particle will be x_{i}. The cartesian components f_{ij} of the deformation gradient are defined by
Using the displacement u_{i} of a material particle X
leads to the alternative expression for the deformation gradient
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 coordinates X_{i}. The deformation gradient can be used to define the Green strain tensor
and the symmetric 2nd PiolaKirchhoff stress tensor
both referring to the undeformed configuration. δ_{ij} is the Cuachy stress tensor. The equilibrium equation
being expressed in the deformed configuration, can be transformed [6] to the initial (undeformed) configuration
where b_{i}^{0} and b_{i} denote the body force with respect to the initial and deformed cofiguration, respectively. The relation between the tractions, t_{i}^{0}, which express the force per unit undeformed area dΓ^{0}, and the tractions, t_{i}, in the deformed configuration dΓ is given by
where n_{i}^{0} and n_{i} 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 Tdisplacements, 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
The Tfunctions are taken in coordinates of the initial configuration for infinitesimal displacements, U_{i}, i.e. they will satisfy the LameNavier 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
Substituting the displacement gradients for the deformation gradient from (3) into (10) results in
Using the reciprocity relation in the form [3]
where E_{ij} and Σ_{ij} denote Tstrain and stress components corresponding to the Tdisplacements. Using this relation, the integration by parts and the Gauss' theorem once again, the equation (11) can be written in the form
where T_{i} are Ttractions derived from the Tstress ( T_{i} = Σ_{ij} n_{i}^{0} ). 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 Tpolynomials (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 Tpolynomials. 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 subdomains, a MD BEM formulation is obtained with the subdomain, 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 subdomains (elements) are chosen to be C^{0} compatible, however, the tractions are incompatible between such elements and the interelement equilibrium and the natural boundary conditions will be satisfied in the weak, integral sense only:
where Γ_{i}, Γ_{t} and Γ_{e} are the interelement 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
where ξ is the local coordinate 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
The last two integrals in equation (13) are nonlinear 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 nonlinear terms have to be linearized. The integrand of the first of them will be
or, if the integrand is written in the form ( for the Hook's material)
and the material is isotropic, with
with λ and μ being the Lame material constants, then the linearized form is
The second integral will be linearized for Hook's isotropic material, for which
as follows
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, Nth, 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
or in the matrix form
where u^{e} and t^{e} are vectors of element nodal displacements and tractions, respectively and
ξ^{(j)} and ω^{(j)} are coordinates and weights of the Gauss quadrature points and J is the Jacobian. The lower case index I denotes the Ith 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 nonlinear term T_{iIJ}^{NL} here, but its obtaining is straightforward from the linerized form of the integrands (18), (21) and (23). The discretized form of equation (14) is
or in the matrix form
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 nonlinear equations is obtained in discretized form
or in the short form
The matrix T^{L} in eq. (29) is computed only once at the beginning of iteration process and contains the integration over the element boundaries, whereas the matrix T^{NL} 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 postprocessing 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 Tpolynomials interpolators are used together with boundary conditions. In nonlinear finite displacement problems, the Tfunctions are nonlinear in the stressdisplacement relation and cannot be find exactly. The smoothing property of the interdomain (interelement) 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 nonlinear 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. elementwise) are of lower accuracy, discontinuous between the elements or subregions, 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 Tpolynomials 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
where [P(x)] is the full polynomial of coordinates {X} = (X_{1}, X_{2}, X_{3}) 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 coordinate system with the origin located in the point of interest is chosen) as
The components of the second Piola Kirchhoff stress tensor for isotropic material are given by
or
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
or, those obtained from deformation gradient in shorter form
The tractions are nonlinear 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 nonlinear and has to be solved in an iteration process as follows. The tractions can be expanded into the Taylor series
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 NewtonRaphson (NR) iteration process of numerical realisation of nonlinear 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:
with
Having chosen the local origin of coordinates in the point of interest and realising from the equation (3) that the displacement gradients are
the following NR 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 nth iteration step minimise:
with
where i = 1,2 and {c_{i}}, {d_{i}} and {t} are the subvector of the vector of coefficients {c} corresponding to the ith component of displacements, the vector of the ith component of the nodal displacement and the vector of prescribed tractions, respectively. [B_{i}] are matrices corresponding to the definition (34) and (35). Solution of the problem (38) is given by:
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) 