Estimation of projected surface area of irregularly shaped fragments Elvedin Kljuno*,Alan Catovic

Mechanical Engineering Faculty,University of Sarajevo,Vilsonovo setaliste 9,71000,Sarajevo,Bosnia and Herzegovina

Keyw ords:Projected area Fragments Trajectory

A B S T R A C TThe essence and the main contribution of this paper are consisted of the suggested novel method for estimation of a projected surface area of an irregularly shaped fragment,w hich represents a signi f i cant step tow ard a new method of an aerodynamic force estimation of a fragment motion through a resistive medium.The suggested method is to use a tri-axial ellipsoid that has a continuous surface(given as a mathematical function)to approximate an irregularly shaped fragment so that the fragment trajectory can be estimated faster taking into consideration that the aerodynamic force is proportional to a projected surface area of the fragment.

1.Introduction

When calculating a body trajectory,it is necessary to know the value of the reference area f i guring in the expression for the aerodynamic force.The most commonly used value for this area is the value of the projected surface area of the body(perpendicular to the velocity vector)at any instant.For axisymmetric bodies,this surface is taken as their cross-section area(circle).How ever,for an irregularly shaped body this surface is irregular and stochastic(not a continuous one),and it is impossible to determine it analytically.One of the methods for estimating the projected surface area of an irregular shape is to create 3D models of such bodies by 3D scanning them and then inserting the model in a CAD softw are to determine their projected surface areas in the required directions.This method is impractical since their orientation changes during the motion,and for the calculation of the trajectory,it is necessary to know the value of the projected surface area of a body at any given moment(for any possible orientation of the body).

Representative examples of the irregularly shaped bodies are fragments of HEw arheads,fragments resulting from the fracture of various structures due to the effects of severe storms,fragments created by improvised explosive devices,fragments resulting from the explosion of larger ammunition storages,meteorites that reach Earth's surface,etc.Generally speaking,it is very dif f i cult to estimate the projected surface area(to the plane perpendicular to the velocity vector)of such bodies,yet this parameter is needed for estimation of their trajectory since the aerodynamic force is proportional to projected surface area of a body.

To get an immediate insight into the idea how to estimate the projected area of the exposed surface of an irregularly shaped rigid body(fragment),the procedure of estimation is show n in Table 1,w hich show s the sequence of steps(a f l owchart).

Although the procedure w ill be described thoroughly throughout the paper,the table gives an overview of the process,starting from the input parameters(the main dimensions of the fragment and the angles of the arbitrary projection direction),calculating the boundary of the projection in the given direction and ending in the estimated projected area of the real fragment.

Table 1The f l ow chart of the projected area estimation process.

2.Literature review

Literature review yields a relative de f i cit of papers dealing w ith this topic.As far as fragments,generated by the detonation of high explosive w arheads,are concerned,in the available literature an approximate method of estimating the projected surface area of the fragment is often implicitly assumed to approximate the shape of a fragment w ith some know n geometrical body(e.g.a plate,cube,sphere,parallelepiped,or a cylinder)and determine the basic geometric parameters(volume and surface)of the body that approximates the fragment;and based on these data and the fragment mass data to determine the ratio of projected surface area of a body to the body mass(w ith the prediction of,so-called,equivalent body dimension).In this method,the data on the dimensions of the fragments from the experimental tests can also be used,w hereby,using the regression analysis,results are obtained that are better suited to the shapes of real fragments.Thus,according to the American standard for static testing of the spatial distribution of fragments[1],the ratio of the average projected surface area of the fragment and its mass is determined,w here correlation factor must be approximated using regression analysis for each fragment.This method does not take into account the generalized case w hen the velocity vector of the f l uid(air)is arbitrarily oriented,so the real projected surface area of the body is not determined at every moment of motion,w hich means that this method is not adequate for a more accurate calculation of the real fragment trajectory.

In the literature,there is also a statistical method of estimating the projected surface area of the HEprojectile fragments.Tw isdale[2]states that a probability theory(uniform random orientation)w as used in the calculation of the trajectory of the fragments to estimate the orientation of the tumbling fragments during the f l ight.Moxnes et al.[3]estimated the expected area of the fragment,considering the know n geometric bodies(plate,parallelepiped,cube).

The researchers from the TNO Institute,w ithin the MISDAC program[4],mention the non-dimensional factor of the fragment shape that is de f i ned as the ratio of the average exposed area of the fragment(according to Cauchy formula equal to 1/4 of the value of the body total surface)and the volume of the fragment represented by V2/3(this parameter w as f i rst introduced by McCleskey in 1988).They also attempted to estimate the effect of fragment tumbling using the assumption that the exposed area of fragment varies from a minimum to a maximum value.Using these methods,the real projected surface area of the fragment during the movement also is not taken into account,but rather its expected(statistical)value,assuming as in the f i rst method that the fragments are approximated by know n geometric shapes.

From the aforementioned,it is clear that there is a need to estimate the projected surface area of the fragment for each possible orientation of the fragment during the f l ight,in order to estimate the force,and ultimately the trajectory of the fragments w ith different shape,mass,and velocity.

3.The p hysical m odel

The developed physical model assumes that the fragments are approximated using the tri-axial ellipsoid.An ellipsoid has three pair-w ise perpendicular axes of symmetry w hich intersect at a center of symmetry,called the center of the ellipsoid.The line segments that are delimited on the axes of symmetry by the ellipsoid are called the principal axes.If the three axes have different lengths,the ellipsoid is said to be tri-axial,and the axes are uniquely de f i ned.

Fig.1 show s a tri-axial ellipsoid approximating the fragment(modeled in 3D Studio Max).Semi-axes of the ellipsoid,a,b,c are half the length of the principal axes.They correspond to the semimajor axis and semi-minor axis of an ellipse.Dimension a is the largest,and c the smallest.

Fig.2 gives a schematic view of the convex surface and plane(normal to the velocity vector)on w hich the elemental surface d A is projected.The unit vector of the velocity direction is denoted by e→vand the unit vector of the orthogonal direction onto the surface element is denoted by n→.

Fig.1.Approximation of the fragment with a tri-axial ellipsoid.

Fig.2.Projection of surface elements(schematic view).

The projection Apof the exposed part of the surface(Fig.2)can be presented as:

The condition on the angle means that only one(upper)part of the surface area is considered(Aexp).The upper exposed surface area hascosφ＞0,the back of the fragment has cosφ＜0 and on the dividing line cosφ=0,since the orthogonal direction on the surface is perpendicular to the velocity direction v→,i.e.n→⊥v→(n→·v→=0).

Although the enclosing the integral in expression(1)is possible(over the closed surface A of the fragment),w hich allow s the conversion to a volume integral using the Gauss-Ostrogradsky formula,it is not a convenient method since di v（e→v）=0 because e→v≠f（x，y，z）,i.e.the projections from both sides are w ith opposite signs and the integral over the closed surface cancels out.Since the upper and low er part have the same projected area Ap,but the projections are w ith opposite signs,then

Fig.3 show s ellipsoidal surface and projection of its element d A on the plane perpendicular to the velocity vector v→.Vector e→vis a unit vector in the line of the velocity vector(but in the opposite orientation),and the vector n→is a unit vector of the orthogonal direction to the surface d A.

Generally,the unit vector of the normal can be obtained as

The term grad→fin(3)is the norm of gradient of f,w here the gradient is in the direction of normal to the surface f.The gradient can be obtained based on the function f,

in the follow ing form

The unit vector of the velocity direction:

Fig.3.Projection of ellipsoidal surface element dA on the plane perpendicular to the velocity.

w here anglesαv，βvandγvare angles betw een the velocity vector and coordinate axes.Based on(7)and(8):Angleφ(Fig.3)is de f i ned asan angle betw een unit vector n→and velocity vector,so the projection of the surface element on the plane perpendicular to the velocity vector is de f i ned using angle cosine:

Based on(5):

w hereαn，βnandγnare angles betw een the unit vector of normal and coordinate axes.Based on(11):

Using(1),(2),(10),(11)and(12):

Fig.4 show s schematically the method for determination of the curve c⊥that separates exposed surface(upper part of an ellipsoid or a f i rst part of the incoming surface in the direction of vector velocity)from the rest of the ellipsoid.Velocity vector is perpendicular to curve c⊥in its every point,w hich means that unit vector n→is also perpendicular to the velocity vector.

Fig.4.Determination of the curve that separates exposed surface from the rest of the ellipsoid.

The idea here is to f i nd the surface limited by curve c⊥(w hich lies in a plane that is generally not perpendicular to the velocity vector)and then this surface should be projected on a plane(the planeπin Fig.4)perpendicular to the velocity vector.In this w ay,the required projection of the complete ellipsoid that the“view er”sees from the direction of the velocity vector is obtained.

Curve cπrepresents an intersection of the planeπand ellipsoid.It can be w ritten as:

In(14)n→πis the unit vector of normal on the planeπ,and r→is the radius vector in the planeπ.From(14):

Along w ith(15),the curve cπalso satis f i es the equation of the ellipsoid:

For curve c⊥,the follow ing applies:

Based on(17):

Fig.5 show s rotation of coordinate system xyz in order to obtain coordinate systemξηζw here one axis(ζ)is perpendicular to the plane w here the curve c⊥is,and other tw o axes(ξandη)belong to this plane.

For c⊥(Fig.5),similarly as in Fig.4,the follow ing applies:

Fig.5.Rotation of coordinate system xyz in order to obtain coordinate systemξηζ.

Based on(21),a vector of the normal on the plane w here curve c⊥is located can be de f i ned as:

Based on(22),unit vector of normal e→ζis obtained(in this case it represents the unit vector of axisζ)and this vector can be determined by dividing(22)w ith its magnitude:

For unit vector e→ζof coordinate axisζgenerally applies:

w here cosαζ==and cosγζ=

In order to f i nd rotation angleφ,unit vector eζcan be divided into tw o vectors e→ζxyand e→ζzas:

w here e→ζxyis the projection of a vector e→ζon xy plane,and vector e→ζzis a component of a unit vector e→ζin the direction of z axis.Components of a vector e→ζcan be obtained using projections of a vector e→ζxyon the axes.This way w e can set up relations from w hich w e can obtain angleφ:

w here eζ=e→ζ

=1,Kζis determined using(22),and intensity of vector e→ζxycan be obtained from(26):

Based on(27)and(28)w e get:

Angleγζbetw een axisζand z axis can be determined using the component of a unit vector e→ζin the direction of z axis.The cosine of this angle represents component of a unit vector e→ζand this component is also given w ith(26),so:After determination of rotation angles of a coordinate system ξηζin relation to old coordinate system xyz,it is possible to transform coordinates,i.e.express old coordinates using new coordinates:The goal is that in expressions w e have the new coordinate of a rotated coordinate system.After the rotation around z axis for angle

φ,w e get temporary coordinate systemξ′η′ζ′w hereζ′≡z,because w e have rotation around z axis.The relation betw een coordinates can be expressed in matrix form:

w here matrix represents rotation matrix for the case of rotation around z axis.Next rotation is rotation of temporary coordinate systemξ′η′ζ′around axisη′for angleθ=γζ,so appropriate connection betw een coordinates is:

The resulting transformation of coordinates after the tw o rotations is given as:

w here the matrix of the tw o transformations of coordinates from ξηζinto xyz is:

Based on(35)w e get:

For the curve c⊥expressed inside the coordinate systemξηζ expressions(36),(3)and(20)apply:

When w e substitute(38)into(20)w e obtain:

w hich can be expressed as:

In order to eliminate the part of(40)next to2ξη,w e need to rotate coordinate system(Fig.6).

For the coordinate system from Fig.6 follow ing applies:

By combining(41-46),w e get:

Fig.6.Rotation of the coordinate system in order to eliminate the part of(40)next to the term2ξη.

We can introduce follow ing substitutions in(48):

In order to eliminate the part of expression(48)next to x y,w e can w rite:

Area enclosed by the curve c⊥(w hich separates the exposed part of the surface from the back part of the fragmental surface)can be w ritten as:

w hich represents a projection of the fragment surface onto the plane containing the curve c⊥and the f i nal formula for projected area on the plane perpendicular to the velocity vector is:w here the angleφvζis the angle betw een the velocity vector direction and the direction perpendicular to the plane of the dividing line c⊥.

In expression(56),follow ing parameters are involved:

In(57,58),the angle is de f i ned by:

Parameters d,e and g are given by(44-46),and anglesφandθ are given by(29)and(30),respectively.

4.Validation of the m odel

The physical model is validated in tw o ways,using the CAD system and analytically.Four examples of the analytical calculation of projected area of an ellipsoid are given in Appendix A.

Regarding model validation using CAD tools,using the software(Ansys System)w e determined the projected area(perpendicular to the velocity vector)of an ellipsoid w ith follow ing dimensions of semi-axes:a=34 mm,b=8.65 mm and c=6 mm,for different ellipsoid orientations and for one full ellipsoid rotation(angle increment 15°,velocity vector was in the direction of y axis,rotation axis w as x axis,Fig.7).The results are compared w ith those obtained using the developed model.

Fig.8 show s the validation of the developed model w ith data obtained from the softw are(Ansys System allow s the determination of the projection of the body exposed surface normal to the coordinate axes)for the tri-axial ellipsoid.

As can be seen in Fig.8,the matching of the data is excellent because the data obtained using our model and using CAD technique practically coincides for all angles(one full ellipsoid rotation,a 15°increment),i.e.for all orientations of the body.

5.Application of a model to an irregularly shaped fragment and analysis of the results

The large percentage of fragments,formed after the explosion of a high-explosive w arhead,has an elongated shape(Fig.9)that can be approximated w ith a tri-axial ellipsoid.

It is assumed that the largest fragment dimension corresponds to the dimension a of the ellipsoid(the largest of the ellipsoid semiaxes),w hile the dimensions b and c are perpendicular to dimension a,w here c is the smallest dimension(the smallest of the tri-axial ellipsoid semi-axes),and the dimension b is larger than the dimension c and smaller than the dimension a.

Table 2 providesdata on the ratio of the dimensions a,b,and c of real fragments(generated after the detonation of projectile 130 mm HE M79).From Table 1 it can be seen that the fragments can be approximated w ith tri-axial ellipsoid because there are generally signi f i cant differences in dimensions a,b and c.

In order to determine the projected area of an irregularly shaped fragment,an arbitrary fragment(Fig.10)w as used(this fragment is similar to one generated by the detonation of the 130 mm HEM79 projectile(show n in Fig.9).To accurately digitize this fragment,Autodesk AUTOCAD softw are w as used.

The three-dimensional(CAD)fragment model(Fig.10)w as made by modeling the real fragment in three projections,and then by softw are manipulation:extruding in the direction of three coordinate axes,joining the extruded projections and de f i ning their intersection as the f i nal 3D fragment model.

Fig.7.CAD model of an ellipsoid used for the validation.

In addition to this procedure,a fragment can be draw n in a 3D graphics program(i.e.3DStudio Max,as show n in Fig.11),but then the 3D model of a fragment has more contours and it is generally more dif f i cult to w ork w ith(also,the modeling process is more complex and takes longer.The third option is to use a 3Dscanner to create a 3D body model and export a model to a CAD system.

Fig.8.Validation of a model with the data obtained using software for ellipsoid.

In order to use the fragment,show n in Fig.10,to calculate its projected area using a developed model,it is necessary to know its principal dimensions in three directions a,b and c,w here the largest dimension is a,and c is the smallest(b and c are perpendicular to a).These dimensions are obtained by measuring the dimensions of the fragment.

Fig.12 show s a fragment w ith an irregular shape(modeled in the CAD system)approximated by a tri-axial ellipsoid.The maximum dimension of the fragment in three directions,a,b and c,corresponds to the semi-axes of the tri-axial ellipsoid.Half of the maximum dimension represents semi-axis a,and tw o dimensions in the plane perpendicular to the semi-axis a and are divided into maximum and minimum.Half of this maximum dimension is b,and half of the minimum dimension is c.

Fig.13 shows the comparison(veri f i cation of results)of the values of fragment projected area(for fragment presented in Figs.10 and 12),obtained using Ansys System(for the 3D CAD model of the fragment),and using a developed model for analytical calculation of projected area.

The velocity vector in the f i rst case(Fig.13,above)w as in the direction of x axis(show n in Fig.12)and rotation of the fragment w as performed around z axis(w ith angular increments of 15°),w hile in the second case(in Fig.13)the velocity vector w as in the direction of axis z(show n in Fig.12)and the rotation of the fragment w as performed around y axis(also w ith angular increments of 15°).

Table 2Ratios of the dimensions a,b,and c of real fragments formed after the detonation of projectile 130 mm HEM79.

In Table 3 are show n the values of projected area of the fragment in both cases(vector velocity in the direction of x and z axis)for different position(orientation)of the fragment,and the relative difference betw een the results obtained by the software(CAD)and using a developed model in w hich the fragment isapproximated by an ellipsoid.

By analyzing the resultsfrom Table 3,w e see that in the f i rst case(velocity vector in the direction of x axis,Figs.12 and 13 above),relative differences are very small for the orientation of the fragment w here the fragment is exposed w ith a larger surface to the velocity vector.Somew hat larger relative differences are observed for angles w here the fragment is exposed w ith a smaller surface to the velocity vector(orientation of 90°and 270°relative to the initial orientation show n in Fig.12).

In the second case,the agreement betw een the results is excellent(velocity vector in the direction of z axis,Figs.12 and 13 below),and the relative differences are less than 6%for all orientation of the fragment.This suggests that the approximation of the fragment by ellipsoid is justi f i ed.

5.Conclusion

Fig.9.Real fragment(formed after the detonation of an artillery projectile 130 mm HEM79)presented in several projections.

The paper show ed a method to estimate the projected surface area of a rigid body(fragment)w ith irregular shape.During the f l ight,the fragment rotates and it can take any orientation w ith respect to the relative velocity vector of the f l uid(air).Since the aerodynamic force is proportional to the projected surface area in the direction of the air relative velocity vector,it is necessary to have a method to estimate the projected area in an arbitrary direction.The w ay this estimation w as done in the paper is by establishing an ellipsoidal surface around the irregular fragment.The ellipsoidal shape has three parameters,the three semi-axes,that can be adjusted such that the fragment f i lls-out the ellipsoidal body in the best w ay.The f i rst semi-axis is obtained as a half of the maximum dimension that can be measured on the irregular fragment.The second and the third semi-axes are obtained in the plane perpendicular to the f i rst semi-axis,as a half of the maximum and the minimum size,respectively,in the plane perpendicular to the f i rst semi-axis.

Fig.10.Digitized 3D model of a real fragment in different projections.

Fig.11.Rendered 3D models of fragments w e created in 3D Studio Max.

Fig.12.Irregularly shaped fragment approximated by a tri-axial ellipsoid.

Fig.13.The comparison of the projected area values for the fragment obtained using the CAD tools and developed model.

Table 3Values of a projected area of fragment in cases w hen vector velocity w as in the direction of x and z axis for different orientation of the fragment,and the relative difference betw een the results obtained by the softw are(CAD)and using a developed model.

The derivation show s that the exposed surface of the fragment model is distinguished from the rest of the fragment(the back part)by a planar curve,the ellipse.By projecting this area enclosed by the ellipse onto the plane w hich is perpendicular to the velocity vector,the projected area of the fragment model is calculated,w hich represents an estimation of the real fragment projected area in an arbitrary given direction.

The estimation was done using a continuous surface given as a mathematical function,w hich is divided into in f i nitesimal elements.The projected area is then obtained as an integral of projected surface elements.

The obtained results w ere tested for several projection directions using shape of a real fragment,w hich was scanned and digitized.The information about real fragment shape w as inserted into Ansys softw are and the projected area w as obtained for speci f i ed directions.The comparison betw een estimated projected area and the area obtained by the digitized model show ed relatively high agreement betw een data w ith relative error below 6%in one case.This analysis show ed that the w ay to estimate the projected area of a fragment is justi f i ed and can be used as a basis for an aerodynamic force estimation.

The contribution and advantage of this approach is re f l ected in the fact that it represents an important step tow ard a dynamic modeling that does not require a CFD result to estimate the aerodynamic force.Further,the fragment trajectory can be estimated in a simpli f i ed w ay,since the resistance force can be estimated at an arbitrary time w ith an arbitrary orientation of the fragment in the direction of motion(in the direction of the velocity vector of the fragment's center of gravity).

Therefore,the future w ork that is going to follow this paper is the aerodynamic(resistance)force estimation,as w ell as the aerodynamic moment estimation.These tw o estimations w ill provide a w ay to calculate the“source term”in the tw o vectordifferential equations:the law of center of mass motion and the law of the momentum change about the center of mass.

Another bene f i t of estimation of the projected surface area is the possibility to calculate the depth of the fragment penetration after the collision w ith an obstacle in a similar way asthe resistance force is calculated for the fragment motion through a resistive medium.The penetration depth estimation of a fragment w ill be one of the future w ork subjects,as w ell.

Appendix A.Examples of projected area calculation and comparison w ith numerical results

If w e w ant to determine,using the model developed in our paper,the projected surface of the ellipsoid in the simplest casew hen the velocity vector is in the direction of the individual axes(x or y or z,Fig.14)then this projected surface can be checked analytically as w ell,because in this case the projected surface is always equal to the product of the ellipsoid semi-axes andπ.a)Let semi-axes of the ellipsoid be:a=0.05 m,b=0.1 m and c=0.02 m.In the case w hen the velocity vector is in the direction of the axis x(Fig.14,left):

Fig.14.Schematic view of projected surfaces of ellipsoid w hen the velocity vector is in the direction of coordinate axis.

Then,according to our model,follow ing applies:

Using know n formula for an ellipse,this projected area has the same value(as obtained here using our model)because the obtained values a1and b1are actually an ellipse semi-axes,perpendicular to the velocity vector.

b)In the case w hen the velocity vector is in the direction of the axis y(Fig.14,center),then the velocity vector is:

so follow ing applies:

Using formula for an ellipse,projected area also has the same value as obtained using our model.

c)In the case w hen the velocity vector is in the direction of the axis z(Fig.14,right),then the velocity vector is:

v→=k→， (A26)

Then,it can be w ritten:

As in previous cases,using formula for an ellipse,projected area also has the same value as obtained using our model.

d)As an addition,w e w ill present general case for the same ellipsoid(semi-axes:a=0.05 m,b=0.1 m i c=0.02 m),but in this case the velocity vector w ill be:

so it can be w ritten(according to our model):

Finally,the semi-axes and the projected area are: