- Email: [email protected]
Ship structure steel plate failure under near-field air-blast loading: Numerical simulations vs experiment
International Journal of Impact Engineering 62 (2013) 88e98
Contents lists available at SciVerse ScienceDirect
International Journal of Impact Engineering journal homepage: www.elsevier.com/locate/ijimpeng
Ship structure steel plate failure under near-field air-blast loading: Numerical simulations vs experiment Patrice Longère a, *, Anne-Gaëlle Geffroy-Grèze b, Bruno Leblé c, André Dragon d a Université de Toulouse, Institut Supérieur de l’Aéronautique et de l’Espace, Institut Clément Ader (EA 814), 10 avenue E. Belin, BP 54032, 31055 Toulouse cedex 4, France b DCNS NA Ingénierie, Rue de Choiseul, 56311 Lorient, France c DCNS Research, 44620 la Montagne, France d CNRS-Institut Pprime (UPR 3346), Ecole Nationale Supérieure de Mécanique et d’Aérotechnique, Université de Poitiers, 1 avenue C. Ader, BP 40109, 86961 Futuroscope - Chasseneuil du Poitou, France
a r t i c l e i n f o
a b s t r a c t
Article history: Received 7 December 2012 Received in revised form 12 June 2013 Accepted 14 June 2013 Available online 29 June 2013
This paper deals with the numerical simulation of the dynamic failure of a ship structure steel plate under near-field air-blast loading. Various energetic levels of air-blast loading, involving variable explosive mass and charge-plate distance, were tested leading to the bulge of the loaded plate for the lowest energy level and to the failure of the loaded plate for the highest level. A modified version of the GursoneTvergaardeNeedleman potential was used to reproduce the response of the material along the damage-plasticity process at stake. The 3D constitutive equations were implemented as user material in the engineering finite element computation code ABAQUSÒ, and numerical simulations were conducted and compared with experiments considering air-blast loaded plates. Several crucial numerical issues are addressed concerning notably the use of ABAQUSÒ conwep function, the hourglass control and the influence of the model constants. Numerical results clearly show the interest of the adopted modelling for the description of salient stages of dynamic structural failure. Ó 2013 Elsevier Ltd. All rights reserved.
Keywords: Air-blasted plate Ductile fracture Viscoplastic material Parametric study Numerical simulation
1. Introduction This paper addresses the numerical simulation of the dynamic failure of a ship structure steel under near-field air-blast loading. Air-blast tests belong to the experiments devoted to study the vunerability of ship hulls regarding explosion loading. Depending on the energetic conditions e induced by both explosive mass and charge-plate distance e one can typically observe the bulge of the structure for the lowest and the failure of the structure for the highest. Many experimental and numerical studies are available in literature aiming to reproduce the deflection of a plate (bulge) under low to moderate air-blast loading induced energy, see e.g. Refs. [1e4]. Conversely, works regarding numerical modelling attempting to describe all the air-blast post-loaded states, i.e. passing from the bulge to the ultimate failure, are scarce. This work aims at reproducing numerically the failure mechanisms of a ship structure constitutive material when submitted to air-blast loading.
* Corresponding author. E-mail addresses: [email protected] (P. Longère), anne-gaelle.greze@ dcnsgroup.com (A.-G. Geffroy-Grèze), [email protected] (B. Leblé), [email protected] (A. Dragon). 0734-743X/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijimpeng.2013.06.009
To qualify the high-purity, ferriticepearlitic mild steel employed as structural material in panels of naval structures under consideration regarding explosion loading, air-blast experiments were carried out. The samples were machined in the form of square thin plates, see Fig. 1a. The steel plate was held down by two frames fixed to the underlying structure. The plastic explosive charge was spherical, hung on a post and braces, see Fig. 1b. The mass of the charge and the distance between the charge and the plate are controlled parameters. The reader may refer to Ref. [5] for an extensive experimental investigation of the air-blast test. Depending on the mass C of the charge, and on the distance D between the charge and the plate, one can typically observe three states resulting from the explosion loading: the deflection of the plate (Fig. 2.1), the macrocracks incipience and growth (Fig. 2.2e 2.3), and the propagation of macrockacks resulting in the so-called petalling failure (Fig. 2.4). The maximum residual deflection of the plate is reported in Fig. 3 for variable charge mass and variable charge/plate distance. One can coarsely distinguish two domains: the first domain, covering long distance, for which the plate is able to consume the explosion induced energy by plastic deformation; the second domain, covering short distance, for which the explosion leads to the catastrophic failure of the plate. It must be noted that the notion of long distance and associated ‘far’ field is very
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
89
Fig. 1. Air-blast experimental set-up. (a) Schematic view. (b) Real view.
relative, as shown in Fig. 1b. The distance values on both sides of the transition, corresponding to the cracking initiation in Fig. 3, separating the domains of whole consumption of the explosion induced energy by plastic deformation (bulge) and catastrophic failure, are consequently very close. Predicting loading range and mode of ultimate structural failure implies identifying the underlying damage mechanisms at stake and describing their consequences. In the material at stake, damage has been seen to proceed from void initiation, growth and coalescence and to occur late in the deformation process. The
consequences of void growth induced damage are generally double: a progressive degradation of the overall properties of the bulk material and the appearance, in addition to the isochoric plastic deformation due to dislocation glide in the matrix material, of an inelastic dilatancy due to void growth. Gurson [6] microporous model is a micromechanics based model widely used for dealing with ductile fracture, see, e.g. Refs. [7,8]. Recently, Longere et al. [9] proposed an extended version of Gurson’s model in order to reproduce the delay experimentally observed between plastic deformation occurrence and hole nucleation, on one hand, and hole
Fig. 2. Final state of the plate after air-blast loading for various loading configurations (DCNS-DGA): constant charge mass C and decreasing distance D from 1 (far) to 4 (close).
Maximum residual deflection (mm)
90
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
Cracking initiation
Petalling
C2 = 2C1
and microporous metal plasticity and the definition of conditions for the microvoid induced damage initiation.
C1 (i) Yield stress: The yield stress sy includes a rate independent contribution sy and a strain rate induced overstress svp :
sy ¼ sy þ svp
Bulge
Close field
‘Far ’ field
Charge centre / plate distance (mm) Fig. 3. Maximum residual deflection as a function of the charge/plate distance for two charge mass values. Experimental results.
growth under null stress triaxiality, on the other hand, see also Longere and Dragon [10]. The three-dimensional constitutive equations corresponding to this model were implemented as user material in the engineering finite element computation code ABAQUSÒ, and numerical simulations were conducted and compared with experiments considering air-blast loaded ship structures. The strongly nonlinear three-dimensional model, accounting for the effects of strain rate, temperature and stress triaxiality on void initiation and growth, is outlined in Section 2. Section 3 is devoted to the numerical aspects of air-blast loaded structures. Several crucial numerical issues are addressed concerning notably the use of ABAQUSÒ conwep function, the hourglass control and the influence of the model constants. 2. Material behaviour constitutive modelling The material of the present work is a ferriticepearlitic mild steel (DH36) widely used in the automotive and railway industries as well as in ship building. An important experimental campaign including laboratory thermomechanical tests allowed for characterizing the behaviour of the ship structure constitutive material at stake under high strain rate conditions. This experimental study was used to develop a physically motivated modelling approach aiming at describing the combined effects of strain hardening, thermal softening, viscoplasticity and void growth induced damage until ultimate failure. In the sequel, only the most significant features of this strongly nonlinear three-dimensional model are presented. The reader may refer to Longere et al. [9] for a detailed presentation of the guiding concepts of the model, notably the concept of kinematic pressure shift related to microporosity growth. 2.1. Constitutive modelling In the present approach, as observed experimentally, defect nucleation is supposed to require a certain amount of plastic deformation. This implies distinguishing at least two steps in the process of ductile damage: the sound (undamaged) material behaves in elasticeviscoplastic manner and its response may be described using the J2-flow theory, i.e. via a HubereMises type yield criterion and potential; once defects nucleate and begin to grow, the (damaged) material may be considered as microporous and its response becomes pressure dependent. Such a scenario implies a transition between conventional dense metal plasticity
(1)
Based on the experimental campaign (not detailed here, see [9]), the rate independent contribution sy in (1) incorporates the combined effects of strain hardening, via a VOCE type law, and thermal softening, via a power law:
n
sy ¼ R0 þ RN ½1 expðkkÞb
o
1
T
m
Tmelt
(2)
where T represents the absolute temperature, k the isotropic strain hardening variable (or equivalent plastic strain), and where ðR0 ; RN ; k; bÞ are isotropic hardening related constants and ðTmelt ; mÞ thermal softening related constants, with Tmelt the melting point. The slight tension/compression asymmetry observed experimentally is here considered as resulting from a thermally activated mechanism involving the mean stress, see e.g. Graff et al. [11]. This approach leads to the following expression of the strain rate induced overstress svp in (1):
svp ¼ Y k_ exp
Va pm kB T
1=n (3)
where (Y,n) are viscosity related constants and (Va, kB) behaviour asymmetry related constants, with Va ¼ Vh [3 where Vh is a constant and [ Burgers vector magnitude ([ ¼ 2.5Å), and with kB Boltzmann constant (kB ¼ 1.3804 x 1023 J/K). (ii) Plastic potentials and evolution laws: In the following, we are considering the decomposition s ¼ pm d þ ð2=3Þseq n, with pm the pressure,pm ¼ Tr s=3, d the second order identity tensor, seq the equivalent stress, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi seq ¼ ð3=2Þs : s, s being the deviatoric part of the Cauchy stress tensor ðTr s ¼ 0Þ, and n the direction of the isochoric plastic yielding, n ¼ ð3=2Þs=seq . e Sound (undamaged) material: The viscoplastic yielding of the material in the undamaged state is supposed to be well described by the HubereMises type plastic potential
F0 ¼ s~2eq 1 ¼ 0; s~eq ¼
seq sy
(4)
The sound material satisfies the conditions for standard inelastic behaviour, i.e. is governed by the generalized normality rule in the framework of the irreversible thermodynamics internal variables formalism.
dp ¼ L
s~eq vF0 vF0 pD pD ¼ 3_ 0 n; 3_ 0 ¼ L ¼ 2L ; sy vs vseq rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 p p k_ ¼ d : d ¼ 3_ pD 0 ; L0 3
(5)
where L represents the viscoplastic multiplier. Finally, heating during any adiabatic processes is supposed to proceed
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
91
predominantly from dissipation Dsm, see Refs. [12,13] for further details, yielding
rC T_ ¼ Dsm ¼ seq r k_ 0
(6)
where r and C are the mass density and specific heat, respectively, and where r represents the isotropic hardening force b
r ¼ RN ½1 expðkkÞ
1
T
Tmelt
m :
e Microporous (damaged) material: The viscoplastic yielding of the material in the damaged state is supposed to obey the following modified version of the GTN model: Fig. 5. History of the extra energy due to ABAQUSÒ hourglass control procedure.
FG ~m p
3 ~m þ p ~ r Þ 1 þ q3 f 2 ¼ 0; ~2e0 q þ 2q1 f cosh q2 ðp ¼ s 2 pm pr ~r ¼ ¼ ; p
sy
3_
sy
pD G
(7)
s~e0 q pM vFG vF ¼ 2L ; 3_ G ¼ L G sy vseq vpm 3 ~ ~ sinh 2 q2 ðpm þ pr Þ ¼ 3q1 q2 f L ¼ L
sy
pr ¼ blnðq1 f Þ
(8)
where f represents the void volume fraction, (q1, q2, q3, b) are positive material constants, and pr represents the kinematic pressure shift related to microporosity growth, see [10]. This pressure shift is a specific modification introduced to the GTN model allowing to describe void growth near shear loading. When the kinematic pressure shift term is neglected in the pressure-dependent term of (7), the Gurson-GTN generalized form, corresponding to the hollow sphere model, including several simplifying approximations, is retrieved, see e.g. [10]. After Berg [14], the normality rule applies to the damaged material:
vF vFG 1 vFG d dp ¼ L G ¼ L n 3 vpm vs vseq
!
1 pM _ d ¼ 3_ pD G n þ 3 3G
(9)
According to (102), dilatation may occur even in absence of pressure, i.e. for null stress triaxiality values, as encountered in simple shearing. The evolution law of the isotropic hardening variable k is deduced from the equality of the macroscopic plastic work rate with the microscopic one, see Gurson [6]:
k_ ¼
_ pM seq 3_ pD G pm 3 G ð1 f Þsy
(11)
According to the aforementioned considerations, adiabatic heating is evaluated from
_ pM rC T_ ¼ seq 3_ pD k_ G r pm 3 G
(12)
The porosity rate f_ is decomposed into a contribution due to growth of existing defects and a contribution due to the formation of new defects, see (131). The former, namely f_ g, is deduced from the classical hypothesis of matrix incompressibility, see (132) below, whereas the latter, namely f_ n, will be expressed further, see (161).
Maximum residual deflection (mm)
_ pM where the distortional and dilatational parts, namely 3_ pD G and 3 G , p respectively, of the inelastic strain rate d are given by
(10)
EXP : C1 - D NUM : C1* - D - M=1.0
Charge centre / plate distance (mm)
Fig. 4. History of the air-blast induced pressure using the empirical conwep approximation for fluid/structure interaction in ABAQUSÒ.
Fig. 6. Maximum residual deflection as a function of the charge/plate distance for a constant charge mass C1, with C1* standing for TNT equivalent charge mass. Experimental vs. numerical results.
92
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
EXP : C1 - D Maximum residual deflection (mm)
Maximum residual deflection (mm)
EXP : C1 - D NUM : C1* - D - M=1.0 NUM : C1* - 1.62xD - M=1.0
NUM : C1* - D - M=1.0 NUM : C1* - D - M=0.7 NUM : C1* - D - M=0.6 NUM : C1* - D - M=0.5
Charge centre / plate distance (mm) Charge centre / plate distance (mm) Fig. 7. Maximum residual deflection as a function of the charge/plate distance for a constant charge mass C1, with C1* standing for TNT equivalent charge mass. Influence of the distance multiplication factor. Experimental vs. numerical results.
pM f_ ¼ f_ g þ f_ n ; f_ g ¼ ð1 f ÞTr dp ¼ ð1 f Þ_3 G ; fg ð0Þ ¼ f0
(13)
* ‘Primary’ hole nucleation criterion: To ensure the instantaneous transition between dense metal plasticity and microporous metal plasticity, the hole nucleation criterion and the microporous metal potential must have possibly the similar aspect and corresponding expressions. As a modified version of the GTN derived plastic potential is employed in this paper, see (7), the hole nucleation criterion FI0 is proposed in a close form:
3 2 pm pr ¼ ; b pr ¼
FI0 ¼ bs 2 þ 2q1 f0 cosh q2 ð bp m þ bp r Þ 1 þ q3 f02 ¼ 0; b s ¼
sy ; b pm sc
sc
sc (14)
In a first approximation, the dependence of the critical stress sc on the strain rate k_ and on the temperature T is postulated in a form close to (1). The critical stress sc is thus tentatively postulated in the simple form:
Fig. 9. Maximum residual deflection as a function of the charge/plate distance for a constant charge mass C1, with C1* standing for TNT equivalent charge mass. Influence of the pressure multiplication factor. Experimental vs. numerical results.
sc ¼ sI þ svp ; sI ¼ aðR0 þ RN Þ
(15)
A limited quantity f0 of microporosity is supposed to spurt instantaneously when the criterion (14) is satisfied. In such a phenomenological way is managed the (instantaneous) transition between dense metal plasticity and microporous metal plasticity, see [9]. As the form (14) involves both the yield stress and the pressure pm, hole nucleation is controlled by the stress triaxiality. Moreover, to describe the delaying effects of the strain rate and the temperature on the nucleation the critical stress sc should increase with increasing strain rate k_ and temperature T. As a first approximation the dependence of sc on k_ and T, see (15), is postulated in a form close to (1). Setting 0 < a < 1 favors the convergence of the dense metal potential with the microporous metal potential; tighter limits can be possibly set. The expression of the equivalent plastic strain at ‘primary’ nucleation k0 ¼ kðf0 Þ can be deduced from (2) and (14), see its physical soundness in [9]. It must be noted that there is some affinity with the Chu and Needleman viewpoint of the plastic strain controlled nucleation, see Ref. [15]. * ‘Secondary’ hole formation kinetics law: The ‘secondary’ voids include microvoids nucleating later than the initial, ‘primary’ holes. They include the nano-voids germinating between microvoids and favoring ultimately the coalescence by
NUM : (5/6)xC1* - D - M=1.0 NUM : (2/3)xC1* - D - M=1.0
Maximum residual deflection (mm)
Maximum residual deflection (mm)
EXP : C1 - D NUM : C1* - D - M=1.0
Petalling
D4
Cracking initiation
C2 = 2C1
D2
D3
D1
Bulge
Charge centre / plate distance (mm)
Charge centre / plate distance (mm) Fig. 8. Maximum residual deflection as a function of the charge/plate distance for a constant charge mass C1, with C1* standing for TNT equivalent charge mass. Influence of the charge mass multiplication factor. Experimental vs. numerical results.
Fig. 10. Maximum residual deflection as a function of the charge/plate distance for the charge mass C2. Experimental results.
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
localized micro-shearing. Consequently, the secondary voids nucleation kinetics is assumed to be controlled by the rate of hardening, see (16)1. This is consistent with the fact that hole germination (here, the secondary stage of it) necessitates a certain quota of plastic straining. As shown in [9], the secondary nucleation rate is zero at primary nucleation, then is progressively increasing until a maximum, and later decreases to zero. Based on the works by Molinari and Wright [16], the kinetic law for secondary nucleation is assumed to be well described by a Weibull type distribution function:
g g1 g FI0 f_ n ¼ B s_ y ; B ¼ fsup exp FI0 ; fn ð0Þ ¼ 0
sc
93
(16)
where g is a constant (g ¼ 2, as suggested by Molinari and Wright [16]) and where fsup represents the upper bound of the nucleated ‘secondary’ void volume fraction. h i represents McCaulay brackets. In a way consistent with the present modelling approach, the magnitude of the critical plastic strain at hole nucleation depends on strain rate, temperature and stress triaxiality.
Fig. 11. Equivalent plastic strain map for constant charge mass C2 and pressure multiplication factor M (M ¼ 1.0) and for variable charge/plate distance D. Influence of the ratio a with fsup ¼ 0.05 and fr ¼ 0.03. Rear views.
94
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
2.2. User material subroutine The constitutive model outlined in Subsection 2.1 was implemented as user subroutine (vumat) in the engineering finite element computation code ABAQUSÒ. The numerical integration is conducted in the GreeneNaghdi rotating frame using the classical return mapping procedure combined with the NewtoneRaphson solving algorithm, see Refs. [17,18] for further details. The thermal dilatation is supposed to be negligible in the present approach. Adiabatic conditions are furthermore assumed to be valid for plastic equivalent strain rate k_ greater than 1 s1. In addition, failure is supposed to occur as soon as the porosity reaches the critical value fr, leading numerically to the deletion of the concerned finite element. 3. Air-blast loading numerical simulation As mentioned in Subsection 2.2, the model outlined in Subsection 2.1 was implemented as user material in the engineering finite element computation code ABAQUSÒ in order to simulate numerically the air-blast loaded experiments. The 3D plate was meshed using 8-node bricks with reduced integration C3D8R, with five elements through the thickness, and the frames with rigid bodies, see the schematic view of the set-up in Fig. 1a. The surfaces common to the plate and the frames are tied and the vertical
translation of the lower frame is constrained. The numerical time integration scheme is explicit. 3.1. Fluid/structure interaction For air-blast loading numerical simulations, fluid/structure interaction is here numerically described using the card ‘conwep charge property’ available in ABAQUSÒ, considering an equivalent mass of TNT for the explosive charge. In the conwep approximation, the explosion is supposed to generate a spherical pressure wave whose history is represented in Fig. 4. The reader can refer to the works of Zakrisson et al. [4] for a comparison of the empirical blast loading conwep function with Eulerian computational analysis. It must be noted that the release of burnt gases is not permitted with ABAQUSÒ (6.10 version) conwep card, implying that explosion induced pressure keeps on being applied on the plate surface all along the loading duration. 3.2. Global energy balance To avoid zero-energy modes in 1-point integration finite element resulting in hourglass shape meshing deformation, ABAQUSÒ incorporates a small artificial stiffness. As shown in Fig. 5, the extra energy generated by the hourglass control procedure in a representative configuration of our numerical simulations
Fig. 12. Equivalent plastic strain map for constant charge mass C2 and pressure multiplication factor M (M ¼ 1.0) and for variable charge/plate distance D. Influence of the ratio a with fsup ¼ 0.05 and fr ¼ 0.03. Side views.
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
95
The sensitivity of the numerical results to the ABAQUSÒ conwep function related parameters on one hand and to the model constants on the other hand is studied in this subsection.
and numerical (NUM), maximum residual deflection values for variable charge/plate distance is shown in Fig. 6, where the triples (charge mass, charge-plate distance, pressure multiplication factor) are given for the numerical simulations as requested in ABAQUSÒ. One can clearly notice that results match for long distance (‘far’ field air-blast loading, as denoted in Fig. 3) and diverge when the distance becomes shorter (close field air-blast loading, as denoted in Fig. 3).
(i) Conwep function related parameters: A first set of numerical simulations was conducted in order to evaluate the ability of the conwep function to reproduce the fluid/structure interaction in the case of a plastically deformed structure. With this aim in view, the tests with the lightest charge (C1) were considered. The fact that no catastrophic failure was observed for these tests indeed authorized using the constitutive model outlined in Subsection 2.1 in a version without damage. The superposition of the experimental (EXP)
To palliate this well-known weakness of the conwep function to reproduce near-field air-blast loadings we are here searching to obtain a best fitting between experimental and numerical results by modifying the conwep function related parameter values. This aim is reached by increasing artificially the charge/plate distance D, see Fig. 7, reducing artificially the charge mass C* (C* standing for TNT equivalent charge mass), see Fig. 8, and/or reducing the maximum explosion induced pressure (see Fig. 4), see Fig. 9. To keep the initial conditions as close as possible to the experimental
reaches a maximum value of 2.5%, which may be considered as reasonable. 3.3. Air-blast loading numerical simulations: parametric study
Fig. 13. Equivalent plastic strain map for constant charge mass C2 and pressure multiplication factor M (M ¼ 0.7) and for variable charge/plate distance D. Influence of the porosity fsup with a ¼ 0.55 and fr ¼ 0.03. Rear views.
96
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
ones, the adaptation of the explosion induced pressure was retained with a pressure multiplication factor M ¼ 0.7. (ii) Microporous material model related constant: A second set of numerical simulations was conducted using the conwep function related parameters previously identified (a configuration with a pressure multiplication factor equal to unity is presented besides M ¼ 0.7) in order to evaluate the sensitivity of the numerical results to the microporous material model related constants. With this aim in view, the configurations considered (D1eD4) are reported in Fig. 10, involving the heaviest charge (C2 ¼ 2C1). To attempt to reproduce the catastrophic failure experimentally observed for the tests with the shortest distance, the constitutive model in Subsection 2.1 was used in its complete version, i.e. with damage. The parametric study consists herein in investigating the influence of the ratio a, entering the expression of the critical stress sc, see (152), see Figs. 11 and 12, the upper bound fsup of the nucleated ‘secondary’ void volume fraction, see Figs. 13 and 14, and the void volume fraction at failure fr , see Fig. 15. e The quantity a is directly proportional to the size of the hole non-nucleation domain (sc in Eq.15) relative to primary hole nucleation criterion. For primary holes to nucleate, the loading path has to reach/intercept the hole non-nucleation frontier (defined by the function FI0 ). For low values of a, the intersection is reached early and damage occurs early, whereas for larger values of a, the intersection is reached later and damage occurs later. Damage initiation will thus be more or less delayed depending on the value given to a.
e fsup represents the upper bound of the nucleated ‘secondary’ void volume fraction. As mentioned before, these ‘secondary’ voids include microvoids of the same nature of the ‘primary’ voids but nucleating later, as well as nano-voids germinating between macro-voids and being consequently at the origin of the coalescence by localised shearing. Increasing fsup consequently leads to an acceleration of the material damage accumulation and further reaching of the void volume fraction at fracture. e fr represents the void volume fraction at fracture triggering the numerical material deletion. Reducing fr consequently leads to a premature volume element failure and finite element deletion. As shown in Figs. 11 and 12, small values of fsup and fr combined with decreasing values of a lead to the formation and ejection of a cap, in one or several parts. As shown in Figs. 13 and 14, small values of a and fr combined with increasing values of fsup lead to the failure of the plate under petalling. As shown in Fig. 15, large values of a and fsup combined with increasing values of fr also lead to the formation and ejection of a cap, in one or several parts. This numerical investigation, summarized in Table 1, allows concluding that the model parameters play a significant role in the air-blast loaded plate failure, on one hand, and, on the other hand, that small values of a and fr combined with a moderate value of fsup are required to reproduce the plate failure experimentally observed in the form of petalling. A low value of a means that the delay between plastic deformation occurrence and damage initiation is relatively short, or equivalently that the
Fig. 14. Equivalent plastic strain map for constant charge mass C2 and pressure multiplication factor M (M ¼ 0.7) and for variable charge/plate distance D. Influence of the porosity fsup with a ¼ 0.55 and fr ¼ 0.03. Side views.
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
97
Fig. 15. Equivalent plastic strain map for constant charge mass C2, charge/plate distance D4 and pressure multiplication factor M (M ¼ 0.7). Influence of the critical porosity fr with a ¼ 0.65 and fsup ¼ 0.10. Perspective views.
plastic strain at damage initiation (which notably depends on the stress triaxiality) is small, around 0.1 according to Fig. 12 in [9]. Damage accumulation inside the material is insufficiently ensured by the growth of the primary voids and requires the nucleation of a non-negligible quantity of secondary voids, explaining that a moderate value of fsup (controlling the secondary void nucleation kinetics) is necessary. However, the void coalescence (not explicitly described in the present work) and further cracking stop suddenly, i.e. for a low value of void volume fraction at fracture fr, this damage accumulation leads to the ultimate ruin of the structure. 3.4. Air-blast loading numerical simulations: discussion The constant values in Longere et al. [9], see Table 2, i.e. identified from the thermomechanical characterization from laboratory tests, were first used to attempt to reproduce the air-blast loading experimental configurations, without success. The parametric study conducted above aimed at evaluating the sensitivity of the numerical results to the conwep function related parameters and to the model constants. A set of microporous material model related constants, reported in Table 3, was accordingly identified as being capable to reproduce as well as possible the petalling of the plate as observed experimentally, see Fig. 16. However, the initiation of cracking seems clearly to be the most difficult stage to reproduce. Table 3 indicates that the upper bound fsup of the nucleated ‘secondary’ void volume fraction, see (162), is greater than the Table 1 Summary of the parametric study. Figures
fsup
fr
a
Cap ejection
Figs. 11 and 12 Figs. 13 and 14 Fig. 15
Low High High
Low Low High
Low Low Moderate
X
Petalling X
X
porosity at failure fr. This is not paradoxical to the extent that fsup represents an upper bound which may be never reached: void coalescence and further meso-cracking (whose consequences are phenomenologically described via fr) may indeed occur much before the germination limit for secondary voids for a particular class of materials under some loading conditions, whereas the opposite may be observed for another class and/or other loading conditions. The numerical simulations of air-blast loading show that the equivalent plastic strain in the highly deformed area can reach 1.3, the maximum equivalent plastic strain rate 3700 s1 and the temperature value 470 C. These values are greater than the maximum values of strain rate and temperature considered for the material thermomechanical characterization. This may at least partially explain the difference between the set of model constants of Tables 1 and 2. In addition, as mentioned previously, the release of burnt gases is not possible when the ABAQUSÒ conwep card, reproducing numerically approximately the fluid/structure interaction, is used. This implies that the explosion induced pressure keeps on being
Table 2 Set of microporous model related constants for the numerical simulation of the tension loading of notched samples, after Longere et al. [9]. q1 1
q2 1
q3 1
f0 10
4
a
b (MPa)
fsup
fr
0.75
115
0.03
0.2
Table 3 Best set of microporous model related constants for the numerical simulation of the air-blast loading of steel plates. q1
q2
q3
f0
a
b (MPa)
fsup
fr
1
1
1
104
0.55
115
0.10
0.03
98
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
Fig. 16. Equivalent plastic strain map for variable charge mass C and charge/plate distance D and a constant pressure multiplication factor M (M ¼ 0.7). Set of model constants as reported in Table 2. Side views.
applied on the plate surface all along the loading duration, provoking a petal reversal more important than in experiment. Consequently, for a better description of the air-blast induced failure of the plate, a more accurate reproduction of the fluidstructure interaction is needed. Finally, even if the viscous feature of the model tends to regularize the initial-boundary value problem in the heating and damage induced softening regime, a mesh size dependence of the numerical results was observed e where the size of the smallest elements (corresponding to the largest mesh number) matches more or less the size of the elements used for the notched specimen meshing in the model constant identification in Longere et al. [9]. Works devoted to attenuate these meshing effects are currently in progress, by using notably enriched finite elements. 4. Concluding remarks The present work was focused on numerical parametric study and further simulation of the dynamic failure process of air-blast loaded steel plates under variable loading parameters. The mechanisms of damage and cracking for the material at stake were previously identified, see Geffroy et al. [5], from micrographic analyses of partially and totally fractured specimens. A multisurface based constitutive model able to describe the transition between dense metal plasticity and microporous metal plasticity has been summarised here. This model, detailed elsewhere (see Refs. [9,10]), reproduces the accelerating effects of stress triaxiality and the delaying effects of temperature and strain rate on damage evolution, and permits to describe the damage growth under shear loading. A parametric study was conducted aiming at evaluating the sensitivity of the numerical results to the conwep function related parameters (regarding fluid/structure interaction for plastically deformed structure) and to the model constants. A set of microporous material model related constants was accordingly identified as being capable to reproduce the petalling of the plate as observed experimentally for severe loading condition. However, the initiation of cracking seems clearly to be the most difficult stage to reproduce. The corresponding results are encouraging when compared qualitatively with experimental data. A better description of the fluid/structure interaction (in particular regarding gas release), numerically reproduced in the present work via the ABAQUSÒ Conwep card, see e.g. Ref. [19], as well as the attenuation of the mesh size dependence of the numerical results in the softening regime, could improve performance of the numerical model.
Acknowledgment The authors would like to acknowledge the financial contribution of the French Association for Research and Technology (ANRT). The authors express also gratitude for the contribution of DGA Naval Systems (DGA-tn), Toulon, France, to the air-blast tests. References [1] Jacinto AC, Ambrosini RD, Danesi RF. Experimental and computational analysis of plates under air blast loading. Int J Impact Eng 2001;25:927e47. [2] Neuberger A, Peles S, Rittel D. Scaling the response of circular plates subjected to large and close-range spherical explosions. Part I: air-blast loading. Int J Impact Eng 2007;34:859e73. [3] Neuberger A, Peles S, Rittel D. Springback of circular clamped armor steel plates subjected to spherical air-blast loading. Int J Impact Eng 2009;36:53e60. [4] Zakrisson B, Wikman B, Häggblad H-A. Numerical simulations of blast loads and structural deformation from near-field explosions in air. Int J Impact Eng 2011;38:597e612. [5] Geffroy A-G, Longère P, Leblé B. Fracture analysis and constitutive modelling of ship structure steel behaviour regarding underwater explosion. Eng Fail Anal 2011;18:670e81. [6] Gurson AL. Continuum theory of ductile rupture by void nucleation and growth: part I e yield criteria and flow rules for porous ductile media. J Eng Mat Tech ASME 1977;99:2e15. [7] Tvergaard V, Needleman A. Analysis of the cup-cone fracture in a round tensile bar. Acta Metall 1984;32(1):157e69. [8] Becker R, Needleman A, Richmond O, Tvergaard V. Void growth and failure in notched bars. J Mech Phys Solids 1988;36:317e51. [9] Longère P, Geffroy A-G, Leblé B, Dragon A. Modelling the transition between dense metal and damaged (micro-porous) metal viscoplasticity. Int J Dam Mech 2012;21:1020e63. [10] Longère P, Dragon A. Description of shear failure in ductile metals via back stress concept linked to damage-microporosity softening. Eng Fract Mech 2013;98:92e108. [11] Graff S, Forest S, Strudel S, Prioul J-L, Pilvin P, Béchade J-L. Strain localization phenomena associated with static and dynamic strain ageing in notched specimens: experiments and finite element simulations. Mat Sci Eng A 2004: 387e9. p. 181e5. [12] Longère P, Dragon A. Evaluation of the inelastic heat fraction in the context of microstructure supported dynamic plasticity modelling. Int J Impact Eng 2008;35(9):992e9. [13] Longère P, Dragon A. Inelastic heat fraction evaluation for engineering problems involving dynamic plastic localization phenomena. J Mech Mat Struct 2009;4(2):319e49. [14] Berg CA. Plastic dilation and void interaction. Proceedings of the Batelle memorial institute symposium on inelastic processes in solids; 1969. p. 171e209. [15] Chu CC, Needleman A. Void nucleation effects in biaxially stretched sheets. J Eng Mat Tech ASME 1980;102:249e56. [16] Molinari A, Wright TW. A physical model for nucleation and early growth of voids in ductile materials under dynamic loading. J Mech Phys Solids 2005;53: 1476e504. [17] Aravas N. On the numerical integration of a class of pressure-dependent plasticity models. Int J Num Meth Eng 1987;24:1395e416. [18] Vadillo G, Zaera R, Fernandez-Saez J. Consistent integration of the constitutive equations of Gurson materials under adiabatic conditions. Comput Methods Appl Mech Eng 2008;197:1280e95. [19] Langrand B, Leconte N, Menegazzi A, Millot T. Submarine hull integrity under blast loading. Int J Impact Eng 2009;36:1070e8.
Contents lists available at SciVerse ScienceDirect
International Journal of Impact Engineering journal homepage: www.elsevier.com/locate/ijimpeng
Ship structure steel plate failure under near-field air-blast loading: Numerical simulations vs experiment Patrice Longère a, *, Anne-Gaëlle Geffroy-Grèze b, Bruno Leblé c, André Dragon d a Université de Toulouse, Institut Supérieur de l’Aéronautique et de l’Espace, Institut Clément Ader (EA 814), 10 avenue E. Belin, BP 54032, 31055 Toulouse cedex 4, France b DCNS NA Ingénierie, Rue de Choiseul, 56311 Lorient, France c DCNS Research, 44620 la Montagne, France d CNRS-Institut Pprime (UPR 3346), Ecole Nationale Supérieure de Mécanique et d’Aérotechnique, Université de Poitiers, 1 avenue C. Ader, BP 40109, 86961 Futuroscope - Chasseneuil du Poitou, France
a r t i c l e i n f o
a b s t r a c t
Article history: Received 7 December 2012 Received in revised form 12 June 2013 Accepted 14 June 2013 Available online 29 June 2013
This paper deals with the numerical simulation of the dynamic failure of a ship structure steel plate under near-field air-blast loading. Various energetic levels of air-blast loading, involving variable explosive mass and charge-plate distance, were tested leading to the bulge of the loaded plate for the lowest energy level and to the failure of the loaded plate for the highest level. A modified version of the GursoneTvergaardeNeedleman potential was used to reproduce the response of the material along the damage-plasticity process at stake. The 3D constitutive equations were implemented as user material in the engineering finite element computation code ABAQUSÒ, and numerical simulations were conducted and compared with experiments considering air-blast loaded plates. Several crucial numerical issues are addressed concerning notably the use of ABAQUSÒ conwep function, the hourglass control and the influence of the model constants. Numerical results clearly show the interest of the adopted modelling for the description of salient stages of dynamic structural failure. Ó 2013 Elsevier Ltd. All rights reserved.
Keywords: Air-blasted plate Ductile fracture Viscoplastic material Parametric study Numerical simulation
1. Introduction This paper addresses the numerical simulation of the dynamic failure of a ship structure steel under near-field air-blast loading. Air-blast tests belong to the experiments devoted to study the vunerability of ship hulls regarding explosion loading. Depending on the energetic conditions e induced by both explosive mass and charge-plate distance e one can typically observe the bulge of the structure for the lowest and the failure of the structure for the highest. Many experimental and numerical studies are available in literature aiming to reproduce the deflection of a plate (bulge) under low to moderate air-blast loading induced energy, see e.g. Refs. [1e4]. Conversely, works regarding numerical modelling attempting to describe all the air-blast post-loaded states, i.e. passing from the bulge to the ultimate failure, are scarce. This work aims at reproducing numerically the failure mechanisms of a ship structure constitutive material when submitted to air-blast loading.
* Corresponding author. E-mail addresses: [email protected] (P. Longère), anne-gaelle.greze@ dcnsgroup.com (A.-G. Geffroy-Grèze), [email protected] (B. Leblé), [email protected] (A. Dragon). 0734-743X/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijimpeng.2013.06.009
To qualify the high-purity, ferriticepearlitic mild steel employed as structural material in panels of naval structures under consideration regarding explosion loading, air-blast experiments were carried out. The samples were machined in the form of square thin plates, see Fig. 1a. The steel plate was held down by two frames fixed to the underlying structure. The plastic explosive charge was spherical, hung on a post and braces, see Fig. 1b. The mass of the charge and the distance between the charge and the plate are controlled parameters. The reader may refer to Ref. [5] for an extensive experimental investigation of the air-blast test. Depending on the mass C of the charge, and on the distance D between the charge and the plate, one can typically observe three states resulting from the explosion loading: the deflection of the plate (Fig. 2.1), the macrocracks incipience and growth (Fig. 2.2e 2.3), and the propagation of macrockacks resulting in the so-called petalling failure (Fig. 2.4). The maximum residual deflection of the plate is reported in Fig. 3 for variable charge mass and variable charge/plate distance. One can coarsely distinguish two domains: the first domain, covering long distance, for which the plate is able to consume the explosion induced energy by plastic deformation; the second domain, covering short distance, for which the explosion leads to the catastrophic failure of the plate. It must be noted that the notion of long distance and associated ‘far’ field is very
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
89
Fig. 1. Air-blast experimental set-up. (a) Schematic view. (b) Real view.
relative, as shown in Fig. 1b. The distance values on both sides of the transition, corresponding to the cracking initiation in Fig. 3, separating the domains of whole consumption of the explosion induced energy by plastic deformation (bulge) and catastrophic failure, are consequently very close. Predicting loading range and mode of ultimate structural failure implies identifying the underlying damage mechanisms at stake and describing their consequences. In the material at stake, damage has been seen to proceed from void initiation, growth and coalescence and to occur late in the deformation process. The
consequences of void growth induced damage are generally double: a progressive degradation of the overall properties of the bulk material and the appearance, in addition to the isochoric plastic deformation due to dislocation glide in the matrix material, of an inelastic dilatancy due to void growth. Gurson [6] microporous model is a micromechanics based model widely used for dealing with ductile fracture, see, e.g. Refs. [7,8]. Recently, Longere et al. [9] proposed an extended version of Gurson’s model in order to reproduce the delay experimentally observed between plastic deformation occurrence and hole nucleation, on one hand, and hole
Fig. 2. Final state of the plate after air-blast loading for various loading configurations (DCNS-DGA): constant charge mass C and decreasing distance D from 1 (far) to 4 (close).
Maximum residual deflection (mm)
90
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
Cracking initiation
Petalling
C2 = 2C1
and microporous metal plasticity and the definition of conditions for the microvoid induced damage initiation.
C1 (i) Yield stress: The yield stress sy includes a rate independent contribution sy and a strain rate induced overstress svp :
sy ¼ sy þ svp
Bulge
Close field
‘Far ’ field
Charge centre / plate distance (mm) Fig. 3. Maximum residual deflection as a function of the charge/plate distance for two charge mass values. Experimental results.
growth under null stress triaxiality, on the other hand, see also Longere and Dragon [10]. The three-dimensional constitutive equations corresponding to this model were implemented as user material in the engineering finite element computation code ABAQUSÒ, and numerical simulations were conducted and compared with experiments considering air-blast loaded ship structures. The strongly nonlinear three-dimensional model, accounting for the effects of strain rate, temperature and stress triaxiality on void initiation and growth, is outlined in Section 2. Section 3 is devoted to the numerical aspects of air-blast loaded structures. Several crucial numerical issues are addressed concerning notably the use of ABAQUSÒ conwep function, the hourglass control and the influence of the model constants. 2. Material behaviour constitutive modelling The material of the present work is a ferriticepearlitic mild steel (DH36) widely used in the automotive and railway industries as well as in ship building. An important experimental campaign including laboratory thermomechanical tests allowed for characterizing the behaviour of the ship structure constitutive material at stake under high strain rate conditions. This experimental study was used to develop a physically motivated modelling approach aiming at describing the combined effects of strain hardening, thermal softening, viscoplasticity and void growth induced damage until ultimate failure. In the sequel, only the most significant features of this strongly nonlinear three-dimensional model are presented. The reader may refer to Longere et al. [9] for a detailed presentation of the guiding concepts of the model, notably the concept of kinematic pressure shift related to microporosity growth. 2.1. Constitutive modelling In the present approach, as observed experimentally, defect nucleation is supposed to require a certain amount of plastic deformation. This implies distinguishing at least two steps in the process of ductile damage: the sound (undamaged) material behaves in elasticeviscoplastic manner and its response may be described using the J2-flow theory, i.e. via a HubereMises type yield criterion and potential; once defects nucleate and begin to grow, the (damaged) material may be considered as microporous and its response becomes pressure dependent. Such a scenario implies a transition between conventional dense metal plasticity
(1)
Based on the experimental campaign (not detailed here, see [9]), the rate independent contribution sy in (1) incorporates the combined effects of strain hardening, via a VOCE type law, and thermal softening, via a power law:
n
sy ¼ R0 þ RN ½1 expðkkÞb
o
1
T
m
Tmelt
(2)
where T represents the absolute temperature, k the isotropic strain hardening variable (or equivalent plastic strain), and where ðR0 ; RN ; k; bÞ are isotropic hardening related constants and ðTmelt ; mÞ thermal softening related constants, with Tmelt the melting point. The slight tension/compression asymmetry observed experimentally is here considered as resulting from a thermally activated mechanism involving the mean stress, see e.g. Graff et al. [11]. This approach leads to the following expression of the strain rate induced overstress svp in (1):
svp ¼ Y k_ exp
Va pm kB T
1=n (3)
where (Y,n) are viscosity related constants and (Va, kB) behaviour asymmetry related constants, with Va ¼ Vh [3 where Vh is a constant and [ Burgers vector magnitude ([ ¼ 2.5Å), and with kB Boltzmann constant (kB ¼ 1.3804 x 1023 J/K). (ii) Plastic potentials and evolution laws: In the following, we are considering the decomposition s ¼ pm d þ ð2=3Þseq n, with pm the pressure,pm ¼ Tr s=3, d the second order identity tensor, seq the equivalent stress, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi seq ¼ ð3=2Þs : s, s being the deviatoric part of the Cauchy stress tensor ðTr s ¼ 0Þ, and n the direction of the isochoric plastic yielding, n ¼ ð3=2Þs=seq . e Sound (undamaged) material: The viscoplastic yielding of the material in the undamaged state is supposed to be well described by the HubereMises type plastic potential
F0 ¼ s~2eq 1 ¼ 0; s~eq ¼
seq sy
(4)
The sound material satisfies the conditions for standard inelastic behaviour, i.e. is governed by the generalized normality rule in the framework of the irreversible thermodynamics internal variables formalism.
dp ¼ L
s~eq vF0 vF0 pD pD ¼ 3_ 0 n; 3_ 0 ¼ L ¼ 2L ; sy vs vseq rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 p p k_ ¼ d : d ¼ 3_ pD 0 ; L0 3
(5)
where L represents the viscoplastic multiplier. Finally, heating during any adiabatic processes is supposed to proceed
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
91
predominantly from dissipation Dsm, see Refs. [12,13] for further details, yielding
rC T_ ¼ Dsm ¼ seq r k_ 0
(6)
where r and C are the mass density and specific heat, respectively, and where r represents the isotropic hardening force b
r ¼ RN ½1 expðkkÞ
1
T
Tmelt
m :
e Microporous (damaged) material: The viscoplastic yielding of the material in the damaged state is supposed to obey the following modified version of the GTN model: Fig. 5. History of the extra energy due to ABAQUSÒ hourglass control procedure.
FG ~m p
3 ~m þ p ~ r Þ 1 þ q3 f 2 ¼ 0; ~2e0 q þ 2q1 f cosh q2 ðp ¼ s 2 pm pr ~r ¼ ¼ ; p
sy
3_
sy
pD G
(7)
s~e0 q pM vFG vF ¼ 2L ; 3_ G ¼ L G sy vseq vpm 3 ~ ~ sinh 2 q2 ðpm þ pr Þ ¼ 3q1 q2 f L ¼ L
sy
pr ¼ blnðq1 f Þ
(8)
where f represents the void volume fraction, (q1, q2, q3, b) are positive material constants, and pr represents the kinematic pressure shift related to microporosity growth, see [10]. This pressure shift is a specific modification introduced to the GTN model allowing to describe void growth near shear loading. When the kinematic pressure shift term is neglected in the pressure-dependent term of (7), the Gurson-GTN generalized form, corresponding to the hollow sphere model, including several simplifying approximations, is retrieved, see e.g. [10]. After Berg [14], the normality rule applies to the damaged material:
vF vFG 1 vFG d dp ¼ L G ¼ L n 3 vpm vs vseq
!
1 pM _ d ¼ 3_ pD G n þ 3 3G
(9)
According to (102), dilatation may occur even in absence of pressure, i.e. for null stress triaxiality values, as encountered in simple shearing. The evolution law of the isotropic hardening variable k is deduced from the equality of the macroscopic plastic work rate with the microscopic one, see Gurson [6]:
k_ ¼
_ pM seq 3_ pD G pm 3 G ð1 f Þsy
(11)
According to the aforementioned considerations, adiabatic heating is evaluated from
_ pM rC T_ ¼ seq 3_ pD k_ G r pm 3 G
(12)
The porosity rate f_ is decomposed into a contribution due to growth of existing defects and a contribution due to the formation of new defects, see (131). The former, namely f_ g, is deduced from the classical hypothesis of matrix incompressibility, see (132) below, whereas the latter, namely f_ n, will be expressed further, see (161).
Maximum residual deflection (mm)
_ pM where the distortional and dilatational parts, namely 3_ pD G and 3 G , p respectively, of the inelastic strain rate d are given by
(10)
EXP : C1 - D NUM : C1* - D - M=1.0
Charge centre / plate distance (mm)
Fig. 4. History of the air-blast induced pressure using the empirical conwep approximation for fluid/structure interaction in ABAQUSÒ.
Fig. 6. Maximum residual deflection as a function of the charge/plate distance for a constant charge mass C1, with C1* standing for TNT equivalent charge mass. Experimental vs. numerical results.
92
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
EXP : C1 - D Maximum residual deflection (mm)
Maximum residual deflection (mm)
EXP : C1 - D NUM : C1* - D - M=1.0 NUM : C1* - 1.62xD - M=1.0
NUM : C1* - D - M=1.0 NUM : C1* - D - M=0.7 NUM : C1* - D - M=0.6 NUM : C1* - D - M=0.5
Charge centre / plate distance (mm) Charge centre / plate distance (mm) Fig. 7. Maximum residual deflection as a function of the charge/plate distance for a constant charge mass C1, with C1* standing for TNT equivalent charge mass. Influence of the distance multiplication factor. Experimental vs. numerical results.
pM f_ ¼ f_ g þ f_ n ; f_ g ¼ ð1 f ÞTr dp ¼ ð1 f Þ_3 G ; fg ð0Þ ¼ f0
(13)
* ‘Primary’ hole nucleation criterion: To ensure the instantaneous transition between dense metal plasticity and microporous metal plasticity, the hole nucleation criterion and the microporous metal potential must have possibly the similar aspect and corresponding expressions. As a modified version of the GTN derived plastic potential is employed in this paper, see (7), the hole nucleation criterion FI0 is proposed in a close form:
3 2 pm pr ¼ ; b pr ¼
FI0 ¼ bs 2 þ 2q1 f0 cosh q2 ð bp m þ bp r Þ 1 þ q3 f02 ¼ 0; b s ¼
sy ; b pm sc
sc
sc (14)
In a first approximation, the dependence of the critical stress sc on the strain rate k_ and on the temperature T is postulated in a form close to (1). The critical stress sc is thus tentatively postulated in the simple form:
Fig. 9. Maximum residual deflection as a function of the charge/plate distance for a constant charge mass C1, with C1* standing for TNT equivalent charge mass. Influence of the pressure multiplication factor. Experimental vs. numerical results.
sc ¼ sI þ svp ; sI ¼ aðR0 þ RN Þ
(15)
A limited quantity f0 of microporosity is supposed to spurt instantaneously when the criterion (14) is satisfied. In such a phenomenological way is managed the (instantaneous) transition between dense metal plasticity and microporous metal plasticity, see [9]. As the form (14) involves both the yield stress and the pressure pm, hole nucleation is controlled by the stress triaxiality. Moreover, to describe the delaying effects of the strain rate and the temperature on the nucleation the critical stress sc should increase with increasing strain rate k_ and temperature T. As a first approximation the dependence of sc on k_ and T, see (15), is postulated in a form close to (1). Setting 0 < a < 1 favors the convergence of the dense metal potential with the microporous metal potential; tighter limits can be possibly set. The expression of the equivalent plastic strain at ‘primary’ nucleation k0 ¼ kðf0 Þ can be deduced from (2) and (14), see its physical soundness in [9]. It must be noted that there is some affinity with the Chu and Needleman viewpoint of the plastic strain controlled nucleation, see Ref. [15]. * ‘Secondary’ hole formation kinetics law: The ‘secondary’ voids include microvoids nucleating later than the initial, ‘primary’ holes. They include the nano-voids germinating between microvoids and favoring ultimately the coalescence by
NUM : (5/6)xC1* - D - M=1.0 NUM : (2/3)xC1* - D - M=1.0
Maximum residual deflection (mm)
Maximum residual deflection (mm)
EXP : C1 - D NUM : C1* - D - M=1.0
Petalling
D4
Cracking initiation
C2 = 2C1
D2
D3
D1
Bulge
Charge centre / plate distance (mm)
Charge centre / plate distance (mm) Fig. 8. Maximum residual deflection as a function of the charge/plate distance for a constant charge mass C1, with C1* standing for TNT equivalent charge mass. Influence of the charge mass multiplication factor. Experimental vs. numerical results.
Fig. 10. Maximum residual deflection as a function of the charge/plate distance for the charge mass C2. Experimental results.
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
localized micro-shearing. Consequently, the secondary voids nucleation kinetics is assumed to be controlled by the rate of hardening, see (16)1. This is consistent with the fact that hole germination (here, the secondary stage of it) necessitates a certain quota of plastic straining. As shown in [9], the secondary nucleation rate is zero at primary nucleation, then is progressively increasing until a maximum, and later decreases to zero. Based on the works by Molinari and Wright [16], the kinetic law for secondary nucleation is assumed to be well described by a Weibull type distribution function:
g g1 g FI0 f_ n ¼ B s_ y ; B ¼ fsup exp FI0 ; fn ð0Þ ¼ 0
sc
93
(16)
where g is a constant (g ¼ 2, as suggested by Molinari and Wright [16]) and where fsup represents the upper bound of the nucleated ‘secondary’ void volume fraction. h i represents McCaulay brackets. In a way consistent with the present modelling approach, the magnitude of the critical plastic strain at hole nucleation depends on strain rate, temperature and stress triaxiality.
Fig. 11. Equivalent plastic strain map for constant charge mass C2 and pressure multiplication factor M (M ¼ 1.0) and for variable charge/plate distance D. Influence of the ratio a with fsup ¼ 0.05 and fr ¼ 0.03. Rear views.
94
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
2.2. User material subroutine The constitutive model outlined in Subsection 2.1 was implemented as user subroutine (vumat) in the engineering finite element computation code ABAQUSÒ. The numerical integration is conducted in the GreeneNaghdi rotating frame using the classical return mapping procedure combined with the NewtoneRaphson solving algorithm, see Refs. [17,18] for further details. The thermal dilatation is supposed to be negligible in the present approach. Adiabatic conditions are furthermore assumed to be valid for plastic equivalent strain rate k_ greater than 1 s1. In addition, failure is supposed to occur as soon as the porosity reaches the critical value fr, leading numerically to the deletion of the concerned finite element. 3. Air-blast loading numerical simulation As mentioned in Subsection 2.2, the model outlined in Subsection 2.1 was implemented as user material in the engineering finite element computation code ABAQUSÒ in order to simulate numerically the air-blast loaded experiments. The 3D plate was meshed using 8-node bricks with reduced integration C3D8R, with five elements through the thickness, and the frames with rigid bodies, see the schematic view of the set-up in Fig. 1a. The surfaces common to the plate and the frames are tied and the vertical
translation of the lower frame is constrained. The numerical time integration scheme is explicit. 3.1. Fluid/structure interaction For air-blast loading numerical simulations, fluid/structure interaction is here numerically described using the card ‘conwep charge property’ available in ABAQUSÒ, considering an equivalent mass of TNT for the explosive charge. In the conwep approximation, the explosion is supposed to generate a spherical pressure wave whose history is represented in Fig. 4. The reader can refer to the works of Zakrisson et al. [4] for a comparison of the empirical blast loading conwep function with Eulerian computational analysis. It must be noted that the release of burnt gases is not permitted with ABAQUSÒ (6.10 version) conwep card, implying that explosion induced pressure keeps on being applied on the plate surface all along the loading duration. 3.2. Global energy balance To avoid zero-energy modes in 1-point integration finite element resulting in hourglass shape meshing deformation, ABAQUSÒ incorporates a small artificial stiffness. As shown in Fig. 5, the extra energy generated by the hourglass control procedure in a representative configuration of our numerical simulations
Fig. 12. Equivalent plastic strain map for constant charge mass C2 and pressure multiplication factor M (M ¼ 1.0) and for variable charge/plate distance D. Influence of the ratio a with fsup ¼ 0.05 and fr ¼ 0.03. Side views.
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
95
The sensitivity of the numerical results to the ABAQUSÒ conwep function related parameters on one hand and to the model constants on the other hand is studied in this subsection.
and numerical (NUM), maximum residual deflection values for variable charge/plate distance is shown in Fig. 6, where the triples (charge mass, charge-plate distance, pressure multiplication factor) are given for the numerical simulations as requested in ABAQUSÒ. One can clearly notice that results match for long distance (‘far’ field air-blast loading, as denoted in Fig. 3) and diverge when the distance becomes shorter (close field air-blast loading, as denoted in Fig. 3).
(i) Conwep function related parameters: A first set of numerical simulations was conducted in order to evaluate the ability of the conwep function to reproduce the fluid/structure interaction in the case of a plastically deformed structure. With this aim in view, the tests with the lightest charge (C1) were considered. The fact that no catastrophic failure was observed for these tests indeed authorized using the constitutive model outlined in Subsection 2.1 in a version without damage. The superposition of the experimental (EXP)
To palliate this well-known weakness of the conwep function to reproduce near-field air-blast loadings we are here searching to obtain a best fitting between experimental and numerical results by modifying the conwep function related parameter values. This aim is reached by increasing artificially the charge/plate distance D, see Fig. 7, reducing artificially the charge mass C* (C* standing for TNT equivalent charge mass), see Fig. 8, and/or reducing the maximum explosion induced pressure (see Fig. 4), see Fig. 9. To keep the initial conditions as close as possible to the experimental
reaches a maximum value of 2.5%, which may be considered as reasonable. 3.3. Air-blast loading numerical simulations: parametric study
Fig. 13. Equivalent plastic strain map for constant charge mass C2 and pressure multiplication factor M (M ¼ 0.7) and for variable charge/plate distance D. Influence of the porosity fsup with a ¼ 0.55 and fr ¼ 0.03. Rear views.
96
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
ones, the adaptation of the explosion induced pressure was retained with a pressure multiplication factor M ¼ 0.7. (ii) Microporous material model related constant: A second set of numerical simulations was conducted using the conwep function related parameters previously identified (a configuration with a pressure multiplication factor equal to unity is presented besides M ¼ 0.7) in order to evaluate the sensitivity of the numerical results to the microporous material model related constants. With this aim in view, the configurations considered (D1eD4) are reported in Fig. 10, involving the heaviest charge (C2 ¼ 2C1). To attempt to reproduce the catastrophic failure experimentally observed for the tests with the shortest distance, the constitutive model in Subsection 2.1 was used in its complete version, i.e. with damage. The parametric study consists herein in investigating the influence of the ratio a, entering the expression of the critical stress sc, see (152), see Figs. 11 and 12, the upper bound fsup of the nucleated ‘secondary’ void volume fraction, see Figs. 13 and 14, and the void volume fraction at failure fr , see Fig. 15. e The quantity a is directly proportional to the size of the hole non-nucleation domain (sc in Eq.15) relative to primary hole nucleation criterion. For primary holes to nucleate, the loading path has to reach/intercept the hole non-nucleation frontier (defined by the function FI0 ). For low values of a, the intersection is reached early and damage occurs early, whereas for larger values of a, the intersection is reached later and damage occurs later. Damage initiation will thus be more or less delayed depending on the value given to a.
e fsup represents the upper bound of the nucleated ‘secondary’ void volume fraction. As mentioned before, these ‘secondary’ voids include microvoids of the same nature of the ‘primary’ voids but nucleating later, as well as nano-voids germinating between macro-voids and being consequently at the origin of the coalescence by localised shearing. Increasing fsup consequently leads to an acceleration of the material damage accumulation and further reaching of the void volume fraction at fracture. e fr represents the void volume fraction at fracture triggering the numerical material deletion. Reducing fr consequently leads to a premature volume element failure and finite element deletion. As shown in Figs. 11 and 12, small values of fsup and fr combined with decreasing values of a lead to the formation and ejection of a cap, in one or several parts. As shown in Figs. 13 and 14, small values of a and fr combined with increasing values of fsup lead to the failure of the plate under petalling. As shown in Fig. 15, large values of a and fsup combined with increasing values of fr also lead to the formation and ejection of a cap, in one or several parts. This numerical investigation, summarized in Table 1, allows concluding that the model parameters play a significant role in the air-blast loaded plate failure, on one hand, and, on the other hand, that small values of a and fr combined with a moderate value of fsup are required to reproduce the plate failure experimentally observed in the form of petalling. A low value of a means that the delay between plastic deformation occurrence and damage initiation is relatively short, or equivalently that the
Fig. 14. Equivalent plastic strain map for constant charge mass C2 and pressure multiplication factor M (M ¼ 0.7) and for variable charge/plate distance D. Influence of the porosity fsup with a ¼ 0.55 and fr ¼ 0.03. Side views.
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
97
Fig. 15. Equivalent plastic strain map for constant charge mass C2, charge/plate distance D4 and pressure multiplication factor M (M ¼ 0.7). Influence of the critical porosity fr with a ¼ 0.65 and fsup ¼ 0.10. Perspective views.
plastic strain at damage initiation (which notably depends on the stress triaxiality) is small, around 0.1 according to Fig. 12 in [9]. Damage accumulation inside the material is insufficiently ensured by the growth of the primary voids and requires the nucleation of a non-negligible quantity of secondary voids, explaining that a moderate value of fsup (controlling the secondary void nucleation kinetics) is necessary. However, the void coalescence (not explicitly described in the present work) and further cracking stop suddenly, i.e. for a low value of void volume fraction at fracture fr, this damage accumulation leads to the ultimate ruin of the structure. 3.4. Air-blast loading numerical simulations: discussion The constant values in Longere et al. [9], see Table 2, i.e. identified from the thermomechanical characterization from laboratory tests, were first used to attempt to reproduce the air-blast loading experimental configurations, without success. The parametric study conducted above aimed at evaluating the sensitivity of the numerical results to the conwep function related parameters and to the model constants. A set of microporous material model related constants, reported in Table 3, was accordingly identified as being capable to reproduce as well as possible the petalling of the plate as observed experimentally, see Fig. 16. However, the initiation of cracking seems clearly to be the most difficult stage to reproduce. Table 3 indicates that the upper bound fsup of the nucleated ‘secondary’ void volume fraction, see (162), is greater than the Table 1 Summary of the parametric study. Figures
fsup
fr
a
Cap ejection
Figs. 11 and 12 Figs. 13 and 14 Fig. 15
Low High High
Low Low High
Low Low Moderate
X
Petalling X
X
porosity at failure fr. This is not paradoxical to the extent that fsup represents an upper bound which may be never reached: void coalescence and further meso-cracking (whose consequences are phenomenologically described via fr) may indeed occur much before the germination limit for secondary voids for a particular class of materials under some loading conditions, whereas the opposite may be observed for another class and/or other loading conditions. The numerical simulations of air-blast loading show that the equivalent plastic strain in the highly deformed area can reach 1.3, the maximum equivalent plastic strain rate 3700 s1 and the temperature value 470 C. These values are greater than the maximum values of strain rate and temperature considered for the material thermomechanical characterization. This may at least partially explain the difference between the set of model constants of Tables 1 and 2. In addition, as mentioned previously, the release of burnt gases is not possible when the ABAQUSÒ conwep card, reproducing numerically approximately the fluid/structure interaction, is used. This implies that the explosion induced pressure keeps on being
Table 2 Set of microporous model related constants for the numerical simulation of the tension loading of notched samples, after Longere et al. [9]. q1 1
q2 1
q3 1
f0 10
4
a
b (MPa)
fsup
fr
0.75
115
0.03
0.2
Table 3 Best set of microporous model related constants for the numerical simulation of the air-blast loading of steel plates. q1
q2
q3
f0
a
b (MPa)
fsup
fr
1
1
1
104
0.55
115
0.10
0.03
98
P. Longère et al. / International Journal of Impact Engineering 62 (2013) 88e98
Fig. 16. Equivalent plastic strain map for variable charge mass C and charge/plate distance D and a constant pressure multiplication factor M (M ¼ 0.7). Set of model constants as reported in Table 2. Side views.
applied on the plate surface all along the loading duration, provoking a petal reversal more important than in experiment. Consequently, for a better description of the air-blast induced failure of the plate, a more accurate reproduction of the fluidstructure interaction is needed. Finally, even if the viscous feature of the model tends to regularize the initial-boundary value problem in the heating and damage induced softening regime, a mesh size dependence of the numerical results was observed e where the size of the smallest elements (corresponding to the largest mesh number) matches more or less the size of the elements used for the notched specimen meshing in the model constant identification in Longere et al. [9]. Works devoted to attenuate these meshing effects are currently in progress, by using notably enriched finite elements. 4. Concluding remarks The present work was focused on numerical parametric study and further simulation of the dynamic failure process of air-blast loaded steel plates under variable loading parameters. The mechanisms of damage and cracking for the material at stake were previously identified, see Geffroy et al. [5], from micrographic analyses of partially and totally fractured specimens. A multisurface based constitutive model able to describe the transition between dense metal plasticity and microporous metal plasticity has been summarised here. This model, detailed elsewhere (see Refs. [9,10]), reproduces the accelerating effects of stress triaxiality and the delaying effects of temperature and strain rate on damage evolution, and permits to describe the damage growth under shear loading. A parametric study was conducted aiming at evaluating the sensitivity of the numerical results to the conwep function related parameters (regarding fluid/structure interaction for plastically deformed structure) and to the model constants. A set of microporous material model related constants was accordingly identified as being capable to reproduce the petalling of the plate as observed experimentally for severe loading condition. However, the initiation of cracking seems clearly to be the most difficult stage to reproduce. The corresponding results are encouraging when compared qualitatively with experimental data. A better description of the fluid/structure interaction (in particular regarding gas release), numerically reproduced in the present work via the ABAQUSÒ Conwep card, see e.g. Ref. [19], as well as the attenuation of the mesh size dependence of the numerical results in the softening regime, could improve performance of the numerical model.
Acknowledgment The authors would like to acknowledge the financial contribution of the French Association for Research and Technology (ANRT). The authors express also gratitude for the contribution of DGA Naval Systems (DGA-tn), Toulon, France, to the air-blast tests. References [1] Jacinto AC, Ambrosini RD, Danesi RF. Experimental and computational analysis of plates under air blast loading. Int J Impact Eng 2001;25:927e47. [2] Neuberger A, Peles S, Rittel D. Scaling the response of circular plates subjected to large and close-range spherical explosions. Part I: air-blast loading. Int J Impact Eng 2007;34:859e73. [3] Neuberger A, Peles S, Rittel D. Springback of circular clamped armor steel plates subjected to spherical air-blast loading. Int J Impact Eng 2009;36:53e60. [4] Zakrisson B, Wikman B, Häggblad H-A. Numerical simulations of blast loads and structural deformation from near-field explosions in air. Int J Impact Eng 2011;38:597e612. [5] Geffroy A-G, Longère P, Leblé B. Fracture analysis and constitutive modelling of ship structure steel behaviour regarding underwater explosion. Eng Fail Anal 2011;18:670e81. [6] Gurson AL. Continuum theory of ductile rupture by void nucleation and growth: part I e yield criteria and flow rules for porous ductile media. J Eng Mat Tech ASME 1977;99:2e15. [7] Tvergaard V, Needleman A. Analysis of the cup-cone fracture in a round tensile bar. Acta Metall 1984;32(1):157e69. [8] Becker R, Needleman A, Richmond O, Tvergaard V. Void growth and failure in notched bars. J Mech Phys Solids 1988;36:317e51. [9] Longère P, Geffroy A-G, Leblé B, Dragon A. Modelling the transition between dense metal and damaged (micro-porous) metal viscoplasticity. Int J Dam Mech 2012;21:1020e63. [10] Longère P, Dragon A. Description of shear failure in ductile metals via back stress concept linked to damage-microporosity softening. Eng Fract Mech 2013;98:92e108. [11] Graff S, Forest S, Strudel S, Prioul J-L, Pilvin P, Béchade J-L. Strain localization phenomena associated with static and dynamic strain ageing in notched specimens: experiments and finite element simulations. Mat Sci Eng A 2004: 387e9. p. 181e5. [12] Longère P, Dragon A. Evaluation of the inelastic heat fraction in the context of microstructure supported dynamic plasticity modelling. Int J Impact Eng 2008;35(9):992e9. [13] Longère P, Dragon A. Inelastic heat fraction evaluation for engineering problems involving dynamic plastic localization phenomena. J Mech Mat Struct 2009;4(2):319e49. [14] Berg CA. Plastic dilation and void interaction. Proceedings of the Batelle memorial institute symposium on inelastic processes in solids; 1969. p. 171e209. [15] Chu CC, Needleman A. Void nucleation effects in biaxially stretched sheets. J Eng Mat Tech ASME 1980;102:249e56. [16] Molinari A, Wright TW. A physical model for nucleation and early growth of voids in ductile materials under dynamic loading. J Mech Phys Solids 2005;53: 1476e504. [17] Aravas N. On the numerical integration of a class of pressure-dependent plasticity models. Int J Num Meth Eng 1987;24:1395e416. [18] Vadillo G, Zaera R, Fernandez-Saez J. Consistent integration of the constitutive equations of Gurson materials under adiabatic conditions. Comput Methods Appl Mech Eng 2008;197:1280e95. [19] Langrand B, Leconte N, Menegazzi A, Millot T. Submarine hull integrity under blast loading. Int J Impact Eng 2009;36:1070e8.