mercredi 30 avril 2008
Numerical simulation of gun powder and explosives
NUMERICAL SIMULATION OF THE DYNAMICAL MECHANICAL PROPERTIES OF HOMOGENEOUS AND GRANULAR ENERGETIC MATERIALS.
B. Schaeffer
25th Int. Annual Conference of ICT, june 28- July 1, 1994, Karlsruhe.
Abstract
Energetic materials are usually highly viscoelastic or viscoplastic. They are brittle at low temperatures or under impact. Some are composite materials and their crystalline filler is responsible for a dilatant behaviour.
The products made with such materials are subjected to static and dynamical mechanical loads during or after fabrication. In a gun, for example, acceleration loads are very high in the gun powder as well as in the explosive charge which may rub against the casement, be ignited and explode prematurely. Cracks increase, sometimes catastrophically, the combustion area. The breaking of the bonding between the propellant and the liner and decohesion between the binder and the crystalline filler in composite propellants or explosives contain are two other types of fracture that may influence combustion.
A microcomputer software (Deform2D) working on PC’s (Windows) and Macintosh has been developed for simulating the dynamical mechanical behaviour of non-linear materials, particularly viscoplasticity and fracture. The computing method is based on finite differences with automatic meshing. Bimaterial structures may be studied: adhesive joints and composite materials. Composite materials are defined very easily with the help of graphic patterns: one pixel is associated with one mesh. The mesh is made of one material or the other depending on the pixel. Anisotropy is the result of an inhomogeneous structure, anisotropically distributed. The bonding strength between the elastomer binder and the particles is approximated by the strength of the solid phase.
Numerical simulations of tension and compression tests on composite materials showing dilatancy have been realised. The propagation of elastic waves has been simulated in a bar and a powder bed. The computed results are compared with experiment.
1. INTRODUCTION
Numerical simulation of small propellant and explosive tests can help to understand the behaviour of full-scale weapon systems [1 ]. Data for code inputs are often scarce, particularly at high rates of straining. The time-temperature shift may be used to extrapolate low-rate tests, but the embrittlement under impact is also due to the finite speed of the stress waves producing a stress concentration at the point of impact [ ]. It is often assumed that wave propagation may be neglected in small specimens [3], but at high strain rates, for deformation speeds near or above 10 m/s, the strain is no more homogeneous [4]. The software presented here is particularly suited to simulate tests in the laboratory, static or dynamical, and shows effects like dilatancy related to the composite nature of many energetic materials. It will be applied to a few simple experiments.
2. NUMERICAL MODEL
Deform2D is a finite differences software where the specimen to be studied is divided into elementary parts, e.g. quadrilaterals (figure 1) each made of a given material. The material may be different from one mesh to the other, but the choice is between two materials only. In a mesh, stresses and strains are constant. The boundary conditions may be of given speed or given pressure. Unilateral contact with friction at rigid contours is modelled.
The finite strain tensor may be computed from the lengths of the sides of a triangle built on the diagonals of the quadrilateral meshes used to discretise the specimen.
The stresses are true stresses (relative to the instantaneous geometry) and not the nominal or engineering stresses (relative to the initial geometry), the fundamental laws of dynamics having to be applied to the actual, not the initial geometry.
The stresses are obtained from the strains through the constitutive law, a combination of hypoelastic and viscous behaviour. The resulting effort on an element (built by joining the four immediate neighbours of a node) is then calculated. Newton's law gives the acceleration, integrated twice to obtain the new coordinates of the nodes.
The computing cycle is repeated for each node and each time step. The Courant-Friedrich-Lewy criterion states that, for stable computation, the propagation of the elastic wave has to be smaller than one mesh during one time step. All calculations are here in plane strain. More details may be found elsewhere [5].
The numerical experiment, for example a tensile test on a sample hold in a fixed grip at the bottom and in a moving grip at the top (figure 1). Sample and grips are drawn with the mouse. The contours (sample, fixed, moving grips and also transducers, applied pressure and symmetry regions) are analysed by Deform2D after clicking inside them with the mouse (fig. 1a). The sample is meshed automatically. The numerical values are introduced through a dialog window (figure 2).
Figure 1 - At left, the drawing used to define the geometry; the contour of the part to be meshed is obtained simply by clicking in it with the mouse, the software analyses it and meshes it automatically. The mobile and fixed clamps are analysed in the same fashion after a slight modification of the drawing to make the contours simply connected. At right, numerical model of a tension test (schematic).
3. DILATANCY IN HIGHLY LOADED ELASTOMERS
When composite materials are deformed, there is debonding (localised fracture) between the binder and the reinforcement producing an increase in volume of the material. This will be simulated in tension and in compression.
3.1. Tension test
The test is schematised on figure 1. The speed of the moving grip is 1 m/s, as shown on figure 2. At this speed the test may be considered as quasi-static. Of course, a lower speed could be chosen, but the computing time being inversely proportional to the speed, it would take too much time. One night is necessary on a Macintosh II (500 nodes, 5000 computing cycles, result shown on figure 3).
Figure 2 - Legend on the next page.
Legend of figure 2
The first series of datas are material constants (maximum two constituants). Some numerical constants are not entered directly. For example, Poisson’s ratio is calculated from the bulk modulus and the shear modulus. If one knows Poisson’s ratio and Young’s modulus, it is necessary to vary the bulk and shear moduli until obtaining the desired values for Poisson’s ratio and Young’s modulus.
The upper and lower yield stresses have very high values: elastomers are usually not plastic.
The ultimate tensile strength is the true strength at fracture, not the engineering strength.
The tensile strength of the binder is taken very high so that no fracture occurs in the binder. On the contrary, the crystalline phase is supposed to break. This may look as a strange hypothesis, but it was necessary to use it for simulating the fracture at the binder and filler interface. When a crack is formed, there is an increase in volume. It may be detected from the volume change in a part of the specimen. If the whole specimen is chosen, it simulates a dilatometer (for example a Farris dilatometer).
The strain at fracture is the maximum local strain in tension, not the mean strain across the specimen (engineering strain). Usually, these values have to be adjusted by simulating a tensile test and verifying that the maximum load and extension correspond to the experimental values. The strain at fracture has no effect on the numerical results when there is no plasticity, the behaviour being supposed to be entirely elastic until fracture for both the binder and the filler. The strain at fracture is used to compute the (linear) strain-hardening coefficient.
The desired compressive strength will be obtained by adjusting the value of the friction coefficient, according to the Coulomb criterion of fracture. The friction coefficient between the specimen and the platens (in compression) is assumed to be the same as the internal friction coefficient.
The viscosity is rarely known from experiments. It has to be chosen empirically in order to have a realistic damping. The viscosity coefficient allows the simulation of flow at constant stress but not of relaxation experiments.
The speed of the longitudinal and transverse waves is calculated from the bulk modulus, shear modulus and specific mass, according to the classical formulae.
The next parameters are related to the numerical experiment itself.
The number of nodes determines the precision of the computation, but the computing times increases proportionally (in fact as the 2/3 power of the number of nodes).
The size of the specimen is given by the next three values: overall height, width and depth.
The interval between views is the time between graphic displays. Usually, the results of the calculations are visualised for one out of ten computing cycles.
The time step is the time corresponding to a computing cycle. It depends on the bulk modulus, the dimensions and the number of nodes.
The duration of the sollicitation is the time during which the speed of the moving part or the pressure is applied.
The maximum duration is the time at which the calculation stops. This is not the duration of the computation measured on a clock, but related to the time of the physical phenomenon simulated.
The safety factor is the inverse of the Courant-Friedrich-Levy number govering the stability of the calculation. It should never be larger than one.
The vertical acceleration is usually the the acceleration of gravity.
The ambiant pressure is usually the atmospheric pressure.
The applied pressure is the absolute pressure, if any, applied on the specimen, locally or all around.
For composite materials, a structure may be chosen with the help of a graphic pattern in an other dialog window. The density, the percentage of reinforcement in volume and weight are automatically calculated from the pattern.
Figure 3 -
Fractured tensile specimen. In black, the crystalline filler and in white, the binder. Fracture occurred at the centre of the specimen where high distortions are to be seen. Stress and dilatation as a function of mean strain. The behaviour is linear until fracture. The volume change is always positive and increases abruptly at decohesion, here immediately followed by fracture.
3.2. Compression test
A material stressed in tension under pressure or in compression should break at stresses and strains larger than in simple tension. This will be checked in this numerical compression experiment. A simple algorithm is used to simulate the contact between specimen and the rigid platens. The deformed specimen at complete fracture is shown on figure 4. The stress and the dilatation are shown on figure 5 in function of the strain.
Figure 4- Specimen shape at complete fracture. It is highly distorted as in a real material. One may distinguish the contours of the fractured regions. The compression speed of 10 m/s makes the deformation larger at the top of the specimen. The number of nodes is 200. The height of the sample is 0.02 m for a width of 0.01 m. All other data are the same as in tension above (figure 2).
Figure 5 - The compressive stress is increasing almost linearly in the elastic range. The plateau corresponds to the beginning of cracking (decohesion). The volume of the specimen decreases first, due to elastic compression, then increases with decohesion takes place.
3.3. Comparison with experiment
The agreement with the experimental curves [6], [7] is not too bad. In tension, the volume increases all the time, first linearly, as may be expected from Poisson’s ratio, then much more rapidly than observed experimentally.
The maximum strain is much higher in compression than in tension as in most materials. In compression, there is a decrease of the volume, due to elastic compression, followed by an increase, due to decohesion, as observed experimentally.
4. WAVE PROPAGATION IN AN EXPLOSIVE SIMULANT
4.1. Numerical model
A 150 mm length and 10 mm thickness bar, made of an explosive simulant (the same as above but entered in the software as an homogeneous material for the sake of simplicity) is impacted at one end at a speed of 1 m/s during 100 µs. The variation of the longitudinal stress at the other end is shown on figure 6.
Figure 6 - Oscillations in a bar impacted at one end. The computed curve (in black) is compared with the experimental curve (in gray)
4.2. Comparison with experiment
4.2.1. Experimental set-up
A 150 x 10 x 10 mm explosive simulant bar was equipped with an accelerometer connected to an Apple II microcomputer with a data acquisition board. The bar was impacted at one side with a pen and the signal from the accelerometer at the opposite side was recorded.
4.2.2. Viscosity
The data from the numerical simulation (gray curve) agree with the experimental results (black curve). They have been obtained with a viscosity of 1000 Pa.s. Hence the viscosity for a solid composite explosive should be around this value. The maximum measurable value during polymerization of composite solid propellants was found to be of 106 Poises (105 Pa.s) [ ] with a Brookfield viscosimeter. It is hundred times smaller; the reason for the discrepancy is not clear. Unfortunately, practically no information can be found on this subject in the literature. There are many papers about viscoelasticity, but the results are expressed in terms of damping or relaxation times. The logarithmic decrement is found to be about 0.3 from the experimental data.
4.2.3. Speed of sound
The longitudinal wave speed is found to be of 270 m/s for assumed speeds of 810 and 140 m/s respectively for longitudinal and transverse waves (the material is taken as homogeneous here). With the classical formula for bar wave speed one finds 235 m/s.
5. WAVE PROPAGATION IN GUN POWDER
5.1. Numerical model
In the numerical experiment, two transducers are simulated in the gunpowder bed at its top and at its bottom, on order to “measure” the speed of an elastic wave.
The powder bed has two constituents, nitrocellulose (reinforcement) and air (matrix). The elastic modulus of the solid phase is taken as 1.5 GPa, with a Poisson’s ratio of 0.45 and a density of 1600 kg/m3. The global density is 850 kg/m3. The height of powder is 0.05 m. The impact speed is 0.1 m/s, for a duration of 20 µs. The geometry of the numerical experiment is shown on figure 7. The non-linearities due to contact between the grains are not taken into account.
The computed recordings of the longitudinal stress at top and bottom of the powder bed are given on figure 8. From these one obtains the propagation time of the pulse from top to bottom.
Figure 7 - Numerical experiment for simulating a wave propagation in a powder bed.
The powder is contained in a container and compressed by a moving anvil. Two regions (dotted line) at the top and bottom simulate transducers for recording the stress as a function of time.
The stresses in the “transducers” are the mean stress in the regions defined by the contours that may be seen near the top and bottom of the powder bed.
Figure 8 - The propagation time of the elastic wave is obtained from the delay between the computed stress curves in the two “transducers”, e.g. 0.09 ms, corresponding to a wave speed of 400 m/s.
5.2. Comparison with experiment
Compression tests (figure 9) are performed in brass cylinders containing the powder. The compression anvil receives an impact from a Charpy pendulum, generating a pulse propagating through the powder bed. A quartz transducer detects the instant of impact. A strain-gage transducer at the bottom detects the incoming wave pulse.
The signals from the two transducers are recorded with an oscilloscope. From the delay between the records given by the two transducers the speed of the elastic wave is obtained.
Figure 9 - Experimental set-up for measuring wave speed in a powder bed. An impact of known energy is obtained with a Charpy pendulum. The effort is recorded on top with a quartz transducer and on bottom with a strain-gage transducer. Both signals are recorded on an oscilloscope. The wave speed is obtained from the delai between the two signals.
The speed of the wave varies between 179 and 354 m/s, depending on the weight and falling height of the Charpy impacter, the prestress, the height of the powder bed (figure 10). The speed given by the numerical simulation, 400 m/s, is higher than all the experimental values, which may be due to the contact non-linearities between the grains, not taken into account in the simulation.
Using the Schubert formula [9] one finds 240 m/s, with a Young’s modulus of 1.5 GPa and a Poisson’s ratio of 0.45 for the nitrocellulose, the porosity of the bed being 40 %. The assumptions are the same used by Schubert, e.g. radius of curvature at the contact point equal to the grain radius and a tension of 1 kPa.
Figure 10 - The speed of the wave increases when the energy of the impactor (weight and falling height) increases and when the height of powder and prestress decreases. The software is not sophisticated enough to predict such effects.
6. CONCLUSION
The computed curves of the deformation (stress and dilatancy) are qualitatively in accord with experiment for highly filled elastomers. Volume is always increasing in tension and, in compression, a decrease is followed by an increase due to decohesion.
The software shows how oscillations are generated in a composite explosive simulant bar by an impact. The study of the corresponding wave propagation shows the influence of viscosity, a parameter often cited in rheological studies but to which a numerical value is rarely given. The damping of the oscillations corresponds to a viscosity 100 times smaller than measured directly towards the end of polymerization of a composite propellant. The propagation speed of an elastic wave in a powder bed is reasonably in accord with experimental results although the contact phenomena were not modelled.
Deform2D is a simple tool to understand the static and dynamical behaviour of energetic materials. It may be useful to optimize an experimental program or a more sophisticated numerical simulation with a standard software. It appears the necessity of knowing the mechanical properties of not only the composites but also of their components (binder and crystalline phase). The same is valid for gun powder where one should determine the mechanical properties (static and dynamic) of the powder material and of the powder bed.
7. REFERENCES
vendredi 30 novembre 2007
Equilibrium shape of a tree branch
Abstract
The shape of tree branches is computed using the theory of flexure of beams. Numerical computation is used and takes into account large displacements and growth. The coupling between elastic bending and growth results in a branch shape having an inflection point connected with a permanent deformation.
Growth stresses
The study of stresses in living materials has been the subject of many papers. For example the stresses due to anisotropic growth [1], influence of stresses on growth [2], and even viscoelasticity [3] or plasticity. The influence of gravity on growth is known under the names of gravimorphism [4] or geotropism [5]. The proportions of trees have been found to be limited by elastic criteria to stand under their own weight [6]. Growth stresses and particularly those due to reaction wood play an important role in the mechanics of trees [8], but this paper is focused essentially on the physical influence of gravity. Growth stresses are predominant in the trunk of a vertical tree, where gravitational forces creating nearly hydrostatic pressures may be neglected [8]. In branches, where bending is dominant, stresses due to gravitational forces are much larger than in the trunk. In this paper, growth stresses will be considered as a distinct phenomenon that will not be taken into account in a first approximation.
Upright growth of trees
A forest tree that has been shifted from its normal upright position during a storm will probably grow straight again. Reaction wood may force the tree to an upright position again, but bending strain was not found to affect radial xylem growth in Douglas-fir [5]. Old branches have an inflection point and their curvature remaining almost unchanged when cut, their deformation is permanent. It is often not possible to straighten these branches without breaking them.
Bending of beams and of branches
In the absence of gravity a branch would probably grow straight in the direction of light. With growth rings of constant thickness and a constant yearly increase in length, the shape of a branch would be a cone. Trees may be considered as structures made of beams (dead branches) subjected to gravity. Classical Strength of Materials considers a beam with a given shape, applies a load and calculates the resulting shape, slightly different from the original one for thick beams. Thin beams may have large deflections resulting in non linear displacement, though still elastic. To apply elasticity to a tree branch, it is necessary to know its shape before any calculation. The theory shows that the curvature of the bent beam has a constant sign, there is no inflection point. Elasticity is adequate to calculate the deflection of a branch when the changes in thickness and length are negligible as for a branch temporarily loaded with snow or fruits.
Irreversibility of the bend
When the load is permanent, viscous or plastic deformation may occur, but the coupling between the simultaneous increase of the load and the thickening of the branch is more important. A branch bends continuously with time, even if it thickens. Even with a constant load there would be no decrease of the deflection when the branch thickens. On the contrary, a decrease of the thickness of the branch, would also increase the deflection. The process of growth being time-dependent, the shape of a branch is a function of time.
Method to calculate the branch bending
Dividing time into infinitely small intervals, it is possible to divide each time step in a few more steps: growth, loading and bending. Growth produces an increase in length and thickness of the branch and therefore a small change in geometry. Using the new geometry, the increase in load being small, linear elasticity may be used to compute the deflection, giving another change in geometry of the branch. The growth cycle may then be repeated for the next time step and so on. Because of the changes in thickness and length, the principle of superposition does not apply directly. Therefore, the stress distribution through the branch is no more linear and a permanent deformation occurs even if the incremental stress distribution is linear, e.g. elastic behaviour of the wood.
Numerical method
Many methods have been devised by engineers to analyse the mechanical behaviour of structures. Numerous configurations have been solved with analytical formulae, but they are valid only for simple geometries. With the advent of computers, numerical methods, such as finite elements and finite differences have been developed to calculate complicated structures. None of them (at our knowledge), takes into account the growth phenomenon (occurring in the same manner in the construction of buildings as for trees). In order to calculate the shape of a branch, a simple finite difference method has been used to integrate the differential equation of flexure step by step along the branch and iterated in a time marching process. At each time step, thickness, length and deflection of the branch are adjusted to take care of growth. Few input data are necessary: the thickness of the annual growth rings, the annual length increase of the branches, the growth angle, the density, the longitudinal elastic modulus of wood and the acceleration of gravity.
Numerical results
The results are synthesised in a picture of the whole tree: the young branches are at the top, almost straight and pointing in the direction of growth. The old branches, at the bottom of the tree, are curved and deflected towards the bottom. Figures 1 to 3 show the influence of the growth angle. On figure 3 b, after touching the ground, the branches grow horizontally, which is not realistic. On figure 4, the branches continue to grow with a support on the ground. Figure 6 shows a real tree where the branches have touched the ground.
The shapes of the branches are characterised by an inflection point, not predictable with the elastic criteria of MacMahon [6].
Fig 1: Small growth angle. Growth direction at 80°.
Fig. 2 : Growth direction at 45°.
Fig. 3 a : 90° growth angle without ground.
Fig. 3 b : 90° growth angle with ground.
Fig 4 : 45° growth angle.
Fig. 4 : Growth direction at 45°, growth speed twice slower than on the other figures. This shape is not geometrically similar to figure 2 (scale effect).
Fig. 5 : Old tree with branches growing up again after touching earth.
Fig. 6 a : Old cypress from the Bois de Boulogne in Paris
Fig 6 b : Another old branch having grown after touching the ground.
Fig 7 : An araucaria araucana looking like the numerical model.
Here is the original paper accompanied with the listing of the software : Branche