A Le Van

P Le Grognec

Modeling and Numerical Computation of Necking in Round Bars Using a Total Lagrangian Elastoplastic Formulation

Keywords: elastoplasticity, necking bifurcation, total Lagrangian formulation, finite element analysis

Necking is a bifurcation phenomenon observed in round bars under tensile loading and has been investigated in numbers of papers. In the present work, it is modeled within the framework of finite rateindependent plasticity. The theory is based on thermodynamic foundations developed for standard materials and results in a total Lagrangian formulation for finite plasticity, where the total strain is decomposed additively according to [Green and Nagdhi 1965)] and the hardening is characterized by a nonlinear isotropic hardening law of the saturation type. The discretization and consistent linearization of the elastic-plastic equation set using the standard finite element procedure lead to a low-cost algorithm, robust enough to deal with necking problems. The numerical computations of necking are performed on cylindrical bars with various boundary condition types and the corresponding results compared with those obtained in the literature.

Introduction

Necking is an bifurcation phenomenon observed in circular cylindrical bars under tensile loading and has been investigated in numbers of works for its theoretical modeling and numerical computation. The pioneering papers dated back to the 70's and were due to [START_REF] Needleman | A numerical study of necking in circular cylindrical bars[END_REF]], [START_REF] Argyris | On the large strain inelastic analysis in natural formulation. Part I : quasistatic problems[END_REF]]. A lot of significant results were obtained in these works, although the computed necked configurations did not resemble experimental observations. Later, [START_REF] Besdo | Total Lagrangian strain-spacerepresentation of the elasto-plasticity of metals[END_REF]] was among the first to give a realistic necking shape. More recently, studies using other formulations have been developed by [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. Part I : Continuum formulation[END_REF]], [START_REF] Simo | Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory[END_REF]], [START_REF] Simo | Computational Inelasticity[END_REF]] and [Brtinig (1998)], providing further results on necking. The numerical solution by Simo necessitates mixed finite elements, moreover the line search method is used to get across the limit point.

Whereas most papers considered a cylindrical circular bar with both ends assumed to be shear free, [START_REF] Needleman | A numerical study of necking in circular cylindrical bars[END_REF]] dealt with two types of boundary conditions, shear free ends and cemented to rigid grips ones. [Brtinig (1998)] went further by also considering the whole tensile specimen with shear free stiff heads and then compared the results corresponding to the three cases.

Most of the previous investigations were performed within the framework of finite rate-independent plasticity, assuming a yield criterion in the stress space and an isotropic hardening. Outside this context, analyses using the strain space representation of the finite plasticity were proposed by [START_REF] Besdo | Total Lagrangian strain-spacerepresentation of the elasto-plasticity of metals[END_REF]] and [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. Part I : Continuum formulation[END_REF]]. Other nontraditional approaches can be used to successfully solve the necking problem, as shown by [START_REF] Atluri | New general and complementary energy theorems, finite strain, rate sensitive inelasticity and finite elements : some computational studies[END_REF]], using a modified complementary energy principle instead of the conventional principle of virtual work. Also, the integral equation method was shown to be capable of solving bifurcation and post-bifurcation problems, see [START_REF] Okada | A full tangent stiffness field-boundary-element formulation for geometric and material non-linear problems of solid mechanics[END_REF]].

On the other hand, thermal aspects in necking were taken into account by [START_REF] Kleiber | On solving problems of advanced plastic flow[END_REF]] who presented a simple solution scheme for elastoplasticity as well as thermoelastoplasticity. Also, one can use none of finite plasticity but a viscoplastic type constitutive law in small strains in order to describe necking behavior in viscoplasticity, e.g. see [START_REF] Rappaz | Modlisation Numrique en Science et Gnie des Matriaux. Volume 10 of the collection "Trait des Matriaux[END_REF]].

In the present work, necking is studied within the framework of finite rate-independent plasticity. The theory . is based on thermodynamic foundations established for generalized standard materials [START_REF] Halphen | Sur les materiaux standard generalises[END_REF]], [Germain, [START_REF] Nguyen | Continuum thermodynamics[END_REF]] andresults in a total Lagrangian formulation for finite plas-ticity. The total strain is decomposed additively according to [START_REF] Green | A general theory of an elastic-plastic continuum[END_REF]]. The yield function is of von Mises type and depends on the symmetric Fiola-Kirchhoff stress tensor. The hardening in the material is characterized by a nonlinear isotropic hardening law of the saturation type, as in [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. Part I : Continuum formulation[END_REF], Part II] and [START_REF] Simo | Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory[END_REF]]. However, here it relates Lagrangian (material) variables, namely the Lagrangian hardening variable and the Lagrangian effective plastic strain.

The numerical implementation using standard finite element analysis is simple and yet robust enough to accurately reproduce the necking phenomenon. The local integration of the elastic-plastic equation set is carried out using the fully implicit integration scheme (backward Euler difference scheme), which ensures both the stability and the symmetry of the consistent tangent modulus. The numerical computations of necking are carried out both on shear free end bars including an initial small geometrical defect and geometrically perfect bars with gripped ends. The corresponding results are then discussed in some detail and compared with those obtained in the literature.

Theoretical background

symmetric too, see [START_REF] Green | A general theory of an elastic-plastic continuum[END_REF]], [START_REF] Casey | A remark on the use of the decomposition F=FFP in plasticity[END_REF]].

In a quasi-static process, the equilibrium equation reads (3) where F is the deformation gradient, I: the second Fiola-Kirchhoff stress tensor (symmetric), Po the reference mass density and f 0 the body force. The product II= is the first Piola-Kirchhoff stress tensor, it is not symmet---T ric whereas II.F is.

The total Lagrangian formulation presented herein is based on thermodynamic foundations developed for the so-called generalized standard materials [START_REF] Halphen | Sur les materiaux standard generalises[END_REF]], [Germain, [START_REF] Nguyen | Continuum thermodynamics[END_REF]]. Within this framework, we consider that the state of any particle and at any time is characterized by the following state variables : the total strain E, fue plastic strain EP and another internal variable to account for the deformation history, represented by the scalar hardening variable denoted by a.

To complete the equation set for the time-independent elastic-plastic problem, we shall provide the state laws together with the ~volution laws. The state laws give the symmetric stress 'f: and the thermodynamic conjugate A, dual of the hardening variable a, by differentiating the free energy per unit volume w Let us consider a three-dimensional continuum body which at initial time t 0 occupies a region no bounded by surface S 0 and assume that the body undergoes a quasistatic elastic-plastic transformation. The initial configuration no will be chosen as the reference one and the r; equations of the problem expressed using Lagrangian variables relative to this configuration. As usual, the stored-energy w is assumed to be the sum of the elastic energy we -function of the elastic strain The motion of the body in the Lagrangian description is defined by the displacement field denoted by U. The deformation inside the body is measured by the Green strain tensor defined in terms of the displacement According to [START_REF] Green | A general theory of an elastic-plastic continuum[END_REF]], [Green A. E.;[START_REF] Green | Some remarks on elastic-plastic deformation at finite strain[END_REF]], the strain tensor E is decomposed additively as where the so-called plastic strain EP is assumed to be =e symmetric and objective, so that the elastic strain E is =e E only -and the hardening energy wa -function of the hardening variable a

(5)

For metals, we can assume small elastic strains and choose the elastic energy as a quadratic function of the elastic strain (6) where f> = 211i + AJ denotes the isotropic elasticity tensor (i is the fourth-order identity tensor, l = ei 0 e j 0 e j 0 ei ' J = I 01. I the second order identity tensor). In order to describe nonlinear hardening effects, the hardening energy is chosen as a function containing exponential terms 1 2 ( exp(-c3a))

wa.(a) =coa+ 2 cza +(ct-co) a+ c 3 (7) where the coefficients ci are material-dependant constants.

The state laws then yield the stress-strain relation and the nonlinear relation between hardening variables

Numerical implementation

In order to numerically solve the above stated elastoplastic problem in the context of the finite element method, it is more convenient to replace the equilibrium equation (3) by the virtual power principle \fU* ' =T :=__,.

_, _,. where Sis the deviator of stress I: and 0' 0 the initial yield stress.

The evolution laws result from the generalized maximum plastic dissipation rule and provide the plastic strain and hardening variable rates by differentiating the yield function with respect to its variables ~p= af -a=}..! (10) where the pl~tic multiplier }.. satisfies }.. ;:::: 0 if f(I:,A) = 0 and : I: 2:: 0, otherwise ~0. Of course, there is no relation between the Lame constant A and the plastic multiplier }... Relation (1 Oa) is the plastic flow rule and relation (lOb) the hardening law. They identically fulfill the second law of thermodynamics for continua, written in terms of objective Lagrangian quantities.

Since the plastic strain is a Lagrangian variable, its rate is obtained by ordinary differentiation with respect to kinematic time t. There follows from relations ( 9)-( 10) that the hardening rate is equal to the Lagrangian equivalent plastic strain rate

• Jz=P =P

a=P= 3 E :E (11)
Relations ( 1)-( 3), ( 8)-( 10) form a coupled system of 30 scalar partial differential equations with 30 scalar unknowns U, Ee, EP, I:, A, a, and i.

The discrete formulation is obtained using the follow-~g vector and matrix notations. A second-order tensor T (symmetric or not) is represented by a column vector {T} of nine components as follows {Tf {Tt, Tz, T3, T4, Ts, T6, T7, Ts, T9}r = {Tu, Tzt, T31, T12, Tzz, T3z, T13, Tz3, T33}r

Thus, the ~omponents of vector {T} are related to those of tensor T by {T}k = T pr H k=3(r-1 )+p. It should be noted that second order tensors are represented by vectors with nine components instead of six as usually done.

In fact, this representation mode, which is necessary for unsymmetric matrices such as the displacement gradient H or the stress TI, does not entail extra computational cost, as we shall see below. Similarly, in matrix notation, a fourth-order tensor Tis represented by a (9x9) matrix In the above relation, vector {H} has 9 components, {UYis the vector of element nodal displacements with 3nne components, matrix [G] is of dimension (9x3nne) and defined by where Ni are the usual shape functions and l3 denotes the unit matrix of the 3rd order.

Note that matrix [G] does not depend on the nodal displacements. Moreover, it has the advantage of being sparse and regularly populated. This allows us to compute the internal force vector as well as the tangent matrix in a simple and efficient way, as shown in proposition 2 below.

The discretization and consistent linearization of the elastic-plastic equation set lead to a nonlinear matrix system which can be solved by a incremental procedure. The convergence at every time step requires a Newton type iterative scheme. Within each iteration, the local integration of the elastoplastic equations ~et (8)-(1 0) is carried out in order to compute the stress II and the consistent tangent ten~or arr;a'F corresponding to a given strain increment L\E. Then we can build up the element internal force vector and the element tangent matrix defined as

{'PY = 1raY {IT}dno (17) [Kr = ~~~F = 1 raYr~~][G]dno (18)
The local integration procedure is formally the same as in small strains, the rate equations are discretized using the fully implicit integration scheme (backward Euler scheme) and we get the correction for the plastic multiplier at the k-th local iteration [START_REF] Simo | Consistent tangent operator for rate-independent plasticity[END_REF]]

(3 =El l y 2 S -O'o -An-I o/Jk) = -L\/Jk-1) + (19) 3p+h =E
where h=a 2 wa/aa 2 , J1 is the shear modulus, S the deviator of the elastic trial stress, =f = 2p ( dev (E) -E:_ 1 ) •

The subscript (n-1) denotes known quantities at the previous time step. The plastic multiplier increment L\A. obtained at convergence allows us to ~p~ate the unsymmetric stress II by using relation II F.E and (20) On the other hand, the solution convergence at each iteration requires computing the consistent tangent modulus ai:;a:E Eventually, the consistent tangent tensor arrjaF entering in the construction of global stru~ur~ tangent matrix is derived from the above tensor ai:;a:E by means of the following proposition.

=E E = "!-21l{z"x ~~~~~ + ~ (3A+ 2p) (trE) l

Proposition 1

By using an appropriate change of variables, the following symmetry cond~on ~an always be satisfied for the expression of stress E: VT (not necessarily symmetric), E (T, state at tn-t) = E (TT ,state at tn-1)

(22) Then one has an a'F ( -=)T =a'!: =T LE +F.aE.F (23)
It can be shown that the above discretization procedure leads to existence of a discrete potential w(F, stateattn-t}_su:_hthatiTT =aw(F) jaF. We deduce that tensor arr;a'F possesses themirror symmetry:

( arr;a'F) = (arr;a'F) and that its matrix represenijkl lkji tation [aiijaF] is symmetric. From (18) there follows the symmetry of the tangent consistent modulus.

To end this section, let us note that one can take advantage of the sparse and regular nature of matrix [G], Relation ( 16), in order to prove the following proposition and where implicit summations are made over r, sE{1,2,3}, whereas integers a, bE{l, ... ,3nne}, and p, qE{l,2,3}are expressed in terms of integers i, j by a=Int(i~2 ), b=Int(jj 2 ), p=i-3(a-1), q=j-3(b-l) (Int=lnteger part of)

Relations ( 24) and ( 25) show that computing the element internal force vector and the element tangent stiffness matrix is not so tremendous a task.

4 Numerical results

Geometry and boundary conditions

The numerical computations are performed on a circular cylindrical bar of initial length 2£ 0 =48mm and of initial radius ro=4mm.

Concerning the boundary conditions, the bar is modeled in two ways:

{i) The first one is most commonly adopted in the literature. During the loading process, the ends of the bar remain shear free so that the radial displacement thereat is not prevented. To trigger the initiation of necking, an axisymmetric geometric imperfection is introduced in shape of a full cosine wave with an amplitude of 1% of the initial radius.

(ii) Another way to compute the necking behavior is to consider a geometrically perfect bar with its ends cemented to rigid grips, as in [START_REF] Needleman | A numerical study of necking in circular cylindrical bars[END_REF]] and [Brtinig (1998)]. Note that in this case there is no need of a geometrical defect for the onset of necking.

Material properties

The Young modulus is E=2.10 11 Pa, the Poisson ratio v=0.3 and the initial yield stress cr 0 =400MPa. The coefficients of the hardening law (8b) are so adjusted that a. Meshno. 1 b. Mesh no. 2 ( 4 x 125 elements) ( 4 x 400 elements)

Figure 1 : Perfect bar with gripped ends. Initial and final meshes.

the numerically derived load-displacement curve best approximates the experimental curve of a given standard steel (XC48C) : co=O (no initial hardening), c 1 =220MPa, c2=-560MPa and c3=15.

Meshes

Because of obvious symmetry, we only have to consider one quarter of the specimen the length of which is fo=24mm. For each boundary condition type, the numerical computations are carried out using two axisymmetric meshes (Figure 1) :

(i) A coarse mesh with 125 axisymmetric eight-node quadrilaterals.

(ii) In order to check the sensitivity of the numerical results to mesh refinement, another computation is performed over a finer mesh with 400 elements, figure lb.

Numerical results in the case of the shear free end bar with imperfection

The bar is subjected to axial end displacements corresponding to a maximum relative elongation of 23% attained after 25 time steps and Figure 1 shows the uncleformed and final configurations.

With the 125-and 400-element meshes, on average 4 iterations are necessary to satisfy the convergence tolerance of l.e-6, with a peak of 8 iterations just after going through the limit point on the load-displacement curve shown in Figure 2. As a matter of fact, the plots corresponding to the coarse and the finer meshes are found to be indistinguishable, thus validating the insensitivity of the numerical results to mesh refinement.

The load-displacement curve representing the tensile load F versus the axial end displacements U B is shown in Figure 2. The tensile load, after reaching its maximum value equal to 1.48 times the elastic yield force, decreases when the necking takes place. For a relative elongation of 23%, it falls down to 76% of its maximum value.

For comparison, we also perform the computations over a perfectly cylindrical bar. As expected, this bar remains cylindrical during the whole deformation process and the necking never takes place. The corresponding tensile force plotted in Figure 2 attains its maximum value later than the imperfect bar and very slightly decreases afterwards. We also find that bifurcation occurs about the maximum load point, as in [START_REF] Needleman | A numerical study of necking in circular cylindrical bars[END_REF]] and [START_REF] Briinig | Numerical analysis and modeling of large deformation and necking behavior of tensile specimens[END_REF]], although the precise location of bifurcation point is not known. It is noteworthy that the maximum load point (1.48 times the elastic yield force) and the bifurcation point (engineering strain of 16% approximately) are very close to those obtained by [START_REF] Okada | A full tangent stiffness field-boundary-element formulation for geometric and material non-linear problems of solid mechanics[END_REF]] for the necking of a tensile plate in plane strain. However, beyond a certain elongation the radial displacement becomes suddenly more pronounced in the central region of the bar, and eventually the high concentration of strain there leads to necking phenomena.

The curve in Figure 3 giving the central cross-section radius versus elongation is tangent as expected to that of the geometrically perfect bar, which is virtually a straight line.

Beyond an engineering axial strain of 15% the former rapidly deviates from the latter, and at an elongation of 23% the current radius drops to about 36% of its initial value. This behaviour is very close to the numerical and experimental results depicted by [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. Part I : Continuum formulation[END_REF]] although the material considered and particularly the formulations adopted in the previous reference and the present work are quite different. The curve shown in Figure 3 can also be compared with that computed by [START_REF] Briinig | Numerical analysis and modeling of large deformation and necking behavior of tensile specimens[END_REF]], yet noting that the last Reference gives the change in crosssection area instead of the change in cross-section radius.

Figure 4 shows the evolution of the plastic region at different deformation stages of the imperfect shear free end bar. The plastic deformation begins in the central region when the axial elongation UB/£ 0 is equal to ojE=0.2%

(Figure 4a), and it rapidly sweeps through the whole bar afterwards. Up to an elongation of 15% -i.e. until the maximum load level (see Figure 1) the bar still remains entirely plastic (Figure 4b). Yet immediately after-at 16% elongation -there suddenly occurs a remarkable de- crease of the plastic zone (Figure 4c ). For larger extensions of the bar, bifurcation has occurred and unloading proceeds rapidly. At the final elongation of23% only the central necked region of the bar remains plastic whereas all other parts are unloaded elastically (Figure 4d).

The changes in the shape of the plastic zone observed above agree qualitatively well with those depicted in [START_REF] Briinig | Numerical analysis and modeling of large deformation and necking behavior of tensile specimens[END_REF]]. However, in the last Reference the plastic zone decreases more rapidly, for instance at 20% elongation of the bar it is virtually concentrated in the central 1/1 Oth of the bar length. Elastically unloaded regions are also detected in the necking of a tensile plate in plane strain, as shown in [START_REF] Atluri | New general and complementary energy theorems, finite strain, rate sensitive inelasticity and finite elements : some computational studies[END_REF]], and their spread is quite similar to that depicted in Figure 4.

The maximum value of the Lagrangian equivalent plastic strain P (see Equation ( 11)) at the final elongation of 23% amounts approximately to 87% and is found on the lateral surface, in the immediate vicinity of the central cross-section of the bar. At the same time, the minimum value of Pis about 13%, at a point on the circumference of the end sections.

Numerical results in the case of the gripped end perfect bar

So far the ends of the bar are shear free and the onset of necking is due to an initial small geometrical defect, as most commonly modeled in the literature. It should be noted that [START_REF] Needleman | A numerical study of necking in circular cylindrical bars[END_REF]] and [Bmig (1998)] also considered another manner to simulate the necking, by taking a geometrically perfect bar with its ends cemented to rigid grips. For comparison, we take up again the numerical computations with these new geometry and boundary conditions.

With the coarse mesh (125 elements) the maximal engineering strain of 23% requires 25 time steps, whereas with the finer mesh (400 elements) it requires one more step, just to get across the limit point when the load reaches its maximum value. Apart from this slight difference in computational effort, all the curves resulting from the two meshes are quite identical.

Figure 5 shows the undeformed and final configurations.

Here the gripped end boundary conditions act as an alternative geometrical imperfection and allows the necking bifurcation from the uniform tensile stress state. Despite the evident difference in the deformed shapes, the loaddisplacement curves shown in Figure 2 are almost the same for the gripped and free end bars.

In contrast, Figure 3 shows that the change in the crosssectional radius in the necked region is rather different : in the case of the perfect bar with gripped ends, the change in the central cross-section radius is slower on the end part of the curve. At an elongation of 23%, the radius is equal to 44% of its initial value instead of 36% as in the case of shear free ends.

The necking radius curves in Figure 3 are different from those given in [START_REF] Needleman | A numerical study of necking in circular cylindrical bars[END_REF]]. However, they compare very well with those in [Brilnig (1998)]. Note that in the last Reference another geometry -more realistic -was also studied, which is a whole cylindrical tensile specimen with two shear free thick ends. This geometry is not considered in the present work, but according to [Brilnig (1998)] the results of the griped end perfect bar must be identical to those of the whole specimen.

Figure 6 shows that the evolution of the plastic zones for the gripped end perfect bar differs notably from that of the shear free imperfect bar, especially at the beginning.

When the axial elongation DB/fo is equal to croiE=0.2%, there exist separate plastic zones of different sizes in the gripped end bar. The small plastic zones, hardly visible in Figure 6a, are situated on the end cross-sections and near the lateral surface. The much larger plastic zones lie around the axis, between the central cross-section and the ends of the bar. As the overall extension increases, these plastic zones grow up together and eventually becomes connected. When the maximum load is reached -i.e. at about 15% elongation (see Figure 1) -the entire bar is plastically loaded, except for two small regions near the ends that will remain elastic throughout the whole range of straining (Figure 6b ). Just after, with a 1% additional elongation only, unloading takes place and then regularly spreads out (Figure 6c). At the final elongation of 23%, the plastic zone is reduced to the central region of the bar, over a volume slightly smaller than that of the shear free end imperfect bar (Figure 6d).

The development of plastic zones shown above is very close to that described in [START_REF] Needleman | A numerical study of necking in circular cylindrical bars[END_REF]]. In particular, as in the quoted Reference, we also observe separate plastic zones which become connected, and the existence of small regions near the ends which remain continuously elastic. In the case of the gripped end perfect bar considered here, less details are given in [Brilnig ( 1998)] for the beginning of the plastic zones evolution, however for the final stage all the results are again in quite good comparison.

The maximum value of the Lagrangian equivalent plastic strain P at the final elongation of 23% amounts approximately to 83% and is found on the lateral surface, a little bit farther from the central cross-section of the bar than in the case of the shear free end imperfect bar. Of course, the equivalent plastic strain P remains zero at the few points near the ends that remain elastic throughout.

Conclusions

It is the purpose of this work to show the ability of the proposed formulation of finite elastoplasticity to deal with necking bifurcation problems. Below are summarized the main features of the theoretical and numerical formulations. Besides, further comments are also added on the results obtained in the preceding sections.

(i) The proposed theory is based on thermodynamic foundations established for generalized standard materials and results in a total Lagrangian formulation for finite plasticity. Since all relations are expressed in terms of Lagrangian variables, the material frame indifference is automatically satisfied within this formulation. The deformation history in the material has been characterized by a nonlinear relation of the saturation type between the isotropic hardening variable and its conjugate, Equation (8b).

(ii) The numerical implementation using standard finite element analysis is simple and yet robust enough to accurately reproduce the necking phe~omenon. Written in terms of the Green strain tensor E and the symmetric Piola-Kirchhoff stress tensor 1:, the local integration procedure is formally the same as in small strains. The rate equations are discretized using the fully implicit integration scheme (backward Euler scheme), which ensures both the stability and the symmetry of the consistent tangent modulus. In computing the global tangent matrix, the update of matrix [di1/dF] is made by means of Relations ( 21) and ( 23).

On the other hand, Relations ( 24) and ( 25) show that computing the element internal force vector and the element tangent stiffness matrix is not a highly timeconsuming operation.

(iii) The numerical computations of necking are carried out both on imperfect bars with shear free ends and geometrically perfect bars with rigid gripped ends.

For each boundary condition type, the numerical computations carried out using a coarse and a finer meshes are found to yield indistinguishable results with the same rate of solution convergence, thus validating the insensitivity of the numerical results to mesh refinement.

The robustness of the implemented solution procedure is successfully assessed since 23% elongation can be attained in only 25 time steps, giving sensible necking simulation of the tensile bar.

All the obtained results have been quantitatively and qualitatively analyzed in some detail. Comparisons are also made with other formulations in the literature, showing agreements or differences in the load-displacement relationship, the necking radius versus axial elongation, and the evolution of plastic zones during the straining history.

(iv) It should be noticed that necking is usually referred to as a bifurcation phenomenon in the literature, whereas a true bifurcation computation of bar necking-as done in [START_REF] Needleman | A numerical study of necking in circular cylindrical bars[END_REF]] -is rare. [START_REF] Needleman | A numerical study of necking in circular cylindrical bars[END_REF]] considered a geometrically perfect bar with shear free ends and performed a real bifurcation analysis by discretizing a variational principle of Hill with the finite element method. As commonly done in other works, here we have solved the necking as a standard deformation problem by either introducing an initial imperfection of 1% in the radius or preventing the radial end displacement.

In connection with bifurcation aspects, let us note an interesting issue not often considered, which concerns the sensitivity to geometric or material parameters. In the present work, the numerical computations have been carried out with a unique set of bar dimensions and material values. The sensitivity study requires more extended discussions and is beyond the scope of this work.

(v) The computer program used here can be improved with path-following techniques such as arc-length and branch-switching methods, in order to deal with more general problems in finite displacements and plastic strains, including limit point instabilities and bifurcation points. First results due to these enhancements have been obtained for the pinched cylinder modeled with shell elements in finite plasticity [Le Grognec and Levan (2000)]. Further investigations to treat static plastic buckling of shell-type structures are being developed.

1

  Laboratoire de Mecanique et Materiaux, Division Mecanique des Structures, Ecole Centrale de Nantes. 1, rue de Ia Noe, BP 92101 -44321 Nantes cedex 3, France

  fn.o TI : gradV*dQo-fn.o f 0 V*dQo -fso U* .TI.NdS0 = 0 (12) (Sa) The nominal stress vector T = TI.Nrelated to the reference normal Nis known on the part of boundary where (Sb) dead loads are prescribed. The elastic and plastic ranges at any particle in the body are defined by the yield function f _?f von Mises type, depending on the symmetric stress I: and the hardening variable

  components of matrix [1r] are related to those of tensor T by Tkt=Tprqs H k=3(r-l)+p, £=3(sl)+q. Within an element with nne nodes (nne = node number of the element), the displacement gradient H = gradUis approximated by {H} = [G]{UY (15)

  remarkably reduce the operation number for the computation of the element internal force vector {'PY and the element tangent stiffness matrix[KY , Relations ( 17

Figure 2

 2 Figure 2 : Load-displacement curves

Figure 3 :

 3 Figure 3 : Necking radius vs. axial elongation

Figure 4 :

 4 Figure 4 : Geometrically imperfect bar with shear free Figure 5 : Perfect bar with gripped ends. Initial and final ends. Evolution of plastic regions (shaded area) meshes.

Figure 6 :

 6 Figure 6 : Perfect bar with gripped ends. Evolution of plastic regions (shaded area)


Valérie Gabelica
email: valerie.gabelica@inserm.fr

A Pilgrim's Guide to G-quadruplex Nucleic Acid Folding

Every day, when I bike back home a little southwards of Bordeaux, I pass by a discrete post symbolizing "The Way of St James". If I would bike, from post to post, for a thousand more kilometers, I would finally arrive in Santiago de Compostela in Spain, going down a pathway followed by pilgrims from medieval times to this day. Translated to French, "The Way of St James" becomes "Les Chemins de St Jacques": it's a shift from the singular to the plural. And the routes to Compostela are indeed many, as pilgrims were coming from all over Europe to converge to their final destination. It is, among other significations, the convergence of these multiple routes towards a final goal that is symbolized by the shell of St James. Of course, depending on where you were coming from, there are preferred routes, nodes where several routes converge, and recommended landmarks to stop over for a rest, some of which being located just before major difficulties such as the crossing of the mountains between France and Spain.

Biomolecule folding routes are a lot like The Way of St James. The goal -the final state -is the folded structure. But the map -the free energy landscape -shows many different routes, converging to the final goal like a funnel. Protein folding pathways [START_REF] Onuchic | Protein folding funnels: the nature of the transition state ensemble[END_REF], and now RNA folding pathways [START_REF] Thirumalai | RNA and protein folding: Common themes and variations[END_REF], are nowadays usually described in this framework of folding funnels. The folding funnel concept is not to be opposed with the concept of "intermediates". However, for biopolymers, intermediates should not be viewed as single states (as a chemist would depict chemical reaction intermediates), but as ensembles of states sharing structural and kinetic properties. Intermediates are regions of the free energy landscape where a fraction of the molecular population gathers for a while before crossing a potential energy barrier.

Over the years, DNA folding had attracted less attention than protein and RNA folding, until other folding motifs that the Watson-Crick double helix came under focus. For example, G-quadruplex structures are formed by guanine-rich strands thanks to hydrogen bonding between four guanines that form base tetrads (the G-quartets), and thanks to cation-assisted stacking of successive Gquartets. Importantly, G-quadruplexes can interfere with replication, transcription and genome maintenance. However, for a given G-rich sequence, it remains difficult to predict which guanines form G-quartets with which ones and, with a few exceptions [START_REF] Da Silva | Design of a G-quadruplex topology through glycosidic bond angles[END_REF], to predict the resulting strand topology.

In particular, the human telomeric sequence, consisting of (TTAGGG) repeats, one of the earliest and most studied one given its importance for cell division and cancer, turned out to be one of the most compound DNA folding problems, even for oligonucleotides containing barely two dozens of nucleotides. To this day, no less than six intramolecular folding topologies (as distinguished by relative strand orientation, loop placement and syn/anti conformation of the guanines, see Figure 1, south of the map) have been solved at high resolution. Which of these structures is the most stable in long telomeric repeats in vivo is still debated. Nor is it known whether the thermodynamically most stable structure is the relevant one in the cellular environment, or whether kinetically controlled structures could prevail. As a consequence, it is important to explore not only the structural and thermodynamic, but also the kinetic aspects of G-quadruplex folding. In other words, it is important to learn something about all plausible forms of G-rich sequences on their folding pathways towards a G-quadruplex. We need to learn more about the North part of the map as well. The six high resolution stable structures of 4-repeat human telomeric G-quadruplexes reported to date. Each guanine base is represented by a rectangle, with dark red/light blue color code representing the syn/anti conformation, respectively. Four bases in a square plane form a G-quartet. Structures A-E have three G-quartets and structure F has two G-quartets. North of the map: some possible stopovers explored by unbiased molecular dynamics: slipped structures were found plausible points of interest for the fully parallel, all-anti structure A whereas they were found impossible for hybrid structures due to steric clash between bases (see for example strand slippage from structure C) [START_REF] Stadlbauer | Structural dynamics of possible late-stage intermediates in folding of quadruplex DNA studied by molecular simulations[END_REF]. In contrast, triplex structures were found plausible for structures containing both syn and anti bases (see for example the triplex of C), but unstable for the parallel, all-anti structure derived from structure A [START_REF] Stadlbauer | Triplex intermediates in folding of human telomeric quadruplexes probed by microsecond-scale molecular dynamics simulations[END_REF]. The map of the Way of St James and the shell-shaped logo were reproduced from Wikimedia Commons [START_REF] Zentgraf | Map of the way of St James In Europe[END_REF].

There are several ways to explore the folding landscape and find points of interest, where the molecules would stop for different amounts of time (a pleasant meadow in a basin, a safe shelter before a difficult mountain pass, a detour to a dead-end maybe), or definitely (the final destination). The approach that made possible the characterization of six folded structures was to modulate slightly the base sequence or the environment (nature of monovalent cations, crowding conditions for X-ray cristallography), so as to slightly tilt the free energy landscape and drive nearly 100% of the population to the same end point. In such conditions, atomic structural investigation by NMR or Xray crystallography is possible. Of course, whether these tilted landscapes represent the landscape really at stake in vivo is a matter of debate. Now let's explore the folding routes at some further distance from the final destination. A first possibility is to start from the final destination and then tilt the landscape a little more drastically to see where the population ends up. This is done in denaturation or unfolding experiments: the free energy landscape is changed by increasing the temperature, by adding chemical denaturing agents, or by applying a force to separate the two ends of the molecule. It's important to realize that the route taken will depend on the method used to drag the molecule away from its final state, and that it is by exploring a sufficiently wide variety of roads that the map will be built.

A second possibility is to start from an unfolded ensemble in a non-native environment, and then initiate the folding by changing the environment to a native-like one. This is the equivalent of suddenly transforming the map from flat to 3D and let the pilgrims find their way. The molecular population coming from all over the free energy map is then monitored, and intermediates are distinguished based both on structural properties (region of the conformational landscape) and kinetic properties (one is more likely to identify as intermediates those kinetic bottlenecks on the folding pathway, like are like the valleys at the bottom of mountain passes.

The most challenging aspect is to characterize the intermediates that have been distinguished based on kinetic and/or physical (i.e. spectral, hydrodynamic, or distance) properties. For structural ensembles that are neither totally unfolded nor totally stable, it's almost impossible to obtain structural information with atomistic detail. Among the five topologies fully characterized by NMR to date, the one having two G-quartets instead of three attracts attention as a possible intermediate: not each and every guanine has found its partners, and one G-tract is slipped by one step (see Figure 1F). However strand slippage involves breaking and reforming G-quartets. Also, not all slippage mechanisms are sterically allowed, because the guanine syn-anti patterns in the quartets before and after the strand slippage should remain unchanged, but must remain compatible mutually with other bases of the new quartets formed. As a result, strand slippage is kinetically hampered for all structures except the fully parallel, all-anti structure A. Therefore, such stable structures with slipped strands can also be dead ends on the folding route, and are not necessarily possible to reach on an unfolding route.

A faster folding/unfolding route would be the stepwise assembly of quartets strand per strand, a route on which partial dimers containing G-G base pairs and partial trimers containing G-triplets would be plausible intermediates. Unfortunately the structural interpretation of physical properties monitored in folding/unfolding studies is not straightforward. Force-extension distributions in singlemolecule unfolding are compatible with an intermediate having one G-tract removed from the core, which would remain a so-called G-triplex [START_REF] Koirala | Intramolecular folding in three tandem guanine repeats of human telomeric DNA[END_REF][START_REF] Mashimo | Folding pathways of human telomeric type-1 and type-2 G-quadruplex structures[END_REF]. Intermediates having a circular dichroism spectrum with a maximum at 265 nm and a minimum at 242 nm were identified both upon thermal unfolding [START_REF] Buscaglia | Thermodynamic characterization of human telomere quadruplex unfolding[END_REF] and K + -jump folding [START_REF] Gray | Folding and unfolding pathways of the human telomeric Gquadruplex[END_REF] were also assigned as G-triplexes, based on the similarity with the CD spectrum of a 3-repeat telomeric sequence [START_REF] Koirala | Intramolecular folding in three tandem guanine repeats of human telomeric DNA[END_REF]. Intriguingly, these CD spectra are identical to those expected for fully parallel G-quadruplexes with homo-stacking between successive guanines [START_REF] Karsisiotis | Topological characterization of nucleic acid G-quadruplexes by UV absorption and circular dichroism[END_REF]. This intriguing G-triplex therefore raises the following questions: What is (are) their structure(s)? (e.g., are they parallel-stranded or antiparallel-stranded?) Are they stable for a sufficient amount of time to be called intermediates? In other words, are G-triplexes local minima on the folding routes of telomeric DNA and if so, where are these points exactly on the conformational map?

To answer these questions, molecular modeling must come into play. A state-of-the-art molecular mechanics (MM) force field is like a faithful local 3D map for the molecule, which will wander around valleys and hills, with more or less energy depending on its temperature. Also, force fields describe the molecular structure with atomistic resolution, providing exquisite details on the structure and how it changes over time. Molecular dynamics (MD) simulations allow one to sample a trajectory of the molecule in its force field. The drawback is that limitations in computer power do limit the time given for the molecule to explore its landscape, and hence the explored area.

There are ways to expand the sampling (the explored zone on the landscape), using for example metadynamics or replica exchange between multiple trajectories at different temperatures. The first atomistic structure of a G-triplex was provided by a 80-ns metadynamics MD simulation of the folding of the thrombin binding aptamer sequence, a short DNA sequence forming an antiparallel structure with two G-quartets: in one of the intermediates the G-tract at the 3'-end protruded out in the solvent [START_REF] Limongelli | The G-triplex DNA[END_REF]. The stability of the rest of the structure, the G-triplex, was experimentally confirmed by NMR, CD, and thermal unfolding studies on the truncated sequence missing the 3'-end. A drawback of enhanced sampling methods is that the sampling is biased: it is like if the landscape was tilted so as to drive exploration further away from the starting point.

To allow the molecule to explore the natural landscape, unbiased sampling is required, and as a consequence, long simulations are needed to allow the molecule to explore the landscape in a meaningful way. In a previous work, Stadlbauer et al. investigated strand-slipped structures as latestage folding intermediates [START_REF] Stadlbauer | Structural dynamics of possible late-stage intermediates in folding of quadruplex DNA studied by molecular simulations[END_REF]. In a feature paper of the present issue [START_REF] Stadlbauer | Triplex intermediates in folding of human telomeric quadruplexes probed by microsecond-scale molecular dynamics simulations[END_REF], the authors report the longest to date unbiased sampling MD simulations of G-triplexes derived from the human telomeric sequence, with three tracts of three guanines each.

Importantly for unbiased sampling, the choice of the starting structures is crucial, and must be driven by a scientific question. Here nine different starting structures of three-triad telomeric sequences were tried out, corresponding to different loop placements, with the underlying question: which of these triplexes are stable on the microsecond timescale and which are not? The MD simulations revealed which proposed starting structures remain around their starting point, or on the contrary whether they move away and how quick they do so. This must be tested for a long enough time scale and repeated a few times for each starting point to check consistency. The present study reveals that several antiparallel triplexes (for example, C(triplex) on Figure 1) remained "stable" on the microsecond time scale (i.e. the postulated triplex starting structure was indeed close to a local minimum on the landscape, and molecules tended to stay around), whereas fully parallel triplexes (e.g., A(triplex) on Figure 1) were unstable on that time scale. The lessons learned are many. Molecular dynamics simulations help to validate intermediates on the folding energy landscape: the simulations determine whether they are points of interest by observing whether molecules placed in that state remain in this local minimum for a significant amount of time. Importantly, they also provide refined details at the atomistic level both on their structure and dynamics within these regions of interest. Also, interplay between theory and experiment is the key to untangle folding problems. Without modeling, structural interpretation of physical measurements is only speculation. And in current state of the art and hardware, without thoughtful speculations simulations would wander about the landscape with no sense of direction. Mutual understanding between specialists of both communities is therefore crucial, and gaining understanding in the advantages and limitations of state-of-the-art MD simulations is another very good reason for experimentalists to read the present feature article from Stadlbauer et al. [START_REF] Stadlbauer | Triplex intermediates in folding of human telomeric quadruplexes probed by microsecond-scale molecular dynamics simulations[END_REF].

Figure 1 :

 1 Figure 1: The folding pathways of G-quadruplexes are like the Way of St James: it's an extremely multi-pathways map with multiple possible stop-overs. South of the map (A-F): The six high resolution stable structures of 4-repeat human telomeric G-quadruplexes reported to date. Each guanine base is represented by a rectangle, with dark red/light blue color code representing the syn/anti conformation, respectively. Four bases in a square plane form a G-quartet. Structures A-E have three G-quartets and structure F has two G-quartets.North of the map: some possible stopovers explored by unbiased molecular dynamics: slipped structures were found plausible points of interest for the fully parallel, all-anti structure A whereas they were found impossible for hybrid structures due to steric clash between bases (see for example strand slippage from structure C)[START_REF] Stadlbauer | Structural dynamics of possible late-stage intermediates in folding of quadruplex DNA studied by molecular simulations[END_REF]. In contrast, triplex structures were found plausible for structures containing both syn and anti bases (see for example the triplex of C), but unstable for the parallel, all-anti structure derived from structure A[START_REF] Stadlbauer | Triplex intermediates in folding of human telomeric quadruplexes probed by microsecond-scale molecular dynamics simulations[END_REF]. The map of the Way of St James and the shell-shaped logo were reproduced from Wikimedia Commons[START_REF] Zentgraf | Map of the way of St James In Europe[END_REF].


A Le Van

J Royer

A NUMERICAL PROCEDURE FOR REGULARIZED INTEGRAL EQUATIONS FOR ELASTOSTATIC PROBLEMS OF FLAT CRACKS IN OPENING MODE

Abdract-The regularid integral equations for flat cracks in opening mode are solved numerically by means of a simple procedure using new interpolation functions. The convergence and the accmncy of the method are tested on problems of circular and elliptical cracks embedded in the infinite medium, which admit analytical solutions. Computational results ensure the validity of this new approach. then they are compared with those already obtained without reguk&ation by the authors.

INTRODUCTION

FEW TIUDIMENSIONAL crack problems admit an analytical solution. Treated in the inWte medium case, they mainly concern cracks of particular geometries and subjected to particular loadings. The most well-known among these problems is that of the penny-shaped crack under an axisymmetrical pressure [l]. More recently, the case of more general loadings is fully studied in refs [2] and [3]. Like the penny-shaped crack, the elliptical crack also attracts many authors' attention because of its particular geometry. The elliptical crack subjected to a uniform pressure is treated in ref. [4], whereas the case of a linear pressure is studied in ref.

[4 and the case of a polynomial pressure in ref. [6]. For elliptical cracks subjected to shear loadings, the main contributions may be found in ref. [7] and more recently in refs [8] and [9]. The rectangular crack, though not very realistic, is another crack geometry amenable to analytical solutions [lO, 111. For more general crack configurations, one often has to turn to numerical computations. Integral equations derived from the potential theory have the common feature of being singular in the sense of the principal value. The numerical treatment of the principal value is difhcult; a recent successful attempt to deal with this singularity in two dimensions can be found in ref. [12]. However, computations without l5articular treatment of the singular part are able to yield highly satisfactory results as shown in ref. [START_REF] Van | Boundary element analysis of an inte8ral equation for threbdhnmsional crack problems[END_REF].

The regularization task aims at lessening the principal value singularity, In theoretical respects, this yields weakly singular integral equations, which must be more tractable than principal value ones. In numerical aspects, the kernel of regularized equations is found to be lengthy and its discretization in the three-dimensional case remains for the time being inextricable. Recently, a numerical procedure based on the least squares global displacement gradient smoothing method was worked out in ref. [START_REF] Polch | solutions for flat cracks[END_REF] for flat cracks in opening mode.

The present paper describes a different approach to problems of flat cracks in opening mode, using new interpolation functions denoted by MC" and J The proposed numerical procedure has the advantage of being simple: similar to the non-regularixation method, the solution of the discretixed equation is carried out by the collocation method and the integration classically by Gaussian points.

DISCRETIZATION OF THE REGULARIZED INTEGRAL EQUATION

Consider a planar crack S embedded in the in6nite medium and subjected to an arbitrary internal pressure (Fig. 1). it is know!@ jl that under these conditions the opening mode is uncoupled from the other ones. The crack surface can be parametrized with Cartesian coordinates: (Implicit sums on Greek indices are made from 1 to 2.) Denoting by S, a neighbourhood of y, one can split up the surface integral in eq. ( 1) as follows: v = r A II, r is the vector tangent to the crack edge 8s.

Let us introduce the following notations for the discretixation of eq. ( 2). In a general way, the symbol ( ) relates to a row-vector, whereas ( } relates to a column-vector, and [] to a matrix.

The geometry transformation will be denoted: e. 3 t? +y E e, where e, is the reference element of the physical element e. The coordinates of point y are expressed in the usual way: Va E { 1,2}, y. = (iV(e)){y,}:, (N(t)) designates the shape function.

The Jacobian matrix of the geometry transformation is denoted by WOI = [ J&> = ~G,lD m.c, *)9 II * .

its determinant by either detJ(t) or D(y)/D(<), and its inverse by

The last three integrals are quite regular and thus require no particular investigation. For instance if use is made of isoparametric elements, one has for any element e lying in S\S,,,:

It has been shown in ref. [START_REF] Levan | Theoretical basis of regularized integral equations for elastostatic crack problems. 8%[END_REF] that the first integral in eq. ( 2) is ordinary improper. So its study is less critical than that of a principal vahae integral integral. However, for numerical purposes an adequate integration scheme remains to be devised in order to take into account the weakly singular behaviour of the type r -*+r, 0 < y' d 1. In this study, a numerical procedure will be proposed which makes it possible to solve this weak singularity. First we define from (N(t)), [.I(<)] and b(r)], the following five functions which will be referred to as "mod&d". The first of them was introduced in ref.

[17J and widely used in I&S [START_REF] Polch | solutions for flat cracks[END_REF] and [18]. For each couple of points (<, &,) in e, we will consider the poiar ~~i~~ (A, (3) of 4, related to centre &, defined by (Fig.

J(t) -det J(G) = AC!&&, 4, @)
The modified inverse of the Jacobian matrix Lf(&,, A, t?)] is such that:

w WI W 3
For simple elements, e.g. six-node triangles and eight-node quadrilaterals, the existence of functions (m) and (M(")), a E { 1,2}, can easily be verified, and their expressions are also readily available. Concerning the other modified functions, the following proposition shows how to determine them: ,1,8)) a E (1,2} In the first or last two terms in eq. ( 4b), t: and to are interchangeable. Clearly, d # det .?. The foregoing proposition, which can easily be shown, enables us to obtain systematically the expressions for [s], a and [I], provided the existence of (M@)), a E { 1,2}. Expressions for the mod&d function (fi) can be found in ref. [18]i for the sake of completeness, we give them in the Appendix, together with the new functions (M@)) and [j]. Now we divide the neighbourhood S,,, of the singular point into elements e 3 yo, and transform the improper integral in eq. ( 2) over such an element as follows:

Pl-OpOSltlOBl If there exist (M@)(&,
From eq. (3a): y, -y, = (W(0) -W(~o)))W = UWO, A, W{Y.zY Furthermore, the difference in the kernel in eq. ( 5) can be transformed, using eqs (3b) and (3e):

A&) -L&o) = ~j^,,<ro, 1, 0 w Lastly, by noting that: D({,, t2)/D(Iz, 0) = 1, the integral (5) can be rewritten hereafter where the square of polar radius A2 is now cancelled out:

This time the integration is to be carried over (A, f?) and eq. ( 7) involves only regular integrals which should yield no dBculties for the numerical purposes. It should be noticed that the common feature between the approaches in refs [START_REF] Polch | solutions for flat cracks[END_REF] and [18] and"the method herein propo&tf is ttiat regular integrals proceed merely from a numerical procedure. The kernel in eq. ( 5) actually remains (weakly) singular of the type r -2+r, 0 < y I; 1. Here the.use of the interpolation type functions has made it possible to force the kernel to behave as a regular function in eq. ( 7). In the following study, the discretixation of eq. ( 5) is carried out considering only two usual elements: the six-node triande and the eight-node quadrilateral. Nevertheless, the principle remains the same for any more general element.

Use of eight-node quadrilaterals

The integration over (A, 13) in the reference element e, is effected by dividing this element into subtriangles and by making another change of variables in each of these subtriangles. This operation will allow us to bring again the integration domain onto the square [ -1, + 11' (Fig. 3).

According to whether the antecedent I& of pointy, is located at a comer or on a side of element e,, the subtriangle number is 2 or 3. Case when to is located at a corner of element eo: Let point to be located at node i E { 1,3,5,7} of element e. (Fig. 4). There are two subtriangles and it can be checked for the expressions of the polar coordinates (A,@ and the Jacobian of the transformation (u,, z+) --) (n,(I):

Subtriangles of type 1:

Subtriangles of type 2: Case when I& is located at a mid-side of e,:

1 1 +vi = cos ( > i(l + 01) A 1 +v, =
The four possible subcases are shown in Fig. 5. The node number of to is i E {2,4,6,8}. Here there are three subtriangles. By putting 0, = tan-l(2), the following relations can easily be established:

Use of six-node triangles

Here use is made of the same technique of dividing element e, into subtriangles and making the change of variables [ -1, + 112 3 v = (v,, vl) + (A, 6). If co is at a comer (resp. at a mid-side) of e,, the number of subtriangles is 1 (resp. 2). For triangular elements, there is no general rule and the six possible cases are to be considered separately. We recall the notation 0, = tan-'(2).

Case of to at node 1 (Fig. 6a) x-00 D(v)

1 _ l/%1+ 4)

Discretized regularized equation

Eventually, the integral in eq. ( 5) assumes the following form: 84,

*(y)-%(yo) ay, II

The integration over the variables (v,, v2) is made in the v-squares [-1, + l]*. Since the kernels involved are regular, the Gaussian scheme for numerical integrations is quite acceptable. Rgure 7 presents the successive steps to pass from the former square. [ -1, + 112 3 (v,, q) to the final physical element in which the point y = (yl, y2) lies. 9

Equation ( 8) together with eq. ( 2) constitutes the discretized regularized equation for a flat crack in opening mode. As remarked above, the weak singularity in eq. ( 2) has been removed with the aid of functions fi, Mea) and J ; of the numerical scheme there remains only <regular integrals.

NUMRRICAL PARTS Convergence test: the circular crack problem

Consider the problem of a circular crack centred at 0, with radius a, imbedded in the infinite medium and subjected to a uniform internal pressure p.

Like any flat crack in the infinite medium, the opening mode is uncoupled from the other modes. Four crack mesh configurations are chosen, containing respectively 85, 161, 209 and 321 interior nodes; the first and the last of them are presented in Fig. 8. Use is made of both six-node triangular and eight-node quadrilateral elements; those adjacent to the crack edge are quarter-point elements [l9,20]. All elements have straight edges; apart from the comer and quarter-point nodes the others are located at the mid-side of the edges.

For elements not containing the collocation point (i.e. lying in S\S,, see notations in eq. ( 2)), the integration is performed with 7 points over triangles and 9 points over quadrilaterals. Integration over elements containing the collocation point must be carried out using the change of variables which leads to a square integration region in the plane (a,, Q) (Fig. 7). The Gaussian quadrature scheme requires 3 x 3 = 9 points. Here the number of Gaussian points is twice as high as that used in the same study without regularization[l3]. This is accounted for by the fact that the original physical element is more "distorted" through the transformation (see Fig. 7) leading to the square element in the v-plane.

As expected, the numerical results show a revolution symmetry for the values of the displacement discontinuity & (Fig. 8). It should be noticed that the Young modulus E datum is unnecessary for computing the stress intensity factor K,. On the contrary, the Poisson ratio v must be specified; here it is set to be equal to 0.3.

Table 1 shows minimal and maximal errors on the displacement discontinuity & as well as errors on the stress intensity factor K,, with respect to the well-known analytical solution[l]: As remarked in the study of notched strips in ref.

[l9], although the stiffness of quarter-point elements is reasonably accurate, the local values of stress and displacement in these elements are poor, thus our error estimations on & do not take into account quarter-point node values. Furthermore, the &-error is not maximal at the crack centre but at some radial distance varying according to the mesh configuration; in present cases it remains about 5% in magnitude. As regards tie &error, it markedly decreases as the mesh becomes denser. On the strength of these results, another investigation concerning an elliptical crack is treated in the next part.

Elliptical cracks

Consider an elliptical crack of semi-axes a, b, imbedded in the infinite medium and subjected to a uniform internal pressure p. Figure 9 shows the meshes for two cracks with semi-axes ratio a/b respectively equal to 1.5 and 2; they include 161 (resp. 113) interior nodes and 64 (resp. 48) elements. Here use is also made of quarter-point elements along the crack edge. As predicted, the values of the stress intensity factor at vertex B are always worse than those at vertex A, since the mesh is not smooth enough at 3, especially in the second mesh. The maximal &-error in the last mesh occurs at a point of the major axis, close to vertex B; it undoubtedly contributes to the large error on &(A). In order to improve the results, more sophisticated mesh configurations must be devised, which result in an adequate refinement at the crack front and at the same time less elements at the crack centre where there is a small displacement gradient.

CLOSURE

The results show that the regularization method allows a high convergence and is quite operational for plane crack problems. Nevertheless, when computing rectangular cracks, for unknown reasons the results obtained by the authors are too poor to be published here. For circular cracks, the errors on K, are quickly reduced to less than 1%. For elliptical cracks, despite a coarse mesh and a bad crack contour representation, &-values remain quite acceptable where they are larger. Obviously, the circular mesh type (with radial rays) leads to too many elements at the crack centre without enough elements along the crack front. Thus more sophisticated mesh generation techniques are required to obtain better results.

Moreover, compared to the classical transformation (c,, &) + (y,, y2), the transformation (tr,, oz)-Q,, yz) entails a greater distortion of the physical element and implies a denser integration scheme as a compensation. The Gaussian point number used here is twice or four times as high as that for the method without regularization (9 or 16 points instead of 4). The above distortion may also account for the faet that the elements seem to be very sensitive to their aspect ratio (ratio of their longer to shorter side): the longer elements are, the worse the results. In connection with the non-regularization method, there appear two questions which call for answers.

Can one avoid using quarter-point elements? Indeed, the use of quarter-point elements condemned us to the range of six-node triangles and eight-node quadrilaterals. The non-regularization method, as shown in ref.

[l3], gave good results without using special elements.

The regularized eq. ( 8) was obtained and used in Cartesian coordinates. Is it possible to obtain a regularized equation in other coordinate systems? With reference to the discretized equation for the non-regularization method, it enabled us to cope with more general coordinates and in particular to compute the stress intensity factor on a local basis.

Under present conditions, the results obtained with regularization are encouraging and should encourage investigations with an improved mesh generation technique, and in a further stage to numerical developments with the view to studying three-dimensional cracks.

1 Fig. 1 .

 11 Fig. 1. Planar crack. y,,: singular point; as: crack edge.

Fig. 2 .

 2 Fig. 2. Polar coordinatea (A, 0) in the nferenoc element e,. Case of a quadrilateral element.

  2): r,=&J,+acose & = C& + 1 sin 8. The modified shape function (i?(&,, R, 0)) is such that: wm -<Wd) = w%, 5 8)) W and: lim (&(&,, A, 0)) < co. (The notation A < 00 signifks that A is bounded if A is a real; that all ~rn~nen~ (resp. elements) of A are bounded if A is a vector (resp. a matrix)). The modified derivated shape functions (A@)(&,, 1, O)), a E (1, 2}, are such that: and: lim (M(*)(&,, A, 6)) < co. The modified Jacobian matrix [3(&,, A, O)] is such that: [X)1 -mill = 43tem A ell and: lim[3go, A, e)] < Q). The modified Jacobian determinant h(&,, A, 0) is such that: and: lim &too, A, e) < ~0.

  det

  , then there also exist the functions [j(&,, rl, e)], d(&,, rl, a), u(&,, A,@] and they ark obtained by: Va, B E {1,2), &9(&, A,@ = W%0, A, W{~P>~ #a) &%,J, 0) = U2cr> +&&(r,) -&<r> -&J&) (4b) where (y, S) is equal to (2,2) if (a, /I) = (1, l), to (1, 1) if (a, 8) = (2,2), and to (a, jl) otherwise.

Fig

  Fig. 3. Change of variables square [-1, + I]* 3 t, = (u,, 4) -+ (A, 0) 6 subtriangle in Q,.

  Fig. 4. & located at a corner of a quadrilateral clcamt.

  Fig. 6. Positions of to io a six-node triangle. Deihition of 2 and 8.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Change from (u,, u2) to (y,, yz). Case when the element is a quadrilateral and ya lies at a comer of element e.

Table 1 .

 1 Circular crack. Errors on ds and 4, with respect to the

	,	analytical solution
	NUUlbCt of d.o.f.	Mipimal errors -4,	Maximal errors on 4,	ElT0l-S on K,
	85	0.28%	5.12%	6.06%
	161	-0.05%	3.05%	2.29%
	209	0.42%	-4.41%	0.78%
	321	-0.02%	-4.12%	0.11%

Table 2 .

 2 Elliptical crack. Errors on 4, and KI

	NUXllber	Minimal errors	Maximal errors	Rrrora on	&rot3 on
	of d.o.f.	on 93	0114,	K,(A)	RI (B)
	161 113	-0.10% 0.16%	-5.36% 8.22%	4.21% 8.71%	-2.25% 3.13%

The mesh generation technique is rudimentary: the crack is divided up with elliptical curves shaped in the crack contour, and radial straight lines issued from the crack centre. This inevitably gives rise to some distorted (elongated) elements which behave badly in numerical purposes. Moreover, this does not allow an accurate representation of the crack edge, whenever the semi-axes ratio becomes large. Figure 9 clearly shows that the mesh is particularly ill-suited to shape the vertex of the major axis. These reasons explain why computations are limited to small semi-axes ratio geometries. They also justify the use of 12 (resp. 16) Gaussian points on triangles (resp. quadrilaterals). Numerical results will be compared with the theoretical solutions [4]:

where E(k) = j&i -k2 second kind. & together with errors on the stress intensity factor K, at the vertices A, B of the semi-axes (Fig. 9).

11

APPENDIX

For the six-node trian8k and the eight-node quadrilateral, we give functions (N) together with the modified ftm&ms (#>, (M@> and [I). Function [_$I can be obtained from relations (4).

Six-node triangular elemenfs

Table 1A. By noting that: 8,' -5;, = -X(cos 6 + sin 8), we obtain Table 3A. Here the Iw(@ are independent of 1.

The modi@ed Jacobian inverse /

Application of relations (4) &es, after performing all calculations:

UG, 1, @)I =

Eight-t&e quadrilaterals

(1 -&I(1 -r:, Table 5A. Function A