A Levan

J Royer

Part-circular surface cracks in round bars under tension, bending and twisting

Circular-fronted cracks in round bars subject to tension, bending and twisting are considered. Numerical expressions are given allowing the calculation of stress intensity factors K., K 1 ., Km at every point on the crack front for a wide range of crack geometries. Comparisons are made with analytical, experimental and numerical results available in the literature. Crack shapes satisfying the iso-K 1 criterion are also determined, making it possible to investigate the problem of crack growth behaviour under tensile or bending fatigue loads.

Introduction

As cylindrical specimens are easily machined, components with a round cross section are commonly used in engineering structures. Bars, shafts, wires, reinforcements, bolts, screws or pins are examples of cylindrically-shaped structural elements. In many applications the loading conditions are quite complex. Under cyclic or repeated loads, fatigue cracks can occur in such components. Experimental works [1][2][START_REF] Nezu | The 25th Japan Conyress 011 Material Research[END_REF][START_REF] Wilhem | Adt•ann:s in Fracture Research • Fracture X 1. Proceedin~Js of the 5th International Conference on Fracture Research (ICF-5)[END_REF][START_REF] Forman | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF] showed that surface cracks created by fatigue have approximately circular fronts. In order to predict the growth of such a crack and the strength of the cracked component, an accurate stress analysis is required. Under linear elastic conditions, this leads us to precisely calculate the stress intensity factors which govern the mechanical state in the structure. In the general case of a combined load, a mixed mode I + II + III situation exists along the crack front. As in any arbitrary three-dimensional configuration, the variation of stress intensity factors Kt. K 11 , Km on the crack front must be taken into account.

Several works have been devoted to the problem of a surface crack in a round bar. Table 1 summarizes the main features of some studies relating to this problem. In Table 1, K 1 A denotes the stress intensity factor at the deepest point of the crack, K 18 the stress intensity factor at the surface terminal point (intersection between the lateral surface and the crack front), K. 1 the average stress intensity factor calculated on the crack front. Additionally, straight-fronted cracks are considered as particular cases of circular or semi-elliptical cracks. The bibliography, although representative, makes no claim to completeness -rather, it is a selection of the many references found in this field.

In this paper, use is made of a boundary integral equation specifically developed for fracture problems [22] to determine the stress intensity factors K 1 , K 11 , K 111 for circular-fronted cracks in a round bar subjected to tension, bending and torsion. Polynomial expressions are provided allowing the calculation of these stress intensity factors at every point on the crack front for a wide range of crack geometries. Crack shapes satisfying the iso-K 1 criterion are computed, and the problem of crack growth behaviour under tensile or bending fatigue loads is discussed. The results are compared with other experimental and numerical results available in the literature Bush (t976) [START_REF] Bush | Proceedings o( the Society .fiJr Experimental Stress Analysis[END_REF] Straight-fronted crack Three-point Compliance K., bending measurements [START_REF] Daoud | [END_REF] [START_REF] Daoud | [END_REF] Straight-fronted crack Tension Compliance K.l measurements Astiz et al. (1981) [START_REF] Astiz | Proceedinys o(the Societ_vfhr Experimental Su•e.~.~• Analysis, Spriny Meeting[END_REF] Straight-fronted crack Tension Photoelasticity KIA Salah el din et al.

Actual shape :::-::circular Tension and Fatigue tests KIA (1981) [2] frontal crack bending Bush (1981) [START_REF] Bush | [END_REF] Straight-fronted crack Tension Compliance K., Kr [START_REF] Wilhem | Adt•ann:s in Fracture Research • Fracture X 1. Proceedin~Js of the 5th International Conference on Fracture Research (ICF-5)[END_REF] fronted crack [START_REF] Nezu | The 25th Japan Conyress 011 Material Research[END_REF] [START_REF] Nezu | The 25th Japan Conyress 011 Material Research[END_REF] Actual shape :::-:: circular Tension Fatigue tests A power Ia w between fronted crack crack propagation rate da fdN and range of stress intensity factor AK Forman et al. ( 1986) [START_REF] Forman | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF] Actual shape :::-:: circular for mode I crack problems. Comparisons are also made with analytical results for some limit configurations.

Formulation

Geometry

Figure 1 represents a round bar of radius R and height h ~ R containing in its median cross section a surface crack. The crack front is part of a circle of radius R'. When R' equals the crack depth a, the crack is a so-called semi-circular crack. On the other hand, when R' tends to infinity, the crack is referred to as 'straight-fronted'. Any intermediate crack geometry between the two above limiting cases can be defined by the crack shape parameter 11.-B 0 B/B 0 B 1 , ex E [0, 1] (Fig. 2). The computation is carried out on 4 crack shapes: 11. = 0 (semi-circular), !, 1 and 1 (straight-fronted). The two intermediate cracks have been determined in such a way that they divide arc B 0 N 1 into three equal sub-arcs. Six relative crack depths are considered: a/ R = 0.04, 0.12, 0.24, 0.40, 0.60, 0.85. Since there are 4 crack shapes for each crack depth, 24 crack geometries are analyzed (Fig. 3). Fig. 3. Computed crack geometries.

Loads

The bar is successively subjected at its ends to three loads (Fig. 4):

• a uniform tensile stress a.

• a linear tensile stress with outer fiber maximal value am .

• a linearly distributed shear with outer fiber maximal value rm, resulting from torques twisting the bar.

As far as stress intensity factor solutions are concerned, these loads are respectively equivalent to a uniform pressure, a linear pressure and an axisymmetricallinear shear applied on the crack faces.

Equations of the problem

The bar is assumed to be made up of a homogeneous, isotropic, linear elastic material characterized by Young's modulus E and Poisson's ratio v. For solving the problem use ts made of the integral equation given in [22]. Let us introduce the following notations:

• x:
field point with co-ordinates (x 1 , x 2 , x 3 ) in a rectilinear rectangular system of axes Oe1 e2 e 3 (Fig. 1)

• x 0 :
source point, r = II x -x 0 II, er = (x -x 0 )/r (Fig. 5 • t(x 0 , nx 0 ): the stress vector at x 0 , related to normal Dxo• • s . cr• the surface of the crack, which is represented by a Cartesian parametrization the cylindrical lateral surface of the bar, which admits the parameterization

Flat: Lllat = [0, 2rr[ X [ -h, h] 3 (8, XJ)-X E slat•
In fact, the lateral surface S 1 at also includes the upper and lower flat bases of the bar, and another surface must be added at the junction of the bases and the cylindrical surface (Fig. 6) in order to avoid any discontinuity of the normal vector on S 1 at• The outer surface of the bar is then sufficiently smooth as required by the theory.

• <pcr:Scr3X-(/Jcr(X)EIR\ <plat:S 1 a 1 3X~<p 1 at(x)EIR 3 : the unknown densities defined respec- tively on Scr and S 1 a 1 , which are shown to be equal to the displacement jumps through the surfaces Scr and S 1 a 1 respectively [23], [START_REF] Kupradze | Three-dimensional Problems of the Mathematical Theory of Elasticity and Thermoelasticity[END_REF].

• <l>cr = (/Jcr o Fen <l>lat = <plat o Flat (compounds of functions F and ({J).

The set of equations of the problem can then be written in the following form, which expresses the stress vector at any point x 0 and related to normal nxo, in terms of the densities <per and (J)Jat:

( I ) Fiq. 6. Typical mesh modelling the cracked bar.

where Ker cr and Ker 1 a 1 are the kernels defined by In relation (2), c:iik is the Levi-Civita symbol. The function F is defined by F: .L\ 3 (u1, u 2 )-+ x, a, /JE{1,2}

For Kercr : (u 1 , u2) = (xt. Xz), F = F er• <I>= <l>cr> Dxo = e3 For Ker 1 at: (u 1 , u 2 ) = (8, x 3 ), F = F 1 at• <I>= <1> 1 at• nxo is the outward normal at the point x 0 to the lateral surface.

The symbol pv before J indicates that the integral is understood in the sense of the principal value.

Partial derivatives of any function f with respect to variables ui is denoted f .i• Implicit summation is made over any repeated index.

The stress vector at the point x 0 is written in (1) as the sum of the surface intergals over the crack and the lateral surface of the bar. Eventually one has to solve the coupled set of equations (1) with unknowns q>cr and ({>tat•

Numerical results

For all numerical purposes, the Poisson's ratio vis taken as equal to 0.3. Both the crack and the lateral surface are discretized into finite elements and Eqns. (1) are solved by the collection method. Figure 6 shows a typical mesh modelling the cracked bar. Eight-node or six-node isoparametric 20 elements are used throughout the structure. Quarter-point elements [START_REF] Barsoum | [END_REF] are specifically used along the crack front. Depending on the surface to which the element belongs, the geometry transformation and the interpolation of the density are performed with different variables [26]:

-for elements belonging to the crack surface, the mapped variables through geometric transformation are Cartesian co-ordinates (x 1 , x 2 )E ~cr and the interpolated functions are Cartesian components of ((>cr. q>~> q> 2 and q> 3 -for elements belonging to the lateral surface of the bar, the mapped variables are cylindrical coordinates (8, x 3 ) E ~Jat and the interpolated functions are here again Cartesian components of q> 131 for simplicity. Cylindrical variables (8, x 3 ) allow us to shape the finite elements into curved elements fitting the cylindrical surface. However, this cannot be seen in Fig. 6 since the sides of the elements are represented by segments.

The mesh covers the crack surface together with the lateral free surface. Of course, if the different loads were treated separately, then symmetries or skew-symmetries could be exploited in order to reduce the problem to the study of one quarter or one half of the bar, provided adequate boundary conditions are added. This approach is not chosen here for two reasons: first, the whole structure is preserved so that several loads can be applied simultaneously, thus requiring only one solution of the algebraic system. Secondly, the matrix of the system being fully populated and moreover non-symmetrical (as is the case with any boundary integral equation method), obtaining symmetrical or skew-symmetrical final results ensures the accuracy of input data.

Eventually, the solution of the discretized equation of (1) provides the densities cpcr and cp 1 a 1 , which in turn yield the complete elastic solution for the problem of the cracked bar, in particular the stress intensity factors can be evaluated along the crack front. Figure 7 shows the deformation of the lateral surface under a tensile load and of the crack under torsion. The latter is strongly emphasized in order to make it visible, explaining why some points on the crack seem to overlap the crack surface. The normalized stress intensity factors at any point located on the crack front are determined as follows, using the notations of Fig. 8:

(a) Case of a uniform pressure on the crack (tensile load): 1 E

~-K_,_ -I E 2 lim {(P3 {br} = 2rr,/Jjj lim { E ,2 cp/ :_}' afo afo8(1 -v )p~o
1 . { )¥n} _ 2nfo 1 . { E .!1!.2_} 2 Im <p 3 - Im 2 ' (Jmfo 8(1 -V ) p .... O p (Jmfo p .... o 16n(l -V ) JP (3.2)
(c) case of a linear shear on the crack (torsion):

1 E 1 . { )¥n} _ 2nfo 1 . { E ~} C -----, 2 =--tm <fJv - C tm 2 r:, , Tmy na 8(1-v )p .... o p Tmy nap .... O 16n(l-v ) V p (3.3) _K-=m= _ 1 E lim {<p {hr_2n} = 2nfo (t _ v) lim { E _<fJ_t} rmfo rmfo 8(1 + v)p .... O t~P rmfo p-o 16n(1-v 2 ) JP ' (3.4)
where <fJv, <p 0 <p 3 are components of <fJcr in the local basis (ev, e" e 3 ) generally varying along the crack front (Fig. 8).

It should be mentioned that relations (3) are taken just as definitions for series intensity factors without any assumption on the elastic state at the point of interest. Neither plane stress nor plane strain states, which are certainly predominant at the neighbourhood of the free surface or the deepest point of the crack, are assumed. The transition from one state to another will not be discussed in this paper.

As shown by (2), solving (1) gives directly E/(16n(l -v 2 ))<fJcrthe term bracketed together in (3)explaining why it is unnecessary to specify the value of Young's modulus E when computing stress intensity factors. Figure 9 shows the stress intensity factors obtained for the geometry corresponding to aj R = 0.4 and a = l, versus the relative abcissa sfsm (Folio wing the notations in Fig. 1, s is the curvilinear abcissa of a point on the crack front, sm is arc length AB,

sfsm = -1 at B', = + 1 at B).
Thus, for each geometry, 4 discrete curves are obtained corresponding respectively to K 1 in tension, K, in bending, K 11 and K 111 in twisting. As predicted, the K 1 value in the bending case is always smaller than that in tensile loading at homologous points on the crack front. As for K 111 values, they are found to be negative for all geometries. The negative sign of K 111 is merely due to the fact that <p 3 equals the normal displacement at the upper crack face minus that at the lower face, and to the choice of the local basis orientation as shown in Fig. 8. Now using the least square method, the set of discrete values of stress intensity factors obtained for the 24 computed geometries are fitted in the following polynomial form:

K 3 3 (a)i ( s )k (a)
r m rra i = 0 j = 0 k = 1' 3 R Sm I I I c~W) ~ ry_i .!... , 3 3 (
)i ( )k i=Oj=O k=0 , 2.4. 6 R Sm
where the Ciik coefficients are given in Table 2. All mean quadratic errors resulting from the fitting expressions (4) are about 1 percent, which is quite acceptable when compared to the accuracy of the finite element method. Expressions (4) should be used with crack depths of less than one-half diameter (0 ~ aj R ~ 0.9) and with lf.E [0, 1]. Figure 10

Stress intensity factors versus crack depth

Here, consideration will be limited to the deepest point where the stress intensity factors are given from (4) by putting s/sm = 0 KA ""

(a)i .

C = L, L, Cijo R ex/,
(Jv nu i j [START_REF] Forman | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF] where K stands forK,, K 11 or K 111 and CJ also stands for CJm or Tm • Figure 11 shows normalized K" K 111 at the deepest point of the crack resulting from the basic loads -tension, bending and torsion-versus the relative crack depth a/ R. Factor K 11 is not plotted since it is identically zero at s = 0. For each load there are four K, or K 111 curves which relate individually to one crack shape, ::x = 0, • !, 1 or 1. It can be seen that the stress intensity factors vary continuously with the crack shape but they are not necessarily monotone functions of the crack depth. For instance in decreases for small crack depths when normalized by a~. It appears from Figs. 11a andb that when a/ R takes the limit value zero, K 1 at the deepest point takes almost the same value for the crack shapes a = 1 and !. This would mean that, despite the notable change in radius of curvature -R' falling off from infinity when r:t. = 1 (straight-fronted crack) to a finite value when a = tthe stress intensity factors remain almost unchanged. In this respect, one can readily establish the limit value for radius R' when aj R tends to zero, the radius R of the bar and crack shape IJ' . being kept constant

(}.2 lim R' = 2 R. a /R--+0 J -'Y. ( 6
)
Relation [START_REF] Bush | Proceedings o( the Society .fiJr Experimental Stress Analysis[END_REF] gives both trivial results -for a = 1, R' tends to infinity and for r:t. = 0, R' tends to zero-and rather non-obvious results: for r:t. = iR' tends to R/ 8 = 0.125R and for 'Y. = fR' tends to 4R/ 5 = 0.8R, as shown in Fig. 12.

Stress intensity factors along the crack front -Application to the crack growth behaviour problem

Figure 13 plots the normalized stress intensity factors versus relative abcissa s/sm along the crack front, for the relative crack depth aj R = 0.4. Figure 13a shows more visibly than Fig. 10 that in the case of a tensile loading, the curvature of the K 1 curve changes in sign when passing from the semi-circular crack (' Y. = 0) to the straight-fronted one (a = 1). Meanwhile, in this instance when 'Y. is approximately -1, a crack shape can be observed such that K 1 remains virtually constant along the crack front. This means that if the iso-K 1 criterion is chosen to predict the propagation of cracks created under mode I fatigue, the actual shape is that corresponding to 'Y. = -t.

Likewise, Fig. 13b shows that under a bending load the iso-K 1 propagation criterion in this instance leads to the shape corresponding to ' Y. = j.

More precisely, for a given crack depth, the crack shape satisfying the iso-K 1 criterion can be computed in the following manner. Bearing in mind that R, a and a are constant, deriving a=l.

---------;::_~_=:J~,_....f-.r-,JC#C::..::=-_•---------R I --) 00 ,__...,.,. ~ expression for K 1 given in (4) with respect to the relative abcissa s/sm yields a K

( a )i ( s )k-1 ( 1 ) = afo~ ~ I kCiik R a/ - a .!__ I 1 k :;;. 1 Sm Sm ( 7
)
As it is clear that the derivative of K, cannot be identically equated to zero for all values of s/sm E [ -1 , 1], the solution crack shape rx is such that it minimizes max II aK tfa(sfsm)ll, i.e. rx is the solution of the min-max condition ( a)i •( s )k-t min max

I I I kCiik -rt. 1 - ae[O, l] s/ sme[ -1 . 1] i j k ;;. l R Sm ( 8
)
In order to outline the crack propagation during mode I fatigue, ( 8) is approximately solved and the results obtained are within 5 percent accuracy. Table 3 lists the solution crack shapes rx for the relative crack depths aj R considered in this paper. Figure 14 compares the crack shapes predicted by the iso-K, criterion with those of the cracks intersecting the free lateral surface at 90° angles. The latter cracks will be referred to, for brevity, 2(R -a) [START_REF] Bush | [END_REF] It can be observed that the difference between iso-K, cracks and 90: angle ones is fairly small in the case of tension (Fig. 14a) whereas this difference becomes nota ble for medium-sized cracks under bending (Fig. 14b). This would mean that in tensile load fatigue tests, cracks verifying the iso-K, criterion are almost 90" angle cracks, while in bending fatigue tests this is true only for cracks with very small or very la rge depths. It should be remembered that this assessment has been made with v equal to 0.3, and that the study of the influence of Poisson's ratio on the shape of cracks created by fatigue, which would not be negligible. is beyond the scope of this work . Figures 15,[START_REF] Raju | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF] show the change in shape as the crack depth increases during the propagation. At the crack initiation, the crack is almost semi-circular. The further it grows up, i.e. the more the crack depth increases, the more it resembles a straight-fronted crack. Also represented in Figs. 15, 16 are the following stress intensity factors for a 90° angle crack: K 1 at the deepest point 1 (in this case it is the minimum K, on the crack front) and the average K 1 denoted by K 1 (10) It is clearly shown in Fig. 15 that for 90° angle cracks under tension, K 1 A is approximately K, (the curves e and fare quite the same) and again one can recognize that 90° angle cracks verify the iso-K 1 criterion. On the other hand, Fig. 16 shows that in the bending case, 90° angle cracks verify the iso-K 1 criterion only for ajR ~ 0 and ajR ~ 0.9.

Figures 13c and13d show the stress intensity factors K 11 and K 111 arising from a twisting moment applied at the ends of the bar. As already mentioned, K 11 varies almost as a linear function of the curvilinear abcissa s. On the other hand, for all crack shapes the absolute value of K m is maximal at the deepest point. In any case, the aspect of K 11 and K 111 curves are little influenced by the crack shape, contrary to what happens to K 1 • Regarding the crack propagation behaviour with the presence of anti-plane shear, experimental investigations in pure mode III by [27] clearly show that the crack no longer grows by extending in its own plane but by "l"TT"T"T"T"T"T"T"T"T'T"'rT"'I"'T'"'r-T"'r-T generating multiple penny-shaped cracks which straddle the original crack front. The study of the crack deviation in a bar under torsion must be even more difficult at the points on the front, where Kn and K 111 take comparable values. Similarly, under combined tensile and anti-plane shear loading (mode I + III), the crack grows by developing multiple lancelike fracture facets surrounding the crack front, as shown in [28]. These observations should be taken into account when the crack instability is studied under general loading conditions. Lastly, it is noteworthy that the provided stress intensity factors values are valid only over about 80 percent of the crack front length. The K values for sfsm approaching + 1 are affected by phenomena extraneous to the present work, such as the vicinity of the surface terminal points B or B' (Fig. 1) that modifies the crack tip singularity. Also, the poor refinement of the finite element mesh around these zones must lower the accuracy of the numerical results. In any event, the results obtained in [29] prove that the crack tip singularity at the surface point depends on the Poisson ratio v and the terminal point incident angle 1/J (see Fig. 1) between the crack front and the surface line BB'. For a given value of v, there exists a limiting value of 1/J for which the stress intensity factor K 1 tends to a non zero finite value. If 1/J is less than this limit value, K 1 falls off to zero and if 1/J is greater, K 1 becomes infinite. In both cases K, classically defined loses its physical meaning.

Comparison with theoretical results

No theoretical results are available for surface cracks in round bars. However, stress intensity factors at the deepest point (point A in F ig. 1) of very shallow straight-fronted cracks can be effectively compared with analytical solutions for the single-edge crack in the half-space. This is due to the likeness between the latter geometry and the crack configuration viewed from point A when aj R tends to zero.

Bar under tension

The theoretical value for an edge crack of length a under a uniform pressure (J is given by Koiter [30] K 1 = 1.1215(J~.

On the other hand, making a/ R tend to zero in [START_REF] Forman | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF] gives

I . KIA -'\' c<It> j Im c_-~ OjOCX, aj R-+0 (JV 1ta j which yields the limit value of normalized K 1 at deepest point A of straight-fronted cracks lim fo ~ L cgyo ~ uz6. aj R-+O,rz= 1 (J na j (11) (12) (13)
This value agrees well with the theoretical value (11) within an error of +0.4 percent. As regards the limit value of K 1 A for semi-circular cracks when aj R tends to zero, one expects it to be greater than the theoretical value 2/n = 0.6366 for the penny-shaped crack embedded in the infinite body, as is easily explained by the presence of the free lateral surface of the bar which must allow a wider crack opening, thus a greater K 1 A• Indeed, here the limit value for semi-circular cracks is found to be 0.668, which is about 5 percent higher than the foregoing analytical value 2/n.

a i a ( K ) ( ) i-1 a(a/R) (Jfo -i~1 ~ rC,i~ R a . ( 14
)
Hence the slopes at aj R = 0 of the normalized K 1 A curves are

I . a ( KIA ) c<H> o 028 -1 20 1m a( /R) C. = 1 oo = . = tan af R -+O ,rz = 0 a (J V na ( 15
) ( 16
)
respectively for semi-circular and straight-fronted cracks. These are quite small values that compare well with the zero slope drawn from relation (11), lim a(Kd(J fo)/a (a/R) = 0.

Bar under bending

From Table 2 one obtains the limit value of normalized K 1 at deepest point A of straight-fronted cracks in the case of bending load [START_REF] Nord | [END_REF] This agrees well with Koiter's theoretical value [30] which is also valid in the bending case, within an error of + l percent. Note that the slopes at a/R = 0 for semi-circular and straight-fronted cracks can also be computed in a similar way as in tensile loading lim

a ( K,~) = C\'gb = -0.471 = tan-1 ( -25 "), a! R -+ 0.2 = 0 o(aj R) amJna ( 18) lim 2 ( K' ' ;-) = "LC\'J6 = -0.961 = tan -1 (-44"), a/ R-+0,7 = l o(a/ R) ()myna J ( 19
)
although no easy comparison with theoretical values is possible to our knowledge.

Bar under torsion

Lastly, consider an edge crack of length a in the semi-infinite space, subjected to a uniform shear r parallel to the crack front. The a nalytical result for the problem is given in [START_REF] Tada | The Stress Analysis of Cracks Handbook[END_REF] r;--:

K, 11 = -Ty'(rra).

(

) 20
Applying the same reasoning to (5) as in proving [START_REF] Fan | [END_REF], one obtains lim

KIIIA = L cg~~ = -0.875. a/R -+ 0 .. 7= 1! f;;z 1 • mv ( 21
)
This value is 13 percent higher tha n the theoretical value -I. The difference can be accounted for by the fact that the shear is in-plane stress. In Fig. 17 difference between the problem of a shallow crack in a round bar and that of an edge crack in the semi-infinite medium. In mode I problems, the applied stresses are normal to this area and their effect can be neglected around the deepest point A of interest. On the contrary, in mode III problems the stresses are parallel to the crack front, hence they may have some influence on K 111 at point A though applied on an area geographically remote.

Comparison with experimental results

Available results in the literature are mainly concerned with mode I and points A, B on the crack fr~ .. ~ (Fig. 1). Sometimes, K 1 is even assumed to be constant on the crack front so that only the average value is considered. Therefore, the somparison with either experimental or numericai results can be made only in the tension and bending cases. Also, as discussed above, stress intensity factor values at the surface terminal point B will be discarded and only the values at the deepest point A-or the average values whenever they are given -are compared. We shall denote the diameter of the bar by D, D = 2R.

Bar under tension

• Straight-fronted cracks

By the compliance method Daoud et al. [START_REF] Daoud | [END_REF] give the average normalized K, for straight-fronted crack

Table 4 gives the normalized K 1 at the deepest point A by photoelasticity determined by Astiz et al. [START_REF] Astiz | Proceedinys o(the Societ_vfhr Experimental Su•e.~.~• Analysis, Spriny Meeting[END_REF]. In fact, the value 1.85 corresponding to aj D = 0.45 in this table will not be retained here as it is unreliable according to [START_REF] Astiz | Proceedinys o(the Societ_vfhr Experimental Su•e.~.~• Analysis, Spriny Meeting[END_REF]. Bush [START_REF] Bush | [END_REF] gives the compliance expression c in terms of the relative crack depth a/D c = 0.0598723 + 0.2680344(;r• + 0.2508381 (;r + 39.43071 (;f. The unity for compliance c in ( 23) is 10-6 injlbf, Young's modulus in ( 24) is 10.6 x 10 6 psi, and D = 3 in. Attention should be drawn to the fact that in [START_REF] Bush | [END_REF] the Poisson ratio is 0.32 and not 0.3 as assumed throughout this paper. However, the results are expected to be close enough to be comparable.

• Cracks intersecting the lateral surface at right angles (t/1 = 90°, Fig. I} A special fatigue marking technique to outline the crack propagating allowed Wilhem et a l. [START_REF] Wilhem | Adt•ann:s in Fracture Research • Fracture X 1. Proceedin~Js of the 5th International Conference on Fracture Research (ICF-5)[END_REF] to express the normalized K, for right angle cracks, constant along the crack front, as K, a (a) 2 (Jv na

D D D ( 25
)
Experiments conducted on fatigue cracks led Forman et al. [START_REF] Forman | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF], [START_REF] Forman | Fatigue crack growth computer program[END_REF] Figure 18 shows a good agreement between different results for small crack depths. Regarding large crack depths, other workers' results relating to 90° angle cracks are rather closer to present K 1 for semi-circular cracks. As for straight-fronted cracks, the discrepancy becomes notable for large crack depths too.

Bar under bending

• Straightjronted cracks By means of compliance measurements Bush [START_REF] Bush | Proceedings o( the Society .fiJr Experimental Stress Analysis[END_REF] provides the following equation for calculating the average stress intensity factor in the case of bending load (27) where the compliance c is determined by either expression Here again, E = I 0.6 x 10 6 psi. and v = 0.32. Relation ( 28) is used with diameter D = 3 in, relation (29) with D = 6 in. In both relations, length I is about 10 in.

• Cracks intersecting the lateral swface at right angles Experiments conducted on fatigue cracks led Forman et al. [START_REF] Forman | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF] to the following approximate expression for the normalized K 1 constant along the crack front K,

2D

with the same accuracy reported above for the bar under tension.

(

) 30
• The comparison with all the available results for the bar under bending are shown in Fig. 19. For 90c angle cracks, a good agreement is observed at small crack depths only. As the crack depth increases, the experimental curve is rather closer to ours computed for semi-circular cracks. This remark agrees with Caspers et al.'s [19] following which 'the theory of perpendicular angles between crack front and shaft circumference seems to be approximately correct for pure bending up to an aj R-ratio of 1.0, but does not seem to be so applicable for tension, especially for increasing a/ R-ratios'.

Concerning the straight-fronted cracks, the computed curve is within the range of Bush's experimental ones. Lastly, it should be mentioned that straight-fronted cracks do not exist naturally and they are actually very difficult to obtain in experimental works. In most cases, a sharp notch was machined to simulate the crack [START_REF] Bush | Proceedings o( the Society .fiJr Experimental Stress Analysis[END_REF][START_REF] Daoud | [END_REF][START_REF] Astiz | Proceedinys o(the Societ_vfhr Experimental Su•e.~.~• Analysis, Spriny Meeting[END_REF][START_REF] Bush | [END_REF]. As the machined notch must be wider for a larger notch depth, it is unlikely to be comparable with a real crack. This should explain why the results agree less for large crack depths.

Comparison with numerical results

Here again, the comparisons are made in tension and bending cases. Moreover, only the K 1 at the deepest point or its mean value is reported.

Bar under tension

Figure 20 shows all the numerical results for the case of a tensile load. Astiz [ 18] considered semi-elliptic cracks with axes a, h. The normalized K 1 factor at the deepest point A is given as a The configurations which are comparable with those considered in this paper correspond to aj b = 1 (semi-circular cracks) and ajb = 0 (straight-fronted cracks).

Investigating part-circular cracks, Caspers et al. [ 19] also gives normalized K 1 at the deepest point in the polynomial form

K lA = L I, CV) az ' ~ J 4 5 ( )"( )" (]Fa i = 0 j = 1 IJ R 1t , (32)
where az =a-R(l -cos 8), (J is the angle at the centre sustended by the arc BB' (Fig. 1). Raju et al. [START_REF] Raju | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF] considered nearly semi-elliptical cracks such that the crack front intersected the lateral surface of the bar at 90° angles. It should be noticed that the so-called crack length in this reference is defined as the arc length measured along the cylindrical surface and not as the major axis of the (nearly) elliptic crack. Thus the results of this reference for ajc = l must be compared to others keeping in mind the difference between the analyzed geometries. The normalized stress intensity factor at point A for ajc = I are listed in Table 5.

• Straight-fronted cracks Daoud et al. [START_REF] Daoud | [END_REF] computed K 1 assumed to be constant along the crack front

K 1 a (a) 2

Bar under bending

Figure 21 shows all the numerical results for the case of a bending load. The normalized K 1 at point A for ajc = I from Raju et al.'s work [START_REF] Raju | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF] is given in Table 5.

Considering a lateral bending load, Caspers et al. [19] express K 1 at the deepest point in the same polynomial form as for tension load (34) 1. [START_REF] Forman | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF] ...-... It should be mentioned that in [19] the linearly distributed stress is zero at the deepest point level (x 1 = -R +a, Fig. 1) and not at the diameter level (x 1 = 0). An adequate combination of ( 32) and ( 34) must be carried out to obtain the results for the bending case.

1.25 G G KaA ..__ ~ 1.0 "'0 Q.) N 0.75 •--~ E --o--o--0 8 ":••:&--~-&--0~ a- 0 • z 0.5 • • KIA • 0.25 Q a/R
• Straight-fronted cracks

• Figures 20 and21 show a good agreement between different numerical results concerning semi-circular cracks, both in tension and bending cases. As regards straight-fronted cracks, the results agree well only for small crack depths. It should be mentioned that KrA values for straight-fronted cracks given in [19] are rather low and thus are closer to K 1 •

• Limiting configuration of very small semi-circular cracks

Let us now consider the limiting case when the radius of the semi-circular crack R' = a tends to zero. The limit value of K 1 at the deepest point A is found to be 0.668 and 0.670, respectively, for the tensile and bending load. On the other hand, since this value must be the same as for a semi-circular crack in the half-space, it can be compared to the limit values obtained by earlier investigators treating the half-space problem. Table 6 shows that the differences are really small, except for Nisitani et al.'s value which amounts to 5 percent above the present result.

Conclusion

The obtained results have brought additional information to the problem of the elastic cracked bar, especially in mode I. A methodical use of integral equations has proved efficient in that it allows the following possibilities:

-the computation of stress intensity factors at any point of the crack front (except for the surface terminating points), -the simultaneous study of modes I, II and III, -the polynomial fitting of numerical values providing the stress intensity factors as functions of three geometrical parameters: the crack depth a, the crack radius R' and the abscissa s measured on the crack front, the results applicable to three basic loads: tension, bending and torsion.

Limiting values of stress intensity factors when the crack depth tends to zero have been found to be in full agreement with the theoretical values for the half-space problem. Furthermore, the numerical values have been compared with both experimental and numerical results available in the literature, which are mainly concerned with the opening mode. A good agreement between the different results has been observed for small crack depths. Nevertheless, for large crack depths for factor K, increases more rapidly than others. The discrepancy is the more important for straight-fronted cracks.

Since no comparison is possible in modes II and III, it is assumed that the computation which is validated in mode I, also gives acceptable results for other modes. In fact, tests conducted on problems having known analytical solutions have shown that the numerical accuracy is higher in mode I problems than in others.

The polynomial fitting of the results has allowed a rapid and accurate calculation of the stress intensity factors as a function of the crack configuration and the applied loads. Also, the knowledge of these factors all along the crack front has made it possible to investigate either elastic fracture criteria or the fatigue crack growth behaviour. As an application of the results, the crack shapes verifying the iso-K 1 criterion have been computed in mode I fatigue problems.

Of course, the obtained results are not influenced by the Young modulus value but they do depend on Poisson's ratio v as shown by the governing equations (I) and (2). The computation has been carried out with v equal to 0.3.

The effect of Poisson's ratio upon the crack shape can be investigated as an extension of the present work. Besides, the proposed integral equation formulation can be applied without major modifications to the problem of thick-walled cylinders containing through-or part-through cracks.

  and K.,

  and K18

Fan

Fig. 1 .

 1 Fig. 1. Geometry.

Fig. 2 .

 2 Fig. 2. Crack shapes defined by parameter :x. Fig. 3. Computed crack geometries.

Fig. 4 .

 4 Fig. 4. Loads. (a) Tension. (b) Bending. (c) Torison.

Fig. 5 .

 5 Fig. 5. Variables of the integral equation set.

Fiy. 7 .

 7 VP a,_jrra p-->O 16rr(l-\ )...;p ( 3.1) a b (al Deformation of the lateral surface under tensile loading. (b) Deformation of the crack under torsion.

Fig. 8 .

 8 Fig. 8. Determination of stress intensity factors. p = Distance to the crack front. e,: in-plane normal to the crack front, outward with respect to the crack. er: in-plane tangent, er = e3 1\ e,..

Fig. 9 .

 9 Fig. 9. Normalized stress intensity factors along the crack front. s = curvilinear abcissa, sm = arc length AB (see Fig. I).

  in the tension case: fo = I I I C~W -ryJ -, in the bending case: fo = I I I Cl~~lry) -, ; •( s )k (c) in the tension case: ):;;_ = I I I C~Wlf. 1 -,

Fiy. 10 .

 10 Stress intensity factors as a functi o n of the c rack shape Y. and the rela ti ve abcissa s;'s.,,. (a ) Normalized K1 in tension. (b) Norma lized K 11 in torsion.

Fig. 11 .

 11 Fig. 11. Normalized stress intensity factors versus the relative crack depth. (a) Tension. (b) Bending. (c) Torsion.

Fiy. 12 .

 12 Limit values of crack radius R' fo r very shallow surface cracks.

IFig. 13 .

 13 Fig. 13. Normalized stress intensity factors versus the relative abcissa on the crack front. (a) K 1 in tension. (b) K 1 in bending. (c) K 11 in torsion. (d) Km in torsion.

Fiq. 14 .

 14 Predict ion of crack growth under mode I fatigue. as 90c (intersecting) angle cracks. For such cracks, 1/J = 90(; (Fig. I) and one can readily prove that the crack radius R' is related to the radius R of the bar and the crack depth a by the relation , a(2R -a) R = -----.

Fig. 15 .

 15 Fig. 15. Normalized K 1 in tension at the deepest point and average normalized K 1 for different crack shapes. Comparison with K1 computed from the iso-K 1 criterion. (a) K 1 at the deepest point for semi-circular cracks. (b) Average K1 for semi-circular cracks. (c) K 1 at the deepest point for straight-fronted cracks. (d) Average K 1 for straight-fronted cracks. (e) K1 at the deepest point for right angle cracks. (f) Average K 1 for right angle cracks. (g) K 1 for cracks verifying the iso-K 1 criterion.

Fig. 16 .

 16 Fig.[START_REF] Raju | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF]. Analogous to Fig.15, the loading is now a bending.

Further

, the hatched area represents the a b Fiq. I 7 .

 b7 Difference between the problem of a crack in a round bar and that of a n edge crack in the semi-infinite body.

Figure 18

 18 Figure 18 compares K, value of the present work with other authors' experimental values which are summarized below.

( 23 )Fig. 18 .

 2318 Fig.18. Comparison with experimental results for K 1 in tension: -8--: Daoud ct al.[START_REF] Daoud | [END_REF]. e : Astiz et al.[START_REF] Astiz | Proceedinys o(the Societ_vfhr Experimental Su•e.~.~• Analysis, Spriny Meeting[END_REF]. -+ -+ Bush[START_REF] Bush | [END_REF]. --G--G--: Wilhem et al.[START_REF] Wilhem | Adt•ann:s in Fracture Research • Fracture X 1. Proceedin~Js of the 5th International Conference on Fracture Research (ICF-5)[END_REF]. ~-: Forman et al.[START_REF] Forman | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF],[START_REF] Forman | Fatigue crack growth computer program[END_REF]. --: Present results (a) Average K 1 for semi-circular cracks. (b) Average K 1 for right angle cracks. (c) Average K 1 for straight-fronted cracks. (d) K 1 at the deepest point for straight-fronted cracks.

(a) 3 C

 3 .= 0.690-0.197 -+ 2.394 -+ 1.965 -, 0.15 < aj D < 0.45.

Fig. 19 .

 19 Fig.19. Comparison with experimental results for K 1 in bending: --8-: Bush[START_REF] Bush | Proceedings o( the Society .fiJr Experimental Stress Analysis[END_REF]. -+ -+ -: Forman et al.[START_REF] Forman | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF].:Present results (a) Average K 1 for semi-circular cracks. (b) Average K 1 for right angle cracks. (c) Average K 1 for straight-fronted cracks.

c = 4 . 2 '

 42 338749 + 23.66921 (~) + 131.0767 (~ r oo-• injlbf).

( 28 )

 28 c = 0.2910587 + 2.54535 (~ )'' oo-• in/lbf).

2 -

 2 -,....--= 0.92 -923 + 0.199(1-sin;~)] cos -•

Fig. 20 .

 20 Fig. 20. Comparison with other numerical results for K 1 in tension: ---8-••• • : Daoud et al. [7]. e : Raju et al. [16].-+ -+ : Astiz[18]. -B-: Caspers eta!.[19]. --:Present results (a) K 1 at the deepest point for semi-circular cracks. (b) Average K 1 for straight-fronted cracks. (c) K 1 at the deepest point for straight-fronted cracks.

(a) 3 (a) 4 -

 34 ----=~= = 1.11-3.59-+ 24.87 --53.39 -+ 57.23 -, 0.06 ~ ajD ~ 0.7.

Fig. 21 .

 21 Fig. 21. Comparison with other numerical results for K 1 in bending: 0

Table 1 .

 1 Summary review

	Reference
	Crack geometry
	Load	Experimental works
	Method
	Main results

  depicts the variation of normalized stress intensity factors K 1 in tension and Kn in torsion as a function of the crack shape If.

and the relative abcissa s/sm when the relative crack depth is 0.4. It is clearly shown that for nearly straight-fronted cracks (If. ~ 1) K 1 is maximum at the deepest point of the crack (s/sm = 0), and for semi-circular cracks (If. ~ 0) the maximum value for K 1 is reached in the neighbourhood of the free surface (ls/sml ~ 1). As for K 11 , it varies almost linearly with either If. or sf sm.

Table 2 .

 2 Ciik coefficients

		for tension	for bending	for torsion			for torsion
	0 0	0.66837E + 00	0.67003£+00	-0.48863 E + 00	0 0		-0.12653E+01
	0 0	-0.12819£+00	-0.11851£+00	0.53272£ + 00	0 0 3	0.54361E+00
	0 0	0.65362£ + 00	0.62139£+00	-0.66724£ -01	0 1 1	0.28415£ +01
	0 0	-0.63476£ +00	-0.60142E+00	-0.12661E +00	0 1 3	-0.89110£+00
	0	0.14917£+01	0.14660£+01	-0.12523£ +01	0 2		-0.29326£+01
	0	-0.15181£+01	-0.14844£ +01	-0.60349£ +00	0 2 3	0.10159£+01
	0	0.17418£+01	0.17943£+01	-0.49469£ +01	0 3 1	0.11584£+01
	0 1	-0.36700£ + 01	-0.36828£ +01	0.70620£+01	0 3 3	-0.66930£ +00
	0 2	-0.17108£ + 01	-0.16117£+01	0.15275£ + 01	1 0 1	0.81442£+00
	0 2	0.34585£ + 01	0.34011 E+01	-0.16732£ +00	0 3	-0.14275£+01
	0 2	-0.11730£ +02	-0.11906£+02	0.17182£ + 02	1	1	-0.13862£+02
	0 2	0.14317£+02	0.14447£+02	-0.22195£+02	1	3	0.55863£+01
	0 3	0.67685£+00	0.61063£+00	-0.66148£ +00	1 2 l	0.12860£ + 02
	0 3	-0.22663£ + 01	-0.22579E + 01	0.53588£ + 00	1 2 3	-0.21359£+01
	0 3	0.85036£ + 01	0.86794E + 01	-0.11410£+02	l 3 I	-0.25711£+01
	0 3	-0.93203£+01	-0.94854£+01	0.14460£+02	1 3 3	-0.10148£ +00
	1 0	0.27839£-01	-0.47133£ +00	0.36699£ + 00	2 0		-0.15475£+01
	0	0.17235£+01	0.20042£ + 01	-0.43199£+00	2 0 3	0.22907£+01
	0	-0.62703£+01	-0.57214£+01	0.18574£ + 00	2 1 t	0.23260£ + 02
	0	0.64590£ +01	0.57645£ + 01	-0.19042£+00	2 1 3	-0.87984£+01
	1	-0.81658£ +01	-0.80954£+01	0.55393£ +01	2 2 I	-0.15307£+02
	1	0.17088E + 02	0.16657E + 02	-0.67371£ +01	2 2 3	-0.76504£+01
	1	-0.47996£ +02	-0.46863E+02	0.48886£ + 02	2 3 1	-0.23390£+01
	1	0.57275£ + 02	0.55311£ +02	-0.56251 E +02	2 3 3	0.10195£ + 02
	2	0.18184£+02	0.16158£+02	-0.10060£ + 02	3 0 1	0.98730E + 00
	2	-0.51812£ +02	-0.50220£+02	0.14305£+02	3 0 3	-0.16671£ + 01
	2	0.18923£+03	0.18476£ +03	-0.13617E +03	3 1 1	-0.13498£ + 02
	2	-0.20458£+03	-0.19806£+03	0.16139£+03	3 1 3	0.35067£ + 01
	3	-0.10094£+02	-0.85514£ +01	0.49413£+01	3 2 1	0.53723£ + 01
	3	0.34890£ + 02	0.33390£ + 02	-0.77503£+01	3 2 3	0.12547£ + 02
	3	-0.13413£ +03	-0.13068£+03	0.82669£ + 02	3 3		0.44763£+01
	3	0.13902£ + 03	0.13465£ + 03	-0.99616£ +02	3 3 3	-0.12247£ +02
	2 0	0.37008£ + 00	0.43077£ + 00	-0.16124£ +00
	2 0	-0.48335£ +01	-0.43 163£+01	0. 73697E-01
	2 0	0.18126£+02	0.15820£ + 02	0.44922£ + 00
	2 0	-0.18218£+02	-0.15788£+02	-0.11729£ +00
	2	0.15716£+02	0.15715£+02	-0.98606£ +01
	2 1	-0.33545E + 02	-0.35751£ + 02	0.20672£ + 02
	2 1	0.94702£+02	0.10330£+03	-0.11358£+03
	2	-0.11449£+03	-0.11960£ +03	0.11944£+03
	2 2	-0.40103£+02	-0.35295E + 02	0.19091£ +02
	2 2	0.11148£+03	0.11409£+03	-0.37027£ +02
	2 2	-0.40124£+03	-0.41823£ +03	0.27773£+03
	2 2	0.43253£ + 03	0.44272£ + 03	-0.30932£ +03
	2 3	0.23819£+02	0.19974£+02	-0.98142£+01
	2 3	-0.77165£+02	-0.77778£ +02	0.16355£+02
	2 3	0.29188£+03	0.30048E + 03	-0.15546£ +03
	2 3	-0.30135£+03	-0.30666£ + 03	0.17910£ +03
	3 0	0.18566£ + 00	-0.80054£-01	0.55516£-01
	3 0	0.40254£ + 01	0.31925£ + 01	0.56762£-01
	3 0	-0.13204£+02	-0.11411£+02	-0.34231 E+OO
	3 0	0.13326£+02	O.ll431E+02	0.97845E-01
	3	-0.76417£+01	-0.84732£+01	0.55027£+01
	3 1	0.16116£+02	0.19953E +02	-0.14710£ +02
	3	-0.37988£ +02	-0.53843£ + 02	0.74788£+02

j k C!'.p IJ C!'.~' IJ C!}~' j k C!W Tahlt! 2. con/ d.

j k

Tahle 3 .

 3 Crack shapes verifying the iso-K 1 criterion. Com pa rison with right angle cracks

		a/R	0.04	0.12	0.24	0.40	0.60	0.1<5
	lso-K 1	Y.	0.0~	0.06	0.14	0.30	0.42	0.70
	criterion in
	TENSION	R',•R	0.048	0.147	0.331	0.694	1.1 72	2.940
	Iso-K,	:X I	0.04	0.10	0.29	0.4 7	(}59	0.70
	criterion in
	BENDING	R':R	0.05 1	0.161<	0.48 1	1.028	1.756	J.o97
	Right
	intersecting	'l.	0.003	0.02	().()7	0.17	0.36	0.73
	angle
	criterion	R'/ R	0.04!	0.1 28	0.278	0.533	1.050	J25X

Table 4 .

 4 Normalized K 1 at the deepest point of straight-fronted cracks[START_REF] Astiz | Proceedinys o(the Societ_vfhr Experimental Su•e.~.~• Analysis, Spriny Meeting[END_REF]

	a/D	0.45	0.21	0.30	0.39	0.41	0.46		0.06	0.12	0.31	0.42
		1.85	1.10	1.45	1.63	1.77	1.91	l st series	1.15	1.02	!.58	1.89
								2nd series	0.95	1.23	1.41	1.69

  Tah/e 5. Normalized stress intensity factor at the deepest point[START_REF] Raju | ASTM STP 905, 17th National Symposium on Fracture Mechanics[END_REF]

	ajD		0.050	0.125	0.200	0.275	0.350
	'2 K I c	tensile load	1.012	1.015	1.038	1.087	1.175
		bending load	0.938	0.836	0.749	0.683	0.629

7T.j X IA / <Jy 7T.a

Table 6 .

 6 Normalized K 1 at the deepest point for a semi-circular crack in the half-space under tensile load. In parentheses: value given in the bending case. Percent difference =(present result/ referenced result -I) x 100

	Reference	Smith et al.	Newman	Nisitani et al.	Newman et al.	lsida et al.	Present
		(1967) [33]	(1973) [34]	(1974) [35]	(1981 ) [36]	(1984) [37]	result
	KIA!(JJ(rra)	0.656	0.656	0.636	0.662	0.659	0.668
					(0.662)	(0.659)	(0.670)
	Difference	2	2	5		1.4
	(percent)

  Daoud et al. [14] computed the average K, over the crack front "m}a = 1.04-3.64 ~ + 16.86(~)' _ 32.59(~)' + 28.41 (~r 0.0625 ~ a/D ~ 0.625.


Maria Predari

Aurélien Esnard

Aurélien Esnard Graph

Graph partitioning techniques for load balancing of coupled simulations

come

Maria Predari & Aurélien Esnard (Inria, University of Bordeaux, France)

	G G r r a a p p h h p p a a r r t t i i t t i i o o n n i i n n g g t t e e c c h h n n i i q q u u e e s s f f o o r r l l o o a a d d
	b b a a l l a a n n c c i i n n g g o o f f c c o o u u p p l l e e d d s s i i m m u u l l a a t t i i o o n n s s

maria.predari@inria.fr & aurelien.esnard@ubordeaux.fr


Annette Calisti
email: annette.calisti@univ-amu.fr

Sandrine Ferri

Caroline Mossé

Bernard Talin

Ionization-potential depression and other dense plasma statistical property studies -Application to spectroscopic diagnostics

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

For various purposes, it is necessary to simulate virtual plasma composed of electrons and ions in different ionization states. Depending on the plasma conditions, e.g. when they are such that the electron-ion interactions are strong or those generated by direct target interaction with intense laser beams, it is required to account, in the simulations, for bound states as well as quantum mechanical processes to create and destroy them. This permits to avoid non-physical Coulomb collapse and to extend the applicability of classical MD simulations, which are ideally suited for fully-ionized, non-degenerate plasma, to a widest range of plasma conditions. Various methods have been developed and proposed by different groups. Among them, one can mention the methods based on the density functional theory (DFT) such as the quantum Molecular Dynamics method (QMD) [1] which is widely used to provide benchmark data but is limited to low temperature or the orbital free density functional theory molecular dynamics (OF-DFT-MD) [2] which has been developed to extend QMD to higher temperatures. Other methods exist which do not use DFT but simulate interacting classical ions and electrons incorporating some quantum mechanical effects through the interaction pseudo-potential between ions and electrons [3,4].

We have also chosen to simulate the ions and electrons as classical particles and to incorporate a minimum of quantum information through a regularized potential allowing to model collisional ionization and recombination processes. The BinGo-TCP code has been designed to deal with neutral mixtures composed of ions of the same atom with different charge states and electrons [5,6]. Within the limits of classical mechanics, all charge-charge interactions are accounted for in the particle motion. In this article, after a brief recall of the modeling basis, the importance of such a model to get a good understanding in correlation effects will be emphasized through a study of the ionization potential depression (IPD) and some other statistical properties in dense plasmas. In particular, some results obtained on aluminum plasmas at various conditions of densities and temperatures, corresponding to different conditions of coupling parameter, will be presented and analyzed.

Modeling

A classical multi-component plasma molecular dynamics (MD) code, the BinGo-TCP code, has been developed to follow the evolution of plasmas involving ions of various charge states and electrons. The overall structure of the BinGo-TCP code is consistent with that of a standard code of classical MD simulation. Non-relativistic point-like charged particles are simulated in a cubic box. The system inside the simulation box is neutral. The periodic boundary conditions are used together with the minimum image convention and the Newton's equations are solved using a Velocity-Verlet algorithm. A molecular dynamics simulation requires the definition of the potential function that governs interactions between the particles inside the simulation box.

Soft potentials

The electron-electron or ion-ion interactions are taken to be screened Coulomb:

V ii,ee = Z 2
i,e e 2 e -r/λ /r,

where the screening length, λ s/2, with s the box size, is chosen large enough not to affect the simulated results. Simulating ions and electrons together can be performed provided that an ionelectron potential finite at short distances avoid electron-ion collapse. Such effective regularized potentials are designed to approximately account for known quantum properties appropriate for the kind of investigated plasma. Here, the ion-electron regularized potential depends on Z i :

V ie (r) = -Z i e 2 e -r/λ (1 -e -r/δ(Z i ) )/r, (2)
with the regularization distance δ chosen to reproduce the ionization energy E i of the unperturbed ion of charge Z i in the ground state when the electron is located at the ion (r = 0).

δ(Z i ) = -Z i e 2 /E(Z i ) (3)
In other words, an electron located at an ion (r = 0) occupies the fundamental state of the ion whose charge is Z with a nucleus charge Z + 1.

Ionization/recombination protocol.

In the simulation, the particle motion is followed with very small time-steps ∼ 10 -20 s appropriate for the description of the micro motion of electrons around ions. These time steps required by the simulations are small compared to typical collision rates, thus any statistics on collisional events are expensive. The setting up of the population of electrons temporary trapped in the ion wells, i.e. the reaching of the required equilibrium state, depends on collisional events between electrons, and thus, is a very slow process. The choice of the previous soft ion-electron potential allows to implement a ionization/recombination protocol to control the plasma ion charge distribution and the trapping of electrons in the ion wells. Usually, the ionization-recombination mechanisms rely on an approximate analysis of collisional events between one ion and 1 or 2 electrons. In dense plasmas, the concept of collisions is not straightforward as the interactions involve all the particles within the screening length. To design a ionization/recombination protocol the definition of a collisional process is crucial but necessarily empirical. An analysis of the local environment of a given ion is performed to evaluate if the ion is susceptible to be ionized or recombined. The location of the two nearest neighbor electrons together with the sign of their total energy is used to evaluate if locally the plasma, at one step of its evolution, is favorable to an ionization (in case of positive energy) or a recombination (in case of negative energy) of the ion. When the local environment is favorable to an ionization, the ion of charge Z is changed to an ion of charge Z + 1 with an electron located on it. This pre-ionization state will give rise to an effective ionization if the interactions of the trapped electron with the surrounding plasma, permit it. A preparation phase of the particle set (into the simulation box) before extracting any sampling from simulations, is necessary. At the end of the preparation phase the system follows a quasi stable evolution with stationary trapped and free electron populations. The ionization/recombination process implemented in the code has two fundamental functions. It allows the evolution of the charge state population towards a stationary state depending on temperature, density and composition of the plasma and it favors the setting up of a population of electrons temporary trapped in the ion wells.

With this model, one gains the ability to describe the ion-electron coupling accounting for mixtures of ions undergoing changes of their ionization stages. It is important to note that the coupling of electrons with radiation is ignored and that the notion of discrete energy for the ionic excited states is here replaced by its continuous equivalent.

Some plasma properties.

Depending on its total energy and its nearest neighbor ion, an electron is either trapped or free. An illustration is given by the electronic total energy distribution functions measured in aluminum plasmas for various conditions. The analysis of these distribution functions together with the knowledge of the electron number in the simulation box, allows to infer the proportion of free electrons and thus the mean charge value. During a simulation run, various charge states of the same ion exist. The distribution function of these charge states is not representative of the plasma equilibrium as some electrons are trapped in the ion wells (revealed by the negative part of the energy distribution functions, cf. Figure 1). Therefore, to estimate the equilibrium charge states distribution function, a detailed examination of the electrons around each charge state has to be performed. The ionic charge state distribution function obtained for an aluminum plasma at T e =T i =20eV and ρ = 0.65g.cm -3 , is shown in Figure 2. It can be seen that the result which accounts for the trapped electrons, compares well with the one obtained with the Flychk code [7].

The ionization potential depression.

Due to its environment, an ion embedded in a dense plasma will undergo a potential significantly different from the one experienced by the same isolated ion. This plasma-dependent potential will lower the energy required to free a bound electron. This mechanism which is known as the ionization potential depression (IPD), will affect the number of bound states and consequently the radiation emitted by the system.

The evaluation of the IPD has important implications for dense plasma physics and in particular for the detailed prediction of the dense plasma equation of state and radiative opacity in stellar interiors, inertial confinement fusion research, or planetary interiors. Recent experiments [8,9] leading to IPD measurements in situ has renewed interest for this issue (see for example [10,[START_REF] Iglesias | High Energy Density Phys[END_REF][START_REF] Son | [END_REF]13,14,6]). Taking advantage of the particular design of the ionization protocol in the BinGo-TCP code, it is possible, when the ion is in a pre-ionization state, to measure the necessary energy to free an electron in the ground-state of a given ion and this, by accounting for all the interactions with the surrounding plasma. Figure 3 shows the distribution of this energy for four charge states existing in an aluminum plasma at T e =T i =50eV and ρ = 2.7g.cm -3 . The curves are not normalized and their relative intensities reflect the number of pre-ionization events for each ion of charge Z. It can be seen that due to a fluctuating local environment of the ions, the ionization energy is characterized by a distribution function. If one compares the mean energies deduced from these distribution functions with the corresponding energies of the isolated ions, it is possible to infer the ionization potential depression due to the interactions with the environment. The results are shown in Figure 4 with the results obtained with two theoretical IPD models, the Stewart-Pyatt (SP) [15] and the Ecker-Kröll (EK) [16] model, which apply across a wide range of densities. In those conditions of plasma, the electron Figure 4. The ionization potential depression obtained with the BinGo-TCP code (down triangles) for an aluminum plasma at T e =T i =50eV and ρ = 2.7g.cm -3 , compared with the IPD calculated with the SP (dot line) and EK (dot-dash line) models for the same conditions. The high density/low temperature limit of the SP model is also plotted (full line) for comparison. coupling plasma parameter is Γ e = 0.35 and the measured mean charge is Z b = 7.4 leading to a ion coupling plasma parameter, Γ i = 9.96. The effects of the charge correlations on the lowering of the necessary energy to free an electron are very strong and cannot be ignored. The comparisons with the theoretical models show that for these conditions, our results are in better agreement with the EK model than with the SP model. Other results obtained for other conditions of plasma (aluminum plasmas at T e =T i =20eV and ρ = 0.65g.cm -3 and T e =T i =550eV and ρ = 6.5g.cm -3 ), not presented here, show a good agreement with the SP model and thus, one cannot conclude on the relevance of one or the other model.

Discussion.

Due to the specific implementation of the ionization/recombination protocol, our classical MD method gives access to the potential surrounding an ion accounting for the influence of the free electrons and neighboring ions. By comparison with the ionization potential of the equivalent isolated ions, the ionization potential depression can be estimated. Our simulations provide data for further discussion on IPD models.

Moreover, the local electric fields produced by the electrons and ions of the plasma, which derived from the force calculations on individual charges, can be measured by accounting for all the interactions between charges. Characterizing properly these electric fields is very important for spectroscopic diagnostics as they interact with the emitters and modify the radiated spectra. The TCP simulation is well designed to obtain time sequences of electric field. These can be used in the resolution of the Schrödinger equation describing the time evolution of the emitter wave functions in the time-dependent field, or to infer, via a statistical analysis, the field distribution and correlation functions, characterizing the contributions of ions or electrons, or both [17]. As Figure 5. Ionic charge and field histories in an aluminum plasma at T e =T i =20eV and ρ = 0.65g.cm -3 . Bare ionic charge -dot dash line -; Ionic charge accounting for trapped electrons (effective charge)-dash line-; Electronic electric field measured at the bare ionic charge -dot line -; Electronic electric field seen by effective charges -full lineit has been mentioned previously, it is possible to identify, during a simulation run, the trapped and free electrons and thus to measure the electric fields at a given charge. Examples of two electronic electric field histories measured at two different bare ionic charges (Z = 6 and Z = 5) are shown in Figure 5(a) and (b), respectively, for the same plasma conditions as Figure 2. The dot line corresponds to the total electronic electric field at the ion and the full line corresponds to the electronic electric fields without the trapped electrons then accounted in the ion charge. It is also shown the value of the bare ionic charge (dot-dash line) and the ionic charge accounting for trapped electrons, the effective charge (dash line). A statistical analysis of the field histories at the effective charges will give access to statistical properties of electric fields, such as field distribution or correlation functions in plasmas composed of various ionic charges and electrons. The same field histories used in the resolution of the emitter evolution equation will permit to infer the influence of correlations on spectral line shapes. And finally, a statistical analysis of the ionization and recombination events will be performed to evaluate the corresponding rates. This work is in progress.

Figure 1 .

 1 Figure 1. Electronic total energy distribution function measured in aluminum plasmas at (a) T e =T i =20eV and ρ = 0.65g.cm -3 , (b) T e =T i =50eV and ρ = 2.7g.cm -3 and (c) T e =T i =550eV and ρ = 6.5g.cm -3 .

Figure 2 .

 2 Figure 2. Charge state distribution function in an aluminum plasma at T e =T i =20eV and ρ = 0.65g.cm -3 measured with (full line and stars) and without (dot-dash line) accounting for trapped electrons compared with the Flychk result (dash line and dots).

Figure 3 .

 3 Figure 3. Distribution of the ionization energy of an electron in the ground state of an ion of charge Z in an aluminum plasma at T e =T i =50eV and ρ = 2.7g.cm -3 .