Fulltext02 PDF

MASTER'S THESIS
Supersonic Artillery Projectile Fin

Deployment Simulation Methodology
Mia Forsblom
2013
Master of Science in Engineering Technology

Mechanical Engineering
Luleå University of Technology

Institutionen för teknikvetenskap och matematik
Acknowledgements
Seeing this master thesis as the end of my five years of university studies, there
are a few persons I would like to express my thanks to. First and foremost my
family, for the energy and support they give but even more for reminding me
that there are more important things in life than writing reports. I would also
like to thank Prof. Roland Larsson at LTU for encouraging me to specialise in
Engineering Mechanics. I have not regretted it a single minute.
I would of course like to thank my supervisor Prof. Pär Jonsén, always helpful
and never too busy to spare a few moments, for his friendly ways and for
sharing his knowledge. Not forgetting the rest of the Division of Mechanics of
Solid Materials at LTU where everyone was very friendly and helping me with
both practical problems as well as modeling ones.
Furthermore I very much appreciate the quick and thorough response form
David Karlsson and Daniel Hilding at Dynamore Nordic regarding my
questions about LS-DYNA, which were both quite a few and not always well
defined.
I am also grateful to Thomas Pettersson, Christer Thuman and Björn Edström at
BAE Systems Bofors who always took their time for good discussions
whenever I ran in to problems and wondered in which direction to focus my
work. Thanks also to Thomas for the simulations made in CFD++.
Last but not least, I would like to thank my supervisor at BAE Systems Bofors
Sven Strömberg. He has not only taught me a lot about artillery shells and
guided me throughout my work but first and foremost made me feel very
welcome in Karlskoga and is a true inspiration when I now start my career.
I would like to encourage any student thinking about where to do their master
thesis to put BAE Systems Bofors on top of their list.
Mia Forsblom
Karlskoga, June 2013
i
ii
Abstract
Sophisticated artillery munitions with guiding capabilities might be fin
stabilised and thus deployment time of the fins is then critical to achieve
aerodynamic stability. As the travelling velocity of the projectile during fin
deployment can be well in to the supersonic range, the wind loads are large
enough to affect the deployment.
When developing fin-stabilised projectiles numerical simulations can be a tool
in the design process, in order to reduce both development time and cost. The
dynamic event require a simulation method that can both capture the
aerodynamic flow around the projectile and the deployment of the fins, and
most important how they interact with each other.
The goal of this master thesis is to evaluate whether a Fluid-Structure
Interaction simulation using Arbitrary Lagrangian-Eulerian formulation in LS-
DYNA is suitable to simulate fin deployment on a projectile travelling at
supersonic speed.
Several simulations were made in LS-DYNA starting with a risk assessment
evaluating different important principles and modeling alternatives. The air
flow around a projectile with static fins was then compared with a traditional
CFD simulation. A fin deployment was finally simulated, both with and without
wind loads.
The main findings are firstly that LS-DYNA capture the important aspects of
the air flow around a projectile, however further validation is needed regarding
the accuracy of the pressure and velocity fields. Secondly, this type of
simulation requires capable computers, both regarding calculation capacity and
memory. Finally, the fin deployment times are in the range found from gun
firing tests of fin-stabilised artillery munitions.
iii
iv
List of abbreviations
ALE – Arbitrary Lagrangian-Eulerian
AoA – Angle of attack
CE/SE – Conservation Elements / Solution Elements
CFD – Computational Fluid Dynamics
FE – Finite Element
FSI – Fluid Structure Interaction
MMALE – Multi-material Arbitrary Lagrangian-Eulerian
v
vi
Table of contents
1 Introduction .................................... 1
1.1 Purpose of the project ............................................................................. 3
1.2 Background ............................................................................................ 3
1.3 Problem description ................................................................................ 4
1.4 LS-DYNA .............................................................................................. 6
2 Theory ............................................. 7
2.1 Inviscid compressible air flow ................................................................ 7
2.2 Isotropic linear elastic solids................................................................. 12
2.3 Finite Element analysis ......................................................................... 14
2.4 Reference frames .................................................................................. 16
2.5 Advection algorithms for ALE ............................................................. 18
2.6 Fluid-structure interaction modelling.................................................... 21
2.7 Artificial bulk viscosity ........................................................................ 23
2.8 Hourglass viscosity............................................................................... 23
3 Method .......................................... 25
3.1 Material modelling ............................................................................... 25
3.2 Risk assessment .................................................................................... 27
3.3 Static comparison with CFD ................................................................. 31
3.4 Deployment of one fin .......................................................................... 34
3.5 Full 3D model of fin deployment.......................................................... 38
3.6 Alternative model ................................................................................. 41
4 Results .......................................... 43
4.1 Risk assessment .................................................................................... 43
vii
5 Analysis and discussion ............. 65
5.1 Risk assessment .................................................................................... 65
5.6 Future work .......................................................................................... 71
6 Conclusions.................................. 73
7 Bibliography ................................ 75
Appendix; Meshing of the ALE domain 79
viii
1 Introduction
This master thesis is a part of the degree project in Mechanical Engineering

specialising in Engineering Mechanics at Luleå University of Technology
(LTU). It is carried out in collaboration with BAE Systems Bofors AB and the
Division of Mechanics of Solid Materials at LTU.
BAE Systems is a global defence, aerospace and security company with close to
100 000 employees worldwide. The Company delivers a full range of products
and services for air, land and naval forces, as well as advanced electronics,
security, information technology solutions and customer support services. BAE
Systems Bofors AB is a part of the Weapon Systems & Support business unit
developing and manufacturing artillery systems, direct fire naval guns, remote
weapon stations as well as both conventional and intelligent types of munitions
for the gun systems.
Traditionally, artillery munitions are spin-stabilised. When fired, rifles in the
barrel engage in a driving band on the projectile forcing it to rotate. Leaving the
muzzle it has an angular velocity, spin, in the range of several hundred Hz [1]
and the gyroscopic effect help to stabilise the projectile. Such munitions travel
in a ballistic trajectory determined partly by the elevation of the barrel and the
charge used.
More sophisticated artillery munitions, such as GPS guided, often require some
sort of fins to deploy at some position in the trajectory. The fins can be used to
generate lift, for guiding capabilities or range extension. A high spin is then
undesired and fins can be used also to increase aerodynamic stability.
For the fins used for aerodynamic stability, the traveling velocity of the
projectile during fin deployment can be almost Mach 3, giving high wind loads.
The angle of attach (AoA), i.e. the inclination of the projectile compared with
the travelling direction, leads to wind loads in more than one direction.
An example of a GPS guided munitions with both fins for guiding at the front
and for stability at the back can be seen in Figure 1.
1
Figure 1 An example of a GPS guided artillery shell, the M982 Excalibur.
http://news.cnet.com/2300-1008_3-6243124-4.html (credits U.S. Army)
Deployment time is an important parameter if fins are used for stability. If they
deploy too slowly, the wind loads on the projectile will be too high and the fins
will not be able to stabilise it. The projectile then tumbles out of control and
will not travel far. If they on the other hand deploy too fast, the projectile has
not completely left the muzzle, and the fins will be damaged by the muzzle
brake.
Another important parameter is the angular velocity of the fins when they reach
their fully deployed state. This, together with the wind loads acting on the fin,
can be used to evaluate their structural integrity.
If a fin deployment can accurately be simulated this can aid during development
of new types of artillery munitions as well as during modification of existing
ones. The amount of testing, which is costly both considering time and money,
could then be reduced.
2
1.1 Purpose of the project
The goal of this master thesis is to evaluate whether Fluid-Structure Interaction

(FSI) simulations using Arbitrary Lagrangian-Euler (ALE) formulation in the
Finite Element (FE) program LS-DYNA is suitable to simulate fin deployment
on a projectile travelling at supersonic speed. The most important variables to
determine are the fin deployment time followed by the velocity of the fins when
they reach their fully deployed state.
If it is, the goal is to develop a suitable methodology for such simulations.
If not, the goal is to analyse and describe the shortcomings of this method and
to see if there are other simulation methods in LS-DYNA that can be used.
1.2 Background
Previously, fin deployment simulations have been done at BAE Systems Bofors
using ANSYS. The wind loads on static fins where simulated using CFD and
the pressure on the fins transferred to loads in a fin deployment simulation
using implicit FEM. The results were acceptable, but the actual interaction
between the air and the fins during deployment is not captured.
The air flow around projectiles traveling at supersonic speeds, above Mach 1,
involves shock waves. Blast induced shock waves and their interaction with
steel plates has been successfully modeled using ALE and FSI formulations in
LS-DYNA.
Zakrisson use LS-DYNA for most of the simulations in his doctoral thesis
Numerical Simulations of blast loaded steel plates for improved vehicle
protection. The second paper focus on simulations of air blast loading
Zakrisson et al. [2] note that reasonable numerical results using reasonable
model sizes can be achieved from near-field explosions in air. A penalty based
coupling algorithm was there used for FSI and the van Leer advection algorithm
for advecting ALE material.
Chafi et al. [3] study blast-induced wave propagation and the interaction with
structures using penalty based FSI. They note that one can correctly predict the
relevant aspects of the blast–structure interaction problem, including: the
3
propagation of the blast wave in the medium and the structure response
subjected to the blast loading.
Štiavnický et al. [4] use ALE and FSI formulation in LS-DYNA to simulate a
gunshot from a rifle with caliber 5.56 mm. They used both constraint based and
penalty based FSI and noted that a good bullet flight simulation was achieved.
Van Loon et al. [5] compare various FSI methods for deformable bodies and
noted that the ALE method is easy to implement and an accurate method with
low computational cost. However, if re-meshing is required the method will be
both more time consuming and less accurate.
1.3 Problem description
A hypothetical artillery projectile of caliber 155 mm is analysed. It has fins

attached to its base which at gun launch are folded around the base and covered
by a hood. When the projectile leaves the barrel, the hood is pushed off by an
internal pressure and the folded fins deploy due to their elasticity. This can
schematically be seen in Figure 2.
Figure 2 Illustration of a projectile before and after fin deployment.
4
1.3.1 Delimitations
Some simplifications of the problem description have inevitably been made.
The focus is evaluating the capabilities of the FSI and ALE formulation in LS-
DYNA to simulate a fin deployment and develop a simple methodology for
such simulations. It is not to make a realistic, full scale model of a real world
scenario.
The models will be very simple and conceptual in the beginning and made more
complex as time goes.
Some general simplifications assumed are:
 There are no effects from the propellant gases
 The air is treated as an inviscid fluid and a perfect gas
 The fins are treated as isotropic and linear elastic
 The hood removal time is not calculated
 Simplified models of hypothetical projectiles and fins are used
 Simulation time will be until the fins are fully developed
The propellant gases flowing out of the muzzle brake influence the flow around
the projectile, however mainly before and in the early stages of the fin
deployment.
Supersonic air flows can often be approximated as inviscid [6]. The Reynolds
number is high and the important parameter is neither the drag forces acting on
the projectile nor temperature effects but the wind loads on the fin.
Air at moderate temperatures and pressures are often treated as an ideal gas. For
example, Zakrisson et al [7], Chafi et al. [3] and Slavik [8] all chose to model
air as an ideal gas during blast simulations, which also contains shock waves.
Štiavnický et al. [4] also model air as an ideal gas for gunshot simulations.
The hood removal time is generally well known and requires a different type of
simulation. A chosen hood removal time can easily be applied as a boundary
condition on the hood.
If a good simulation methodology can be developed for simple, hypothetical,
yet realistic geometries, the extension to real geometries should be
straightforward.
5
When the fins reach their fully deployed state, the deployment time is obviously
known, and the angular velocity can then be used in further simulations to
calculate the stresses and strains on the fins.
1.4 LS-DYNA
LS-DYNA is a general purpose FE program developed by Livermore Software

Technology Corporation [9]. The code's origins lie in highly nonlinear, transient
dynamic finite element analysis using explicit time integration. Some of LS-
DYNA’s analysis capabilities are
 Nonlinear dynamics
 Arbitrary Lagrangian-Eulerian formulation
 Fluid-Structure Interaction
The software includes a wide range of element formulations, material models
and contact algorithms.
6
2 Theory
Some relevant theory is briefly described here, starting with governing

equations for inviscid compressible air flow and isotropic linear elastic solids.
These are followed by an introduction to finite element analysis and a
discussion of three different reference frames. Two advection algorithms used
in LS-DYNA are then described together with two ways to model fluid-
structure interaction. The section ends with some words on artificial bulk
viscosity and hourglass viscosity.
2.1 Inviscid compressible air flow
The general governing equations for a fluid flow are the well-known non-linear
partial differential Navier-Stokes equations [10]. Analytic solutions are not
known but for a few simple cases, and different simplifications are often made,
depending on the characteristics of the flow, to get more applicable equations.
2.1.1 Characterization of air flow
Aerodynamic flows can be categorized in several different ways depending on
their properties and which aspects of the flow one consider important. If the
viscosity of the air is considered negligible, the compressibility considered
important and the flow velocity is higher than the speed of sound one is dealing
with inviscid, compressible supersonic air flow.
Inviscid
An inviscid fluid is considered to have zero viscosity, i.e. no resistance to shear
stress. As a result, inviscid flow exhibits neither friction nor thermal
conductivity or diffusion. No real world fluid is purely inviscid, but it can be a
good enough approximation if such effects are unimportant enough to neglect.
Theoretically, one approach an inviscid flow if the Reynolds number, ation 1
(1)
goes to infinity. Here and are the density and dynamic viscosity of the fluid
respectively, the flow velocity and a characteristic length.
7
The Reynolds number is thus a measure of the ratio of internal forces to viscous
forces. A projectile with a diameter m, air at atmospheric pressure
3
kg/m , Pas and a free stream velocity of m/s
in Equation (1) give a Reynolds number of . The inertial forces are
thus in the order of 10 million times higher than the viscous forces and the flow
is essentially inviscid [11].
Supersonic
If an air flow is considered supersonic is determined by the ratio of the flow
velocity and the speed of sound in the media . This is the definition of the
Mach number, 2
(2)
and the flow is considered subsonic if , sonic if , and finally
supersonic if In a supersonic air flow, disturbances created in some
point of the flow cannot travel upstream which will lead to shock waves around
objects [6].
A free stream velocity of m/s and a speed of sound in air of
m/s in the definition of the Mach number (2) gives M= thus a
supersonic flow.
Compressible
The compressibility of a fluid is the change in density of a fluid due to a
change in pressure and defined as 3
(3)
Generally, gases need to be considered compressible except for low speed flows
and a rule of thumb for air is that when the air is to be considered
compressible [6].
2.1.2 Governing equations for inviscid compressible flow
The governing equations for inviscid compressible flow are three conservation
equations, for mass, momentum and energy respectively, followed by an
equation of state and an equation for internal energy. A thorough derivation and
analysis of the governing equations can be found in an Aerodynamics or Fluid
Mechanics textbook, e.g. [6] or [10].
8
Conservation of mass
The continuity equation, which is an equation for conservation of mass, builds
on the principle that mass can neither be created nor destroyed. It is given by 4
(4)
and states that decrease of mass at a point equals the net flow out of that point.
Here is the material derivative and
Conservation of momentum
The equation of motion, which is an equation for conservation of momentum,
builds on Newton’s first law. It states that the net force acting on a body equals
the time rate of change of momentum of the body. Neglecting the viscosity
and assuming no body forces , such as gravity and electromagnetic forces,
gives 5
(5)
Conservation of energy
The equation for conservation of energy builds on the principle that energy can
neither be created nor destroyed, it can only change form. Assuming inviscid,
adiabatic flow without body forces gives 6
[ ( )] ( ) (6)
where is the specific internal energy.

Equation of state and internal energy
The three conservation equations contains five unknowns, pressure velocity
density , temperature and specific internal energy Two more equations
are thus needed, namely an equation of state and an equation for the internal
energy.
If the intermolecular forces between atoms in are neglected the pressure is
given by the ideal gas law 7
9
(7)
where is the specific gas constant. This relation can also be expressed in
terms of and the specific heat at constant pressure and volume
( ) . (8)
respectively according to 8
If the specific heats are considered constant, which is a reasonable assumption
for air at temperatures below 1000 K [6], the specific internal energy is given
by 9
. (9)
2.1.3 Shock waves
Shock waves occur at supersonic flows as disturbances introduced at some
point in the flow cannot travel upstream. The oncoming flow does thus not
know it will hit an object until it is too late and there is no time to direct the
motion around it. A shock wave is formed, as opposed to subsonic flows which
give smoothly varying streamlines around an object as illustrated in Figure 3.
a)
b)
Figure 3 Illustration of a) subsonic flow with smooth streamlines and b)

supersonic flow with a shock wave.
A shock wave is a very thin region of drastic changes in flow properties; the
pressure increase is almost discontinuous. A shock wave can be either normal
or oblique to the incoming flow. At the nose of a projectile travelling at
supersonic speed, which can be rather blunt, the shock wave is detached from
10
the body. A picture of a detached shock wave in front of a blunt body can be
seen in Figure 4. There is a then a small area at the center where the shock wave
is normal to the flow while the remaining part is oblique.
Figure 4 Detached shock wave in front of a blunt body.

http://www.hq.nasa.gov/office/pao/History/SP-440/ch6-2.htm
Often the change in properties across the shock wave is of interest, and across
the shock wave mass, momentum and energy are conserved. Denoting the state
before and after the shock 1 and 2 respectively the continuity, energy, and
momentum equations for a normal shock wave are 10
(10)
11
(11)
and 12
(12)
where is the enthalpy . Together with the equation of state (7) the
specific energy for an ideal gas (9) and the definition of the Mach number (2)
this leads to 13
(13)
11
14
( )
( ) (14)
and 15
( ) (15)
where the ratio between the specific heats . This means that the change in
properties across a normal shock wave can be described by the upstream Mach
number only. The pressure and density will increase across a normal shock
wave while velocity will decrease to subsonic levels.
Across an oblique shock wave, the normal components of the flow will behave
as across a normal shock wave while the tangential components are unchanged.
The flow downstream of an oblique shock wave can thus still be supersonic as
only the normal component of the flow velocity is reduced.
2.2 Isotropic linear elastic solids
Treating a body as linear elastic means that the following three assumptions are
made; [12]
1. The body is continuous and remains continuous under the action of
external forces. In other words, the atomic structure of the body is
disregarded.
2. The body obeys Hooke’s law. If the body is acted on by the forces
the resulting displacement is given by
where are constants independent on the magnitude of the

forces.
3. There is a unique unstressed state of the body to which it returns
whenever all the external forces are removed.
Furthermore, assuming an isotropic material means having the same properties
in all directions. These assumptions are reasonable for many steel structures
when subjected to moderate temperatures and stresses below the yield stress of
the material [13].
12
2.2.1 Equations and relations for isotropic linear elastic solids
The governing equations for linear elastic solids are the conservation equations
for mass, linear momentum and angular momentum respectively. These are
used together with a constitutive equation relating stresses and strains.
Thorough derivations and analysis of the equations can be found in a solid
mechanics textbook such as [12].
Conservation of mass
The continuity equation (4) is general and applicable to solids as well as fluids
and repeated here for clarity
Conservation of momentum
The equation of motion, which balance the linear momentum, for a linear
elastic solid written in terms of the Cauchy stress tensor is 16
(16)
where are the body forces per unit volume. The conservation of angular
momentum requires that the stress tensor is symmetric; .
Constitutive equation
The generalised Hooke’s law relates stresses and strains for a general elastic
solid and if thermal effects are ignore it is given by 17
(17)
where is the Elastic modulus tensor and the infinitesimal strain tensor.
The latter is in terms of displacements given by 18
( ). (18)
For an isotropic linear elastic material this reduces to the more convenient 19
(19)
where and are the so called Lamé constants. They can be described in terms
of the more well-known material parameter; the elastic modulus , shear
modulus and Poisson’s ratio according to ( )( )
and
13
2.3 Finite Element analysis
FE analysis is a numerical technique extensively used in engineering analysis

for obtaining approximate solutions to wide range of problems [14].
Using a continuum description, trying to find the solution of a field variable e.g.
pressure or displacement, gives a problem with an infinite number of
unknowns. The FE method will reduce the problem by dividing the continuum
in to a finite number of interconnected subregions, elements. Each element has
a set of nodes, and the nodal values of the field variable are instead a finite
number of unknowns to solve for.
Shape functions, also called interpolation functions, defined in terms of the
field variables at the nodes are then used to describe the variation within an
element. The nodal values and shape functions thus give a complete description
of the field variables.
Huebner et al. [15] describes the FE method as a seven step-procedure:

1. Discretize the continuum
I.e. divide it in to a finite number of elements. There are a wide range of
elements types, with varying complexity and applicability. A
comprehensive description is can be found in many FE textbooks e.g.
[14] or [15].
2. Select shape functions
The shape functions should preferably be continuous and have continuous
derivative across element boundaries. Polynomials are often used as they
fulfill these criteria and are easy to integrate [15].
3. Find the element properties
The relevant properties such as stiffness, mass and damping are described
in matrix form.
4. Assemble element properties to obtain the system equations
A matrix notation on the same form as for each element is used, simply
with more terms. They are assembled on the principle that two elements
sharing a node have the same value of the field variable at that node.
5. Impose boundary conditions
I.e. define known nodal values or nodal loads.
14
6. Solve the equations
If the problem is unsteady the nodal unknowns are functions of time and
differential equations needs to be solved using time integration.
7. Make additional computations
Often other than field variables solved for are the desired, such as stresses
and strains which can be calculated from the displacements.
FE analysis is often described as an art, [15] requiring experience and know-
how in order to achieve accurate and reliable models. One must bear in mind
that even though FE analysis is a useful tool for solving complex mathematical
models of physical problems, it is just only models that are solved. The real
physics of a problem can never be fully captured and even the best FE analysis
is never better than the chosen mathematical models.
2.3.1 Time integration for dynamic problems
In broad terms, the problems solved with FE-analysis can be classified in to
three categories; steady state, propagating or eigenvalue problems [14].
In a steady state problem, the response of the system does not change with time
and time is thus not included in the equilibrium equations. Propagating
problems, on the other hand, do depend on time and can be further classified in
to two subcategories. Firstly pseudo steady state, if the time effect on the
equilibrium relations is negligible but load functions are time dependent and
secondly actual propagating, or dynamic, if the equilibrium equations are time
dependent. The third category, eigenvalue problems, is where there is no unique
solution to the response of the system, but various possible solutions are sought
[14].
The equations of motion governing the linear dynamic response of a system of
finite elements can be written in matrix form as 20
̈ ̇ (20)
where and are the mass, damping and stiffness matrices respectively
[16]. ̈ ̇ and are the acceleration, velocity and displacement vectors and
finally is a vector with the external loads. This equation describes a system
that in static equilibrium at time and time integration is needed to find the
solution at a later time .
15
Central difference explicit time integration
The explicit time integration scheme used in LS-DYNA is the central difference
scheme [17]. The acceleration and velocity at time is then given by 21
̈ ( ) (21)
and 22
̇ ( ). (22)
The displacements are calculated using the equation of motion (20), thus
a known solution at time is used to calculate the solution a . An implicit
time integration scheme on the other hand use the equation of motion at .
Further information on the difference between implicit and explicit time
integration together with a number of different schemes can be found in a FE or
Numerical Analysis textbook such as [14] or [18].
Stability condition
A numerical solution is said to be unstable if errors grow exponentially for a
problem which have a bounded solution [18]. The central difference scheme is
conditionally stabile, requiring a time step of 23
(23)
to be stable [9]. is the highest eigenfrequency of the system.
2.4 Reference frames
Different kinematic descriptions, reference frames, can be used to describe the

kinematics of a continuum thus the relation between the deforming continuum
and the computational mesh. Two classical are the Lagrangian and the Eulerian
formulations and the ALE formulation has been developed trying to combine
the advantages of the two.
16
2.4.1 Lagrangian formulation
The Lagrangian formulation is the dominant in FE deformation analysis [3].
The computational mesh is then attached to the material and deforms with it so
that each Lagrangian element contain the same subset of material throughout
the simulation. An advantage with this reference frame is that the interface
between interacting parts are precisely defined and tracked. It can also handle
materials with history-dependent constitutive relations [19]. Large
deformations, material movements, can however lead to severe distortion of the
computational mesh. This in turn leads to both less accurate solutions and
decreasing time steps hence longer simulation times.
2.4.2 Eulerian formulation
The Eulerian formulation is widely used in fluid dynamics [19]. The
computational mesh is then fixed in space and the material moves with respect
to it, which means that large material movements can be handled. However, the
relative motion of the material and the computational mesh leads to numerical
difficulties due to advection of material. The Eulerian reference can be
computationally expensive as the computational domain needs to be large
enough to capture the relevant material movements at the same time as a small
element size is needed to achieve high accuracy. Furthermore, there are no
well-defined interfaces to couple with structures [20].
2.4.3 Arbitrary Lagrangian-Euler formulation
The ALE formulation makes use of a computational mesh which can be given
an arbitrary movement. It thus includes both the Lagrangian formulation, if the
mesh completely follows the material, and the Eulerian formulation, if the mesh
if fixed. As the computational mesh is allowed to move it is possible to have
less advection of material between elements than using pure Eulerian and less
mesh distortion than using pure Lagrangian [21].
ALE in LS-DYNA
LS-DYNA use the Operator split technique for ALE formulation which means
that first a pure Lagrangian step is taken and secondly an advection step.
The advection step includes;
1) Deciding which nodes to move
2) Moving boundary nodes
3) Moving interior nodes
4) Calculating the transport of the element centered variables
5) Calculating the moment transport and updating velocities.
17
LSTC notes that the Lagrangian step is usually less costly than the advection
step and during advection most time is spent on material transportation
calculations [17].
The alternative is a fully coupled formulation where material deformation and
convective effects are coupled in the same equations. Gadala notes that the
Operator split technique is a computational convenient technique, the main
advantages being less costly to implementation in current Lagrangian codes
and that the decoupling of the Lagrangian and Euler steps result in simpler
equations to be solved [22].
The Multi Material ALE (MMALE) formulation in LS-DYNA allows more
than one material to be defined in an element. This increase in flexibility
however comes to the cost of computational time. When more than one material
is allowed in each element additional data is needed to specify the materials in
each element and this must be updated by remap algorithms [17].
2.5 Advection algorithms for ALE
Using an ALE formulation, there is advection of material due to the relative

motion of the computational mesh and the material. The element centered
variables, density, internal energy, the stress tensor and the history variables,
are only known at the center of an element but needs to be approximated at the
element faces in order to determine the amount to be advected. Two alternative
advection algorithms implemented in LS-DYNA for advection of element
centered variables are the Donor Cell and the van Leer algorithms. Furthermore
there are node centered variables, i.e. three velocities and temperature (if
applicable) where the Half Index Shift (HIS) algorithm is used [17]. This is also
used for advection of momentum, which is product of the node centered
velocities and the element centered density [23].
According to Olovsson a good advection scheme should be monotonic,
conservative and as little dissipative and dispersive as possible [20]. An
advection scheme is monotonic if it does not introduce new minimum or
maximum values of the solution variable and conservative if the total properties
such as mass and momentum of the system are conserved. A dissipative
advection scheme smears out the variable fields and dispersive make high
18
spatial frequencies in the solution variable field travel slower than the mass
flow velocity.
Both the Donor Cell and van Leer are designed to conserve the total momentum
in the advection, on the expense of kinetic energy. This leads to a loss of kinetic
energy, thus total energy, with time [7].
2.5.1 Donor Cell algorithm
The Donor Cell algorithm, also called Upwind algorithm is monotonic,
conservative, stabile and also relatively simple. It is however only spatially first
order accurate and strongly dissipative. It is also dispersive but that is masked
out by the dissipation as high frequencies that travel slow are quickly damped
out [20].
The value of a variable at the face of an element is simply taken as the average
value in the element from which the material is flowing. The change of a
variable at the center of element, in between the nodes j and j+1, during a
time step is thus given by 24
( ) (24)
where
| |
( ) ( )
and is the velocity at node [17].
2.5.2 Van Leer algorithm

The van Leer algorithm is also monotonic and conservative but requires more
computational time compared with the Donor cell. It is spatially second order
accurate on a rectangular mesh but distorted elements introduce errors making it
only first order accurate [20].
The van Leer algorithm use a piecewise linear solution field variable instead
of assuming that is constant over an element. This gives 25
| |
( ) ( ) (25)
19
where
( )
( )
and an approximation of the slope. [17]
2.5.3 Critical time step regarding advection

The critical time step for the central time integration scheme is given by
Equation (23) and in LS-DYNA approximated as 26
[ ] (26)
where is the number of elements, a characteristic size of an element

and is the speed of sound. The latter is for a two dimensional element given
by the ratio of the Youngs modulus and the density
(7)
√
The advection imposes further restrictions on the time step to ensure that a
particle does not flow across more than half an element during a time step [20].
This gives 27
[ ] (28)
where is defined by the difference between the material velocity ̇ and

the mesh velocity
‖ ̇ ‖. (29)
20
2.6 Fluid-structure interaction modelling
When a fluid and a structure are in contact, there is continuous interaction

between them. As an example, the pressure of the fluid deforms the structure
which in turn changes the flow of the fluid, which will change the pressure
acting on the structure. In order capture the interaction between a fluid,
modeled in an ALE reference frame, and a structure, modeled in a Lagrangian
reference frame, a coupling method is needed. Two different ways to model this
FSI in LS-DYNA is by constraint based or penalty based coupling. Using either
method, the ALE and Lagrangian mesh should overlap but does not have to
coincide.
2.6.1 Constraint based FSI
The constraint based FSI is a method that conserves momentum but not energy.
Fluid nodes that comes in to contact with a Lagrangian structure are projected
on to the structure and the fluid mass, momentum and force are lumped to it.
The projected fluid nodes are then forced to follow the motion of the structure.
Consider a fluid with an initial velocity impacting a stationary solid, as
illustrated in Figure 5. If both have the same mass the momentum and
energy before the impact are 28
(30)
and 29
(31)
After the impact both have the same velocity, which is in order to conserve
the momentum . The energy is then reduced to
21
𝑢 𝑢
𝑢 𝑢
𝑢 𝑢
Before impact After impact

Figure 5 Illustration of a fluid impacting on a stationary solid using constraint
based FSI.
2.6.2 Penalty based FSI
The penalty coupling method used for FSI is similar to the one for standard
Lagrangian contacts in LS-DYNA. When a Lagrange element penetrates into an
ALE material the location of the coupling points at the interface are tracked.
Coupling forces are then applied to both elements, proportional to the relative
displacement between the two, as illustrated in Figure 6.
By not including damping in the contact formulation, the coupling force can be
defined as 30
(32)
where is the contact stiffness, the penetration distance and the unit
surface normal at the contact point. Alternatively, a coupling pressure can be
defined as a function of penetration distance. The penalty based coupling
conserves both momentum and energy, but is not as stable as the constraint
based.
ALE Coupling
point at Coupling forces
material
interface
ALE reference Lagrangian Penetration

mesh element distance
Figure 6 Illustration of penalty based FSI.
22
2.7 Artificial bulk viscosity
In order to handle the almost discontinuous changes in flow properties in a

shock wave many FE-programs make use of an artificial viscosity. The method
was first introduced by Von Neumann and Richtmeyer and will smear out the
shock front [24].
In LS-DYNA the artificial bulk viscosity is given by 31
( ) (33)
where is the density, is a measurement of the element size, is
the trace of the strain rate tensor, is the speed of sound in the media and
are dimensionless constants. The artificial bulk viscosity damps out the highest
frequencies in the solution and the ripple around the pressure front can be
eliminated [9].
2.8 Hourglass viscosity
When having under integrated elements, there is a need to control so called

hourglass modes having zero energy. These are unphysical modes of
deformation producing neither stress nor strain [25]. Two ways to avoid these
are by adding viscous damping or a small elastic stiffness and several different
algorithms are implemented in LS-DYNA [17].
Type 8 – full projection warping stiffness.
The type 16 element in LS-DYNA is a fully integrated shell element which
lightens problems with locking and enhances in-plane bending behaviour [17].
As these elements are fully integrated they do not exhibit traditional
hourglassing, but instead warping of the elements might degrade the solution.
Hourglass type 8 can then be used in LS-DYNA to avoid this by including a
warping stiffness [9].
23
24
3 Method
Starting with very simple models, the complexity was increased towards a full
3D model of a fin deployment with find loads. The work was divided into the
following steps;
1. Firstly a risk assessment was performed trying to identify possible
hazards and predict if LS-DYNA is a suitable tool for the desired
simulations.
2. Secondly a comparison was made between LS-DYNA and the CFD
software CFD++ on a half-symmetry model of a projectile with two rigid
fins. Two models
a. Without AoA – visual comparison of velocity and pressure fields
b. With 2° AoA – comparison of forces acting on the fin
3. At the same time a model of a fin deployment, excluding wind loads, was
developed. This was in turn done in two steps
a. Fold the fin
b. Deploy the fin
so that one folding only is needed and can be used for several
deployments.
4. Finally a full 3D model of a fin deployment on a projectile with eight fins
traveling at 950m/s was developed including spin of the projectile.
The latest solver available, R7.0.0 was used with double precision for all
simulations but the full 3D model. The earlier version R6.1.1 with double
precision was used for the full 3D-model due to errors in the ALE mapping
using R7.0.0.
3.1 Material modelling
Four material models have been used throughout the simulations; an elastic and
a rigid to model the Lagrangian parts together with two null material models for
the ALE air.
25
3.1.1 Lagrangian parts
Throughout the work, the structural parts were either modeled with a rigid or an
elastic material model with the properties given in Table 1. The elastic
parameters of the rigid material are used in the contact definition between
different parts.
Table 1 Material properties used in rigid and elastic material models.
Material Density Elastic modulus Poisson’s
model [kg/m3] [GPa] ratio
Rigid 7850 210 0.30
Elastic 7850 193 0.25
3.1.2 ALE parts

The air was modeled with a null material model together with an Ideal gas
equation of state with the properties given in Table 2 and Table 3. A dummy
vacuum material was also defined to avoid warnings using the MMALE
formulation.
Table 2 Material properties for the null and vacuum material models.
Dynamic
Density
Material model viscosity
[kg/m3]
[Pas]
Null 1.2 0
Vacuum 1.2·10-3 0
Table 3 Material properties for the ideal gas equation of state.

Initial Initial
Equation of Cv Cp
temperature relative
state [J/kgK] [J/kgK]
[K] volume
Ideal gas 717 1004 293 1
Throughout, artificial bulk viscosity was used with standard parameters, i.e.
and together with standard hourglass viscosity with the
recommended coefficient [9] .
26
3.2 Risk assessment
A risk assessment was made in order to identify possible problems at an early

stage and evaluate their effect on the simulation methodology development.
Four basic questions were identified which capture the essence of the desired
simulation;
1. Can the FSI between a fin and air be captured?
2. Can velocities of Mach 3 be applied?
3. Can the elastic springback of a fin be captured?
4. Can all of the above be simulated at the same time?
To answer them, simulations on greatly simplified geometries, not focusing on
mesh quality, were performed.
3.2.1 FSI and applying Mach 3 velocities
A solid fin, 10 mm x 10 mm x 1 mm, surrounded by a 200mm x 100mm x 100
mm volume of air was used to evaluate the risk assessment questions 1 and 2
and 4. The domain can be seen below in Figure 7 with the fin in black, and the
ambient section to the left in a). Four different models were used, which are
summarised in Table 4.
Table 4 Summary of the four different models used to model the fluid structure
interaction between a fin and air.
Advection
Model nr Fin elements Wind loads
scheme
1 Solid Air flow Donor Cell

2 Shell Air flow Donor Cell
3 Shell Air flow van Leer
4 Shell Moving fin Donor Cell
27
Air
flow
a) b)
Figure 7 The fin and the air domain used in the risk analysis, viewed a) from
the side and b) from the inlet.
Modelling the air
The air was modelled as an ALE domain consisting of 150 000 hexagonal solid
MMALE elements with an element length of 2 mm. The lower boundary of the
air domain was defined as a slip wall while the remaining was left open. A
reference pressure of 100.9 kPa was applied to all open boundaries. Both the
Donor Cell and the van Leer advection algorithms were used.
Modelling the fin
The fin was modelled as a Lagrangian structure either using solid hexagonal
elements or shell elements, both types having an element length of 1 mm.
Either -constant stress solid element or Belytschko-Tsay shell element was
used, the later having 1mm reference thickness.
Both the rigid material model, for evaluating the flow field variables, and the
elastic, for evaluating stresses in the fin, were used. The fin was constrained in
all directions for model 1, 2 and 3 but free to move in the x-direction for the 4th.
Applying the wind loads
The wind loads on the fin was applied in two ways, either by letting air flow
over the fin or by moving the fin through the air. The difference in velocity
between the fin and the air was in either case 1000 m/s, about Mach 2.9.
Air flowing over the fin
Air flowing over a stationary fin was modelled by defining an ambient inlet
section, which can be seen as the grey section to the left in Figure 7. This
section keeps the thermodynamic state defined by equation of state throughout
the simulation. The left end nodes were then given a prescribed motion in the x-
direction.
28
Fin moving through the air
When moving the fin through the air, the fin was given a prescribed motion in
the negative x-direction. The ALE mesh was then defined to follow the
movement of the fin and no ambient section was defined.
FSI modelling
The interaction between the fin and the air was modelled using penalty coupling
with a penalty factor of 0.1. The coupling was applied in the normal direction in
both compression and tension. A coupling grid of 4x4 nodes per Lagrangian
element was used as the element size of the Lagrangian and Eulerian mesh are
equal. When using shell elements, both the edge of the shell and the reference
thickness were considered in the contact definition.
3.2.2 Elastic springback of a fin
A 50x25mm elastic fin rotating about a rigid end was used to answer question
3. The fin with the elastic elements in black and the rigid end to the left in grey
can be seen in Figure 8.
Figure 8 Showing fin for elastic springback with elastic elements in black and
the rigid end to the left in grey.
Modelling the fin
The fin was modelled as a Lagrangian structure using 5 mm x 5 mm
Belytschko-Tsay shell element with 1 mm reference thickness. The rigid
material model was used for the end which the fin rotates about while the
elastic material model was used for the remaining part.
29
Applying loads and constraints
The rigid end of the fin was constrained from translating in any direction and to
only rotate about the y-axis. The rotation about the y-axis was furthermore
constrained to be between -2 and 0 rad, where 0 rad is when the fin lies in the
yz-plane. The rigid end will then rotate about an axis parallel to the y-axis going
through the center of mass of the rigid end, and the elastic part will follow it.
The nodes at the right end of the fin were loaded in the x-direction as seen in
Figure 9. The load was increased from 0 to 100 N during 1 ms, hold at 100 N
during 0.1 ms and then removed during 0.1 ms.
Axis of rotation 100 N load

X
Figure 9 Showing a loaded fin together with the axis of rotation.
3.2.3 Doing all at once

In order to evaluate the fourth risk assessment question model nr 3 in Table 2
was used. The fin was then translated with a constant velocity of 10 m/s in the
-y-direction, perpendicular to the air flow, as shown in Figure 10.
Wind
direction
Fin
movement
Figure 10 Showing the air domain from above and the translation direction of
the fin.
30
3.3 Static comparison with CFD
In order to evaluate the possibilities of simulating the air flow in LS-DYNA,

comparisons were made with a CFD simulation made in CFD++. A half
symmetry model of a projectile with two static, fully deployed fins was used
and wind loads of Mach 2.7. The model, cut through the middle, can be seen in
Figure 11. Firstly the pressure, density and velocity fields at the symmetry
plane were compared for 0° AoA and secondly the net force in the y-direction
acting on the fin for 2° AoA.
Figure 11 Outline of the projectile, a fin and the surrounding air domain.
3.3.1 Modelling the projectile
A simple model of a hypothetical projectile was used having a total length of
970 mm. The diameter at the back end is 155 mm and at the front end 30 mm.
The projectile is defined by the outline of the ALE domain and is thus modelled
as a boundary condition on the ALE elements.
3.3.2 Modelling the fin
A 100 mm high, 140 mm long and 1 mm wide fin was modelled in the zx-plane
at the back end of the projectile. Belytschko-Tsay shell elements were used with
an element size of 0.5 mm x 0.5 mm and 1mm reference thickness. The wing
was modelled as rigid and constrained in all directions.
3.3.3 Modelling the air
An air volume extending 75mm in front of projectile and 750mm behind, i.e.
having a total length of 1795 mm was modelled using hexagonal MMALE
elements. The area in the yz-plane is increasing from 300 mm x 150 mm at the
31
front to 1500 mm x 750 mm at the back. The air volume including the outline
of the projectile and the fin, viewed from above, can be seen in Figure 12.
Two models have been used, one finer with ~1 320 000 element and one
coarser with ~550 000 elements. The mesh close to the fin is 1 mm x 1 mm x 1
mm or 2 mm x 2 mm x 2 mm for the fine and coarse model respectively. Both
van Leer and Donor Cell advection methods were used for the coarser mesh
while only the Donor Cell advection was used for the finer. Further details
about the meshing of the ALE domain can be found in Appendix II.
Figure 12 Showing the air domain form above, including the outline of the
projectile and the fin.
Boundary conditions
A free slip boundary condition was applied to the bottom area i.e. the symmetry
plane including the outline of the projectile. All other ALE boundaries where
defined as open with a reference pressure of 100.9 kPa.
3.3.4 Applying wind loads
The front, top and sides of the air domain were defined as ambient thus held at
constant properties throughout the simulations. The ambient section can be seen
as the dark region in Figure 13. The wind loads were applied by defining an
initial velocity at the nodes of the ambient section. An inflow velocity of Mach
2.7 was used with either 0° or 2° AoA giving the x- and y-component of the
velocities in Table 5.
32
To reduce the simulation time before a stabile flow pattern is reached, most of
the remaining ALE elements were given an initial velocity equal to the inflow
velocity. The elements next to the fin and in the wake were not given any initial
velocity. Figure 14 shows a cut in the xz-plane in the middle of the domain
were elements in dark have been given an initial velocity.
Table 5 Showing inflow and initial velocities for different angles of attack, used
for comparison with CFD++.
Angle of attack x-velocity [m/s] y-velocity [m/s]
0° -926.6 0
2° -926.1 32.3
Figure 13 Showing the air domain surrounding the projectile. The darker
elements define the ambient section.
Figure 14 Showing initial velocities at an xz-plane through the middle of the

domain. The dark region has initial velocity Mach 2.7, the bright region has
zero initial velocity.
33
3.3.5 FSI modelling
Penalty based FSI was used to couple the air and the fin. The coupling was
applied in the normal direction, in both compression and tension, and both the
edge of the shell and the reference thickness were considered in the contact
definition. A coupling grid with 2x2 nodes per Lagrangian element was used as
the Lagrangian element length is half of the ALE element length next to the fin.
The coupling pressure was given as a linear function of penetration distance so
that zero penetration give zero coupling pressure and 0.5 mm penetration give
coupling pressure of 4 MPa. The maximum pressure next to the fin is
approximately 0.2 MPa thus giving 25 µm penetration.
3.3.6 Simulation in CFD++
The simulation in CFD++ was performed by Thomas Pettersson at BAE
Systems Bofors. The finer mesh was used with a plane of elements removed to
define the outline of the fin. All boundaries but the symmetry plane and the
projectile including the fin were given a velocity far field boundary condition
where static temperature and static velocity were defined. The symmetry plane
was given a symmetry boundary condition and the projectile and fin was
defined as a free slip wall.
The air was modelled as inviscid and the simulation was initialized with zero
velocity in the wake. The convergence criteria used was a reduction of the
residuals by three orders of magnitude and that the forces and moments acting
on the projectile and fin had been stabilised.
3.4 Deployment of one fin
A fin deployment simulation of one fin was performed in two steps, firstly the
fin was folded under a hood and secondly the hood removed and the fin
deployed due to its elasticity. This was done in separate simulations and the
stress state in the fin from the folding simulation was used as an initial
condition in the deployment simulation. Only one fin folding is thus needed and
can be used for several deployment simulations.
When folding the fin, it was constrained between two rigid half-cylinders to
imitate the location of the fin in between the base of the projectile and the hood
during gun launch, which can be seen in Figure 15. The fin and its stress state
34
after folding, together with the rigid parts, were exported using
*INTERFACE_SPRINGBACK_LSDYNA. The deployment of the fin was
firstly done without any surrounding air and secondly within the ALE domain
used for CFD comparison.
Folding
direction
Fin
Hood
Inner
cylinder
Figure 15 Showing fin folding where the hood is slid over the fin and the inner
cylinder.
3.4.1 Modelling the fin

A fin approximately 110 mm long and 70 mm high was modelled using
2.5 mm x 2.5 mm fully integrated shell elements with a reference thickness
slightly larger than 1 mm. Lobatto integration was used with three through
thickness integration points thus giving one integration point at each shell
surface and one in the middle. Hourglass type 8 was used to activate full
projection warping stiffness using standard settings.
Most of the wing is modelled as elastic, however the last 2.5 mm at the right
end, as seen in Figure 16, is however rigid. The rigid end is free to rotate about
the x-axis going through its centre of mass; all other rotations and translations
are constrained. The elastic part is attached to the rigid end but not further
constrained.
35
Figure 16 Showing the fin with the rigid end which it rotates about, to the right
in grey.
3.4.2 Modelling hood and inner cylinder

The inner diameter of the hood is 146 mm and the outer diameter of the inner
cylinder is 128 mm. Both have a width of 120 mm, thus slightly wider than the
fin. Only a quarter of the hood and inner cylinder are modelled. The centreline
of the inner cylinder is shifted 5 mm in the y-direction and 6.5 mm in the
z-direction, compared with the centreline of the hood. Directions are given in
Figure 15.
The rigid material model and 2.5 mm x 2.5 mm fully integrated shell elements
were used.
3.4.3 Folding the fin
Three different folding procedures were investigated. Firstly either the hood
was translated against the inner cylinder or vice versa as seen in Figure 17. The
hood was then given translation in the –z-direction and the inner cylinder was
constrained in all directions or the inner cylinder was given a translation in the
z-direction and the hood was constrained in all directions. In both cases, the
magnitude of the translation was 75 mm in 1 s. The rigid end of the fin was free
to rotate about the x-axis and given the same translation as the inner cylinder.
36
a) b)
Figure 17 Showing find folding procedures where a) the hood is translated
towards the inner cylinder b) the inner cylinder is translated towards the hood.
Secondly a folding procedure where the hood is slid over the fin thus folding it
against the inner cylinder was used, which can be seen in
Figure 15. The end of the hood first coming in to contact with the fin was then
bent upwards with a radius of 20 mm to give a smoother folding. The inner
cylinder was constrained in all directions while the hood was set to rotate 1 rad
in 1 s about the x-axis going through its centreline. The rigid end of the fin was
free to rotate about the x-axis.
Contact algorithms
Two different contact algorithms were investigated when rotating the hood;
firstly the *AUTOMATIC_SURFACE_TO_SURFACE and secondly the
*AUTOMATIC_SURFACE_TO_SURFACE_MORTAR. When translating
either the hood or the inner cylinder, only the *AUTOMATIC_SURFACE_
TO_SURFACE contact was used. In all cases, the contact was assumed
frictionless and the curvature of the elements was considered. Contact was
modelled between the fin and the hood and the fin and the inner cylinder,
excluding the rigid end of the fin.
3.4.4 Deployment of the fin
The fin deployment was simulated in a separate simulation by importing the
stress state in the folded fin and rigid parts. The alternative where the hood was
37
slid over the fin using the *AUTOMATIC_SURFACE_TO_SURFACE_
MORTAR contact algorithm was used. The hood was then removed by
translating the hood in the –x direction, using the same contact algorithm as in
the folding simulation. A hood-removal time of 3.4 ms with constant
translational velocity was used.
Deployment with wind loads
The fin deployment simulation with wind loads was performed in the same way
as without wind loads, but inside an air domain. The ALE domain used was the
coarse mesh described in Section 3.3 with Donor Cell advection algorithm. To
test the principle of deploying a fin in the ALE domain without excessive
simulation times, the fin was deployed immediately, i.e. well before the flow
field was fully developed.
FSI modelling
The fluid structure interaction between the fin and the air was modeled using a
penalty based coupling algorithm, with 0.5 mm penetration giving 4 MPa
coupling pressure, as described in Section 3.3.5. The coupling grid was
however increased to 4x4 nodes per element as the mesh of the fin is slightly
coarser than ALE mesh surrounding it. The hood and the inner cylinder were
not included in the FSI modelling.
3.5 Full 3D model of fin deployment
Finally, a full 3D model containing eight fins and including the spin of the
projectile was developed. The fin deployment times and angular velocities of
the fins where evaluated using 1 Hz and 50 Hz spin rates; firstly without wind
loads and secondly with wind loads using 0° and 2° AoA.
3.5.1 Modelling the fins
The model of one fin together with the hood and inner cylinder described in
Section 3.4 was duplicated and rotated about the x-axis to give model with eight
fins at 45° separation. The fins were located at a circle with diameter 155 mm to
fit on the back end of the projectile. As the fins and hoods are slightly
overlapping, each fin and each hood was defined as separate parts. Contact was
defined between a fin and the hood covering it, but not with the other hoods and
not between the different fins.
38
In order to include the spin of the projectile, the rigid ends of the fins were
removed and replaced with constraints on the node set at the new edge. The
constrained node set at the edge of the fin, the hoods and the inner cylinders
were then given an angular velocity about the global x-axis. The projectile is
spinning in the counter-clockwise direction as seen from the front as in
Figure 18. The angular velocities used are given in Table 6.
A hood-removal time of 1.4 ms was used and to ensure that all fins deploy
simultaneously, they were constrained to have equal rotation about the local
x-axis at the constrained node set.
Figure 18 Showing the eight fins in grey and the inner cylinders and hoods in
black
Table 6 Angular velocities for different spin rates.

Angular
Spin rate [Hz] velocity [rad/s]
1 6.28
50 314
39
3.5.2 Modelling the air
The coarser mesh from Section 3.3 was mirrored at the symmetry plane to give
a full 3D model of the air. Only the outline of the projectile was given a slip
boundary condition while the remaining boundaries were open with a reference
pressure of 100.9 kPa. The Donor Cell advection algorithm was used.
An air flow velocity of 950 m/s was used with either 0° or 2° AoA giving the x-
and y-velocity components in Table 7. These velocities were used both as a
boundary condition on the ambient nodes and as an initial condition for all ALE
elements apart from in the wake behind the projectile.
Table 7 Showing inflow and initial velocities for different angles of attack, used
in the full 3D model.
Angle of attack x-velocity [m/s] y-velocity [m/s]

0° -950 0
2° -949.4 33.2
Mapping of air properties

The field variables of the ALE domain was mapped from a model containing
the ALE domain only to a second model also including the fin, hoods and inner
cylinders using *INITIAL_ALE_MAPPING. As one map file can be used for
several fin deployment simulations, this can considerable reduce the total
simulation time. A simulation time of 5 ms was used in the first simulation.
This means that the solution is not yet stabile. However a stabile solution would
require computational times exceeding one week, which was considered not
feasible.
3.5.3 FSI modelling
The fluid structure interaction between the fins and the air was modeled using a
penalty based coupling algorithm, with 0.5 mm penetration giving 4 MPa
coupling pressure, as described in Section 3.3.5. The coupling grid was
however increased to 4x4 nodes per element as the mesh of the fin is slightly
coarser than ALE mesh surrounding it. The hood and the inner cylinder were
not included in the FSI modelling.
40
3.6 Alternative model
An alternative model was considered where an FE-model of a projectile using

solid elements given by Christer Thuman at BAE Systems Bofors was
immersed in a square ALE mesh. The model, which does not include any fin,
can be seen in Figure 19.
Figure 19 Alternative setup showing the outline of the projectile in black and
the ambient section in dark grey.
The ALE domain was 1500 mm wide and 1500 mm high and extending 75 mm
in front of the projectile and 750 mm behind it. Approximately 620 000
hexagonal 1-point integration ALE elements were used to model the air. The
element size was smallest around the nose and the back end of the projectile,
4.5 mm x 5.5 mm at minimum and increasing outwards. The solid elements
used for the projectile was approximately the same size as the smallest elements
in the ALE mesh. The ALE mesh seen from the symmetry plane can be seen in
Figure 20 with the ambient section in darker grey and the outline of the
projectile in black.
41
Figure 20 Showing the ALE mesh and the outline of the projectile at the
symmetry plane of the alternative model
The front, top and sides of the ALE domain were defined as ambient, thus held
constant, and given an inflow velocity of Mach 2.7. The symmetry plane was
defined as a free slip wall and the back end was left open with a reference
pressure of 100.9 kPa. An initial velocity of Mach 2.7 was used for all nodes of
the ALE mesh excluding inside the projectile and in the wake behind it.
Penalty based FSI was used between the projectile and the air in the same ways
as between the fin and the air described in section 3.3.5. However a coupling
grid with 6x6 nodes per ALE element was used.
42
4 Results
The most important results are given in this section, starting with the risk
assessment. The results from comparable simulations in LS-DYNA and CFD++
are then analysed. This is followed by fin deployments, both with and without
wind loads, firstly of one fin and secondly of eight fins including spin. Ending
the section are some words on the results using the alternative model.
4.1 Risk assessment
Four different models of a small fin and surrounding air where used to simulate
air flow of Mach 3. Firstly the FSI between the air and the fin are visualized and
secondly the pressure and velocity in front of the fin are compared together with
total energy and computational times.
The spring-back of an elastic fin without FSI was also simulated and its
deployment time and angular velocity given here.
4.1.1 FSI and applying Mach 3 velocities
A small fin was subjected to wind loads of Mach 3 using four different models
as described in Section 3.2.1 which are summarised in Table 4. The table is
repeated below for clarity.
Model nr Fin elements Wind loads Advection scheme

1 Solid Air flow Donor Cell
2 Shell Air flow Donor Cell
3 Shell Air flow van Leer
4 Shell Moving fin Donor Cell
43
The total energy as a function of simulation time can be seen in Figure 21 for
all four models. When the air is flowing over the wing, the total energy is stable
after approximately 0.7 ms and almost immediately when the wing is translated
through the air. For the latter, the total energy is approximately 30% lower than
for the other.
Figure 21 Plots of total energy as a function of simulation time for the four
different models.
The velocity field around the fin show that there is an interaction between the
fluid and the structure, as the air flows around the fin rather than through it, and
that the shell element reference thickness of the fin is considered. The pressure
and density fields show a typical shock wave pattern with a pressure and
density increase in front of the fin and extending backwards in a V-shape. The
velocity, pressure and density of the air at an xy-plane at the top of the fin can
be seen in Figure 22 for model nr 2 at t=1 ms.
When using the elastic material model, there are stresses in the fins due to the
fluid-structure interaction, which can be seen in Figure 23.
The average resultant velocity and pressure of four elements in front of the
wing, which can be seen in Figure 24, are given in Table 4 for the four different
models at t=1 ms. The total computational time is also given in the same table
using 4 CPU’s on an 8 CPU Intel Core i7 3.2 GHz processor in Windows x64.
44
1025
820
616
411
207
3
a) Velocity in m/s
605
496
387
278
169
60
b) Pressure in kPa
3.03
2.55
2.07
1.59
1.10
0.62
c) Density in kg/m3
Figure 22 Variable fields of a) resultant velocity b) pressure and d) density at
an xy-plane at the top of the fin for model nr 2 at t=1 ms.
45
64.8
52.1
39.4
26.7
14.0
1.4
Figure 23 Showing von Mises stress in MPa in the fin at t=0.5 ms for model 1.
a) b)
Figure 24 Showing the fin in black and the elements used for comparison
between the different models in grey. Viewed from a) the front and b) the side.
Table 8 Comparison of computational time and average values of four elements

in front of the fin for the four different models.
Model nr Resultant Pressure Computational
velocity [kPa] time [min]
[m/s]
1 467 379 12
2 496 356 13
3 504 380 20
4 679 212 18
46
When comparing the resultant velocity and pressure, the three models where air
is flowing over the wing gave similar results. However, when the wing is
translated through the air, the velocities are considerable higher and the
pressure considerable lower.
The van Leer advection algorithm increase computational time compared with
the Donor Cell algorithm and so does translating the wing through the air
compared with air flowing over the wing. However, translating the fin through
the air gives a stable solution much faster thus requiring shorter simulation time
if only a stabile air flow is desired.
4.1.2 Elastic springback of a fin
An elastic fin with a rigid end is loaded with 100N at the elastic end as
described in section 3.2.2. The maximum von Mises stress just before
unloading is slightly larger than 1200 MPa, which can be seen in Figure 25.
When unloaded, the fin rotates 90° about the y-axis in 1.3 ms. The angular
velocity as a function of time can be seen in Figure 26. The fin is unloaded at
between t=1.1 ms and t=1.2 ms and the fin deployment is rather fluttery as can
be seen from the sinusoidal-shaped angular velocity.
1234
988
741
494
247
0
Figure 25 von Mises stress in the fin just before unloading given in MPa.
47
Figure 26 Angular velocity as a function of time of a fin due to unloading
between t=1.1 ms and t=1.2 ms.
4.1.3 Doing all at once
A simulation of a fin translating through an air flow of Mach 3, as described in
section 3.2.3, show that the field variables change as the wing is moved. The
change in velocity when the fin is moved 4.2 cm in the –y direction can be seen
in Figure 27 at an xy-plane at the top of the fin. The movement of the fin
neither affected the time step nor the total simulation time.
t=0.3 ms
1070
856
642
589
428 t=4.5 ms
214
Figure 27 Showing that the velocity field changes as the fin is translated in the
–y direction. The xy-plane is at the top of the fin and velocities are in m/s
48
A comparison is made between simulations performed in CFD++ and LS-

DYNA of the air flow around a half-symmetry model of a projectile with two
static fins as described in Section 3.3. Firstly the projectile is traveling at Mach
2.7 with
0° AoA and the velocity and pressure fields at the symmetry plane are
compared. Secondly the AoA is increased to 2° and the resultant forces acting
on the fin are compared.
4.2.1 Projectile travelling with zero AoA
The resultant velocities and pressures at the symmetry plane can be seen in
Figure 28 and Figure 29 for three different models described in Section 3.3. To
summarise, the three different models are
a) Using CFD++ and a fine mesh
b) Using LS-DYNA, a fine mesh and Donor Cell advection
c) Using LS-DYNA, a coarse mesh and Donor Cell advection
CFD++ finds a steady state solution while the simulation in LS-DYNA is
dynamic and a simulation time of 7ms was used. The variable limits from
CFD++ are used for displaying the LS-DYNA results as well. Note that the
colour bars are slightly different; the top one is for a) and the bottom one for b)
and c).
Visually comparing the different models the results are fairly similar. A shock
wave is formed in front of the projectile, extending backwards and a smaller
one is also formed in front of the fin. The maximum pressure in front of the
nose is however slightly higher for the simulations in LS-DYNA, 1075 kPa
compared with 927.7 kPa in CFD++.
49
976.9
782.2
a) 587.5
392.8
198.1
33.8
976.9
b)
782.2
587.5
392.2
198.1
33.8
c)
Figure 28 Velocity at the symmetry plane in m/s. Comparison between

a) CFD++, fine mesh b) LS-DYNA, fine mesh c) LS-DYNA coarse mesh.
50
972.7
785.5
a) 598.3
411.1
224.0
36.8
972.7
b)
785.5
598.3
411.1
224.0
36.8
c)
Figure 29 Pressure at the symmetry plane in kPa. Comparison between

a) CFD++, fine mesh b) LS-DYNA, fine mesh c) LS-DYNA coarse mesh.
51
4.2.2 Projectile travelling with 2° AoA
A projectile travelling with 2° AoA was simulated by introducing a velocity in
the y-direction as described in Section 3.3. Four different models were used
which can be summarised as;
a) Using CFD++ and a fine mesh
b) Using LS-DYNA, a fine mesh and Donor Cell advection
c) Using LS-DYNA, a coarse mesh and Donor Cell advection
d) Using LS-DYNA, a coarse mesh and vanLeer advection
The net force acting on the fin in the y-direction for the four models can be
found in Table 9. CFD++ finds a steady state solution while the simulation in
LS-DYNA is dynamic and the results are after a simulation time of 7 ms. The
simulation time for the b) is 113 h while for the coarser c) only 27 h, using 4
CPU’s on an 8 CPU Intel Core i7 3.2 GHz processor in Windows x64.
Table 9 Total force in y-direction on fin from simulations using CFD++ and at
t=7 ms using LS-DYNA.
Model Y-force on fin
[N]
a) 498.9
b) 686.0
c) 663.5
d) 685.1
Comparing the force on the fin as a function of time for the three models in
LS-DYNA b), c) and d), which can be seen in Figure 30 however show that the
solution is not stable after 7ms as the net force is increasing with time. This can
also be seen by looking at how the kinetic energies change with time, given in
Figure 31, which are still decreasing after 7 ms simulation.
52
Figure 30 Showing net force acting on the fin in the y-direction as a function of
time for models b) c) and d) using 2°AoA.
Figure 31 Showing kinetic energy as a function of time for models b) c) and d)

using 2°AoA.
A longer simulation of 60 ms was performed with model c), showing the
shortest computational time per ms simulation. The y-force on the fin and the
kinetic energy as functions of time can be seen in Figure 32 and Figure 33. The
force level out at 751.6 N after 27 ms while the kinetic energy still have a slight
increase but is fairly stable after 30 ms. Compared with the simulation in
CFD++, the y-force is thus 50% larger.
53
Figure 32 Showing net force acting on the fin in the y-direction as a function of
time for model c) using 2°AoA.
Figure 33 Showing kinetic energy as a function of time for model c)

using 2°AoA.
A fin deployment simulation was done in two steps. Firstly the fin was folded
and the stresses in the fins using three folding procedures and two contact
algorithms were investigated. Secondly the fin was deployed, with and without
wind loads, and the deployment time and angular velocities of the fin
compared.
4.3.1 Folding the fin
Three different procedures for folding the fin were investigated, translating the
hood towards the inner cylinder, translating the inner cylinder towards the hood
or sliding the hood over the fin.
54
Translating either the hood or the inner cylinder towards the other did not give a
smooth folding of the fin. The fin was given an unnatural curvature, which can
be seen in Figure 34. The resulting stress in the fin was also higher than
anticipated and the maximum stress found too far away from the rigid end.
Unnatural curvature
Figure 34 Folded fin by translating the inner cylinder towards the hood giving
an unnatural curvature of the fin.
Sliding the hood over the fin gave a smoother folding of the fin and the highest
stresses were found closer to the rigid end as anticipated. The von Mises stress
in the folded fin using the two different contact algorithms can be seen in
Figure 36. The MORTAR contact gave lower maximum stress, 1724 MPa
compared with 1894 MPa, and also lower contact force between the hood and
the fin which can be seen in Figure 35. Comparing computational time, using
the MORTAR contact resulted in an increase from 2 h to 4 h using 4 CPU’s on
an 8 CPU Intel Core i7 3.2 GHz processor in Windows x64.
55
Figure 35 Contact force between the hood and the fin when sliding the hub over
the fin using the two different contact algorithms.
1894
1519
1143
767
392
16
a) b)
Figure 36 Showing von Mises stress in MPa in the folded fin using
a)*AUTOMATIC_SURFACE_TO_SURFACE
b) *AUTOMATIC_SURFACE_TO_SURFACE_MORTAR.
56
4.3.2 Deployment of the fin
The fins were deployed by removing the hood in 3.4 ms, with and without wind
loads, as described in Section 3.4.4. The resulting fin deployment times and
angular velocities at the fully deployed state are given in Table 10. The angular
velocities of the fins as functions of time are also given in Figure 37.
The fin deployment is rather fluttery as can be seen from the sinusoidal-shaped
angular velocity. The wind loads increase the amplitude of the fluctuations, and
also the mean value resulting in a shorter deployment time. Without wind loads
the angular velocity is fluctuating about 260 rad/s and with wind loads about
300 rad/s.
The deployment of the fin with wind loads seen at a zx-plane through the
middle of the fin can be seen in Figure 38 together with the velocity field
surrounding it.
Table 10 Deployment time and angular velocity at the fully deployed state.
Deployment Angular
Wind loads
time [ms] velocity [rad/s]
No 10.3 278
Yes 8.9 225
 With wind loads  Without wind loads
Figure 37 Angular velocities of the fins during deployment. The dashed line
show when the hood is fully removed.
57
976.9
782.2
587.5
392.2
t=3.4ms t=6.0ms t=8.9ms 198.1
33.8
Figure 38 Showing the velocity field around the fin at a zy-plane through the
middle of the fin during deployment. Velocities given in m/s.
The deployment times and angular velocities without wind loads, using the
model described in Section 3.5, are presented for three different spin rates.
Some errors encountered when mapping the properties of the ALE domain are
then described followed by a fin deployment sequence with wind loads where
the fins however did not deploy simultaneously.
4.4.1 Deployment without wind loads
The fins deployment times and angular velocity at their fully deployed state are
given in Table 9 for different spin rates. The hood removal time is 1.4 ms and
included in the fin deployment time. A fin deployment sequence with zero spin
rates can be seen in Figure 39. The angular velocity is given for the fin located
at the top at the beginning of the simulation. The angular velocity of this fin as
function of time can also be seen in Figure 40.
Comparing the deployment times show that increasing the spin decrease the
deployment time. The angular velocities as functions of time show, just as for
the deployment of one fin, that the deployment is fluttery. These also show that
even if the angular velocity at the fully deployed state is higher at 0 Hz spin
than 1 Hz, the mean velocities are almost identical. At a spin rate of 50 Hz, the
angular velocity increase throughout the simulation.
58
Table 11 Fin deployment times and angular velocities at the fully deployed state
for different wind loads and spin rates.
Wind loads Spin rate Deployment Angular velocity
[Hz] time [ms] [rad/s]
0 6.8 338
None 1 6.7 319
50 4.3 733
t=0ms t=1ms t=3ms
t=5ms t=6.8ms
Figure 39 Fin deployment sequence at zero spin rate without wind loads.
59
Figure 40 Angular velocities of a fin after hood removal for different spin rates.
The dashed line show when the hood is fully removed.
4.4.2 Deployment with wind loads
To simulate a fin deployment with FSI, the variable fields of an air domain was
mapped from a model only containing the air to a model also including the fins
hoods and inner cylinders as described in Section 3.5.2.
The mapping of the velocity field was not correct at a few nodes of the ambient
section which were given zero velocity. The resultant velocity of the air using
2°AoA immediately after mapping can be seen in Figure 41, viewed from
above. As these nodes with zero velocity are a part of the ambient section, with
a prescribed motion and thermodynamic state, the velocities were correct at the
next time step.
There were an accurate mapping of the density, but the pressure field after
mapping was not right at all. All elements were given pressures between 99 kPa
and 103 kPa. However, the pressure is given by the equation of state if density
and specific heats are known, so after a time step it was recalculated to correct
values. The pressure field at the symmetry plane immediately after mapping
and then 1ms later, using 2° AoA, can be seen in Figure 42.
In the Figure 42 b) one can see how the FSI with the wings have changed the
pressure field at the back end of the projectile. The deployment of the fins with
wind loads was however not successful. They were not constrained correctly to
deploy simultaneously, resulting in varying deployment times, which can be
seen for 2° AoA and 1 Hz spin in Figure 43.
60
950.0
760.0
570.0
380.0
190.0
0.0
Figure 41 Resultant velocity of the air directly after mapping, seen from above.
1015.0
a)
812.5
609.7
406.9
204.1
1.3
b)
Figure 42 Pressure in kPa at symmetry plane a) directly after mapping

b) 1 ms later.
61
Figure 43 Fins and inner cylinders after 5 ms deployment using 2° AoA and
1 Hz spin, seen from the font.
The pressure and velocities at the symmetry plane using the alternative model,
as described in Section 3.6, can be seen in Figure 44 at a simulation time of 7
ms. This show that there is a pressure rise at the nose, extending backwards,
and that the velocity field is bent around the projectile. There is however
leakage through the Lagrangian structure, as there is both a pressure and a
velocity increase at some ALE elements and nodes located inside the projectile.
62
976.9
782.2
587.5
a)
392.8
198.1
33.8
972.7
785.5
598.3
b)
411.1
224.0
36.8
Figure 44 a) velocity in m/s and b) pressure in kPa at the symmetry plane.
63
64
5 Analysis and discussion
The results from the four main parts; the risk assessment, the comparison with
CFD, the fin deployment of one fin and the full 3D fin deployment, are
analysed here. This is followed by some words on the results using the
alternative setup and the section ends with suggestions on future work.
5.1 Risk assessment
The four basic questions to be answered with the risk assessment where;
1. Can the FSI between a fin and air be captured?
2. Can velocities of Mach 3 be applied?
3. Can the elastic springback of a fin be captured?
4. Can all of the above be simulated at the same time?
The short answer to all four questions is simply; yes. There is clearly a fluid-
structure interaction between the fin and the air as the air flows around the fin
rather than through it, either using solid or shell elements to model the fin. The
FSI can also be seen as there are stresses in the fins due to the wind loads. The
compressibility of the air is also captured, and the supersonic velocity of the air
gives a shock wave extending from the wing as anticipated. The mesh used is
however too coarse to properly resolve it.
Using Donor Cell or vanLeer advection algorithms gave approximately the
same average pressure and velocity when looking at four elements in front of
the wing. Modelling the wing with either shell or solid elements also rendered
similar results.
However whether the air is flowing over the wing or the wing is moved through
the air made a large difference. When translating the wing through the air, the
resultant velocity was 40% higher and the pressure 40% lower and the total
energy of the system was also reduced by 30%. The fact that no ambient section
is defined for this model might be the reason for the differences.
One of the questions with an explicit simulation of the air which will not be
encountered when using CFD-simulation is how to determine the length of the
65
simulation. One does not know on beforehand when the variable fields are
stable, if they ever are.
In these models, the addition of hourglass viscosity for the ALE domain is not
necessary. The mesh does not deform and will thus not encounter hourglass
problems. However, using a too high hourglass parameter, e.g. the standard for
steel 0.1 will introduce unwanted viscosity. Specifying a low hourglass
parameter for the air will avoid any problems using default hourglass viscosity
and furthermore avoid warning messages.
Hexagonal elements have to be used when modelling the ALE domain, which is
rather inflexible.
The simulation of springback of the fin gave reasonable values of deployment
time and stresses in the fin due to the given load. Finally the movement of a fin
orthogonal to the flow direction could be simulated, without problems and
without increasing computational time, so a fin deployment simulation
including FSI should be possible.
When making the ALE mesh for comparison with CFD++, the outline of the
projectile is defined by the outline of the mesh. This has both pros and cons
compared with immersing a solid mesh for the projectile in a square shaped
ALE domain;
+ The same mesh can be used in the CFD simulation
+ No elements are wasted on meshing the solid projectile or the ALE mesh
inside the projectile
+ No FSI-coupling is needed between the projectile and the air, and there is
no risk of leakage
— The need to follow the shape of the projectile makes meshing more
difficult
— The hexagonal elements will be more distorted then when using a square
domain
— The forces acting on the projectile will not be captured
66
The comparison of the variable fields at the symmetry plane between
simulations in LS-DYNA and CFD++ show reasonable similarities. By
assuming an inviscid fluid, the wake behind the projectile can not be properly
simulated. Furthermore, as shear stresses in the fluid are disregarded, there will
be no drag forces; however pressure can be accurately predicted. For the
deployment, the pressure acting on the fins is larger, and more important for the
fin deployment, than the frictional forces between the air and the fins. The
assumption of an inviscid fluid is thereby reasonable.
The computational times using the finer mesh is inacceptable long. A slow fin
deployment takes approximately 10 ms, which would result in simulation times
of just less than a week, for a half-symmetry model. And that does not include
stabilisation the ALE variable fields before mapping. A 10ms simulation using
the coarser mesh would still take more than 1.5 days, which is still long. Apart
from waiting a long time to get the results, long simulation times make losses in
network connections or power cuts an issue.
When comparing the y-force acting on the fin for 2°AoA the resulting force
using LS-DYNA is approximately 50% higher than using CFD++. The results
form CFD++ are more likely to be accurate, as CFD simulations has been more
widely used to model such air flows and CFD++ has given correct results in
similar simulations validated with measurements at BAE Systems Bofors.
One possible source of error is the distorted mesh used for the ALE domain.
Even if the same mesh is used in both LS-DYNA and CFD++, CFD++ might be
better on handling the skewed elements with varying mesh size than LS-
DYNA. CFD++ uses the finite volume method and several types of elements
can be used, while the ALE formulation in LS-DYNA requires hexagonal
elements. The mesh is furthermore rather coarse in some areas, including the
shock wave extending away from the nose of the projectile. Mesh refinement in
some areas might be needed.
There could also be errors originating in the contact definition between the fluid
and the structure. For the coarse mesh and Donor Cell advection algorithm, a
simulation was made using 1/10 of the contact pressure so that 0.5 mm
penetration gave 0.4 MPa contact pressure, with identical results. To investigate
this further, one could look at a contour plot of the contact pressure at the wing
and how it develops with time. However, as such a fine mesh was used to
67
model the wing the files for plotting the pressure contours were too large to
open.
Validation between CFD++ or analytic solutions and LS-DYNA for a few
simple models, such as the flow around a thin plate or cylinder, is suggested to
see if the correct forces acting on the structure are found.
The Donor Cell advection is chosen for fin deployment simulations as the
vanLeer advection increase computational times. The vanLeer advection is
furthermore only second order accurate on rectangular grids, and the ALE mesh
used here is rather distorted, apart from around the fin.
The multi-material formulation recommended even for problems including only
one ALE material. This also makes it easy to check whether there is leakage
through the structures by switching material group if a fluid element passes the
structure. Including a dummy vacuum material define different material groups
for the air and the vacuum avoids warning messages when running the
simulation. The alternative would be to use a single material formulation and
then define void areas where the Lagrangian parts are initially located, which
can be rather cumbersome.
The fin deployment simulation is done in two steps so that the stress state in the
fins can be initialised only once and used for several deployment simulations.
In reality, the fins are folded before gun launch and only the deployment is of
interest for an FSI simulation.
Assuming that the material of the fin is isotropic and linear elastic is reasonable
as long as the fin is subjected to moderate temperatures and stresses below the
yield stress.
The smoothest curvature and lowest stress in the fin after folding was when
sliding the hood over the fin. Using the *AUTOMATIC_SURFACE_TO
_SURFACE_MORTAR contact gave even lover stress and was thus used as it
was still somewhat higher than anticipated. Even if this contact increased
computational time, it was still not considerable long. In reality, the hood is
not rigid but deforms thus giving a larger radius of curvature leading to lower
68
stress in the fin. The stress in the fin increase somewhat during the first stage of
deployment, reaching a maximum when the hood is removed by 2/3.
Folding the fin by either translating the hood or the inner cylinder towards the
other resulted in an unnatural curvature of the fin and too high stress. This can
be due to the simplifications made when modelling the rotational joint between
the inner cylinder and the hood and or the location of the rigid parts relative to
each other. It can also be due to dynamic effects when pushing the fin towards
the rigid part. An implicit simulation was done with contact problems giving an
elongated fin and very high stress. A small time step was furthermore needed
for the contact algorithm thus not reducing simulation time compared with an
explicit simulation.
The deployment time of the fin is in the anticipated order of magnitude when
comparing with earlier simulations and test shootings. In reality, the rigid end
of the fin has a greater mass and moment of inertia than the shell elements used
and there is also friction in the rotation, both leading do longer deployment
times.
The fluttery nature of the deployment was also expected and so was the fact that
wind loads decrease the deployment time. As the hood leaves the front end of
the fins first, they start deploying slightly which leads to wind loads in the
deployment direction.
When making the full 3D model with eight fins, they fins were located on a
circle with diameter 155mm. When in immersed in the ALE domain they will
then fit on the back end of the simplified model of the projectile. If placed
outside of the ALE domain, there might take some time for the FSI contact
algorithms to notice that they have entered it. This can also lead to large or
noisy contact forces when the structure is recognised inside the ALE domain.
The spin rate of the projectile increase the angular velocity of the fins and
reduce the deployment time. This was expected, as the gyroscopic effects help
the deployment.
69
The spin rate of 50 Hz gave too high stress in the fins, even if held folded in
between the hood and the inner cylinder. They were not stationary, but
vibrating up and down, so some damping might be needed.
The mapping of the air properties using *INITIAL_ALE_MAPPING was rather
troublesome. It did not work at all in the latest R7.0.0 solver but the earlier
version R6.1.1 had to be used. Any changes made in between the two versions
important for this work were not found. The mapping of the pressure and
velocity were furthermore not correct, and even if it was put right in the next
time step, it does not feel dependable.
One of the advantages of *INITIAL_ALE_MAPPING is that it can map
solutions from one mesh to another, e.g. from 2D to 3D or from a fine mesh to a
coarser, which is not what is done here. A restart analysis could be used instead,
where the fins are added in the restart. However, adding parts requires a full
restart analysis, which is currently not available for the Shared Memory Parallel
licence used.
The full 3D model containing 1.1 million ALE elements and 30k shell elements
gave small restart dump files of 5 GB, i.e. too large to handle in Windows. This
was rather troublesome as the possibility of checking the results at regular
intervals and then restarting the simulation is convenient, especially when
having long simulation times.
The full 3D model did only increase the simulation time by 5% compared with
the half-symmetry model with one fin. The major issue is thus not the number
of elements, but the time step size. The smallest characteristic element length,
which can be smaller than the smallest element length for a skewed hexagonal
element, is determining the time step. Using a rectangular mesh for the ALE
domain would make it easier to control the time step size.
The fin deployment with wind loads was not successful as the constraint to
make all fins deploy simultaneously was not correct. However, the fact that
they did not deploy simultaneously is due to the wind loads, so there is fluid-
structure interaction.
70
An alternative model is used where a solid mesh for the projectile is immense in
a square shaped ALE mesh. One can then make use of existing models for the
projectile, and the meshing of the ALE domain is simplified.
The setup of this model was much faster than the one used for CFD comparison
and fin deployment, both due to the easy meshing of a square ALE domain and
as it was much easier to select appropriate elements and nodes for applying
boundary conditions and defining the ambient section.
The velocity and pressure fields at the symmetry plane are not unrealistic, there
is however leakage through the Lagrangian structure. Further work is needed to
if one wants to compare the two model types in a fin deployment simulation.
5.6 Future work
The work in this report is focused on evaluating the capabilities of ALE and FSI
simulations of a fin deployment in LS-DYNA. A start for future work could be
to continue with the full 3D model and accurately constrain the wings to deploy
simultaneously with wind loads. The deployment times using different angles
of attack and spin rates can then be compared.
The difference in the force acting on the fin compared with CFD++ could also
be investigated further or a comparison made with measurements. Suggestions
are;
 Look closer into the velocity and pressure fields surrounding the fin and
see if they are similar.
 Use a coarser mesh for the fin or select a few of the elements and
examine the development of the pressure contours.
 Make mesh-refinements in some critical areas, e.g. in the shock wave
front
 Investigate different settings of the advection and FSI algorithms and the
effect of bulk and hourglass viscosity
But also to model some very simple cases, such as the flow over a flat plate or a
cylinder, and compare with CFD++ and/or analytic solutions to see if the
71
correct forces are found. One path to follow could then be to increase the
complexity of the model by for example;
 Making a more realistic model of the projectile, the fin and the hood
 Include friction in the rotational joint between the fin and the hood
 Model the air as a viscid fluid
 Include temperature effects
A realistic model could then be compared with fin deployment times from high
speed films of test firings of relevant projectiles.
Another path would be to aim to reduce the simulation time and the size of the
models. One could investigate much the wind loads affect the fin deployments,
and how realistic the wind loads must be to give correct fin deployment times.
As the deployment times are realistic even with a non-stabilised flow field this
suggests that the wind loads doesn’t have to be modelled with high accuracy to
give correct deployment times. A very coarse mesh, not resolving the shock
wave or giving accurate variable fields would be good enough. This could be
combined with a more complex model of the fins.
The alternative setup could also be further developed, including fin deployment.
This type of model would probably be more suitable for a mesh convergence
study, as doubling or halving the element size is more straightforward.
The new CE/SE solver in LS-DYNA is also an interesting area of investigation.
This is partly developed to handle compressible flows with shock waves and
utilise a completely different way of calculating. Furthermore, other than
hexagonal, e.g. tetrahedral, elements can be used for meshing which is
convenient.
72
6 Conclusions
A fin deployment on an artillery projectile travelling at supersonic speed can be

simulated using an FSI simulation in LS-DYNA using ALE formulation. The
fin deployment times are in the interval found from gun firing test of relevant
projectiles.
More work is however needed before a suitable methodology is found and it
has to be developed to suit the available computational capacity. The models
used in this work require fast processors, a lot of memory and an operating
system being able to handle large files.
Further validation of the computation of the wind loads at supersonic speeds in
LS-DYNA is needed. Furthermore, the problems with the initial ALE mapping
and full restart analyses needs to be solved in order to have convenient ways of
modelling were simulations can be done in several steps.
An FSI simulation with ALE formulation differs in many ways from a pure
CFD or a pure structural FE simulation. The dynamic interaction between fluids
and structures is captured and a range of new phenomenon can be analysed. A
numerical model that captures the dynamics of a fin deployment and the forces
acting on the structure will give better understanding of the deployment
sequence and open up new possibilities of optimisation.
73
74
7 Bibliography
[1] R. L. McCoy, Modern Exterior Ballistics, The Launch and Flight

Dynamics of Symmetric Projectiles, Atglen, PA: Schiffer Publishing Ltd.,
1999.
[2] B. Zakrisson, B. Wikman and H. -Å. Häggblad, “Numerical Simulations of

Blast Loads and Structural Deformation from Near-Field Explosions in
Air,” International Journal of Impact Engineering, pp. 38 (7) 597-612,
2011.
[3] M. S. Chafi, G. Karami and M. Ziejewski, “Numerical analysis of blast-

induced wave propagation using FSI and ALE multi-material
formulations,” International Journal of Impact Engineering , pp. 36,1269-
1275, 2009.
[4] M. Štiavnický and P. Lisý, “Gunshot Effects Simulations,” Science and

Military, p. 1, 2011.
[5] R. van Loon, P. D. Anderson, F. N. van de Vosse and S. J. Shervin,

“Comparison of various fluid-structure interaction methods for deformable
bodies,” Computers and Structures, pp. 83 ,833-843, 2007.
[6] J. D. Anderson, Fundamentals of Aerodynamics, New York: McGraw-Hill

Companies, Inc., 2011.
[7] B. Zakrisson, “Numerical simulations of blast loaded steel plates for

improved vehicle protection,” Universitetstryckeriet, Luleå, 2013.
[8] T. P. Slavik, “A Coupling of Empirical Explosive Blast Loads to ALE Air

Domains in LS-DYNA,” in 7th European LS-DYNA Conference ,
Livermore, California, U.S.A, 2009.
[9] LSTC, LS-DYNA Keyword User's Manual Vol I, Version 971 R6.1.0,
Livermore: Livermore Software Technology Corporation, 2012.
75
[10] J. A. Çengel, J. M. Cimbala and M. Kanoglu, Fluid mechanics :
fundamentals and applications, New York: McGraw-Hill, 2012.
[11] “NASA,” [Online]. Available: http://www.grc.nasa.gov/WWW/k-

12/airplane/reynolds.html. [Accessed 14 May 2013].
[12] Y. Fung and P. Tong, Classical and Computational Solid Mechanics,

Vol.1., World Scientific Publishing Co. Pte. Ltd., 2001.
[13] J. William D. Callister, Materials Science and Engineering: An

Introduction, John Wiley & Sons, Inc., 2007.
[14] K.-J. Bathe, Finite Element Procedures, Upper Saddle River, New Jersey:
Prentice-Hall, Inc. , 1996.
[15] K. H. e. a. Huebner, The finite element method for engineers, 4th edition,
John Wiley & Sons, Inc., 2001.
[16] D. Hutton, Fundamentals of Finit Element Analysis, McGraw Hill, 2004.
[17] C. b. J. O. Hallquist, LS-DYNA Theory Manual, Livermore, California:

Livermore Software Technology Corporation, 2006.
[18] S. C. Chapra, Applied Numerical Methods with MATLAB, for Engineers

and Scientists., McGraw Hill, 2006.
[19] J. Donea, A. Huerta, J.-P. Ponthot and A. Rodriguez-Ferran, Arbitrary

Lagrangian-Eulerian Methods.
[20] L. Olovsson, LS-DYNA Training class in ALE and fluid-structure

interaction, Livermore Software Technology Corporation, 2006.
[21] I. Do and J. Day, "ALE Overview LS-DYNA®," 28 09 2005. [Online].

Available: http://awg.lstc.com/tiki/tiki-index.php?page=LSTC+Tutorials+
and+Topic+Presentations. [Accessed 21 03 2013].
[22] M. S. Gadala, “Recent trends in ALE formulation and its applications in

solid mechanics,” Computer methods in applied mechanics and
engineering, pp. 193, 4247-4275, 2004.
76
[23] D. J. Benson, “Momentum advection on a staggered mesh,” Journal of
Computational Physics, pp. 143-162 , 1992.
[24] J. VonNeumann and R. D. Richtmyer, “A Method for the Numerical

Calulation of Hydrodynamic Shocks,” Journal of Applied Physics, pp. 21,
232-237, 1950.
[25] “LS-DYNA Support, Hourglass,” LSTC Inc, [Online]. Available:

http://www.dynasupport.com/howtos/element/hourglass. [Accessed 10
June 2013].
[26] B. Zakrisson, H.-Å. Häggblad and P. Jonsén, “Modelling and simulation of

explosions in soil interacting with defomable structures,” Centeral
Eurpoean Journal of Engineering, pp. 2 (4), 532-550 , 2012.
77
78
Appendix; Meshing of the ALE domain
Three different meshes was used for the ALE domain representing the air
around a projectile, were the outline of the projectile was defined by the mesh.
One coarse and one fine for the half-symmetry model and one coarse for the
full 3D model. Figures with the coarse half-symmetry model are given in this
appendix. The fine half-symmetry mesh is similar, but having approximately
35% shorter element lengths and the full 3D model is simply the coarse half-
symmetry mesh mirrored at the symmetry plane.
The numbers of hexagonal solid elements used are:
1. Coarse half-symmetry: 551 554
2. Fine half-symmetry: 1 320 966
3. Coarse full 3D: 1 103 108
The main focus when meshing the air domain was the area around the fin. The
fine half-symmetry model was used in CFD++ as well, with a plane of elements
removed to define the outline of the fin. 1x1x1mm elements were used to model
this area. A fairly square volume surrounds the fin where the aspect ratio is
below 2.5 and the squish, i.e. deviation from orthogonality, less than 0.1.
Secondly the region in front of the nose was considered as there are large
gradients in pressure, density and velocity in that region. The smallest element
lengths are found here and the aspect ratio is below 3.1 and the aspect ratio less
than 0.49. Thirdly the wake behind the projectile is a region of complex flow
and here the aspect ratio is below 5 and the squish less than 0.49. Smaller and
more equal sized elements are used close to the projectile and larger elements
with worse aspect ratios further away.
79
1795
902 735
300
1040
Figure 45. Coarse half-symmetry mesh seen from the top. Measurements in mm.
370
150
Figure 46. Coarse half-symmetry mesh seen from the side.

Measurements in mm.
80
1500
750
Figure 47. Coarse half-symmetry mesh seen from the front. Measurements
in mm.
Figure 48 Coarse half-symmetry mesh seen from the back
81
A
C C
A
Figure 49. Coarse half symmetry mesh seen from the bottom
Figure 50. Mesh at symmetry plane, close up at back end of

projectile
82
Figure 52. Mesh at symmetry plane, close up at nose of projectile
Figure 51. Mesh at cut A-A
83
Figure 53. Mesh at cut A-A, close-up .
Figure 54. Mesh at cut B-B
84
Figure 55. Mesh at cut B-B, close up.
Figure 56. Mesh at cut C-C
85
Figure 57. Mesh at cut C-C, close up at back end of projectile
Figure 58. Mesh at cut C-C, close up at nose of projectile
86

Fulltext02 PDF

Uploaded by

Copyright:

Available Formats

Fulltext02 PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Fulltext02 PDF

Uploaded by

Copyright:

Available Formats

MASTER'S THESIS

Supersonic Artillery Projectile Fin

Master of Science in Engineering Technology

Luleå University of Technology

Appendix; Meshing of the ALE domain 79

This master thesis is a part of the degree project in Mechanical Engineering

The goal of this master thesis is to evaluate whether Fluid-Structure Interaction

1.3 Problem description

A hypothetical artillery projectile of caliber 155 mm is analysed. It has fins

Figure 2 Illustration of a projectile before and after fin deployment.

LS-DYNA is a general purpose FE program developed by Livermore Software

Some relevant theory is briefly described here, starting with governing

2.1 Inviscid compressible air flow

where is the specific internal energy.

Figure 3 Illustration of a) subsonic flow with smooth streamlines and b)

Figure 4 Detached shock wave in front of a blunt body.

2.2 Isotropic linear elastic solids

where are constants independent on the magnitude of the

FE analysis is a numerical technique extensively used in engineering analysis

Huebner et al. [15] describes the FE method as a seven step-procedure:

to be stable [9]. is the highest eigenfrequency of the system.

2.4 Reference frames

Different kinematic descriptions, reference frames, can be used to describe the

2.5 Advection algorithms for ALE

Using an ALE formulation, there is advection of material due to the relative

and is the velocity at node [17].

2.5.2 Van Leer algorithm

and an approximation of the slope. [17]

2.5.3 Critical time step regarding advection

where is the number of elements, a characteristic size of an element

where is defined by the difference between the material velocity ̇ and

When a fluid and a structure are in contact, there is continuous interaction

Before impact After impact

ALE reference Lagrangian Penetration

In order to handle the almost discontinuous changes in flow properties in a

2.8 Hourglass viscosity

When having under integrated elements, there is a need to control so called

3.1 Material modelling

3.1.2 ALE parts

Table 3 Material properties for the ideal gas equation of state.

A risk assessment was made in order to identify possible problems at an early

1 Solid Air flow Donor Cell

Axis of rotation 100 N load

Figure 9 Showing a loaded fin together with the axis of rotation.

3.2.3 Doing all at once

In order to evaluate the possibilities of simulating the air flow in LS-DYNA,

Figure 14 Showing initial velocities at an xz-plane through the middle of the

3.4 Deployment of one fin

3.4.1 Modelling the fin

3.4.2 Modelling hood and inner cylinder

3.5 Full 3D model of fin deployment

Table 6 Angular velocities for different spin rates.

Angle of attack x-velocity [m/s] y-velocity [m/s]

Mapping of air properties

An alternative model was considered where an FE-model of a projectile using

Model nr Fin elements Wind loads Advection scheme

Table 8 Comparison of computational time and average values of four elements

A comparison is made between simulations performed in CFD++ and LS-

Figure 28 Velocity at the symmetry plane in m/s. Comparison between