Effect of Boundary Constraints in the Formulation of the Partition of Unity Method : One-dimensional Setting

The paper examines an effect of boundary constraints applied to the enhanced degrees of freedom of partition of unity based discontinuous elements. To highlight the present issue the problem is studied in a one-dimensional setting. In particular, an example of a one-dimensional bar element crossed by a set of discontinuities having a finite elastic stiffness clearly shows a need for proper approximation of the displacement field within a discontinuous element in order to correctly represent the structural response. While the discontinuous elements with boundary constraints applied to the enhanced degrees of freedom display an unrealistic dependence of the global response on the locations of the discontinuities, the discontinuous elements with complete approximation of the discontinuous part of the displacement field provide the expected global response independent of the locations of the discontinuities.


Introduction
The partition of unity concept [2], which allows a local enrichment of the standard finite element basis by special functions has been widely used to model the displacement discontinuity in a number of applications, e.g., the quasi-brittle failure of natural stones such as Massangis limestone [6] or continuous-discontinuous modeling of failure in high performance fiber-reinforced cement composites [7].In this framework, the discontinuity in the displacement field is introduced by enriching the standard finite element polynomial basis with the Heaviside function [8].This enrichment, however, results in additional degrees of freedom (enhanced degrees of freedom) in the nodes that belong to the domain affected by this enrichment.As these degrees of freedom are found in a set of displacement degrees of freedom the proper constraints must be applied to those found on the domain boundary in order to maintain the regularity of the resulting system of equations.Although not immediately evident, this step in certain applications may significantly pollute the correct solution of a given boundary value problem.
To introduce the subject, recall the problem of localization of the inelastic deformation in problems free of initial stress concentrators.This problem has been addressed, e.g., in the habilitation thesis of Brocca [5] and recently revisited in [1] using the concept of the partition of unity method, which allows the necessary splitting of the total displacement field into elastic and inelastic displacements associated with the crack opening.To test the ability of the latter approach to provide the desired results, and also in order to gain a clear insight into the problem formulation, we used one-dimensional setting, for which the exact solution is available.The presented numerical examples revealed several drawbacks associated with this approach.Among others, the study showed a possible depreciation of the results when an element crossed by a discontinuity containing a boundary node that had to be eliminated by the boundary constraints.
Motivated by the above result this paper attempts to shed a more detailed light on this problem and to clearly illustrate the need for a complete approximation of the discontinuous part of the displacement field in order to arrive at the correct results.To keep the analysis simple, attention is again limited to a one-dimensional bar element crossed by a set of discontinuities with finite elastic stiffness assigned to each of the predefined discontinuities.
The paper is organized as follows.Section 2 outlines the derivation of the linearized weak form of the governing equations.Application to a one-dimensional problem is then discussed in Section 3 and compared to the analytical solutions provided by the conventional chain of spring elements.

Strong discontinuity problem
This section reviews general steps in the formulation of the problem of embedded discontinuities based on the partitions of unity method.In this framework, the discontinuous modes are introduced through the Heaviside step function directly in the kinematic relations.The standard principle of virtual work is then used to arrive at the discrete system of linearized governing equations.).Suppose that the displacement field can be split into discontinuous and continuous parts

Kinematics of a displacement jump
where ) and $ u and ũ are continuous functions on W. Note that the discontinuity is introduced by the Heaviside function x at the discontinuity surface G d and that the magnitude of the displacement jump [[u]] at the discontinuity surface is given by ũ.For small displacements, the strain field assumes the form where the operator matrix ¶ can be found in [4].

Governing equations
The displacement field can be interpolated over the body W using the concept of the partition of unity method.For the purpose of the present work, it is sufficient to define a partition of unity as a collection of functions j i which satisfy (see, e.g., [2] for more details) where n is the number of discrete points (nodes).The displacement field will be interpolated in terms of discrete nodal values by ( ) where j i is a partition of unity function, a i is the discrete nodal value and b i is the 'enhanced' nodal value with respect to 'enhanced' basis g i .Note that the standard finite element shape functions N i also form a partition of unity, since In the standard finite element method, the partition of unity functions are shape functions and the enhanced basis is empty.When adopting the general scheme (4), the discretized form of the displacement field becomes ( ) where N is the matrix of standard nodal shape functions to interpolate regular nodal degrees of freedom and vector N N ( ) g b serves to introduce certain specific features of the displacement field u using so called enhanced degrees of freedom stored in vector b.Introducing Eq. (4) into Eq.( 2) gives the strain field in the form where B N = ¶ T and B N N A suitable choice of N g may considerably improve the description of the displacement field for a specific class of problems [3].When solving, e.g., the localized damage problem, the discontinuous displacement field can be easily modeled after replacing the matrix N g by a scalar Heaviside function H d G multiplied by a Boolean matrix H (a matrix with entry 1 if the corresponding degree of freedom is enhanced and zero otherwise).Eqs. ( 6) and ( 7) then become where Eqs. ( 8)-( 9) hold for x ÎW G \ d .Thus the constitutive equations for the stress s at a point x ÎW G \ d and tractions developed on the discontinuity surface Employing the principle of virtual work, the resulting discrete system of linear equations receives the traditional form, see [1] for more details, K u f = , (12) where u represents the vector of nodal displacements consisting of standard and enhanced degrees of freedom and f lists externally applied forces where and finally K represents the enhanced stiffness matrix where individual sub-matrices are defined as

Numerical analysis of a one-dimensional problem
The general formulation presented in the previous section will now be given in the context of a one-dimensional discontinuous bar element.In particular, we will consider the elements with one, two or an arbitrary number of discontinuities with a constant elastic stiffness assigned to each discontinuity.The effect of the problem setup in terms of number of elements, location of the crack element with Acta Polytechnica Vol.44 No. 4/2004 Fig. 2: Simple chain model respect to the prescribed boundary conditions and a discontinuity position within an element is of primary interest.

Simple chain model
First, assume a simple chain model consisting of a spring and a set of discontinuities, as shown in Fig. 2, in which k i is the spring constant and h i represents an elastic stiffness of the discontinuity i that relates the force transmitted across the discontinuity to the discontinuity opening displacement.The assumed arrangement of the individual elements in the chain model suggests Substituting from Eq. ( 23) into Eq.( 22) gives and then using Eq. ( 21) provides the effective stiffness k in the form (25) Simple generalization to m springs and n discontinuities yields 1 1 1 Note that in the previous derivation no assumption about the location of the discontinuity is required.It is therefore expected that, if addressing the same problem in the framework of PUM-based discontinuous elements, the jumps across individual discontinuities should be independent of the discontinuity location and should depend solely on the assigned discontinuity stiffness.The latter condition arises from the fact that the tensile stress in the structure should remain constant and equal to s = F A, where F is the applied force and A is the element cross-sectional area, recall Eqs. ( 22)-(25).Fulfillment of the above requirements will now be explored for several configurations.

PUM-based discontinuous elements
Three different configurations will be examined.First, we consider the most simple structure consisting of a spring and a single discontinuity.An element with two discontinuities is studied next, and finally we provide general results for an element with n discontinuities.

PUM-based element with one discontinuity
Two representatives of the possible numerical models appear in Figs.3(a),(b).However, before commenting on individual configurations we present the derivation of the element stiffnesses for typical elements in Figs.3(a),(b).To that end, we introduce the following notation where E, A, a, b are the Young modulus, cross-sectional area and lengths of individual elements, respectively, k a and k b then represent in analogy with Fig. 2 the corresponding spring constants and h is reserved for the discontinuity elastic stiffness.To proceed, consider an element in Fig. 3(a).By analogy with Eq. ( 13) the element degrees of freedom are ordered as Note that a one-dimensional bar element crossed by a discontinuity has two degrees of freedom in each node, one standard and one enhanced.With reference to Fig. 3 Assuming standard linear interpolation functions for a one-dimensional bar element given by and employing the notation in Eq. ( 28) we arrive at the following element stiffness matrix Finally, after imposing the boundary constraints to both standard and enhanced degrees of freedom and introducing the applied loading we get for the configuration displayed in Fig. 3  ) Eq. ( 38) clearly shows that the solution of the first configuration violates the basic requirement of being independent of the discontinuity location.This can be attributed to the fact that in this case the discontinuous element displacement field is not well approximated, as one of the two enhanced degrees of freedom is constrained.Consequently, the above solution when introduced into Eq.(10) gives a linear variation of the stress over the element, which is in direct contradiction with the results summarized in Section 3.1.
On the contrary, rather different results are discovered for the configuration of Fig. 3

(b). In this case the global system of equations reads
Acta Polytechnica Vol.44 No. 4

/2004
Similarly we derive the stiffness matrix for the configuration plotted in Fig. 3(b).On the structural level the vector of unknown degrees of freedom assumes the form The element stiffness matrix for element #2 is identical to that given by Eq. ( 34) oncereplacing k a by k b and a by b.It thus remains to determine the element stiffness matrix for element #1.Note that this element contains a node whose support (element #2) is crossed by a discontinuity.Also note that node #1 does not necessarily have to be enhanced, as the support of the associated nodal base function is not crossed by any discontinuity.Here, the node enhanced degree of freedom b 1 is preserved for the sake of simplicity in the der-ivation of the element stiffness matrix and will be eliminated via the boundary constraints.Since the entire element is contained in the domain W + = ( ) a d and is discontinuity-free the element stiffness matrix provided by Eq. (34) reduces to After assembly the global stiffness matrix of the configuration in Fig. 3(b) becomes Fig. 4: Variation of the displacement field for a single crack element for two different configurations of Fig. 3 Graphical representation of the above results derived using the material setting from Table 1 is plotted in Figs.4(a), (b).The figure shows the variation of the displacement field for two different crack locations.Note the expected constant distribution of the tensile stress found for the second configuration and plotted in Fig. 4(b).The same results, however, are not obtained for the first configuration.See Fig. 4(a) suggesting an unrealistic jump in the tensile stresses at the discontinuity location.A similar conclusion can be drawn for the problem of an element with two discontinuities studied below.

PUM-based element with two discontinuities
For the case of two discontinuities placed within an element the two possible configurations are plotted in Figs.5(a Moving in the footsteps of the previous section we first derive the element stiffness matrix.Owing to the presence of two discontinuities there are two enhanced degrees of freedom in each node.The two degrees of freedom in the first node of the second configuration in Fig. 5(b) are, however, inactive since the support of the node base function is not crossed by a discontinuity.For the solution of the underlying problem they will again be eliminated by the boundary constraints.
In order to derive the element stiffness matrix suppose that the enhanced degrees of freedom are ordered consecutively with respect to the individual discontinuities according to Figs. 5(a where individual submatrices are provided by is clearly independent of the discontinuity location d.In addition, the variation of the discontinuous part of the displacement field, recall Eq. ( 39), is constant in the discontinuous element.As expected, the solution in Eq. (47) depends, for the same reasons as already pointed out, on the locations of the two discontinuities and must be disqualified.
In contrast to the first configuration, the solution for the second configuration in Fig. 5(b), Eq. (48), does not suffer from this drawback.The correctness of this solution is again supported by Fig. 6(b), which shows a constant variation of the tensile stress along the bar unlike the plot in Fig. 6(a) derived for the first configuration.Also note the constant variation of the discontinuous part of the displacement field for both discontinuities.

PUM-based element with n discontinuities -general case
To complete our discussion we also present the derivation of the element stiffness matrix for the case of n discontinuities.Keeping the same ordering of the enhanced degrees of freedom as in the previous section, see also Fig. 7, we get where d ij is assumed to represent an identity matrix for i = j and a zero matrix for i ¹ j.By analogy with Eq. ( 48), the solution of the problem plotted in Fig. 7 As in the problems discussed in the previous section the solution of the two problems in Fig. 5 requires the introduction of the boundary constraints and loading.In particular, removing all the degrees of freedom in node #1 then gives after some algebra the solution of the first problem in the form

Conclusions
A simple one-dimensional example was given to demonstrate the essential drawback of the PUM-based discontinuous elements associated with constraining the enhanced degrees of freedom.It was shown that for the correct results to be independent of the locations of discontinuities the discontinuous part of the displacement field must be fully approximated.This can be accomplished by placing the discontinuous element away from the domain boundary.When the discontinuous element, however, contains a boundary node that must be constrained, the free degree of freedom in the other node is not sufficient to provide a correct representation of the discontinuous part of the displacement field resulting in an erroneous response that depends on the discontinuity location.Although the present results cannot be directly transplanted to the general case they suggest possible problems when applying the fixed kinematic boundary conditions to the enhanced degrees of freedom in higher dimensions as well, as typically done in [2].
Consider a body W bounded by a surface G and crossed by a discontinuity G d , Fig. 1.G u represents a portion of G with prescribed displacements u while tractions t are prescribed on ( ) The internal discontinu-ity surface G d divides the body into two subdomains, W + and W -

Fig. 3 :
Fig. 3: One discontinuity model (a) the following global system of equations degrees of freedom then yields

(
), (b), where d 1 and d 2 represent two arbitrary locations of the element discontinuities.
),(b).Thus the degrees of freedom (b 1 , b 2 , b 3 ) correspond to discontinuity #1, whereas the degrees of freedom (b 4 , b 5 , b 6 ) are linked to discontinuity #2.The element stiffness matrix then receives the following form

Fig. 6 :
Fig. 6: Variation of the displacement field for a single crack element for two different configurations of Fig. 5

Table 1 :
Material, geometrical and loading parameters