AP04_4web.vp 1 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-brit- tle 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 polyno- mial 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 cor- rect 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 pro- vide the desired results, and also in order to gain a clear insight into the problem formulation, we used one-dimen- sional 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 disconti- nuities 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 solu- tions provided by the conventional chain of spring elements. 2 Strong discontinuity problem This section reviews general steps in the formulation of the problem of embedded discontinuities based on the parti- tions of unity method. In this framework, the discontinuous modes are introduced through the Heaviside step function di- rectly in the kinematic relations. The standard principle of virtual work is then used to arrive at the discrete system of linearized governing equations. 2.1 Kinematics of a displacement jump Consider a body � bounded by a surface � and crossed by a discontinuity �d, Fig. 1. �u represents a portion of � with prescribed displacements u while tractions t are prescribed on � �� � � �t u t d� � � 0 . The internal discontinu- 32 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 4/2004 Effect of Boundary Constraints in the Formulation of the Partition of Unity Method: One-dimensional Setting M. Audy, M. Šejnoha 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. Keywords: Strong discontinuity approach, Partition of unity method (PUM), Boundary constraints. Fig. 1: Body � crossed by discontinuity �d ity surface �d divides the body into two subdomains, � � and � � (� � �� � �� ). Suppose that the displacement field can be split into discontinuous and continuous parts u u u( , ) �( , ) ( ) ~( , )x x x xt t H t d � � � , (1) where H d� ( )x is the Heaviside function centered at the dis- continuity surface �d (H d + � �( ) ,x x� � �1 and H d� ( )x � 0, � � �x � ) and �u and ~u are continuous functions on �. Note that the discontinuity is introduced by the Heaviside function H d� ( )x at the discontinuity surface �d and that the magnitude of the displacement jump [[u]] at the discontinuity surface is given by ~u . For small displacements, the strain field assumes the form � � � �� �T T� ~ \u u+H d d� � �, x (2) where the operator matrix � can be found in [4]. 2.2 Governing equations The displacement field can be interpolated over the body � using the concept of the partition of unity method. For the purpose of the present work, it is sufficient to define a parti- tion of unity as a collection of functions �i which satisfy (see, e.g., [2] for more details) �i i n ( )x x� � � � 1 1 � , (3) where n is the number of discrete points (nodes). The dis- placement field will be interpolated in terms of discrete nodal values by � �u a bi i i i i n ( ) ( ) ( )x x x� � � � � 1 , (4) where �i is a partition of unity function, ai is the discrete nodal value and bi is the ’enhanced’ nodal value with respect to ‘enhanced’ basis �i. Note that the standard finite element shape functions Ni also form a partition of unity, since Ni i n ( )x x � � � � 1 1 � . (5) 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 dis- cretized form of the displacement field becomes � �u a b( ) = ( ) ( ) ( )x x x xN N N+ � , (6) where N is the matrix of standard nodal shape functions to interpolate regular nodal degrees of freedom and vector N N( )� 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 � ( ) ( ) ( )x x x� �B Ba b� (7) where B N� � T and B NN� �� � T( ). A suitable choice of N� may considerably improve the description of the displacement field for a specific class of problems [3]. When solving, e.g., the localized damage prob- lem, the discontinuous displacement field can be easily modeled after replacing the matrix N� by a scalar Heaviside function H d� 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 u a b( ) = ( ) ) ( )x x x xN N H+H d� ( , (8) � ( ) = ( ) ) ( )x x x xB B Ha b+H d� ( , (9) where Eqs. (8) – (9) hold for x �� �\ d . Thus the constitutive equations for the stress � at a point x �� �\ d and tractions developed on the discontinuity surface �d are �( ) �x � D � � �� ( ( ) ( ) )D B B Hx xa bH d� , (10) t b( ) ~ ( )x x� D N H . (11) Employing the principle of virtual work, the resulting dis- crete 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 consist- ing of standard and enhanced degrees of freedom �u a b� T, (13) and f lists externally applied forces �f f f� a b T, (14) where f ta t � � N T d� � , (15) f tb � �� N T d� �� , (16) and finally K represents the enhanced stiffness matrix K K K K K � � � � � � aa ab ba bb , (17) where individual sub-matrices are defined as K B DBaa � � T � d� � , (18) K K B DBab ba� ��T T � d� � , (19) K B DB N DNbb d � � �� � T T� ~d d� � � � . (20) 3 Numerical analysis of a one-dimensional problem The general formulation presented in the previous sec- tion 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 dis- continuities 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 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 33 Acta Polytechnica Vol. 44 No. 4/2004 Fig. 2: Simple chain model respect to the prescribed boundary conditions and a disconti- nuity position within an element is of primary interest. 3.1 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 ki is the spring constant and hi 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 F F F F Fk k h h� � � �1 2 1 2 , (21) u u u u uk k h h� � � �1 2 1 2 , (22) u F k � , (23) Substituting from Eq. (23) into Eq. (22) gives F k F k F k F h F h k k h h� � � �1 1 2 2 1 1 2 2 (24) and then using Eq. (21) provides the effective stiffness k in the form F k F k F k F h F h � � � � 1 2 1 2 (25) 1 1 1 1 1 1 2 1 2k k k h h � � � � (26) Simple generalization to m springs and n discontinuities yields 1 1 1 1 1 k k hii m ii n � � � � (27) 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 frame- work of PUM-based discontinuous elements, the jumps across individual discontinuities should be independent of the dis- continuity 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 con- stant and equal to � � F A, where F is the applied force and A is the element cross-sectional area, recall Eqs. (22)–(25). Ful- fillment of the above requirements will now be explored for several configurations. 3.2 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 ele- ment with n discontinuities. 3.2.1 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 k EA a a � , k EA b b � , (28) where E, A, a, b are the Young modulus, cross-sectional area and lengths of individual elements, respectively, ka and kb 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 �u � a a b b1 2 1 2 T. (29) Note that a one-dimensional bar element crossed by a dis- continuity has two degrees of freedom in each node, one standard and one enhanced. With reference to Fig. 3 (a), Eqs. (18) – (20) now become K B Baa a EA x� � T d 0 (30) K B Bab d EA x� � T d 0 (31) K B B N Nbb d EA x h� �� T Td 0 (32) Assuming standard linear interpolation functions for a one-dimensional bar element given by N � � �� � �� a x a x a , so that B � � �� � �� 1 1 a a , (33) and employing the notation in Eq. (28) we arrive at the fol- lowing element stiffness matrix 34 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 4/2004 Fig. 3: One discontinuity model K � � � � � � � �� k k k d a k d a k k k d a k d a k d a k d a k d a h d a a a a a a a a a a a a 1 � � � � � � � � � � � � � � � � � �� � � 2 1 1 k d a h d a d a k d a k d a k d a h d a d a a a a a � � � � � � � � � � � � � � � � � � � � � � � � � � � � �k d a h d aa 2 (34) 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 (a) the following global system of equations k k d a k d a k d a h d a a b Fa a a a � � � � � � � � � � � � � � � � � � � � � � � � �2 2 2 0� � � � . Solving for free degrees of freedom then yields � � u � � � � � � � � � � � � � dh ak k d h k ak F a dh ak dk Fa a a a a a( ) , T . (38) Eq. (38) clearly shows that the solution of the first configu- ration 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 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 35 Acta Polytechnica Vol. 44 No. 4/2004 Similarly we derive the stiffness matrix for the configura- tion plotted in Fig. 3 (b). On the structural level the vector of unknown degrees of freedom assumes the form �u � a a a b b b1 2 3 1 2 3 T. (35) The element stiffness matrix for element #2 is identical to that given by Eq. (34) oncereplacing ka by kb 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 b1 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 con- tained in the domain �� �( )a d and is discontinuity-free the element stiffness matrix provided by Eq. (34) reduces to K � � � � � � � � � � � � � � � � k k k k k k k k k k k k k k k k a a a a a a a a a a a a a a a a � � � � . (36) After assembly the global stiffness matrix of the configura- tion in Fig. 3 (b) becomes K � � � � � � � � � � � k k k k k k k k k k k d b k d b k k k a a a a a a b b a a b b b b b 0 0 0 0 d b k d b k k k k k k k d b k d b k k k d b h d b b a a a a a a b b a a b � � � � � � � �� � 0 0 1� � � � � � � � � � � � � � � � �� � 2 1 0 0 1 k d b h d b d b k d b k d b k d b h d b d b b b b b � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � k d b h d bb 2 . (37) 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 configu- ration and plotted in Fig. 4 (b). The same results, however, are not obtained for the first configuration. See Fig. 4 (a) sug- gesting 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. 3.2.2 PUM-based element with two discontinuities For the case of two discontinuities placed within an ele- ment the two possible configurations are plotted in Figs. 5 (a), (b), where d1 and d2 represent two arbitrary locations of the element discontinuities. 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 free- dom 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 consecu- tively with respect to the individual discontinuities according to Figs. 5 (a), (b). Thus the degrees of freedom (b1, b2, b3) cor- respond to discontinuity #1, whereas the degrees of freedom (b4, b5, b6) are linked to discontinuity #2. The element stiff- ness matrix then receives the following form K K K K K K K K K K � � � � � � � � � � aa ab ac ba bb bc ca cb cc , (40) where individual submatrices are provided by K B Baa a EA x� � T d 0 , (41) K K B Bab ba d EA x� � � T d 0 1 , (42) K K B Bac ca d EA x� � �T T d 0 2 , (43) K B B N Nbb d EA x h� �� T Td 0 1 1 , (44) K K B Bbc cb d EA x� � �T T d 0 1 , (45) K B B N Ncc d EA x h� �� T Td 0 2 2 . (46) 36 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 4/2004 k k k k k d b k d b k k k d b k d b k k d b k d b k k a b b a b b b b b b a b b a � � � � � � � � � b b b b b d b h d b k d b h d b d b k d b k d b k d b � �� � � � � � � � � � � � � � � � � � 1 1 2 h d b d b k d b h d b a b1 2 2 �� � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � a b b F3 2 3 0 0 0 � � � � � � � � � � � � � � � � � � � � � � � � � � � � � , and its solution listed in (39) u � � � � �� � � �� � � � � �� � � �� � � � � � � � 1 1 1 1 1 k h F k k h F F h F ha a b , , , � T , (39) 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. k EA aa � [N/m] k EA bb � [N/m] h [N/m] a [m] b [m] F [N] 100 50 50 1 2 100 Table 1: Material, geometrical and loading parameters Fig. 5: Two-discontinuities model 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. �u � � � � � � � �� � � �� � � � a a b b b b k h h F k k h a a b 2 3 2 3 5 6 1 2 1 1 1 1 1 1 1 1 T h F F h F h F h F h 2 1 1 2 2 � � �� � � �� � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � . (48) 3.2.3 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 discontinu- ities. Keeping the same ordering of the enhanced degrees of freedom as in the previous section, see also Fig. 7, we get K B Baa a EA x� � T d 0 , (49) K K B B N Nij ji d d ij iEA x h i j � � ��T T Td 0 min( , ) � , (50) where dij 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 reads �u � �a a b b b b b bn n2 4 2 3 5 6 3 1 3� T (51) © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 37 Acta Polytechnica Vol. 44 No. 4/2004 Fig. 6: Variation of the displacement field for a single crack element for two different configurations of Fig. 5 �u � � � � � � a b b d d h ad d d hk a d d k k d d a a a 2 2 4 1 2 2 2 1 1 2 2 1 2 2 1 2 T ( ) ( ) � �2 2 2 1 2 1 2 1 2 2 2 2 2 h d d d a d d hk a d d a d k F ad a a� � � � � � � � � ( ( )) ( )( ) h d d h d d d a d d hk a d d a d ka a1 2 2 2 2 1 2 1 2 1 2 2 22� � � � � � � �( ( )) ( )( ) F a d d h a d d k d d h d d d a d d hk a( ( ) ) ( ( )) � � � � � � � 1 2 1 2 1 2 2 2 2 1 2 1 22 a aa d d a d k F � � � � � � � � �� � � � � � � � � � �� � � � � � ( )( )1 2 2 2 . (47) 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 Fig. 7: n-discontinuities model u � � � � � � � � � � � � � � � � � � � � � � � 1 1 1 1 1 0 0 k h F k k h F F a ii n a b ii n h F h F h F h F h F h n n 1 1 2 2 � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � (52) revealing again a constant distribution of the discontinu- ous part of the displacement field within the discontinuous element. 4 Conclusions A simple one-dimensional example was given to demonstrate the essential drawback of the PUM-based discontinuous elements associated with constraining the en- hanced 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 dis- continuous 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 repre- sentation of the discontinuous part of the displacement field resulting in an erroneous response that depends on the dis- continuity location. Although the present results cannot be directly transplanted to the general case they suggest possible problems when applying the fixed kinematic boundary con- ditions to the enhanced degrees of freedom in higher dimen- sions as well, as typically done in [2]. 5 Acknowledgment This work was sponsored by research projects MSM 210000001,3. References [1] Audy M.: Localization of inelastic deformation in problems free of initial stress concentrators. Master thesis, Czech Techni- cal University in Prague, Faculty of Civil Engineering, Prague, 2003. [2] Babuška I., Melenk J. M.: “The partition of unity method.” International Journal for Numerical Methods in Engineering, Vol. 40 (1997), p. 727–758. [3] Babuška I., Banerjee U., Osborn J. E.: “On principles for the selection of shape functions for the generalized finite element method.” Computer Methods in Applied Mechanics and Engineering, Vol. 191 (2002), p. 49 – 50. [4] Bittnar Z., Šejnoha M.: Numerical methods in structural en- gineering. ASCE Press, 1996. [5] Brocca M.: Analysis of cracking localization and crack growth based on thermomechanical theory of localization. Ph.D. the- sis, University of Tokyo, 1997. [6] De Proft K.: Combined experimental-computational study to discrete fracture of brittle materials. Ph.D. thesis, Vrije Uni- versiteit Brussel, Brussel, 2003. [7] Simone A.: Continuous-Discontinuous Modelling of Failure. Ph.D. thesis, Delf University, Delft, 2003. [8] Moes N., Belitschko T.: “Extended finite element method for cohesive crack growth.” Engineering Fracture Mechanics, Vol. 69 (2002), p. 813 – 833. Ing. Miroslav Audy phone: +420 224 354 472 fax: +420 224 310 775 e-mail:miroslav.audy@fsv.cvut.cz Doc. Ing. Michal Sejnoha, Ph.D. phone: +420 224 354 494 fax:+420 224 310 775 e-mail:sejnom@fsv.cvut.cz Department of Structural Mechanics Czech Technical University in Prague Faculty of Civil Engineering Thákurova 7 166 29 Prague 6, Czech Republic 38 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 4/2004