Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

53, 2, pp. 487-500, Warsaw 2015
DOI: 10.15632/jtam-pl.53.2.487

APPLICATION OF VARIATIONAL AND FEM METHODS TO

THE MODELLING AND NUMERICAL ANALYSIS OF THE SLITTING

PROCESS FOR GEOMETRICAL AND PHYSICAL NONLINEARITY

Łukasz Bohdal, Leon Kukiełka

Koszalin University of Technology, Koszalin, Poland

e-mail: bohdall@interia.pl; leon.kukielka@tu.koszalin.pl

Finite elementmodelling provides a great deal of support in the understanding of technolo-
gical processes. However, there are few studies of the slitting process, and those that exist
are simplified for use only in the calculation of steady states of such processes. This paper
proposes the application of variational and finite elementmethods for the analysis of slitting
and the nonlinearities of this process. Physical andmathematical models of the process and
a new thermo-elastic/thermo-visco-plasticmaterial model are elaborated. The procedure is
implemented in the finite element code ANSYS/LS-DYNA and themodel is validated com-
paring the numerical and experimental results. The influence of various process conditions
on the strain and stress states and the quality of the final product are analysed. The re-
sults lay the groundwork for further study regarding the numerical analysis of spring-back
behaviour and the effect of tool elasticity on the quality of the final workpiece.

Keywords: slitting, constitutive laws, modelling, numerical simulation, variational formula-
tion, FEM

1. Introduction

Slitting is a sheetmetal-,metal alloy- ormagnetic tapes-cuttingprocess that uses circular knives.
Industries use slitting to splitwide coiled sheets into narrowerwidths or for the edge trimmingof
rolled sheets (Lu et al., 2001; Liu et al., 2005).This process is schematically drawn inFig. 1a.The
sheet is first plastically deformed (A-B). Continued deformation leads to ductile fracture (B-C)
and complete separation of the parts (D) (Wisselink, 2000; Gałęzia et al., 2012). The quality of
the part after slitting depends onmanyprocess parameters such as clearance between the knives
in the horizontal and vertical direction, lubrication method, shearing velocity, knife geometry
and themethod of clamping the sheet plates (Meehan andBurns, 1998; Ma et al., 2006). These
parameters are key factors affecting energy consumption and physical phenomena occurring in
the cutting process such as the state of material displacement, visco-plastic strain and stress.
In turn, certain physical phenomena affect the mechanisms of separation of the material in
the cutting zone, such as visco-plastic shear or rupture, because the mechanism of separation
depends on the quality of both the sheared edge and the product. The amount of adjustable
process parameters and the fact that the influence of these parameters on the process are not
fully understood,makes it difficult to control the slitting process. In practice, the right setup for
the slitters ismostly foundby a trial and errormethod combinedwith the experience.Therefore,
the final product frequently has serious defects and drawbacks such as large deformation and
defects on the sheared edge, such as burrs and edge wave (Fig. 1b).
Computational models, such as the finite element method (FEM), are valuable in reducing

the number of trial-and-error experiments required to predict the state ofmaterial displacement,
residual stresses, strains, material fracture, sheet deformations and quality of the sheared edge.



488 Ł. Bohdal, L. Kukiełka

Fig. 1. (a) Scheme of the slitting process; (b) typical defects in the sheet sheared edge

Some of the researches have studied the process using FEM. Wisselink and Huetink (2004)
developed a finite element model of slitting to calculate the steady state of such a process
using theArbitrary Lagrangian Eulerian formulation. An isotropic plastic constitutive equation
combined with a simple uncoupled damagemodel was used to describe thematerial behaviour.
This model was based on stress triaxiality, which was defined as the ratio of the hydrostatic
stress to the equivalent stress. The simulations performed as a part of the study completely
neglected the effect of damage on the material elastic and plastic behaviour and the effects
of strain rate and temperature-dependent material properties. The developed FE model was
used to simulate residual sheet stresses and deformation, though. Meehan and Burns (1998)
studied the slitting using photoelastic micrographs and made observations of two-dimensional
stress distributions. Aggarwal et al. (2005) conducted a finite element analysis of the magnetic
tape slitting process. The polymeric tape substrate was modeled using a pressure dependent
yield function and the material separation was incorporated using the shear failure criterion.
Saanouni (2006) coupled Continuum Damage Mechanics with an FE model describing elasto-
plastic behaviour to numerically simulate the slitting of a metallic sheet. An implicit Newton-
Raphson scheme, together with a dynamic explicit resolution strategy, was used to analyse the
effects of process configurations on the state of stress in the sheet.
Analysis of the current literature suggests that the main difficulty in the modelling of the

slitting is that only a limitednumberofFEMmodels are currently capable of describing the com-
plete shearing process, including the complete separation of the material parts through ductile
fracturing. A model describing the primary physical phenomena and characterising the mecha-
nical behaviour of the metal sheet is required to numerically simulate this process using finite
element analysis. If the modelling used the updated Lagrangian formulation, FEM considering
processnonlinearity (geometrical, physical, and thermal), largedeformations, strain rate, friction
and nonlinear material characteristics could be used to analyse the thermo-visco-plastic flow of
the material and its stress and strain fields at any time during the slitting process, successfully
numerically reproducing the operation. The slitting differs from the other shearing operations
such as blanking and punching as it involves the rotation of blades and, hence, the material is
cut in two directions simultaneously instead of one. A thorough understanding of the process
of slitting would thus require a three-dimensional analysis of the process to be conducted. The
three-dimensional analysis would allow a complete study of the process enabling the considera-
tion of parameters such as blade geometry, knife radius, shearing velocity, etc., which cannot be
studied using a simplified two-dimensional model.
In this paper, the slitting process is considered as an initial andboundaryvalue problemwith

geometrical, physical and thermal nonlinearity andwith only partial knowledge of the boundary
conditions. The boundary conditions in the tools-sheet contact zones are not known. The paper



Application of variational and FEM methods... 489

focuses on the development of a methodology to study the phenomena that occur during and
after the cutting process (e.g., the change from the steady state to the unsteady state, sheet
deformation, stress and strain states and the sheared edge characteristic features). Physical
andmathematical models of the slitting process and a new thermo-elastic/thermo-visco-plastic
material model are elaborated. The updated Lagrange description is used to describe nonlinear
phenomena on a typical incremental step using a stepwise and co-rotational coordinate system.
The strain and strain rate states are described using Green-Lagrange strain and strain rate
tensors, respectively, and are nonlinear dependences without any linearisation. The procedure
is employed in an application created by the authors usingANSYS/LS-DYNAprogram and the
finite element method. The study presents sample 3D simulations and experimentally verifies
and analyses the impact of various process conditions on the material displacements, stresses,
strains and quality of the final product.

2. Thermo-elastic and thermo-visco-plastic material model

The goal of this Section is to present a material model that considers the combined effects
of thermo-elasticity (in the reversible zone) and thermo-visco-plasticity (in the non-reversible
zone). The model accounts for the strain, strain rate and material temperature histories. The
formulation of the model assumes that the typical small incremental components ∆εij of the

strain tensor T∆ε can be expressed as the sum of thermo-elastic ∆ε
(TE)
ij , visco-plastic ∆ε

(VP)
ij

and thermal∆ε
(TH)
ij incremental strains

∆εij =∆ε
(TE)
ij +∆ε

(VP)
ij +∆ε

(TH)
ij (2.1)

where:∆εij,∆ε
(TE)
ij ,∆ε

(VP)
ij and∆ε

(TH)
ij are the incremental components of the total, thermo-

elastic, visco-plastic and thermal strain tensors, respectively.
The constitutive law for an isotropic, thermo-elastic material with a temperature-dependent

modulus is

∆σij =C
(TE)
ijkl

(
∆εkl−∆ε

(VP)
kl
−∆ε

(TH)
kl

)
+∆C

(TE)
ijkl
ε
(E)
kl

(2.2)

where: ∆σij are the incremental components of the second Piola-Kirchhoff stress, ∆εij are

the incremental components of the Green-Lagrange strain tensor, C
(TE)
ijkl
(T) = λ(T)δijδkl +

µ(T)(δikδjl+δilδjk) is the component of the temperature-dependent elastic constitutive tensor,

∆C
(TE)
ijkl
(∆T)= [∂C

(TE)
ijkl
(T)/∂T ]∆T are its increments

λ(T)=
E(T)ν(T)

[1+ν(T)][1−2ν(T)]
µ(T)=

E(T)
2[1+ν(T)]

E(T) and ν(T) are temperature-dependent Young’s modulus and Poisson’s ratio, respectively,
δij is the Kronecker delta and ε

(E)
kl
are the accumulated components of the elastic strain tensor

at the time t. The components of the thermal increment strain tensor are calculated from the
equation

∆ε
(TH)
ij
∼=αm(T)∆Tδij (2.3)

where αm is the mean coefficient of thermal expansion. The visco-plastic incremental strain is
defined as (Kleiber, 1985)

∆ε
(VP)
ij =∆λ

∂F

∂S̃ij
(2.4)



490 Ł. Bohdal, L. Kukiełka

where ∆λ is the positive scalar variable or Lagrange’s coefficient and S̃ij = Sij −αij is the
component of the reduced stress deviator D̃σ.
The calculation of∆λ requires a hardening rule to be selected. The hardening rule describes

the change in the yield surface with continuing plastic deformation.
For mixed hardening and an isotropic body, the rule is defined by

∆λ=
S̃ijC

(TE)
ijkl

(
∆ε
(TE)
kl
+∆ε

(VP)
kl

)
+ S̃ij∆C

(TE)
ijkl
ε
(E)
kl
−A

S̃ijC
(TE)
ijkl
S̃kl+

2
3
σ2
Y

(
C̃+ 2

3
ET
) (2.5)

where

A=
2
3
σY (ε

(VP)
i , ε̇

(VP)
i ,T)(ĖT∆ε̇

(VP)
i +TT∆T) (2.6)

is the positive scalar variable, tanα = ET(T) = ∂σY/∂ε
(VP)
i is the hardening modulus at T ,

tanβ= ĖT = ∂σY/∂ε̇
(VP)
i is the change in temporary yield stress with a change in the intensity

strain rate, tanγ = TT = ∂σY/∂T is the change of temporary yield stress with a change in
temperature, C̃(T) = [2E(T)ET(T)]/[3E(T)−ET(T)] is the kinematic hardening parameter
and E is Young’s modulus at the temperature T . To further evaluate the above expression
for∆λ, it is necessary to obtain the partial derivatives ∂σY/∂ε

(VP)
i , ∂σY/∂ε̇

(VP)
i and ∂σY/∂T .

It is assumed that the relationship between σY , ε
(VP)
i , ε̇

(VP)
i and T can be derived from the

data obtained in a series of tensile tests at different temperatures and strain rates using virgin
material specimens.
Substituting equation (2.5) into equation (2.4) and using ∂F/∂S̃ij = S̃ij, we obtain

∆ε
(VP)
ij = S̃

∗∗∆εij +∆ε
∗∗

ij (2.7)

where

S̃∗∗ = S̃∗ijC
(TE)
ijmnS̃mn (2.8)

is a positive scalar variable

∆ε∗∗ij = S̃
∗

ij

(
−C
(TE)
ijmnS̃mnαm∆Tδij + S̃mn∆C

(TE)
mnkl
ε
(E)
kl
−A

)
(2.9)

is the component of the incremental strain tensor, and

S̃∗ij =
S̃ij

S̃ijC
(TE)
ijkl
S̃kl+

2
3
σ2Y

(
C̃+ 2

3
ET
) (2.10)

is the component of the stress tensor.
Substitutingequations (2.3) and (2.7) into equation (2.2),weobtain theconstitutive equation

of mixed hardening for an isotropic material that includes the combined effects of thermo-
elasticity and thermo-visco-plasticity

∆σij =C
(TE)
ijkl
(1− S̃∗∗)∆εkl−C

(TE)
ijkl
(∆ε∗∗kl +αm∆Tδkl)+∆C

(TE)
ijkl
ε
(E)
kl
=C

(TE)∗
ijkl

∆εkl+∆σ
∗∗

ij

(2.11)

where

C
(TE)∗
ijkl

=C
(TE)
ijkl
(1− S̃∗∗)

∆σ∗∗ij =−C
(TE)
ijkl
(∆ε∗∗kl +αm∆Tδkl)+∆C

(TE)
ijkl
ε
(E)
kl

(2.12)



Application of variational and FEM methods... 491

3. Damping models

Avery essential and difficult problemduring themodelling of technological processes is conside-
ration of the phenomenon of energy dissipation (Gontarz andRadkowski, 2012). The dissipation
of energy influences the decrease of displacement anddynamic stress in the object. Thedifficulty
depends on selection of the proper physical modelmainly as well as qualification of suitable va-
lues of its parameters. At present (Bathe, 1982), usually theRayleigh damping is used, in which
the damping matrix C is calculated from the equation: C=α1M+α2K, whereM and K are
mass and stiffness matrices, α1 and α2 are constants to be determined from two given damping
ratios that correspond to two unequal frequencies of vibration. In the case of the modelling of
dynamical processes, one should apply dampingmodels dependent on velocity, strain and strain
rate. For this reason, we propose to calculate the dissipation of energy, for a typical step time,
from the following expression

∆Ed ∼=α1
∫

V

ρ∆u̇i∆ui dV +α1
∫

V

ρ∆u̇iui dV +α1
∫

V

ρu̇i∆ui dV +
1
2
α2

∫

V

C
(TE)
ijkl
∆ε̇ij∆εkl dV

+
1
2
α2

∫

V

C
(TE)
ijkl
∆ε̇ijεkl dV +

1
2
α2

∫

V

C
(TE)
ijkl
ε̇ij∆εkl dV

(3.1)

where:∆ui,∆u̇i are the i-th increment components of displacement and velocity vectors at the
typical step time ∆t, respectively, ui is the accumulated i-th component of the displacement
vector at the time t.
Constitutive models (2.11) and (3.1) can be used for practical engineering analysis. In the

next Section, thesemodels are used for variational formulation of a nonlinear equation ofmotion
for the slitting operation.

4. Variational formulation of the equation of motion

In this Section, we develop the equation of motion for a three-dimensional body in the global
Cartesian coordinates {z} and the updated Lagrange incremental formulation. Variational me-
thods are used in continuum mechanics to formulate the equations of motion (Kleiber, 1985;
Orłowska and Czapp, 2012). However, there is no method of accounting for the incremental
formulation. The use of the variational method to formulate the equation ofmotion for a typical
incremental step in nonlinear problems in technological processes was proposed by Kulakowska
and Kukielka (2008), Kukielka et al. (2012). In this paper, this method is used to model the
slitting process. Assuming that a numerical solution has been obtained at discrete time points
∆t,2∆t,. . ., the solution for t+∆t is desired. In this case, the functional increment is formu-
lated for the incremental displacement ∆F(∆üi,∆u̇i,∆ui) = ∆F(·), where ∆ui, ∆u̇i, ∆üi are
the incremental components of the displacement, velocity and acceleration vectors, respectively.
Using the conditions for the functional∆F(·) being stationary (Kukielka et al., 2012)

δ[∆F(·)] =
∂[∆F(·)]
∂(∆ui)

δ(∆ui)=
∂[∆F(·)]
∂(∆ui)

∆ui (4.1)

where the underline (·) denotes “variation in”, and because∆ui is the only variable, we obtain

δ[∆F(·)] =
∫

V

ρ(üi+∆üi)∆ui dV +α1
∫

V

ρ∆u̇i∆ui dV −2ω
∫

V

ρ∆u̇iΩij∆uj dV

+
1
2
α2

∫

V

∆ ˙̃εijC
(TE)
ijkl
∆ε̇kl dV +

1
2
α2

∫

V

∆ ˙̃εijC
(TE)
ijkl
∆εkl dV +

1
2
α2

∫

V

∆ε̇ijC
(TE)
ijkl
∆ε̃kl dV



492 Ł. Bohdal, L. Kukiełka

+
1
2
α2

∫

V

∆ε̇ijC
(TE)
ijkl
∆εkl dV +

1
2
α2

∫

V

∆ ˙̃εijC
(TE)
ijkl
∆ε̃kl dV +

1
2
α2

∫

V

∆ ˙̃εijC
(TE)
ijkl
∆εkl dV

+
1
2
α2

∫

V

∆ε̇ijC
(TE)
ijkl
∆ε̃kl dV +

1
2
α2

∫

V

∆ε̇ijC
(TE)
ijkl
∆εkl dV +

∫

tV

∆ε̃ijC
(TE)∗
ijkl

∆ε̃kl dV (4.2)

+
∫

V

∆ε̃ijC
(TE)∗
ijkl

∆εkl dV +
∫

V

∆εijC
(TE)∗
ijkl

∆ε̃kl dV +
∫

V

∆εijC
(TE)∗
ijkl

∆εkl dV

+
1
2

∫

V

(Tij +∆σ
∗∗

ij )∆ε̃ij dV +
1
2

∫

V

(Tij +∆σ
∗∗

ij )∆εij dV −ω
2
∫

V

ρ∆uiΩijΩij∆ui dV

−ω2
∫

V

ρriΩijΩij∆ui dV −

∫

V

ρ(fi+∆fi)∆ui dV −
∫

Σk

(p̂i+∆p̂i)∆ui dΣk =0

where: Tij is the component of Cauchy’s stress tensor, ∆ε̇ij and ∆ ˙̃εij are the linear and
nonlinear incremental components of the Green-Lagrange strain rate tensor, respectively
(∆ε̇ij = ∂∆εij/∂t; ∆ ˙̃εij = ∂∆ε̃ij/∂t, ∆εij = (∆ui,j + ∆uj,i + ul,i∆ul,j +∆ul,iul,j)/2 and
∆ε̃ij = (∆ui,k∆uj,k)/2), ρ(T) is the temperature-dependent mass density at the time t and
εij are the accumulated components of the total strain tensors at the time t (depending on
the history of deformation and temperature). fi, ∆fi are the components of the internal force
and incremental force vectors, respectively,p̂i,∆p̂i are the components of the externally applied
surface force and surface incremental force vectors in the contact body zones, respectively, and
Ωij are the components of the gyro tensor. However, in the typical updated Lagrange formula-
tion∆εij =(∆ui,j+∆uj,i)/2. The integrations are performed over the volume V and surfaceΣ
of the body, and equation (4.2) provides the basis for the finite element discretisation required
to obtain the solution.

5. Numerical analysis algorithm

Assuming that the temporary step time∆t is very small, it is possible to remove the incremental
stiffness matrix (∆K≈0) and the internal incremental force vector (∆F≈0). Then, using the
principle of incremental decomposition, approximating ṙ and r̈ in terms of r in equation (4.2)
using the central difference method (DEM), and using the following approximate (Bathe, 1982)

t
ṙ=

1
2∆t
(τr−t−∆tr) tr̈=

1
∆t2
(τr−2tr+t−∆tr) (5.1)

we obtain the solution of τtr

M̃
τ
tr=Q (5.2)

where M̃ is the effective mass matrix and Q̃ the effective load vector

M̃= a0M+a1C

Q̃=F+R+a0M(2
t
r−t−∆tr)+a1C

t−∆t
r

(5.3)

with the integration constants a0 =1/(∆t2), a1 =1/(2∆t).

5.1. Step-by-step solution using the Central Difference Method

In this Section, the Central Difference Method (Bathe, 1982; Patyk and Kukielka, 2008) is
adopted to solve Eq. (5.2) with initial and boundary conditions.



Application of variational and FEM methods... 493

• Initial calculations:

1. Form the massM and dampingCmatrices.

2. Initialise 0r, 0ṙ and 0r̈.

3. Select a time step ∆t < ∆tcr = TN/π and calculate the integration constants a0
and a1.

4. Calculate −∆tr.

5. Form the effective stiffness matrix M̃ from equation (5.3).

• For each time step:

1. Calculate the effective load vector Q̃ from equation (5.4).

2. Partitioning the equilibrium equation for two blocks, write the problem in form
[
M̃
n×n
11 M̃

n×w
12

M̃
w×n
21 M̃

w×w
22

]{
τ
tr
n×1
1

τ
tr
w×1
2

}
=

{
Q̃
n×1
1

Q̃
w×1
2

}

where the vectors rw×12 , Q̃
n×1
1 are known, and the vectors r

n×1
1 , Q̃

w×1
2 are unknown.

3. Solve for the displacement vector τtr
n×1
1

τ
tr
n×1
1 =(M̃

n×n
11 )

−1[Q̃1−M̃
n×w
12

τ
tr
w×1
2 ]

4. Substitute the vector τtr
n×1
1 into the equation

Q̃
w×1
2 = M̃

w×n
21

τ
tr
n×1
1 +M̃

w×w
22

τ
tr
w×1
2

and solve it to obtain the load vector Q̃w×12 in the tool-object contact zone.

5. Calculate the vectors tṙ and tr̈ from equation (5.1).

6. Sample solution

A three-dimensional finite element model of the slitting is presented in Fig. 2a. The model is
created based on the experimental configuration and the geometrical parameters of our test
stand shown in Fig. 2b. The slitting machine (KSE 10/10) consists of two rotary knives driven
by the engine.
The contact between the knives and sheet is considered a non-sliding contact that uses a

polyurethane roll, whichmoves the sheet in the horizontal direction. The sheet is clamped by a
sheet holder, and a rake angle (α= 7◦) is applied to the upper knife. The radius of the upper
knife carry out 14mm, the radius of the lower knife is 12mm.Thewidth of the knives carry out
15mm.Numerical calculations are performed for the 3D state of strain and 3D state of stress in
the FEmodel. It is assumed that the cutting process takes place under a constant temperature
(∆T = 0). The models in (2.11) and (4.2) reduce then to the elastic (in the reversible zone)
and visco-plastic (in the non-reversible zone) strains. Moreover: λ(T)=λ, µ(T)=µ, ρ(T)= ρ,
∆ε
(TH)
ij = 0, ∆C

(TE)
ijkl
(∆T) = 0, C(TE)

ijkl
= C(E)

ijkl
and ∆σ∗∗ij = −C

(TE)
ijkl
∆ε∗∗kl . Aluminum alloy

AA6111-T4, which is often employed for exterior panels in the automotive industry, is used
to simulate typical production conditions. An element with length of l = 50mm, width of
wi =30mmand thickness of t=1.5mm is analysed. The kinematics of different components is
as follows: first the upper knife and the lower onemoves vertically in opposite directionswith the
same constant velocity v1 =200mm/s in order to cut the sheet. Then, the sheetmoves along the



494 Ł. Bohdal, L. Kukiełka

Fig. 2. (a) Scheme of the 3D FEmodel of the slitting process; (b) schematic view of the experimental
slitting equipment: 1 – knives, 2 – sheet holder, 3 – engine, 4 – clearance regulator, 5 – drive pedal,

6 – slitting velocity regulator, 7 – force sensor

Z axis with a constant velocity v2 = 500mm/s, and at the same time the upper knife and roll
rotates with a angular velocityw1 =36rad/s and the lower knife rotate in the opposite direction
with the same angular velocity. The value of horizontal clearance carry out c = 0.1mm. The
objects aremeshedwith an 8-nodeSolid164 element typewith reduced integration andhourglass
control, and the mapped mesh with various densities is generated in three areas on the sheet
(Fig. 2a). The size of elements (where a is defined as element length, and h is defined as element
height) in each zone is as follows: in zone 1: a=0.8mm, h=0.37mm, in zone 2: a=0.5mm,
h = 0.37mm, in zone 3: a = h = 0.1mm. The finite element model of the sheet uses 320000
elements, that of the upper knife uses 27648, that of the lower knife uses 11520 and that of
the roll uses 25920. The contact between the tools and the deformable sheet metal is described
usingCoulomb’s frictionmodel, and constant coefficients of static friction µs =0.08 and kinetic
friction µd =0.009 are accepted.
The tools (knives) are considered as rigid bodies, but the sheet is modelled as an isotropic,

elastic/visco-plastic material with nonlinear hardening the temporary yield stress of which is
described with the aid of the Cowper–Symondsmodel (Bohdal andKukielka, 2006)

σY =

(
1+
ε̇
(P)
i

C

)m(
σ0+βEpε

(P)
i

)
(6.1)

where β is the plastic strain hardening parameter, σ0 is the initial, static yield point, ε̇
(P)
i is the

plastic strain rate intensity, C is the material parameter defining the effect of the plastic strain
rate intensity,m=1/n is thematerial constant defining its sensitivity to the plastic strain rate,
ε
(P)
i is the plastic strain intensity andEp =ETE/(E−ET) is thematerial parameter depending

on the modulus of plastic strain hardeningET = ∂σY/∂ε
(P)
i and Young’s modulusE.

TheCowper-Symondsmaterial deformationmodel is frequently used to perform simulations
of dynamic processes using the FE method. Although the majority of material constants are
determined using a classical tension test, problems are associated with the determination of
the Cowper-Symonds constantsC and n. In a previous report (Bohdal andKukielka, 2006), we
established values of theCowper-Symondsconstants at differentdeformation rates. In this study,
the values of the constants for aluminum alloy AA6111-T4 are taken as follows: C = 6500s−1

and n = 4, while ρ = 2720kg/m3, σ0 = 145MPa, E = 76GPa, ν = 0.27, ET = 25MPa, and
β=1.



Application of variational and FEM methods... 495

To simulate the crack initiation and propagation in the material, we use a new constitutive
model proposed by Xue and Wierzbicki (2013). The model covers the full range of plasticity
until the onset of fracture. It is understood that the fracture initiation in uncracked solids is
the ultimate result of a complex damage accumulation process induced by plastic deformations.
Fracture can be predicted for complex loading paths that are not limited to a restricted loading
inwhich the pressure is constant. Themodel also incorporates the coupling between the damage
and the strain hardening function. We determined the model constants experimentally and
presented them in (Bohdal and Kukielka, 2013).

7. Results and discussion

7.1. Stress and strain analysis

Figure 3 shows contours of the normal stresses σxx, σyy and σzz. The stresses σyy and σzz
show high compressive values in the region of the partially slit sheet, which is in contact with
the knives. The compressive stresses occur because of pushing down of the knife surface on the
sheet.

Fig. 3. Contours of the normal stresses; (a) σxx, (b) σyy and (c) σzz (25% step time)

As the sheet slits, it moves tangentially to the blade. This causes the area of contact with
the knife blade on the sheet to be inclined to the horizontal at an angle. The normal compressive
stress is thus split into two components in the direction of the axes Z and Y contributing to
the two normal stresses. In addition, it is observed that σxx is also highly compressive in the
same region of the partially slit sheet. The shear stress shows high values in the region where
the sheet is expected to slit and the values drop down as the knives move away from the region
(Fig. 4). The high shear stress and effective strain is caused by the shearing action of the two
blades on either side of the sheet (Figs. 4 and 5). In the cutting area, material concentration
takes place. After reaching the critical value of the effective strain, the process of integrity loss
occurs (Fig. 5b).The results are close to those obtained by Gąsiorek (2013a) during cutting of
sheet metal bundles using a guillotine.

7.2. Energy balance and slitting force analysis

Figure 6 shows the energy curves obtainedduring theprocess. In shearingprocesses, the total
increase of energy is nearly completely related to the increase of deformations in the integrity
loss area (Gąsiorek, 2013b;Bollen et al., 1989). Theminimal increase of the kinetic energy results
from the process of horizontal movement of a sheet during slitting (Fig. 6).
The evolution of the slitting force and its comparison to experimental measurement versus

the upper knife penetration is shown in Fig. 7 with some deformed configurations of the sheet.



496 Ł. Bohdal, L. Kukiełka

Fig. 4. Values of the maximum shear stresses measured in the shearing region during slitting

Fig. 5. (a) Contours of the effective plastic strain for 75% time step, (b) characteristics of the
plasticization process for a selectedmesh element describing the sheet in the cutting area up to the

moment of decohesion

Fig. 6. Energy balance characteristics [10−6kJ]

Three stages can be clearly observed. The first stage (AB) corresponds to a strong increase
in the slitting force due to the increase of the bending forces in the sheet cross-section when the
knives penetrate the sheet. The second stage extending from points B to C is characterized by
an aquasi-constant slitting force Fsta ≈ 179N which clearly indicates the steady state of this
process stage. The last stage (CD) of the slitting curve corresponds to the drastic decrease of
the slitting force due to unstable propagation of the macroscopic crack along sheet thickness
terminating the slitting operation. The experimental verification of the cutting force at steady
state is shown in Fig. 7b. Values of forces obtained from numerical calculation show a little
difference from the experimental results. This difference does not exceed 12%.



Application of variational and FEM methods... 497

Fig. 7. A comparison of the resultant force from the experiment and the FEM simulation

7.3. Sheared edge contours

In this Section,we illustrate the capability of the present solution procedure to reproduce the
effects of someprocess technological parameters such as slitting velocity andhorizontal clearance
on thequality of the sheared edge. In the industrial practice ahigh slitting sheared edgequality is
required for all good products. In general, good sheared edge quality for cut sheets is indicated
by smooth (burnished), clean, and straight edges with no edge wave (thin sheets), and most
importantly, with no or minimal burr (Fig. 8). The area values (Fig. 8) are measured from
simulations and experiments at different locations over the cut edge in the z-direction, averaged
and presented inFigs. 9 and 10. The results are in good agreementwith those obtained from the
experiments, with an approximate error margin of 13%. The experiments and simulations are
executed for the slitting velocities v = 3, 7, 17, 27, 30, and 32m/min and clearances c= 0.03,
0.05, 0.09, 0.13, and 0.15mm.

Fig. 8. Typical sheared edge profile with marked zones: (a) front view and (b) cross-sectional view

Fig. 9. Influence of the slitting velocity and clearance on the size of: (a) burr and (b) burnish zone



498 Ł. Bohdal, L. Kukiełka

Fig. 10. Characteristic features of chosen sheared edges viewed by optical microscopy “Vision
Engineering” and predicted by FE simulations at (a) v=7m/min, c=0.05mm; (b) v=17m/min,

c=0.09mm, (c) v=32m/min, c=0.09mm

Figure 9a shows experimental results of the effect of clearance and slitting velocity on burr
height. The general trend seen in light alloys during another shearing process (blanking, guillo-
tining) is for the burr height of the workpiece to increase with the increasing cutting clearance
(Khadke et al., 2005). There are differences, though, in the relative rate of burr height increase
with clearance. Increasing theclearance from c=0.03mm(2%of sheet thickness) to c=0.09mm
(6%t) increases the burr height for the analysed slitting velocities. When the clearance is set
at c > 6%t, the burr height decreases. The highest burrs are obtained when the middle range
of velocities is used (v = 17m/min) and when c= 6%t. The height of the burrs in this case is
non-uniform along the line of cutting, with some deep pockets that can serve as stress concen-
trators if stretching is applied parallel to this surface (Fig. 10b). Burrs can become separated
from the cut part if there is a significant gradient of clearance along the shearing line. These
local burrs are subjected to additional forces and torn off from the side of the sheared surface.
When the clearance is set at c=6%t, the burnished zone is reduced (Fig. 9b) and the velocity
effect is reinforced. A high zone value (approximately 66% of the sheet thickness) is obtained
when clearances of c=2%t and c=10%t are used.

8. Conclusions

The objective of this paper is to present the modelling, analysis, and testing of a solution
procedure for the analysis slitting process. Dynamic and damage effects, thermo-mechanical
coupling and contact friction are taken into account. The procedure is implemented into a
finite element code ANSYS LS-DYNA. Themain conclusion is that most trends in the slitting
process are qualitatively well described by the developedmodel. Themodel developedmakes it
possible to analyze the states of stresses and strains at anymoment of the process and to analyse
the causes of defects in the metal sheets after processing (deformation, edge wave, defects of
the sheared edge, e.g., burrs). In the slitting, it is important to apply the required boundary
conditions considering the dimensions of the sheet that is being cut. Special attention should be
paid to the correct modeling of clamping of the sheet. This also means that 3D models cannot
be simply replaced by 2D approximations.
The simulation results agreewellwith experiments, and themodel canbeusedboth to design

the slitting process and as a support to the solution of practical problems. Actual investigations
concern simulations of the slitting in which knives cannot be considered as rigid bodies. Further



Application of variational and FEM methods... 499

work is still to be done in order to introduce some other effects concerning the behavior of rolled
metal sheets with damage induced anisotropy and spring-back effects.

References

1. Aggarwal S., BhushanB., KatsubeN., 2005,Three-dimensional finite element analysis of the
magnetic tape slitting process, Journal of Materials Processing Technology, 170, 71-88

2. Bathe K.J., 1982,Finite Element Procedures In Engineering Analysis, Prentice-Hall, Englewood
Cliffs, N.J.

3. BohdalL.,KukielkaL., 2006,The effect of selectedmaterial parameterson the stress and strain
states in the process of cutting a sheet plate with circular cutters,Task Quarterley, 4, 391-400

4. Bohdal, L., Kukielka L., 2013, The modeling and numerical analysis of the shearing process
with the regard of the geometrical and physical nonlinearity, Ecological Aspects of Applying New
Technologies in Transport, 131-140

5. Bollen D., Deneir J., Aernoudt E., Muylle W., 1989, Shear cutting of PET film, Journal
of Material Science, 24, 2957-2966

6. GałęziaA.,Gontarz S., JasińskiM.,Mączak J., Radkowski S., Seńko J., 2012,Distribu-
ted system for monitoring of the large scale infrastructure structures based on analysis of changes
of its static and dynamic properties,Key Engineering Materials, 518, 106-118

7. Gąsiorek D., 2013a, Modelling And Experimental Investigation Of Dynamic Processes Occur-
ring During Cutting Lithograpic Plates Using Guillotines (in Polish), Wydawnictwo Politechniki
Śląskiej, Gliwice

8. Gąsiorek D., 2013b, The application of the smoothed particle hydrodynamics (SPH) method
and the experimental verification of cutting of sheet metal bundles using a guillotine, Journal of
Theoretical and Applied Mechanics, 51, 4, 1053-1065

9. Gontarz S., Radkowski S., 2012, Impact of various factors on relationships between stress and
eigenmagnetic field in a steel specimen, IEEE Transactions on Magnetics, 48, 3, 1143-1154

10. Khadke A., Ghosh S., Li M., 2005, Numerical simulations and design of shearing process for
aluminum alloys,ASME Journal Manufacturing Science Engineering, 127, 3, 612-621

11. Kleiber M., 1985, Finite Element Method in Non-Linear Solid Mechanics, IPPT PAN, PWN,
Warsaw-Poznań

12. KukielkaL., GeletaK., KukielkaK., 2012,Modelling of initial and boundary problemswith
geometrical and physical nonlinearity and its application in burnishing processes, Steel Research
International, Special Edition, 14th International Conference on Metal Forming, 1375-1378

13. Kulakowska A., Kukielka L., 2008, Numerical analysis and experimental researches of bur-
nishing rolling process with taking into account deviations in the surface asperities outline after
previous treatment, Steel Research International, 2, 42-48

14. LiuC., LuH.,HuangY., 2005,Dynamic steady-state stress field in awebduring slitting,ASME
Journal Applied Mechanics, 72, 157-164

15. Lu H., Wang B., Iqbal J., 2001, Deformation in shear slitting of polymeric webs, Proceedings
of the 6th International Conference on Web Handling, Stillwater, OK, June 10-13, 389-402

16. Ma J., Lu H., Li M., Wang B., 2006, Burr height in shear slitting of aluminum webs, ASME
Journal Manufacturing Science Engineering, 128, 46-55

17. Meehan R. R., Burns S. J., 1998,Mechanics of slitting and cutting webs,Experimental Mecha-
nics, 38, 2, 103-109

18. Orłowska M., Czapp M., 2012, Numerical analysis of heat efficiency of the convective heat
exchanger build with horizontal plates,Annual Set The Environment Protection, 14, 582-586



500 Ł. Bohdal, L. Kukiełka

19. Patyk R., Kukielka L., 2008, Optimization of geometrical parameters of regular triangular
asperities of surface put to smooth burnishing, Steel Research International, 2, 642-647

20. Saanouni K., 2006, Virtual metal forming including the ductile damage occurrence. Actual state
of the art andmain perspectives, Journal of Materials Processing Technology, 177, 19-25

21. Wisselink H.H., 2000, Analysis of guillotining and slitting: finite element simulations, Ph.D.
Thesis, University of Twente, The Netherlands

22. Wisselink, H. H., Huétink J., 2004, 3DFEM simulation of stationarymetal forming processes
with applications to slitting and rolling, Journal of Materials Processing Technology, 148, 328-341

23. Xue L., Wierzbicki T., 2013, Verification of a new fracture criterion using LS-DYNA, 9th In-
ternational LS-DYNA Users Conference, Simulation Technology, 3, 13-22

Manuscript received January 8, 2014; accepted for print December 30, 2014