Simulation and modelling of pellets in large storing systems

(1)

M A S T E R’S T H E S I S

GUSTAF GUSTAFSSON

Simulation and

Modelling of Pellets in Large Storing Systems

MASTER OF SCIENCE PROGRAMME Mechanical Engineering

Luleå University of Technology

Department of Applied Physics and Mechanical Engineering Division of Solid Mechanics

(2)

PREFACE

This thesis was written in request of LKAB at the Division of Solid Mechanics at Luleå University of Technology. It is the final project for the Master of Mechanical Engineering and is written as a part of the Research Trainee program.

Many people have contributed directly or indirectly with the completion of this study. First I would like to thank my supervisor, Professor Hans-Åke Häggblad for his help and support during the work. I would also like to thank Professor Sven Knutsson at the Division of Soil Mechanics for his help with the evaluations of the experimental data. Further I would like to thank Mats Strömsten, Simon Töyrä and Dr. Kent Tano and all the other participants in the project group at LKAB for their ideas and feedback.

I thank the students of the Research Trainee group for their advice and support during the year and finally I thank my family and friends for their support and encouragement.

Gustaf Gustafsson Luleå, June 2006

(3)

ABSTRACT

This thesis is the result of a project done in request of LKAB at the Division of Solid Mechanics at Luleå University of Technology.

Forming and handling of iron ore pellets is an important part in the converting process for LKAB. Knowledge about this sub process is very important for further efficiency progress and increasing product quality. The aim of this work was to create a numerical simulation tool describing the mechanics of a moving bed of iron ore pellets.

A mesh free particle method called Smoothed Particle Hydrodynamics has been used for the numerical simulation. The code for the method was rewritten with equations describing granular material. A continuum material model for pellets describing the yield strength, elastic and plastic parameters has been worked out from experimental data and used in the numerical simulations.

It is possible to simulate pellets flow with the SPH method for simple geometries, is concluded in this work. The axisymmetric model agrees well with 3D simulations except close to the symmetry line. Things also have to be improved around the boundary of the problem domain and more validations are necessary for reliable results.

(4)

1. INTRODUCTION

1.1 Background

Forming and handling of iron ore pellets is an important part in the converting process for LKAB. Knowledge about this sub process is very important for further efficiency progress and increasing product quality. The improvement of the processes will be by computational simulations rather then by costly half or full scale experiments in the future. This work is a pre-study for a future research project in pellets simulation. The final goal is to develop a simple tool for predicting the forces on pellets during different types of flow in the production chain.

1.2 Statement of problem

Today there is no established method for simulating pellets flow when new storage systems are designed. Velocity-, pressure- and density distribution are of interest in a moving bed of pellets. The main problem of this work was to create a numerical simulation model for pellets and evaluate it for simulation of silo flow.

The method used is called Smoothed Particle Hydrodynamics, which is a quite new numerical method used for high velocity impact problems and in fluid mechanics. One part of the work was to make a SPH-formulation for granular material. To speed up the calculations an axisymmetric formulation were necessary. The other part was to work out a continuum material model for pellets from experimental data that could be used in the numerical model.

(6)

2. LKAB AND IRON ORE PELLETS

LKAB is a wholly Swedish state owned company and is one of the world’s leading producers of upgraded iron ore products for the steel industry. The production is situated in the northern part of Sweden, with mines in Kiruna and Malmberget. The iron ore are refined into pellets or sinter fines, where pellets are the largest share of the production. Pellets are sintered, centimetre-sized spheres of ore with high iron content and are produced in two verities: blast furnace pellets (BF) and direct reduction pellets (DR). Most of the iron ore products are sold to the European steel mills but they also produce for the North African, the Middle East and Asian market.

2.1 The handling system of iron ore pellets

The finished products from the processing plants in Kiruna and Malmberget are transported to the customers by rail and by ship via the ports at Narvik and Luleå.

The railway is connecting Narvik in Norway via Kiruna, Malmberget with Luleå at the Swedish coast in the Baltic Sea. Iron ore products from Kiruna are transported to Narvik for customers in the European market and the rest of the world. From Malmberget the products are transported to Luleå for customers in the nearby and countries round the Baltic Sea. 23 million tons per year are transported on the railways to the shipping ports. The cars that are currently in operation on the Ore Railway carry a payload of between 80 and 100 tonnes and each train set consists of 52 cars. LKAB are planning for sets with 68 cars, with a capacity of 100 tonnes each. In the port in Luleå ore products are mainly stockpiled in three silos, with a total capacity of 135,000 tonnes. 10 million tonnes of iron ore are passing this port every year. The harbour in Narvik consists of a terminal for discharging the ore trains, stocks of the various iron ore products, and quays where the vessels dock for loading. The capacity of the port is 25 million tonnes per year. A new ore harbour is about to be built with a storage capacity of 1.5 million tonnes. 12 large storage silos and one smaller silo in the form of rock caverns will be blasted. The silos are cylindrical with a diameter of 40m and a height of 60m. Each silo has a capacity of 110 000ton pellets. Above the storage area, the ore trains will enter a tunnel and bottom-discharge their loads into the silos. A belt-conveyor transports the pellets to the ships via a screening station. It was in the planning of these silos knowledge about pellets, and their behaviour in pellets flow, were found to be investigated.

(7)

3. METHOD AND THEORY

3.1 Shear tests on iron ore pellets

Previous to this work, tests on pellets were done at the Division of Soil Mechanics to measure the crushing of pellets, when an axial load and a cycle of shearing are applied. Experimental data from these tests are used in this work to get a material model, describing the yield strength and the elastic and plastic behaviour of pellets as a continuum. The tests were done on a cylindrical shear cell filled with pellets. A shell made of rubber and wires of steel were surrounding the pellets to avoid the cell to expand in radial direction. An axial load was applied at the top of the cell, and a force in horizontal direction sheared the body at the lower part of the cell. This shear force was applied in a cyclic manner, but in this study only the first cycle is studied. During the test the vertical and horizontal deformation was measured. No measurements were done in radial or tangential direction. It is assumed that the shell totally provides the cell to expand in radial direction, therefore the radial strain is set to zero. This assumption is required to get the relation between stresses and strains (Hooke’s law). A picture and a schematic view of the experiment are seen in Figure 1.

Figure 1: The experimental setup and an explanatory sketch of the test.

The height H of the cell varies between 568mm and 679mm. This is the active height during the compression part, but during the shearing a rigid shell and reinforcement bars prevent the upper and lower parts to shear, see Figure 1. The parameters needed in the material model are the shear modulus G, Poisson’s ratio ν and a plastic yield model, describing the yield stress q as a function of the pressure p.

(8)

There were five experiments done at three different axial loads. The tests were done at screened iron ore pellets to avoid influence of fines. The first step in the tests is the axial loading, where the cell is compressed before the shearing. Figure 2 show the result from the compression part of the tests, where the axial stress, σax

is plotted as a function of the axial strain, εax.

Axial Stress - Axial Strain

y = 3,36E+07x + 1,23E+02

y = 3,14E+07x - 1,89E+04 y = 3,80E+07x - 8,76E+03

y = 4,52E+07x - 2,90E+04

-100000 0 100000 200000 300000 400000 500000 600000 700000

0 0,005 0,01 0,015 0,02 0,025

Axial Strain

Axial Stress [Pa]

S1_P1 (Axial s tres s 46kPa) S2_P6 (Axial s tres s 600kPa) S1_P5 (Axial s tres s 601kPa) S2_P8 (Axial s tres s 627kPa) Linear (S1_P1 (Axial s tres s 46kPa)) Linear (S2_P6 (Axial s tres s 600kPa)) Linear (S2_P8 (Axial s tres s 627kPa)) Linear (S1_P5 (Axial s tres s 601kPa))

Figure 2: Data from the compression part. Axial stress as a function of axial strain.

From this test one of the elastic parameters can be worked out from the slopes of the curves, see calculations below. The slope varies from 31MPa to 45MPa. The mean value has been used for the calculations.

Next step is the shearing, where an increasing force is applied to the lower part of the cell. The result from the shearing can be seen in Figure 3, where the shear stress, τxy is plotted as a function of the shear angle, γxy. The shear modulus is given by the slope of the curves. This result is simplified by taking the slope of the secant from origin to the point where the shear stress is 50% of its maximal value, see Figure 3. This value of G is used in the simulations.

(9)

Shear stress - Shear angle

-50000 0 50000 100000 150000 200000 250000 300000

-0,05 0 0,05 0,1 0,15 0,2 0,25 0,3 0,35

Shear angle [rad]

Shear stress [Pa]

S1_P1 (Axial s tres s 46kPa) S2_P7 (Axial s tres s 327kPa) S2_P6 (Axial s tres s 600kPa) S2_P8 (Axial s tres s 627kPa) Approx curve(S1_P1) Approx curve(S2_P7) Approx curve(S2_P6) Approx curve(S2_P8)

Figure 3: Data from the shearing. Shear stress as a function of the shear angle.

Assuming the radial strain to be zero, the axial stress as a function of the axial strain is given by Hooke’s law:

( )

( )( ) ( )

( )

^ax

ax ax

G

E ε

υ ε υ

υ υ

σ υ

2 1

1 2

1 1

1

−

= −

− +

= − (3.1.1)

Setting the slope of the curves, L from the compression test, Figure 2 equal to the slope given in Equation (3.1.1) and solving for ν:

( )

L G

G L L G

2 2

1 1

−

= −

− ⇔

= − υ

υ

υ (3.1.2)

Poisson’s ratio ν and the shear modulus G for the five different tests are seen in Figure 4 and Figure 5.

(10)

Shear modulus - Pressure

y = 6,30E+00x + 1,70E+06

0 1000000 2000000 3000000 4000000 5000000 6000000 7000000

0 50000 100000 150000 200000 250000 300000 350000 400000 450000 500000 p [Pa ]

G [Pa] Secant values

Shear modulus - Pressure Linear (Shear modulus - Pressure)

Figure 4: Shear modulus G as a function of the pressure.

Figure 4 shows an increasing value of G as the pressure increase. In the simulations described later, the pressure is 150kPa and below. A shear modulus of

has been used in these simulations.

2,0MPa G=

Poisson's ratio - Pressure

0 0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 0,45 0,5

0 50000 100000 150000 200000 250000 300000 350000 400000 450000 500000 p [Pa ]

ny Poisson's ratio- Pressure

Figure 5: Poison’s ratio ν for different pressures.

In Figure 5 a decreasing trend for ν can be seen. With only 5 tests it can not be said by sure if it is pressure dependent but the value is between 0,35 and 0,45. In the simulations it was set toυ =0,4.

(11)

The plastic behaviour of the material is described by Drucker Prager yield function without hardening, see [3]. This model is used for pressure sensitive granular materials for example metal powder, sand and in this case iron ore pellets. The yield function is a function of the pressure and the second stress deviator invariant,f(p,J₂ )=0 and is written:

=0

−

= q

f σ_vm (3.1.3)

Where q is a liner function of p:

2 1p c c

q= + (3.1.4)

σvm can be written in terms of the second stress deviator invariant J2:

2 2

3 2

1 2 3

J J

vm ij

ij ij ij

vm ⇒ =

⎪⎪

⎩

⎪⎪⎨

⎧

′

= ′

′

= ′

σ σ

σ σ σ σ

(3.1.5)

The yield function in the final form is written:

0 3 ₂ − ₁ − ₂ =

= J c p c

f (3.1.6)

When f = 0 the strains have reached the yield surfaces, the material have started to flow and the strains have to be scaled back under the yield surface, see Equation (3.2.27).

As can be seen, the yield function is a function of the von Mieses stress σvm and the pressure p. These identities can be coupled to the measured strains with the following equations, taken from [2], starting with Hooke’s law.

The relation between stresses and strains are given by:

(

_⎥⎦^⎤

⎢⎣⎡ + +

+ −

= + _x _x _y _z

x

E ε ε ε

υ ε ν

σ υ

2 1

1

)

(3.1.7)

(

_⎥⎦^⎤

⎢⎣⎡ + +

+ −

= + _y _x _y _z

y

E ε ε ε

υ ε ν

σ υ

2 1

1

)

^(3.1.8)

(

_⎥⎦^⎤

⎢⎣⎡ + +

+ −

= + _z _x _y _z

z

E ε ε ε

υ ε ν

σ υ

2 1

1

)

^(3.1.9)

(12)

In cylindrical coordinates σ_x =σ_axandσ_y =σ_z =σ_rad. With the assumption

=0

εrad the relation between stresses and strains are given by:

( )

( )( ) ( )

( )

^ax

ax ax

G

E ε

υ ε υ

υ υ

σ υ

2 1

1 2

1 1

1

−

= −

− +

= − (3.1.10)

( )( )

^ax

( )

^ax

rad

G

E ε

υ ε υ

υ υ

σ υ

2 1 2

1

1 = −

−

= + (3.1.11)

To reduce terms the axial and radial stresses are related by:

υ υ σ

κ σ

= −

= _ax 1

rad (3.1.12)

The yield stress q is equal to the von Mieses stress σvm which is given by:

[

x² ²y ²z x y y z z x ³ xy² ³ ²yz ³ zx²

]

¹²

q=σvm = σ +σ +σ −σ σ −σ σ −σ σ + τ + τ + τ (3.1.13) In cylindrical coordinates with τyz = τzx = 0:

[

ax² ² ax² ² ax² ax² ² ax² ax² ³ xy²

]

¹²

q= σ +κ σ +κ σ −κσ −κ σ −κσ + τ (3.1.14)

(

⇔

)

[

ax² ¹ ² ² ³ ²xy

]

¹²

q= σ +κ − κ + τ (3.1.15)

The pressure p is the trace of the stress tensor:

3

z y

p σx +σ +σ

= (3.1.16)

In cylindrical coordinates and with κ as formulated in Equation (3.1.12):

( )

3 2

1 κ

σ +

= ^ax

p (3.1.17)

Inserting the values for ν, σax and τxy from the five tests, q as a function of p can be plotted, see Figure 6:

(13)

Yield stress - Pressure

y = 1,0491x + 28523

0 100000 200000 300000 400000 500000 600000

-100000 0 100000 200000 300000 400000 500000

p [Pa]

q [Pa] Yield s tres s - Pres s ure

Linear (Yield s tres s - Pres s ure)

Figure 6: The Drucker Prager yield function.

In the figure the slope is the constant c1 in Equation (3.1.4) and c2 is given by the intersection of the curve with the y-axis. The slope is called the inner friction of the material and was found to be 1,05, or an angle of 46 degrees. The constant c2

is here set to zero, which means that the pellets continuum is assumed to have no tensile strength.

3.2 Theory of smoothed particle hydrodynamics

This section will explain the basic ideas of SPH, for more detailed information, see [1]. The Smoothed Particle Hydrodynamics (SPH) method was invented independently by Lucy; Gingold and Monaghan, 1977, to solve astrophysical problems in open space. It is a mesh free, particle method for modeling fluid flows, and has also been extended to solve problems with material strength.

Today, the SPH method is being used in many areas such as fluid mechanics (for example: free surface flow, incompressible flow, and compressible flow), solid mechanics (for example: high velocity impact and penetration problems) and high explosive detonation over and under water.

The difference between SPH and grid based methods as Finite Element Method (FEM), is that in the SPH method the problem domain is represented by a set of particles or points instead of a grid. Besides representing the problem domain, the points also act as the computational frame for the field approximation. Each point is given a volume and density, and carries information about spatial coordinate,

(14)

velocity, and internal energy. Other quantities as stresses and strains are derived from constitutive relations. SPH is an adaptive Lagrangian method, which means that in every time step the field approximation is done based on the current local set of distributed points. Another difference from FEM is that the points are free to move from each other in action of the internal and external forces, there are no connections between them. Virtual points are used to describe walls at the boundaries. These points do not have any velocities, but masses and densities equal to the real points. This and the adaptive nature of SPH approximation result in a method that can handle extremely large deformation. This is what makes it attractive for this kind of problems.

3.2.1 The calculation cycle in SPH

The basic idea for a numerical method is to reduce the partial differential equations (PDEs) that describe the field functions (for example density, velocity and internal energy) to a set of ordinary differential equations (ODEs), with respect to time only. These equations can easily be solved with some standard integration routine. With SPH this is done with following key-steps:

1. The problem domain is represented by an arbitrary distributed set of non connected points. (Mesh free)

2. Each field function is described with an integral function. (Kernel approximation)

3. The kernel approximation is then further approximated using the points.

This is called the particle approximation in the SPH method. The integral is replaced with a summation over the neighboring points to the computational point. This is done to each point.

4. The particle approximation is done at every time step, based on the local distribution of points.

5. The particle approximations are performed to all terms of field functions and reduce the PDEs to ODEs with respect to time only. (Lagrangian) 6. The ODEs are solved using an explicit integration algorithm. The time

history for all the field variables for all particles can be saved.

7. Other quantities are derived from the field variables using constitutive relations.

3.2.2 Integral representation of a function

In the first step an arbitrary field function is represented in integral form by a multiplication of the field function and a smoothing kernel function. This is called the kernel approximation. It starts from the following identity:

(3.2.1)

Ω

∫

′

− ′

= x′ x x x

x f d

f( ) ( )δ( )

(15)

Where f is a field function of the three-dimensional position vector x, and δ(x-x´) is the dirac delta function. Ω is the volume of the integral that contains x. So far the integral representation of the function is exact, as long as f(x) is defined and continuous in Ω.

3.2.3 The kernel function

Now the dirac delta function δ(x-x´) is replaced with a smoothing kernel function W(x-x´,h), the integral representation is given by:

(3.2.2)

Ω

∫

′

− ′

>= ′

< f(x) f(x)W(x x,h)dx

In the smoothing function, h is the smoothing length, defining the influence area of the smoothing function W. The integral representation is an approximation of the field function as long as W is not the dirac delta function. This is called the kernel approximation. The kernel function should satisfy some conditions, stated below.

The normalization condition:

(3.2.3) 1

) ,

( − ′ ′=

∫

Ω

x x

x h d

W

The Delta function property:

) ( ) , (

lim0 x−x′ = x−x′

→ W h δ

h (3.2.4)

The compact condition:

0 ) ,

( − ′ h =

W x x When x− x′ >κh (3.2.5)

Where κ defines the support domain (non zero area) of the smoothing function.

Positivity:

0 ) , ( − ′ h >

W x x For any x´ (3.2.6)

Except from these conditions the smoothing function should decrease with the increase of distance from the evaluation point, be an even function and sufficiently smooth.

(16)

The smoothing function is quite important in the SPH method, because it determines the accuracy of kernel and particle approximations. That is why many different smoothing functions have been used in the published literature in the area. In this work a higher order quintic spline kernel is used. This smoothing function is sufficiently smooth even for higher order derivatives and therefore more stable than lower order functions. The equation of the kernel function is seen in Equation (3.2.7).

⎪⎪

⎩

⎪⎪

⎨

⎧

−

− +

−

= 0

) 3 (

) 2 ( 6 ) 3 (

) 1 ( 15 ) 2 ( 6 ) 3 ( )

,

( ₅

5 5

5

R

R R

R h

R

W α_d

3 3 2

2 1

1 0

>

<

≤

<

≤

<

≤

R R R

R

(3.2.7)

Where R is the relative distance between to points x and x´,

R x−hx´

= and αd is a constant equal to120/h, 7/478πh, in one-, two- and three dimensional space respectively.

359 3

/

3 πh

Figure 7: Quintic kernel function.

The integral representation of the derivative of a function ∇f(x) is given by the following equation:

(3.2.8)

Ω

∫

′

− ′

′ ∇

−

>=

∇

< f(x) f(x) W(x x,h)dx

(17)

From the equation, it can be seen that the differential operation on the field function is transmitted to a differential of the kernel function. In other words, instead of deriving the derivative of the function from the derivative of the function itself, the spatial gradient is determined from the values of the function and the derivatives of the smoothing function. This is the first step to reduce the PDEs to ODEs.

3.2.4 Particle approximation

The next step in getting a numerical solution for the field functions is the particle approximation. The entire system is represented by a finite number of points that carry individual mass and density. The continuous integral representation is replaced with a summation over the points in the support domain, to get a representation in discrete form. The infinitesimal volume dx´ is replaced with the finite volume of the point ΔV_j =m_j /ρ_j as follows:

) , ( ) ( )

, ( ) ( )

(

1

h W

m f d

h W

f

f _j

N

j

j j

j x x x

x x x x

x >= ′ − ′ ′≅ −

<

∫ ∑

₌

Ω ρ (3.2.9)

That can be written:

ij N

j

j j j

i m f W

f( ) ( )

∑

1

=

>=

< x x

ρ Where W_ij =W(x_i −x_j,h) (3.2.10) The spatial derivative of a function in discrete form:

ij i N

j

j j j

i m f W

f >= ∇

∇

<

∑

=

) ( )

(

1

x

x ρ Where

ij ij ij

ij ij

ij j i

j i ij

i r

W r r W W

∂

= ∂

∂

−

= −

∇ x

x x

x

x (3.2.11)

It can be shown, see [1] that there is other ways to write the spatial derivative of the field functions at discrete form that is more common to use and which have been used in this thesis:

[

j i i ij

N

j j i

i m f f W

f >= + ∇

∇

<

∑

=

) ( ) 1 (

) (

1

x x

x ρ

]

^(3.2.12)

ij i i

i j

j N

j j i

i f f W

m

f ∇

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡ +

>=

∇

<

∑

= 2 2

1

) ) (

( ( )

(x ρ ρ^x ρ^x (3.2.13)

(18)

Figure 8: Particle approximation with support domain and kernel function.

In Figure 8 the support domain for calculation point i, is shown. Summation of neighbouring particles with influence given by the kernel function gives the searched value for the point i.

3.2.5 Hydrodynamics with material strength

The governing equations describing the density and the internal force in solids with material strength are the conservation equations of continuum mechanics, where the nodal velocities v^α is the primary variables for the calculations:

⎪⎪

⎪

⎩

⎪⎪

⎪

⎨

⎧

=

∂

= ∂

∂

− ∂

=

α α

β αβ α

β β

σ ρ ρ ρ

x v

Dt D

1 (3.2.14)

Where the density ρ, the velocity component v^α and the total stress tensor σ^αβ are the dependent variables. The use of particle approximation converts the equation to discrete form:

(19)

⎪⎪

⎩

⎪⎪

⎨

⎧

=

∂ + ∂

−

=

∂

− ∂

=

∑

=

α α

β αβ αβ

α

β β β

ρ σ ρ σ ρ

i i

N

j i

ij j j i i j i

N

j i

ij j i j i

Dt D

m W Dt

D

m W Dt

D

x v

v x v

1 2 2

1

) (

(3.2.15)

An alternative for the density calculation in Equation (3.2.15), is to use the SPH- formulation directly, called the summation density approach:

(3.2.16)

ij N

j j

i

∑

mW

=

1

ρ

To be able to describe the behavior of the stress tensor σ^αβ, a constitutive model has to be included in the program for the SPH method, see next section.

3.2.6 Constitutive modeling

The constitutive relation described below is taken from [1] and [3]. To be able to describe the behavior of a material, it is necessary to define a material model. In solid mechanics, in general, the stress σ is a function of strain ε and strain rate ε_&

(Hook’s law). The total stress tensor σ^αβ is made up by two parts, one part of the isotropic pressure p and the other part of the shear stress τ.

(3.2.17)

αβ αβ

αβ τ δ

σ = − p

In order to get a frame indifferent strain rate, the Jaumann rate is used for the shear stress calculation, the shear stress rate is given by:

γβ αγ γβ αγ αβ

αβ ε τ τ

τ& = 2G& +R& − R& (3.2.18)

Where G is the shear modulus, ε^&^αβis the strain rate deviator, R&^αβ is the rotation rate tensor. The strain rate tensor ε_&^αβ is a function of the velocity gradient:

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∂ +∂

∂

= ∂ _β^α _α^β εαβ

x v x v 2

& 1 (3.2.19)

(20)

And the rotation rate tensor is given by:

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∂

−∂

∂

= ∂ _β^α _α^β

αβ

x v x v 2

R& 1 (3.2.20)

The traceless part of the stress tensor is called the strain rate deviator, given by:

ε&αβ

γγ αβ αβ

αβ ε δ ε

ε& & &

3

−1

= (3.2.21)

To get the new value of the shear stress in each time step the incremental change in shear stress is multiplied with the time step and added to the previous value.

(3.2.22)

t dt

dt

t⁺ ταβ=ταβ+τ_&αβ

The isotropic pressure p is a function of the volumetric strainε_v multiplied by the bulk modulus.

(3.2.23)

v dt t dt

t⁺ p=−K ⁺ ε

The volumetric strain rate is a summation of the strain rates:

z y x

v ε ε ε

ε_& = _& + _& + _& (3.2.24)

In same way as for the shear stress the volumetric strain rate is calculated:

(3.2.25)

vdt

v t v dt

t⁺ ε =ε +ε_&

The provisional von Mises flow stress J=3J2 is compared with the yield strength q of the material in order to see if the shear stresses have to be scaled back to the yield surface. For steel q is a constant, but for other materials as granular material, that is the case in this study, q is a function of the isotropic pressure, see Section 3.1. Equation (3.2.27) is showing the downscaling of τ^αβ.

αβ αβτ 2τ 1

2 =

J (3.2.26)

2

2 3

3

J q q

J > ⇒τ^αβ =τ^αβ (3.2.27)

(21)

3.2.7 SPH formulation for hydrodynamics with material strength

3.2.7.1 3D FORMULATION

Using the particle approximation the following equations are to be solved in three dimensional Cartesian space.

The governing equations in 3D SPH formulation:

⎪⎪

⎩

⎪⎪

⎨

⎧

=

∂ + ∂

−

=

∂

− ∂

=

∑

=

α α

β αβ αβ

α

β β β

ρ σ ρ σ ρ

i i

N

j i

ij j j i i j i

N

j i

ij j i j i

Dt D

m W Dt

D

m W Dt

D

x v

v x v

1 2 2

1

) (

(3.2.28)

Or with the summation density approach for density calculation:

(3.2.29)

ij N

j j

i

∑

mW

=

1

ρ

The strain rate tensor in 3D SPH formulation:

⎪⎪

⎩

⎪⎪

⎨

⎧

⎥⎦

⎢ ⎤

⎣

⎡

∂

− ∂

∂ +

− ∂

=

⎥⎦

⎢ ⎤

⎣

⎡

∂

− ∂

∂ +

− ∂

=

⎥⎦

⎢ ⎤

⎣

⎡

∂

− ∂

∂ +

− ∂

=

∂

− ∂

=

∂

− ∂

=

∂

− ∂

=

∑

=

N

j i

ij i j i

ij i j j

i zx i xz i

N

j i

ij i j i

ij i j j

i zy i yz i

N

j i

ij i j i

ij i j j

i yx i xy i

i ij i j N

j j i

z i

i ij i j N

j j i

y i

i ij i j N

j j i

x i

x w W z w

u W u m

y w W z w

v W v m

x v W y v

u W u m

z w W w m

y v W v m

x u W u m

1 1 1 1

1 1

) (

) 2 (

1 1

) (

) 2 (

1 1

) (

) 2 (

1 1

) 1 (

ε ρ ε

ε ρ ε ρ ε ρ

&

(3.2.30)

(22)

The rotation rate tensor in 3D SPH formulation:

⎪⎪

⎩

⎪⎪

⎨

⎧

−

=

⎥⎦

⎢ ⎤

⎣

⎡

∂

− ∂

∂ −

− ∂

=

−

=

⎥⎦

⎢ ⎤

⎣

⎡

∂

− ∂

∂ −

− ∂

=

−

=

⎥⎦

⎢ ⎤

⎣

⎡

∂

− ∂

∂ −

− ∂

=

∑

=

xz i zx i

N

j i

ij i j i ij i j j i xz i

yz i zy i

N

j i

ij i j i

ij i j j i yz i

xy i yx i

N

j i

ij i j i ij i j j i xy i

z i y i x i

R R

x w W z w

u W u m R

R R

y w W z w

v W v m R

R R

x v W y v

u W u m R

R R R

&

1 1 1

) (

) 2 (

1 1

) (

) 2 (

1 1

) (

) 2 (

1 1

0

ρ ρ ρ

(3.2.31)

Where u, v and w are the velocity components in x-, y-, z-direction respectively.

Many problems can easily be cast into a three dimensional rectangular Cartesian SPH form. The computational time for such a formulation increases dramatically with the number of approximation points. To get through this some problems can be formulated in a two dimensional axisymmetric SPH formulation. The following axisymmetric formulation has been tested in this thesis, for details see [4] and [5]. The masses remain constant throughout the computation, and are obtained from M=ρ0V0 where ρ0 and V0 represent the initial density of the material and the initial volume represented by the point, Vi = 2πridrdz.

(23)

3.2.7.2 2D AXISYMMETRIC FORMULATION

Governing equations in 2D axisymmetric SPH form:

( )

⎪⎪

⎩

⎪⎪

⎨

⎧

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

∂

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ +

∂ +

∂

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ +

−

∂ =

∂

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

∂

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ +

∂ +

∂

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ +

−

∂ =

∂

− ∂

=

∑

=

N

j i

ij j j

rz j i u

rz i i

ij j j

r j i

u r i j i

i i i

N

j i

ij j j

z j i u

z i i

ij j j

rz j i u

rz i j i

N

j i

ij j i j i i

z W r

r r W r

m r r

t r

z W r

r r W r

m r t

z

m W r Dt

D

1 2 2 2 2

1

2 1 2

1 2

1

ρ σ ρ σ ρ

σ ρ σ π

ρ σ

ρ σ ρ σ ρ

σ ρ σ π

π ρ

θθ

β β β

&

v x v

(3.2.32)

And the summation density approach in axisymmetric SPH form:

ij N

j j i

i m W

r

∑

=

= 2 1

1

ρ π (3.2.33)

The constitutive equations in 2D axisymmetric SPH form:

⎪⎪

⎩

⎪⎪

⎨

⎧

=

−

=

⎥⎦

⎢ ⎤

⎣

⎡

∂

− ∂

∂ −

− ∂

=

⎥⎦

⎢ ⎤

⎣

⎡

∂

− ∂

∂ +

− ∂

=

∂

− ∂

=

∂

− ∂

=

∑

=

0 2

) (

) 2 (

1 1

2 )

( )

2 ( 1 1

2 /

2 )

1 (

2 )

1 (

1 1 1

1 1

z i r i z i r i i z i r i z i r i z i r i

zr i rz i

j N

j i

ij i j i

ij i j j i zr i

j N

j i

ij i j i

ij i j j i zr i rz i

j ij j N j

j j

j i

ij i j N

j j i z i

j i

ij i j N

j j i r i

R R R R R R R R

R

r r v W z v

u W u m R

r r v W z v

u W u m

r r W

u m

z r v W v m

r r u W u m

θ θ θ θ θ θ

θ θ θ θ

ε ε ε ε

ρ π ρ π ε ε

ρ π ε

&

(3.2.34)

(24)

4. RESULT AND DISCUSSION

4.1 Numerical simulations

The simulations were done in Compaq Visual Fortran and the results were revised in Matlab 7.0. A SPH-code from [1], written for fluids were rewritten with equations describing granular material flow, see Section 3.1 and Section 3.2. The material parameters used in the simulations were given from the shear test

experiments. The bulk density was set to ; shear

modulus , Poisson’s ratio

= 2260kg/m3

ρ 2,0MPa

G= υ = 0,4 and failure function

0 05 , 1

3 ₂ − =

= J p

f . To validate the axisymmetric model a three dimensional model with an outlet of 1,8m was created. To find the critical time step the Courant criterion was used:

L E t

t≤Δ _critical = −υ ρ

Δ _min 1 (4.1.1)

In Equation (3.3.1), Lmin is the smallest length between two calculation points and E is Young’s modulusE=2G(1-υ). The computer used for the simulations had a 2,8GHz processor with 504Mb of RAM.

4.1.1 3D silo model

In the three dimensional silo model 27 426 points and 20 040 virtual points were used. This meant 15cm between each node and a time step of 0,00075s. The Quintic kernel function was used and the smoothing length was set twice the distance between two points. A plot of the 3D model with half of the virtual points after 1.0s discharging time is seen in Figure 9.

(25)

Figure 9: 3D silo model after 1,0s discharging.

The time for emptying the silo with this geometry is around 10s and the simulation time was 96 hours. By taking a cut in the middle of the silo the emptying process can be followed in more detail, see Figure 10.

0,0s 1,5s 3,0s

4,5s 6,0s 10,0s

Figure 10: Cross section of the 3D model at 6 different times of discharging.

(26)

4.1.2 2D axisymmetric model

In the axisymmetric two dimensional silo model 4 830points and 302virtial points were used. This meant 5cm between each node and a time step of 0,00025s. The Quintic kernel function was used and the smoothing length was set twice the distance between two points. A plot of the 2D model after 1.0s of discharging is seen in Figure 11.

Figure 11: 2D axisymmetric model after 1,0s discharging.

The simulation time for this geometry was around 50minutes. By mirroring the geometry in the axis of symmetry, the emptying process can be compared with the three dimensional case, see Figure 12. A comparison of Figure 10 and Figure 12 shows similar lapse of discharging. It is in the middle of the silo differences can be seen, se explanation below. The axisymmetric model is based on calculation on the half cross section of the silo, se Figure 11; therefore the closeness of the points is to be higher than the case for the 3D model. This in order to resolve the outlet that is described by half of the size of the 3D model outlet.

(27)

-5 0 5 -9

-8 -7 -6 -5 -4 -3 -2 -1 0

-5 0 5

-9 -8 -7 -6 -5 -4 -3 -2 -1 0

0,0s 1,5s 3,0s

Figure 12: Axisymmetric model at 6 different times of discharging.

As can be seen in Figure 12 above the points in the middle of the silo is thinning out as the simulation steps forward. This happens when points from the sides moves towards the centre. A point in the axisymmetric formulation represents a torus ring, with a volume and a mass that is kept during the simulation. When the point moves from the side to the centre, the radius of the ring will decrease. To keep the volume and mass of the ring constant, the cross section area will increase and the result is that the points moves away from each other, as they moves to the centre, see Figure 13.

4,5s 6,0s 10,0s

Figure 13: Torus rings with different diameters but the same volume.

(28)

4.1.2.1 VELOCITY FIELD

In Figure 14 below the velocity field after 1,0s of discharging is plotted. It is a typical funnel flow where the material at the sides stands still and the silo is emptied from the top and downwards. The magnitude of the velocities is seen in right Figure 14.

-5 -4 -3 -2 -1 0 1 2 3 4 5

-9 -8 -7 -6 -5 -4 -3 -2 -1 0

Figure 14: Left, velocity field. Right, magnitude of the velocities.

4.1.2.2 PRESSURE

In Figure 15 the pressure is plotted at the same time step as for the velocities. It is zero pressure in the free flow and higher pressures at the sides which seems quite reasonable, but the solution is unstable and varies a lot between neighbouring points. The solution is also very oscillating that can be seen in the right Figure 15, where the pressure at the bottom of the silo is plotted while the silo is emptied.

0 10 20 30 40 50 60 70 80 90 100

0 1 2 3 4 5 6 7x 10⁴

Figure 15: Left, pressure distribution. Right, pressure at the bottom as a function of the discharging time.

(29)

4.1.2.3 DENSITY

In the density plot in Figure 16 it can be seen an increased density in the bottom of the silo. Some problems at the boundaries can be observed. At the symmetry line a few point have higher values, and at the rest of the boundary the density is too low. The problem around the symmetry line can be because of the division with the radius in the axisymmetric formulation that cause numerical problems when the radius approach zero. The other density problem is because of the SPH- density calculations, as seen in Equation (3.2.33) the density is calculated by summing the masses of the neighbouring points multiplied with the kernel function. For a fully surrounded node a summation of the kernel is equal to one, see Equation (3.2.3). At the boundary this is not the case and the summation will be less than one (0,5 at the boundary points), resulting in a lower density. The truncation of the kernel function is seen in Figure 17.

Figure 16: Density distribution.

Figure 17: The kernel function and truncation of the boundary.

(30)

4.1.3 Simulation of the experimental silo

For validation of the numerical model, simulations were done on geometry from an experimental silo in Edinburgh. Available data from the silo is the discharging time and the flow pattern. Results and conclusions from the experiments are collected in [6].

4.1.3.1 SILO GEOMETRY

The experimental silo has a cylindrical shape with a filling height of 6,4m and a width of 4,2m. Experiments were done on three different outlets, were the test with a concentric outlet were of interest in this work. It is a circular outlet with a diameter of 480mm. Unfortunately the slide valve was just opened 110mm during the discharging. A sketch of the silo and the outlets are seen in Figure 18.

Figure 18: Experimental silo with dimensions and outlets.

(31)

Because of the small outlet, it would be needed around 10 million points for a three dimensional simulation of the silo which is far too much and an axisymmetric simulation is necessary. For such a simulation the geometry has to be cylindrical symmetric. The silo itself fulfils this requirement, but due to the not fully opened outlet a correct description with an axisymmetric model can not be done. Therefore the simulations were done on a silo with a circular outlet of 480mm (fully opened). Due to the different shapes of the outlets the discharge times of the simulation and experiment can not be compared, but the flow pattern should have a minor effect [6].

4.1.3.2 SIMULATION

The maximal numbers of points available on the computer used was 85 344 real- and 1324 virtual points. With such representation the distance between the points isdx=1,25cm and time stepdt =0,000075s. Figure 19 shows the lapse of the simulation with 6 different level markers.

0,0s 6,0s

Figure 19: Axisymmetric model of the experimental silo with six level markers at two different time steps.

(32)

Figure 20: Flow pattern from the experimental silo.

Figure 20 shows the lapse of the experimental silo discharging. 6 different layers of level markers were located in the silo in order to follow the discharging process. A comparison of Figure 19 and Figure 20 shows similarities in the flow pattern between the experiment and simulation. It would be better to simulate the complete discharging of the silo for a better verification, but after 102 000 time steps (7,5s of discharging) the simulation breaks because the problem domain gets to large. To get through this a box should be placed under the silo in order to collect the departure points. The disadvantage with such a solution is the longer simulation time. With a free outlet the simulation gets faster and faster as the points exit the silo.

The time for the simulation above was 46 hours (6,0s of discharging). From experiment a complete emptying of the silo will proceed in 1,5 hours with a not fully opened outlet. Because of the explicit method, the time steps are too short for simulating such a slow process. Even though the simulation above has a fully opened outlet it would take months to simulate the complete emptying of the silo.

At the outlet about 20 points are needed to resolve the geometry. Due to the small outlet, this result in a lot of points in the total silo because the point density is the same in the whole problem domain. With a larger outlet fewer points would be needed and the process would proceed faster, resulting in a shorter simulation time.

Simulation and modelling of pellets in large storing systems

M A S T E R’S T H E S I S

GUSTAF GUSTAFSSON

Simulation and

Modelling of Pellets in Large Storing Systems

PREFACE

ABSTRACT

TABLE OF CONTENTS

1. INTRODUCTION

2. LKAB AND IRON ORE PELLETS

3. METHOD AND THEORY

( )

( )( ) ( )

( )

( )

( )

(

)

(

)

(

)

( )

( )( ) ( )

( )

( )( )

( )

[

]

[

]

(

)

[

]

( )

∫

∫

∫

∫

∫ ∑

∑

∑

[

∑

]

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

( )

∑

∑

∑

∑

∑

∑

∑

∑

∑

4. RESULT AND DISCUSSION