International Journal of Impact Engineering 109 (2017) 1
13

Contents lists available at ScienceDirect

International Journal of Impact Engineering

journal homepage: www.elsevier.com/locate/ijimpeng

The particle finite element method for the numerical simulation of bird
strike

geD4X X, D5X XR. BomanD6X X, D7X XL. PapeleuxD8X X, D9X XJ.P. PonthotD10X X
TagedPD1X XM.L. CerquagliaD2X X*, D3X XG. Delie
TagedPUniversity of Liege, Alle
e de la De
couverte 9, Liege 4000, Belgium

TAGEDPA R T I C L E I N F O TAGEDPA B S T R A C T

Article History: The Particle Finite Element Method (PFEM) is evaluated in the context of the numerical simulation of bird
Received 9 November 2016 strike events. To assess the possibilities of the method, theoretical analyses are initially performed based on
Revised 10 May 2017 the impact of a water jet on a rigid surface. Then, the influence of the geometry of a more realistic projectile
Accepted 12 May 2017
is analyzed and the capability of the method to take into account separation and fragmentation is
Available online 20 May 2017
highlighted. Finally, the method is tested for impacts against deformable targets, using a partitioned algo-
rithm with dynamic relaxation for the fluid-structure interaction, combining the PFEM for the description
TagedPKeywords:
of the bird with a non-linear Finite Element approach for the description of the impacted structure, which
Impact
Bird strike
can undergo large plastic deformations. To assess the quality of the obtained results a series of numerical
Numerical simulation examples from the literature has been selected and used as a reference throughout the paper. Among the
PFEM studies presented in this work also a novel numerical benchmark for the evaluation of bird impact simula-
Fluid-structure interaction tions is proposed.
© 2017 Elsevier Ltd. All rights reserved

1. Introduction (TagedP SPH) [2] is often employed for the simulation of bird strike [3,4].
SPH is a Lagrangian meshless method that can naturally account for
TagedPIn the last decades the advancements in numerical simulation extremely large deformations. Unfortunately this method is cursed
techniques together with the necessity to improve aircraft perfor- with consistency and stability issues [5,6] that can sometimes com-
mance without excessively increasing costs have made numerical promise the entire simulation. Moreover, it is usually more expen-
methods an essential tool both in design and certification phases of sive than Finite Element approaches (see e.g. [7]). The Particle Finite
aircraft components. When it comes to bird strike, the most difficult Element Method (PFEM), initially introduced in the field of civil engi-
task from a numerical standpoint is the modeling of the bird, which neering [8,9], is a Lagrangian particle method that has proven to be a
undergoes extremely large deformations. Moreover, as well known, powerful tool to solve complex free-surface fluid-structure interac-
when high velocity impact is concerned, the bird is often reduced, tion. The method can account for very large deformations, but pre-
even experimentally, to a surrogate projectile modeled as a weakly serving the robustness and generality of the Finite Element method,
compressible fluid (typically a mixture of water and air), as dis- and thus owning a key advantage over other approaches. The key
cussed by Wilbeck [1]. From a numerical standpoint, the presence of idea of the PFEM is the use of classical Finite Elements combined
a free surface and the strong interaction with the aircraft structure with a very fast remeshing procedure that allows the treatment of
represent a challenge for traditional computational fluid dynamics extremely large deformations, including separation and fragmenta-
methods based on an Eulerian grid. On the other hand, classical tion. In this work a preliminary assessment of this method in the
Lagrangian methods cannot cope with the extremely large deforma- framework of bird strike is proposed. First of all, a short description
tions experienced by the projectile during the impact, and in prac- of the PFEM is given. In the second part some academic tests are con-
tice some artifacts have to be introduced to control the element ducted for impacts on rigid surfaces. First, the case of a water jet
distortion, often leading to a loss of generality and robustness. More- impacting a rigid surface is considered and then the one of a more
over, since explicit time integration is usually employed, the pres- realistic bird impact. Results are compared to available theoretical
ence of distorted elements would lead to a drastic reduction in the and experimental data. The capabilities of the method to take into
time step size, degrading the computational efficiency. As an alter- account separation and fragmentation are also highlighted through
native to Finite Elements, the Smoothed Particle Hydrodynamics a bird impact over a rigid pseudo wing leading edge. Then, the
method is tested on more demanding cases, involving deformable
*
Corresponding author. targets and strong interaction between the bird and the impacted
E-mail address: marcolucio.cerquaglia@ulg.ac.be (M.L. Cerquaglia). structure. In this case, the structure is modeled using a distinct Finite

http://dx.doi.org/10.1016/j.ijimpeng.2017.05.014
0734-743X/© 2017 Elsevier Ltd. All rights reserved.
2 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13

Fig. 1. Schematics of the PFEM.

TagedPElement solver and the coupling is enforced using a partitioned algo- TagedP TagedP1. Define the particle connectivity through a Delaunay triangula-
rithm. For the structural part the software employed in this work is tion;
Metafor, an object-oriented Finite Element solver developed at the TagedP2. identify the domain boundaries using the a-shape technique;
ge.1 It is worth underlining that, to the best of the
University of Lie TagedP3. solve the governing equations making use of linear FE shape
authors’ knowledge, this is the first attempt to use PFEM in combina- functions;
tion with FEM in large plastic deformations in a partitioned way. A TagedP4. use the solution obtained from the previous step to update the
partitioned coupling of PFEM with FEM was already proposed in [10] particle positions.
but it included only elastic material behavior. The method is first val-
idated on a dam break against an elastic obstacle. Then, a pseudo
TagedPThis procedure is schematically described in Fig. 1 for better
bird impact on a clamped beam and a bird impact on a deformable
understanding. In the following, several technical aspects of the
metallic panel are analyzed. Whenever possible, the results are com-
implementation are omitted. For a complete description of the
pared to those available in the literature. It is also to be underlined
method, the reader is referred e.g. to Idelsohn et al. [8] or On eate
that no damage or failure of the structural parts (targets) are taken
et al. [9]. Finally, it is important to note that the present work focuses
into account in any of the examples presented in this work.
on two-dimensional cases.

2.2. Fluid-solid contact in the PFEM

2. The PFEM
TagedPClassically, PFEM employs no-slip conditions for the fluid at fluid-
solid boundaries, that is a fluid that comes into contact with a solid
2.1. PFEM general ideas
surface sticks to it, i.e. its velocity is the same as the solid at the fluid-
solid interface. To detect the contact the method takes advantage of
TagedPIn the PFEM the fluid is discretized using a set of points, hereafter
the remeshing procedure by creating new fluid elements every time
referred to as particles, which actually represent material points of
that a particle approaches a solid boundary, or flows over it, as illus-
the body. In order to evaluate the forces acting on each particle, a
trated in Fig. 2. Once the fluid elements are created, the particle is
new mesh is built at each time step from the entire set of particles
forced to move tangentially to the walls due to the incompressibility
using a Delaunay triangulation combined with the so-called a-shape
constraint: in this way it cannot penetrate the solid boundary. One
technique [11] for the correct identification of the external bound-
drawback of this technique is that some mass is added to the system
aries (see Fig. 1). This allows a fast evaluation of nodal connectivity
even when the total number of particles becomes very large [12].
Classical Finite Element (FE) shape functions can then be defined on
this new mesh to solve the corresponding weak form of the govern-
ing equations (see Section 2.3).
TagedPThe concept of a-shapes is the following: given a Delaunay trian-
gulation, all the triangles whose circumradius is larger than ah, a
being a scalar user-defined parameter and h being the average size
of the elements in the triangulation, are discarded. For most applica-
tions, a good choice of the parameter a is usually between 1.2 and
1.5 (see [13] for a thorough analysis of the role played by this
parameter in the PFEM).
TagedPGiven a set of particles, the main steps of the PFEM within one
time step can be thus summarized as follows:
Fig. 2. Fluid-solid contact detection in PFEM. Newly created fluid elements are in
light blue. (For interpretation of the references to colour in this figure legend, the
1
http://metafor.ltas.ulg.ac.be/dokuwiki/ reader is referred to the web version of this article.)
 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13 3

TagedPby the creation of the new fluid elements. However, the amount of TagedPwhere Nel is the total number of elements. The parameter tePSPG can
created mass is directly proportional to the elements size and can be be computed as
always reduced by refining the discretization, if needed.
t ePSPG ¼ h~
e
TagedPAlternatively, free-slip conditions could be employed, but the ð8Þ
imposition of such conditions is not trivial in the framework of the 2kU ~ k;
PFEM. Recently, the authors have proposed an efficient way to impose where h~ e is defined as the diameter of the circle which is area-equiv-
free-slip conditions in the context of fluid-structure interaction prob- alent to element e and U ~ is a global user-defined scaling velocity
lems solved via the PFEM [14], but for the moment these have not [16].
been applied to the impact problems analyzed in this paper.

2.4. Time integration

2.3. Governing equations and weak form
TagedPTime integration is performed using an unconditionally-stable
T t a given time t, and in the current configuration x, the incom-
agedPA implicit Euler scheme. Positions at time tnþ1 are computed as
pressible Navier
Stokes equations read:
  xnþ1 ¼ xn þ vnþ1 Dt; ð9Þ
Du @p @ @ui @uj
r i ¼
þ m þ þ rbi in VðtÞ; ð1Þ which leads to the system of fully discrete equations
Dt @xi @xj @xj @xi
2 38 9 8 9
1 nþ1 nþ1 Tnþ1 >
> > >
nþ1 > > nþ1 1 nþ1 n >>
@ui 6 Dt M þ K D <
7 v = < f þ M v =
¼ 0 in VðtÞ; ð2Þ 6 7 ¼ Dt : ð10Þ
@xi 4 1 5> > > 1 nþ1 n >
Cnþ1
Dnþ1 Lnþ1 >: pnþ1 >; > : hnþ1 þ C v > ;
where D(¢)/Dt is the Lagrangian, or material, derivative of the quan- Dt Dt
tity (¢ ), V(t) is the volume occupied by the continuum in the current where M is the mass matrix (that can be lumped if needed), K the
configuration, r(x, t) is the current density, m(x, t) is the dynamic matrix containing the viscous terms, D the discrete version of the
viscosity, u ¼ uðx; tÞ is the velocity vector, p(x, t) is the pressure, and divergence operator, and f a vector containing the contribution of
r b(x, t) is a body force, such as gravity for instance. body forces and surface tractions, i.e. the applied forces. C represents
TagedPThese equations have to be complemented with Dirichlet and a dynamic stabilization term and L is the discretized version of the
Neumann boundary conditions: Laplacian operator. A Picard algorithm is used to solve the non linear
uðx; tÞ ¼ uðx; tÞ 8 x 2 GD ðtÞ; ð3Þ system (10). For all the studies presented in this work one or two
non linear iterations were enough to reach convergence.
s ðx; tÞ ¢ n ¼ tðx; tÞ 8 x 2 GN ðtÞ; ð4Þ
where uðx; tÞ and t ðx; tÞ are imposed velocities and surface tractions 3. Theoretical preliminaries: water jet impacting on a rigid flat
respectively, s ¼ s ðx; tÞ is the Cauchy stress tensor, n denotes the wall
unit outward normal to the boundary, GD ðtÞ [ GN ðtÞ ¼ @VðtÞ; where
@V(t) is the boundary of the current volume V(t), and TagedPThe impact of a water jet on a rigid flat wall is analyzed in the
GD ðtÞ \ GN ðtÞ ¼ ⌀. first place. This problem has already been investigated in the past, in
TagedPConsider the following spaces for trial and test functions the stationary case, both from a theoretical and experimental view-
n   o point, for example by Schach [17] and Liu et al. [18], and represents
S ¼ u 2 H1 VðtÞ ju ¼ u on GD ðtÞ ; a natural starting point for our analysis. In this case however, since a
n  
S0 ¼ w 2 H1 VðtÞ jw ¼ 0 on GD ðtÞg; fully Lagrangian approach is employed, the impact of a finite jet,
n  o whose dimensions are given in Fig. 3(a), is considered. Some atten-
Q ¼ q 2 H1 VðtÞ ; tion should be drawn to the way the total impact time, T, and the
where H1(V(t)) is the space of square integrable functions with time step, Dt, are chosen. The total time T is taken of the order of the
square integrable first derivatives over the domain V(t). The weak time employed by the whole jet to collapse against the rigid wall,
form of Eqs. (1) and (2) can be established by multiplying Eq. (1) by a that is:
vector test function w 2 S0 and Eq. (2) by a scalar test function q 2 Q. H
T’ :
Integrating over the current volume and using Green’s theorem, the u0
weak form finally reads:find u 2 S £ [0, T] and p 2 Q £ [0, T] such that
Since no CFL (Courant
Friedrichs
Lewy) condition has to be veri-
Z Z  
R Dui @wi @ui @wi @uj @wi fied, the time step Dt is estimated only based on contact detection
r
V Dt wi d V ¼ p dij d V
m þ dV
V @xj V @xj @xj @xi @xj ð5Þ requirements. In this work it has been chosen as one third of the
R R
þ V r bi wi dV þ GN ti wi dGN 8 w 2 S0 ; time needed to cross one element moving at velocity u0, i.e.:
Z 1 h
@ui Dt ¼ ;
q dV ¼ 0 8 q 2 Q: ð6Þ 3 u0
V @xi
where h is the average element size. These general considerations
T o circumvent the Ladyzhenskaya
Babuska
Brezzi (LBB) condi-
agedPT are valid for all the other studies presented in this work.
tion [15] a Petrov
Galerkin pressure stabilization procedure has
been preferred over a fractional step method as it shows better mass 3.1. Impact at u0 ¼ 10 m=s
conservation properties, as demonstrated by Cremonesi et al. [10].
The stabilized version of Eq. (6) then reads TagedPSome preliminary analyses are conducted at an impact speed u0
Z
@ui ¼ 10 m/s. This value has been chosen arbitrarily. First, a mesh con-

q dV
V @xi vergence study is performed in order to choose the right discretiza-
X Nel Z     tion. The results are shown in Fig. 4(a) and (b), where the pressure
1 @q Du @p @ @ui @uj
þ tePSPG r i þ
m þ
rbi dVe ¼ 0 8 q 2 Q; evolution at the center of the wall is regarded as the quantity of
e¼1 Ve
r @xi Dt @xi @xj @xj @xi
interest. Note that in this analysis the same time step, corresponding
ð7Þ to the finest discretization i.e. Dt ¼ 20 ms, has been used for all the
4 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13

Fig. 3. Impact of a water jet on a rigid flat wall.

TagedPsimulations. A discretization corresponding to 10 particles along the tTagedP he only interesting quantity: the pressure profile along the wall is
half width, R/2, is chosen for further studies. The resulting mesh is also important. In his studies, Schach [17] proposes a theoretical and
shown in Fig. 3(b). However, the pressure at the impact point is not an experimental profile, both reported in Fig. 5(a) together with the

Fig. 4. Impact of a water jet on a rigid flat wall. Convergence analysis at u0 ¼ 10m=s: pressure evolution at the center of the wall. The number of particles refers to the half-
width, R. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13 5

Fig. 5. Impact of a water jet on a rigid flat wall. Analysis of the pressure profile along the wall and stagnation pressure. Discretization: 10 particles along the half-width, R.

TagedPanalytical solution proposed by Liu et al. [18]. These results are com-
pared to the profile obtained with the PFEM, evaluated at the middle
of the simulation, i.e. T/2, this instant being considered as the closest
to the case of a stationary jet. The agreement is fairly good, especially
if one considers that a two-dimensional finite jet is analyzed instead
of a three-dimensional axisymmetric stationary jet and that theoret-
ical results are also derived under some simplifying assumptions.

3.2. Impacts between u0 ¼ 10 m=s and u0 ¼ 400 m=s

TagedPIn this section different impact speeds which fall into the range of
interest for bird impact on aeronautical structures are considered. To
evaluate the quality of the results obtained, the stagnation pressure Fig. 6. Bird impact on a rigid flat wall. Bird models. The two models have the same
at the center of the wall is chosen as a representative quantity. For a mass. The model with flat bottom is discretized with 1115 particles. The model with
hemicylindrical bottom is discretized with 850 particles.
stationary jet, this pressure is equal to pth S ¼ 2 ru0 ; i.e. it grows with
1 2

the square of the impact velocity. From a numerical point of view,

since a steady state is never reached in this case, the stagnation pres-
sure is computed as the average pressure at the center: TagedPmm). Conversely, in the case of flat impact surface the length is com-
RT puted based on the one of the hemicylindrical case, so that the two
pðx ¼ 0; tÞdt projectiles have the same mass. The impact speed is also chosen to
pnum
S ¼ 0 :
T be representative of a real bird strike, i.e. u0 ¼ 117 m/s. The pressure
The results are shown in Fig. 5(b). The numerical solution is in very evolution at the center of the wall is shown in Fig. 7(a). The results
good agreement with the theoretical expectations, especially in the obtained with the two geometries are quite similar except for the
light of the approximation introduced by the analysis of a finite jet. amplitude of the initial peak, which is much higher in the case of a
Indeed, a maximum error of less than 10% has been obtained for an flat impact surface. This fact has already been discussed in the litera-
impact speed of 400 m/s. ture (see e.g. [3]) and could be expected somehow since the impact
is less abrupt in the case of an hemicylindrical impactor. In Fig. 7(b)
4. The PFEM for bird impact simulation the results obtained with a hemicylindrical impact surface are com-
pared to those derived experimentally by Wilbeck [1] for the same
4.1. Bird impact on rigid obstacles impact speed. The agreement is very good, except for the amplitude
of the initial peak. Nonetheless, some remarks should be made.
TagedP4.1.1. Bird impact on a rigid flat wall, u0 ¼ 117 m/s First of all, the correct amplitude of this peak is extremely delicate
TagedPThe method is now tested on the same problem, but for a more to evaluate and a value of 12 times the stagnation pressure is not
realistic geometry, closer to the one of a real bird. The fluid is again unphysical, as it can be found elsewhere in the literature both in
water, which is an acceptable approximation according to the results numerical and experimental results (see again [3]). Moreover, in
of Wilbeck [1]. Two geometries are considered, as shown in Fig. 6. As this case the analysis is limited to two dimensions and to an incom-
it can be observed, the two models differ for the shape of the impact- pressible fluid, while experimental data result, of course, from
ing surface, which is flat in one case (Fig. 6(a)) and hemicylindrical in three-dimensional tests, and from the use of slightly compressible
the other (Fig. 6(b)). The dimensions for the bird have been arbi- projectiles.
trarily chosen in the range proposed in [19], i.e. the half-width of the TagedPAs a further analysis, some oblique impact tests have been con-
body, R, has been fixed to 60 mm. For simplicity, a length of 3R has ducted, for incidence angles of 60°, 45°, and 30°. The time evolution
been chosen in the case of hemicylindrical impact surface (i.e. 180 of the pressure at the center of the wall for the different impact
6 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13

Fig. 7. Bird impact on a rigid flat wall. Pressure at the center of the wall. (For interpretation of the references to colour in this figure legend, the reader is referred to the web
version of this article.)

TagedPangles is shown in Fig. 8(a). As one could expect, as the impact tTagedP arget is still rigid, the only purpose of this study is to prove the
becomes more and more oblique, both the initial peak in the pres- capabilities of the PFEM to take into account extremely large defor-
sure and the total load transmitted to the structure decrease. In par- mations and eventually the complete fragmentation of the impact-
ticular, one would expect the stagnation pressure to be roughly ing body. Indeed, the bird is often cut during the impact, either by
proportional to the sine of the incidence angle. The value of the theo- the leading edge of a wing or by the fan blades in an engine, and this
retical stagnation pressure for the different angles are reported in has of course important consequences on the dynamics of the
dotted lines in Fig. 8(a): the numerical results follow the theoretical impact event. Being capable to take this fact into account is thus of
prediction with reasonable accuracy. For completeness, a compari- paramount importance. The geometry of the problem is described in
son with experimental results [1] for an impact angle of 45° is pro- Fig. 9(a). The front part of the leading edge has been constructed
posed in Fig. 8(b). As for the normal impact, the agreement is again using a B-spline whose control points are A, B and C. The two
very good. It is worth underlining that the negative pressure peak remaining segments (in green) are straight lines. The discretization
appearing in the simulation is a numerical artifact coming from the employed for the simulation is depicted in Fig. 9(b). The initial speed
fact that in an incompressible formulation the pressure is defined up is 100 m/s. The global evolution of the impact event is shown in
to a constant (i.e. a reference pressure). Here, this reference pressure Fig. 10. As it can be observed the method is perfectly capable of
is fixed to zero, meaning that depression regions correspond to neg- accounting for extremely large deformations up to the complete sep-
ative pressures. In this case, the depression is caused by the bird aration of the bird into halves. It has to be noticed that no conver-
flowing almost tangentially to the wall (due to its shape) at the very gence or stability problems arose during the simulation.
beginning of the impact.
4.2. Bird impact on deformable obstacles
TagedP4.1.2. Bird impact on a pseudo leading edge, u0 ¼ 100 m/s
TagedPIn this example, the bird model previously validated is tested for TagedPAfter the preliminary theoretical studies of Section 4.1 and the
the impact on a pseudo-leading edge of an aircraft wing. Since the assessment of the present method for impacts on rigid surfaces, the

Fig. 8. Oblique bird impacts on a rigid flat wall. Pressure at the center of the wall.
 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13 7

Fig. 9. Bird impact on a pseudo leading edge. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

TagedPapproach is evaluated for the case of impacts on deformable targets. sTagedP teps 1
5 are only performed once per time step. In the simplest
This is, of course, one of the necessary steps towards the application case of implicit coupling, on the other hand, the above procedure is
of the present method to real-life cases. Due to the use of a fluid repeated within one time step, until a given convergence criterion is
model for the bird, impacts on deformable surfaces fall into the fulfilled, which can be regarded as a fixed-point algorithm [22]. In
broader field of fluid-structure interaction problems. In this context, this work, an implicit Dirichlet
Neumann procedure is employed
two large families of algorithms can be identified: monolithic and and the convergence criterion is defined as
partitioned approaches.
1
TagedPMonolithic schemes employ the same formulation for both the pffiffiffiffiffiffiffi k rnþ1 k < ɛFSI ; nþ1
riþ1 ¼ d nþ1 nþ1
iþ1
d i ; ð11Þ
neq iþ1
solid and fluid portions of the domain so that no “external” coupling
procedure between the solid and the fluid part is needed: all equa- where n and i represent the time step and iteration numbers
tions are solved at the same time by a single code. This is a big disad- respectively, neq is the number of interface equations, and eFSI is
vantage of the approach since the most efficient codes have been the chosen tolerance, fixed at 10
6 for the examples appearing in
designed to solve either Navier
Stokes or solid mechanics equa- this work.
tions, which require different formulations and specific algorithms. TagedPThe main disadvantage of partitioned algorithms is that the cou-
By contrast, in partitioned approaches, the fluid and solid problems pling procedure may diverge, inducing instabilities or convergence
are solved separately in an iterative way, thus allowing dedicated issues, in the case of explicit or implicit couplings, respectively. This
software to be used in each case. is the case, for instance, when the fluid and solid densities are simi-
TagedPThe most employed partitioned algorithm is the Diri- lar, leading to well-documented problems linked to the so-called
chlet
Neumann procedure, which can be summarized as follows: ‘added mass’ effect (see e.g. [23,24]).
TagedPTo enhance (or even ensure) convergence, a relaxation procedure
is exploited in this work, which means that a weighted displacement
TagedP1. At a given time, tnþ1 ; get a suitable set of displacements, d, and
velocities, u, for the fluid-solid interface (usually coming from d nþ1 nþ1 nþ1 nþ1 nþ1
iþ1 ¼ viþ1 d iþ1 þ ð1
viþ1 Þd i ð12Þ
the previous solution of the solid equations);
is used to impose Dirichlet boundary conditions to the fluid, in step 2
TagedP2. Consider the solid-to-fluid interface as a Dirichlet boundary for
of the above procedure. That is, the predicted solid displacements
the fluid;
are “relaxed” before being applied to the fluid. Since, in the general
TagedP3. Solve the fluid equations using the predicted velocities as Dirich-
case, an optimal value of the relaxation parameter v cannot be
let boundary conditions;
defined a priori, Aitken’s relaxation is employed here, as discussed in
TagedP4. Consider the fluid-to-solid interface as a Neumann boundary for
[22,25]. The relaxation parameter is then computed automatically
the solid;
as:
TagedP5. Solve the solid equations by taking into account the forces com-  
ing from the fluid. rnþ1 T rnþ1
rnþ1
nþ1 nþ1 i iþ1 i
viþ1 ¼
vi   ð13Þ
 nþ1 nþ1 2
r
r 
TagedPThe coupling can be either explicit (see e.g. [20] and references  iþ1 i 
therein) or implicit (see e.g. [21]). When the coupling is explicit,
8 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13

TagedP erify the implementation of the fluid-structure interaction proce-

v
dure previously described. The geometry of the problem is depicted
in Fig. 11(a), where L ¼ 0:1456 m and w ¼ 0:012 m.
TagedPA column of fluid in a reservoir is initially retained by a rigid wall:
when the wall is (instantaneously) removed, the fluid starts flowing
inside the reservoir under the action of gravity and eventually hits
the flexible obstacle, which is fixed to the base of the reservoir. This
benchmark has been presented for the first time by Walhorn et al.
[27], in the context of Space-Time Finite Elements [28,29]. Rafiee
et al. [30] used it to validate a fluid-solid SPH formulation, while
Idelsohn et al. [31] and Marti et al. [32] to assess the quality of mono-
lithic fluid-solid interaction formulations in the framework of the
PFEM. The fluid is water, with density of 1000 kg/m3 and dynamic
viscosity of 0.001 Pa ¢ s. The obstacle is characterized by a density of
2500 kg/m3, a Young’s modulus of 106 Pa and a Poisson’s ratio equal
to zero. The acceleration of gravity is set to 9.81 m/s2. The chosen
discretization is described in Fig. 11(b): the average element size for
the fluid elements is 5 mm, while the panel is discretized with 30
elements along its length and 4 elements across the thickness. A
Chung
Hulbert time integration scheme [33] is employed for the
solid domain. For all the examples proposed in this paper, the
Chung
Hulbert dissipation parameters am and af have been chosen
as am ¼
0:97 and af ¼ 0:01.
TagedPSome significant instants of the simulation are reported in Fig. 12.
The results of Idelsohn et al. [31] are also included for comparison.
The results obtained with the present method are qualitatively very
similar to those of Idelsohn et al. both in terms of the flow evolution
and solid deformation. For a more quantitative comparison, the time
evolution of the horizontal displacement of the top-left corner of the
obstacle is reported in Fig. 13 and compared to others found in the
literature. Though there exist some dispersion between the different
authors, the results are generally in agreement. In particular, a simi-
larity with the results obtained by Idelsohn et al. [31] appears
clearly.

TagedP4.2.2. Bird impact on a clamped beam (LS-Dyna benchmark)

TagedPA 2D bird impact benchmark against a clamped beam coming from
the LS-Dyna software community2 is analyzed. In the benchmark an
ALE formalism is employed to take into account the large deforma-
tions of the bird. The problem geometry is depicted in Fig. 14(a),
where h ¼ 0:03 m, w ¼ 0:001 m, s ¼ 0:0034 m, x0 ¼ 0:0083 m, and
y0 ¼ 0:0243 m. The horizontal component of the initial velocity u0
Fig. 10. Bird impact on a pseudo leading edge. Results at different time instants. (For is 5 m/s while its (downward) vertical component is 1 m/s.
interpretation of the references to colour in this figure legend, the reader is referred TagedPThe beam has a density equal to 1 kg/m3, a Young’s modulus
to the web version of this article.) equal to 104 Pa and a Poisson’s ratio equal to zero. In this study, the
TagedPFor the first iteration, this formula can not be applied, and is replaced elastic beam is discretized using standard fully-integrated 4-nodes
by vnþ1 ¼ maxðvn ; v max Þ; where vn is the last computed value of v quadrangles, while in the LS-Dyna benchmark Hughes
Liu [34,35]
0
in the previous time step and v max is a well-suited user-defined shell elements are employed. The bird is modeled as a square of fluid
value. For the examples contained in this paper, an v max ¼ 0:5 is with a density of 781 kg/m3 and a dynamic viscosity of 0.1 Pa ¢ s. The
employed. Moreover, at each new time step, before entering the cou- fluid is incompressible in the case of PFEM and slightly compressible
pling iterations, an extrapolation is also performed on the computed in the case of the LS-Dyna ALE approach. The discretization is
solid displacements based on previous times, in order to further reported in Fig. 14(b): each side of the square representing the bird
accelerate convergence, as proposed in [26]. Depending on the num- is discretized with 10 elements both in LS-Dyna and in the present
ber of available previous time steps, the prediction can be computed approach. For the beam, 130 elements are employed along the
as: length and 3 elements across the thickness for the PFEM, while 100
shell elements are used in LS-Dyna. A Chung
Hulbert time integra-
for n ¼ 1; Order 1 : d~
nþ1
d~
nþ1
Order 0 : ¼ dn tion scheme [33] is employed for the solid domain in the case of
PFEM while a central difference explicit scheme is used in LS-Dyna.
¼ 2d n
d n
1 for n ¼ 2; Order 2 : d~
nþ1
TagedPSome time instants of the simulation obtained through the PFEM
5 n 1 are showed in Fig. 15. The horizontal displacements of the top (left,
¼ d
2d n
1 þ d n
2 for n3:
2 2 in the case of PFEM) corner of the beam obtained with the two meth-
ods are compared in Fig. 16. The results show very good agreement
TagedP4.2.1. Verification example: dam break against an elastic obstacle up to 0.45 s, when some differences appear. There can be multiple
TagedPBefore analyzing problems involving bird impacts on deformable
2
structures, a dam break against a flexible obstacle is employed to available at http://www.dynaexamples.com/ale/bird/bird-b.
 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13 9

Fig. 11. Dam break against an elastic obstacle. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

TagedPexplanations for these differences, including the different models TagedP ypothesis, all the main features of the phenomenon are taken into
h
used for the beam, but the authors believe that the main reason is account: high energy involved, strong interaction between the bird
the contact interaction between the bird and the beam: in LS-Dyna and the structure, large deformations of the bird, and large strains in
the bird seems to ricochet off the solid surface, while, in the case of the structure that undergoes both elastic and plastic deformations.
PFEM, the fluid sticks to it, since no-slip conditions are imposed at At the authors’ best knowledge, this is the first example in the litera-
the fluid-solid interface. ture of partitioned PFEM-FEM coupling where plasticity is taken into
TagedPIn this sense, the use of free-slip conditions at the fluid-solid account for the solid part and highlighting in this way the interest in
interface (as discussed in [14]) has been investigated. A snapshot of using and optimizing partitioned algorithms, whenever complex
the end of the simulation (t ¼ 0:6s) comparing no-slip and free-slip material laws have to be included in the model. The geometry of the
conditions is reported in Fig. 17. As it can be observed the behavior problem is inspired to the (3D) benchmark presented in [4] and
of the projectile is rather sensitive to the boundary conditions depicted in Fig. 19(a), where L ¼ 0:4 m. The panel thickness is 6.35
that are employed. Nonetheless, the corresponding effects on the mm. The discretization employed is described in Fig. 19(b): the bird
response of the beam are almost negligible, as it can be appreciated is modeled as in previous examples of Section 4.1 of this work, while
in Fig. 18. This is due to the fact that, even if the projectile flows over the panel is discretized with 400 elements along its length and 4 ele-
the beam in the case of free-slip conditions, it still cannot quit the ments across the thickness. To avoid hydrostatic locking in plasticity
surface of the target. Thus, its momentum is transferred to the beam 4-node quadrangular finite elements with selective reduced integra-
in approximately the same amount as when no-slip conditions are tion and constant pressure [36] are employed for the panel. The
employed. ends of the panel are fixed in the horizontal direction and the two
bottom nodes at the two extremities are fixed in both horizontal and
TagedP4.2.3. Bird impact on a metallic panel vertical directions. Typical physical properties of low-alloy steels
TagedPIn this last example an event representative of a real bird impact (see for example [37]) are chosen for the panel in this example. The
on a metallic panel is simulated. Except for the two-dimensional metallic panel has a density of 7800 kg/m3, a Young’s modulus of
10 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13

Fig. 14. Bird impact on a clamped beam. (For interpretation of the references to
colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 12. Dam break against an elastic obstacle. Results at different time instants. Left
column: present method. Right column: Idelsohn et al. [31]. (For interpretation of the
references to colour in this figure legend, the reader is referred to the web version of
this article.)

Fig. 13. Dam break against an elastic obstacle. Time evolution of the horizontal dis- Fig. 15. Bird impact on a clamped beam. Results at different time instants. (For inter-
placement of the obstacle upper-left corner. (For interpretation of the references to pretation of the references to colour in this figure legend, the reader is referred to the
colour in this figure legend, the reader is referred to the web version of this article.) web version of this article.)
 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13 11

Fig. 16. Bird impact on a clamped beam. Time evolution of the beam tip horizontal Fig. 18. Bird impact on a clamped beam. Time evolution of the beam tip horizontal
displacement. displacement using different boundary conditions. (For interpretation of the referen-
ces to colour in this figure legend, the reader is referred to the web version of this
article.)

TagedP210 GPa, and a Poisson’s ratio of 0.3. Moreover, a yield stress of sTagedP tage. At around 1:6  10
2 s the impact can be considered over: the
300 MPa and a linear isotropic hardening with plastic modulus equal bird is completely crushed onto the structure and a great part of its
to 1 GPa are employed for the description of the plastic deformation. kinetic energy has been transferred to the panel. These observations
Again, a Chung
Hulbert time integration scheme [33] is employed are confirmed by the evolution over time of the kinetic energies and
for the solid domain. For the sake of readability a detailed descrip- internal works of the bird and the panel, reported in Fig. 21. Concern-
tion of the solid formulation in the framework of large deformations ing the panel, a stress concentration can be observed close to the fix-
and plasticity is not given here. However, the interested reader can ations. The Von Mises stress there is higher than the yield stress,
refer to Metafor documentation (http://metafor.ltas.ulg.ac.be/doku thus leading to the presence of a plastic region. The stresses are
wiki/), or to [38] for further details. lower in the central part of the panel, which only undergoes very
TagedPThe global evolution of the impact is reported in Fig. 20. In the small plastic strains (see Fig. 20(a) and (b)). During this work the
first phase of the impact the phenomenon is driven by the bird, authors remarked the absence in the literature of 2D benchmarks
that pushes the panel which accelerate, gaining some kinetic for bird strike. Since this example could be employed for further
energy, and undergoes large deformations at the same time, comparison by other authors, the time evolution of the vertical
absorbing a large amount of the initial kinetic energy of the bird. displacement of the central point of the panel is reported in Fig. 22
Starting from about 4  10
3 s, the few kinetic energy left to the for completeness. The result is consistent with the aforementioned
panel is mainly linked to residual vibrations. The mutual interaction evolution of the phenomenon.
between the panel and the bird becomes more significant at this

Fig. 17. Bird impact on a clamped beam. Results at t = 0.6s for no-slip and free-slip Fig. 19. Bird impact on a flexible metallic panel. (For interpretation of the references
conditions. to colour in this figure legend, the reader is referred to the web version of this article.)
12 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13

Fig. 20. Bird impact on a flexible metallic panel. Results at different time instants. (For interpretation of the references to colour in this figure legend, the reader is referred to the
web version of this article.)

Fig. 22. Bird impact on a flexible metallic panel. Time evolution of the vertical dis-
Fig. 21. Bird impact on a flexible metallic panel. Energy balance. placement of the center of the panel.
 M.L. Cerquaglia et al. / International Journal of Impact Engineering 109 (2017) 1
13 13

5. Conclusions TagedP[14] Cerquaglia ML, Delige G, Boman R, Terrapon V, Ponthot J-P. Free-slip boundary
conditions for simulating free-surface incompressible flows through the particle
finite element method. Int J Numer Methods Eng 2017;110(10):921–46
TagedPA first assessment of the Particle Finite Element Method for the Nme.5439. doi: 10.1002/nme.5439.
simulation of bird strike events on aircraft structures has been per- TagedP[15] Brezzi F, Fortin M. Mixed and hybrid finite element method. Springer: Berlin;
formed. The method has been first validated through the analysis of 1991.
TagedP[16] Tezduyar TE, Mittal S, Ray SE, Shih R. Incompressible flow computations with
the impact of a water jet on a rigid flat plate, providing good results stabilized bilinear and linear equal-order-interpolation velocity-pressure ele-
in terms of prediction of the pressure at the impact point, the pres- ments. Comput Methods Appl Mech Eng 1992;95:221–42.
sure profile along the plate and the stagnation pressure for a typical TagedP[17] Schach. W. Umlenkung eines kreisfrmigen flflssigkeitsstrahles an einer ebenen
platte senkrecht ur stre
a{mungsrichtung. Ingenieur-Archiv 1935;6:51–9.
aircraft range of velocities. Then, a geometry closer to a real bird has
TagedP[18] Liu X, Gabour L, Lienhard V J. Stagnation-point heat transfer during impingement
been considered, confirming a good agreement with experimental of laminar liquid jets: analysis including surface tension. J Heat Transfer
results from the literature regarding the pressure evolution at the 1993;115(1):99–105.
TagedP[19] Bureau ATS. Research paper, the hazard posed to aircraft by birds. Tech Rep.
impact point for both normal and oblique impacts. It has also been
Department of Transport and Regional Services; 2002.
shown how the method can easily account for large deformations TagedP[20] Farhat C, van der Zee KG, Geuzaine P. Provably second-order time-accurate
and fragmentation without incurring any numerical stability or con- loosely-coupled solution algorithms for transient nonlinear computational aero-
vergence issue. Finally, the capabilities of the method to take into elasticity. Comput Methods Appl MechEng 2006;195(1718):1973–2001 Fluid-
tructure Interaction. doi: 10.1016/j.cma.2004.11.031.
account the interaction of the bird with deformable targets have TagedP[21] Wall WA, Genkinger S, Ramm E. A strong coupling partitioned approach for flu-
been shown: for this, a partitioned coupling algorithm that combines idstructure interaction with free surfaces. Comput Fluids 2007;36(1):169–83
PFEM and FEM in an original way has been employed. Of course, Challenges and Advances in Flow Simulation and Modeling. doi: 10.1016/j.
compfluid.2005.08.007.
some simplifications were made in this first attempt to model bird TagedP[22] Ku € ttler U, Wall WA. Fixed-point fluid
structure interaction solvers with
impact through the PFEM (namely the use of an incompressible fluid dynamic relaxation. Comput Mech 2008;43(1):61–72. doi: 10.1007/s00466-008-
and the restriction to two-dimensional tests), which make the 0255-5.
TagedP[23] Causin P, Gerbeau J, Nobile F. Added-mass effect in the design of partitioned
method still unsuited for the simulation of real-life bird strike algorithms for fluidstructure problems. Comput Methods Appl MechEng
events. Nonetheless, the preliminary results presented in this work 2005;194(4244):4506–27. doi: 10.1016/j.cma.2004.12.005.
are very promising. In the future, the method will be extended to TagedP[24] Fo € rster C, Wall WA, Ramm E. Artificial added mass instabilities in sequential
staggered coupling of nonlinear structures and incompressible viscous flows.
three dimensions, a compressible material will be employed, and a
Comput Methods Appl MechEng 2007;196(7):1278–93. doi: 10.1016/j.
rigorous analysis in terms of computational costs will be performed, cma.2006.09.002.
in order to assess the PFEM against the other available techniques TagedP[25] Degroote J, Souto-Iglesias A, Paepegem WV, Annerel S, Bruggeman P, Vierendeels
J. Partitioned simulation of the interaction between an elastic structure and free
(FEM, SPH ...).
surface flow. Comput Methods Appl MechEng 2010;199(3336):2085–98. doi:
10.1016/j.cma.2010.02.019.
Acknowledgments TagedP[26] Habchi C, Russeil S, Bougeard D, Harion J-L, Lemenand T, Ghanem A, et al. Parti-
tioned solver for strongly coupled fluidstructure interaction. Comput Fluids
2013;71:306–19. doi: 10.1016/j.compfluid.2012.11.004.
TagedPM.L. Cerquaglia would like to thank Eng. Michele Pittofrati from TagedP[27] Walhorn E, Klke A, Hbner B, Dinkler D. Fluid-structure coupling within a mono-
GDTech (Belgium) for the help received in performing the LS-Dyna lithic model involving free surface flows. Comput Struct 2005;83:2100–11.
simulations. TagedP[28] Tezduyar TE, Behr M, Mittal S. A new strategy for finite element computations
involving moving boundaries and interfaces - the deforming-spatial-domain/
space-time procedure: ii. computation of free-surface flows, two liquid flows,
References and flows with drifting cylinders. Comput Methods Appl Mech Eng
1992;94:353–71.
TagedP [1] Wilbeck J. Impact behavior of low strength projectiles. Technical Report 77-134. TagedP[29] Tezduyar TE, Sathe S, Keedy R, Stein K. Spacetime finite element techniques
Air Force Materials Laboratory, Wright-Patterson Air Force Base; 1978. Ohio for computation of fluidstructure interactions. Comput Methods Appl MechEng
45433. 2006;195(1718):2002–27 Fluid-tructure Interaction. doi: 10.1016/j.
TagedP [2] Monaghan JJ, Gingold RA. Smoothed particle hydrodynamics: theory and appli- cma.2004.09.014.
cation to non-spherical stars. Mon Not R Astron Soc 1977;181:375–89. TagedP[30] Rafiee A, Thiagarajan KP. An {SPH} projection method for simulating fluid-
TagedP [3] Letellier A. Contribution a  la mode
lisation des impacts d’oiseaux sur les aubes hypoelastic structure interaction. Comput Methods Appl MechEng 2009;198
des re
acteurs d’avions. (in French) Universite
d’Evry; 1996. (3336):2785–95. doi: 10.1016/j.cma.2009.04.001.
TagedP[31] Idelsohn SR, Marti J, Limache A, On eate E. Unified lagrangian formulation for elas-
TagedP [4] Smojver I, Ivancevi
c D. Coupled Euler Lagrangian approach using Abaqus/Explicit
in the bird strike aircraft damage analysis. 2010 SIMULIA customer conference; tic solids and incompressible fluids: application to fluidstructure interaction
2010. problems via the PFEM. Comput Methods Appl Mech Eng 2008;197:1762–76.
TagedP [5] Morris JP. A study of the stability properties of smooth particles hydrodynamics. TagedP[32] Marti J, Idelsohn SR, Limache A, Calvo N, D Elia J. A fully coupled particle method
for quasi-incompressible fluid-hypoelastic structure interactions. Meca
nica
Proc Astron SocAust 1996;13:97–102.
TagedP [6] Belytschko T, Guo Y, Liu WK, Xiao SP. A unified stability analysis of meshless par- Computacional 2006;25:809–27.
ticle methods. Int J Numer Methods Eng 2000;48:1359–400. TagedP[33] Chung J, Hulbert GM. A time integration algorithm for structural dynamics with
TagedP [7] Ryabov AA, Romanov VI, Kukanov SS, Shmotin YN, Chupin PV. Fan blade bird improved numerical dissipation: the generalized-a method. J Appl Mech
strike analysis using Lagrangian, SPH and ALE approaches. 6th European LS- 1993;60(2):371–5.
DYNA users conference; 2007. TagedP[34] Hughes TJ, Liu WK. Nonlinear finite element analysis of shells-part ii. two-
TagedP [8] Idelsohn SR, On eate E, Del Pin F. The particle finite element method: a powerful dimensional shells. Comput Methods Appl MechEng 1981;27(2):167–81. doi:
tool to solve incompressible flows with free-surfaces and breaking waves. Int J 10.1016/0045-7825(81)90148-1.
Numer Methods Eng 2004;61:964–89. TagedP[35] Hughes TJ, Liu WK. Nonlinear finite element analysis of shells: part i. three-
TagedP [9] On eate E, Idelsohn SR, Del Pin F, Aubry R. The particle finite element method. An dimensional shells. Comput Methods Appl MechEng 1981;26(3):331–62. doi:
overview. Int J Comput Methods 2004;1(2):267–307. 10.1016/0045-7825(81)90121-3.
TagedP[10] Cremonesi M, Frangi A, Perego U. A lagrangian finite element approach for the TagedP[36] Malkus DS, Hughes TJ. Mixed finite element methods reduced and selective inte-
analysis of fluid-structure interaction problems. Int J Numer Methods Eng gration techniques: a unification of concepts. Comput Methods Appl MechEng
2010;84(5):610–30. 1978;15(1):63–81. doi: 10.1016/0045-7825(78)90005-1.
TagedP[11] Edelsbrunner H, Mu € cke EP. Three-dimensional alpha shapes. ACM Trans TagedP[37] Department of Defense of the United States of America. MIL-HDBK-5J, Metallic
Graphics 1994;13(1):43–72. materials and elements for aerospace vehicle structures. Tech. Rep.. Department
TagedP[12] Idelsohn SR, On eate E. To mesh or not to mesh. That is the question.... Comput of Defense Handbook. Washington, DC.; 2003.
Methods Appl Mech Eng 2006;195:4681–96. TagedP[38] Ponthot J-P. Unified stress update algorithms for the numerical simulation of
TagedP[13] Franci A, Cremonesi M. On the effect of standard pfem remeshing on volume large deformation elasto-plastic and elasto-viscoplastic processes. Int J Plast
conservation in free-surface fluid flow problems. Comput Part Mech 2016: 1–13. 2002;18(1):91–126. doi: 10.1016/S0749-6419(00)00097-8.