...

Svoboda2010-MetalCutting.pdf

by user

on
Category: Documents
4

views

Report

Comments

Transcript

Svoboda2010-MetalCutting.pdf
Home
Search
Collections
Journals
About
Contact us
My IOPscience
Simulation of metal cutting using a physically based plasticity model
This article has been downloaded from IOPscience. Please scroll down to see the full text article.
2010 Modelling Simul. Mater. Sci. Eng. 18 075005
(http://iopscience.iop.org/0965-0393/18/7/075005)
View the table of contents for this issue, or go to the journal homepage for more
Download details:
IP Address: 129.5.228.139
The article was downloaded on 31/10/2012 at 23:10
Please note that terms and conditions apply.
IOP PUBLISHING
MODELLING AND SIMULATION IN MATERIALS SCIENCE AND ENGINEERING
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005 (19pp)
doi:10.1088/0965-0393/18/7/075005
Simulation of metal cutting using a physically based
plasticity model
Ales Svoboda1 , Dan Wedberg2 and Lars-Erik Lindgren1
1
2
Luleå University of Technology, 971 87 Luleå, Sweden
AB Sandvik Coromant, 811 81 Sandviken, Sweden
E-mail: [email protected]
Received 2 February 2010, in final form 8 July 2010
Published 25 August 2010
Online at stacks.iop.org/MSMSE/18/075005
Abstract
Metal cutting is one of the most common metal shaping processes. Specified
geometrical and surface properties are obtained by break-up of the material
removed by the cutting edge into a chip. The chip formation is associated with
a large strain, high strain rate and a locally high temperature due to adiabatic
heating which make the modelling of cutting processes difficult. This study
compares a physically based plasticity model and the Johnson–Cook model.
The latter is commonly used for high strain rate applications. Both material
models are implemented into the finite element software MSC.Marc and
compared with cutting experiments. The deformation behaviour of SANMAC
316L stainless steel during an orthogonal cutting process is studied.
(Some figures in this article are in colour only in the electronic version)
1. Introduction
Cutting tool designs and cutting process parameters may be optimized by experimental
measurements or by numerical analysis. In deformation zones in front of the tool and
beneath the cutting edge, the material is deformed plastically by simultaneous action of large
compressive and shearing stresses. The friction forces and dissipative plastic work generate
high temperatures. This makes direct observations in the cutting zone during the machining
difficult. Alternatively, the chip generation and flow may be observed by freezing the motion
using quick-stop techniques. Nevertheless, these measurements are time consuming and a
limited amount of information can be gained from experiments in comparison with simulations.
The finite element method (FEM) is the most common numerical method used to analyse
metal cutting operations, see Vaz et al (2007). However, the modelling of the special
conditions in the cutting zone requires a robust finite element software including capabilities
such as thermo-mechanical coupling, friction models, material models and effective contact
algorithms. An additional degree of complexity is caused by the necessity to model material
0965-0393/10/075005+19$30.00
© 2010 IOP Publishing Ltd
Printed in the UK & the USA
1
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
removal. The finite element mesh distortion due to large deformations requires a remeshing
technique in the case when a Lagrangian reference frame is used.
An important feature of a model for metal cutting is the material model. It must adequately
represent the deformation behaviour during high rate loading and be able to account for strain
hardening and strain softening, thermal softening as well as large variations in strain rate
and temperature. In the majority of commercial FEM codes, the deformation behaviour
of metals is represented by an empirical relationship, mostly in terms of the power law of
strain and strain rate. These equations, relating the flow stress to the plastic strain, strain
rate and temperature, are not based on any reasoning about underlying physics. The empirical
relations are usually lacking in predictive capabilities beyond the derived range of experimental
conditions at which they were curve-fitted. It is therefore preferable to use models which are
related to the underlying physics of the deformation coupled to the microstructure evolution.
The advantage of using physically based models is an expected larger domain of validity as
compared with phenomenological models. However, the physically based models must still
be feasible for large scale computations.
In this study, which is part of a project aiming at development and validating models
for machining, the physically based dislocation density (DD) model is compared with the
phenomenological Johnson–Cook (JC) plasticity model. The latter is commonly used in
machining simulations. Both models were used to simulate orthogonal cutting of SANMAC
316L stainless steel, a material with a face-centred cubic (FCC) structure and low stacking-fault
energy. The steel has significant amounts of alloying elements which improve the corrosion
resistance and the strength. These properties may be reduced during the machining when large
residual tensile stresses are introduced. Residual stresses in the workpiece may also affect
fatigue life of machined components. Therefore, it is highly technically relevant to be able to
predict a stress field after machining operations in order to optimize the machining parameters.
The measurements of cutting forces and quick-stop tests were performed in order to
evaluate accuracy and performance of the two material models. The simulated cutting forces,
the chip morphology and dimensions are compared with the experimental results.
2. Material modelling
The material models discussed in the following sections were implemented into MSC.Marc
software using user subroutine interface WKSLP. The subroutine makes it possible to introduce
user-defined relations for the yield stress and the corresponding hardening slope as a function of
temperature, equivalent plastic strain and equivalent plastic strain rate. The implementation of
the material models is based on the additive decomposition of the spatial rate of the deformation
tensor. A radial-return type algorithm was used for the integration of constitutive equations.
2.1. Physically based plasticity model
Related to the theory of dislocation mechanics, Zerilli and Armstrong (1987) developed a
constitutive model including different rate controlling mechanisms for body-centred cubic
(BCC) and FCC metals. Jaspers and Dautzenberg (2002) investigated the applicability
of Zerilli–Armstrong (ZA) model to metal cutting simulations in the comparison with the
JC model.
Guo et al (2006) compared predicted chip morphology obtained using the conventional JC
model with the Bammann–Chiesa–Johnson (BCJ) model. The dislocation mechanics based
BCJ model, Bammann et al (1996), incorporates strain rate and temperature sensitivity, as
well as damage, through a yield surface approach. Meyers et al (2002) presented a review of
2
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
physically based models for plasticity due to dislocation glide as well as twinning. They also
discussed high strain rates phenomena. The paper has a particular focus on the ZA model.
The current model is an explicit physically based model which includes a coupled set of
evolution equations for the internal state variables, DD and vacancy concentration. Explicit
means that these internal state variables are used instead of accumulated effective plastic strain
in the previous mentioned (implicit) physically based models. The concept of the DD is the
amount (length) of dislocations for a representative volume element divided by its volume. The
model considers two different densities, a mobile and an immobile DD. The microstructure is
not represented explicitly, but in average sense. Basic equations of the model are shown in the
next section. Details of the model, experimental procedure and parameter optimization are to
be found in Lindgren et al (2008) for low strain rates. The model is based on the dislocation
glide mechanism; see Frost and Ashby (1982).
2.1.1. Flow stress and long-range term contribution. We assume that the flow stress can be
expressed as a combination of the long range and the short-range parts of the resistance to the
motion of dislocations, e.g. Bergström (1970), Estrin (1998), Follansbee and Kocks (1988).
Hence, the yield stress is defined by
σ y = σG + σ ∗ ,
(1)
where σG is due to long-range interactions with the dislocation substructure. It is an athermal
contribution and is related to the immobile DD. One common assumption for the long-range
term is
√
σG = mαGb ρi .
(2)
Here, m is the Taylor orientation factor transforming the effects of resolved shear stress in
different slip systems into effective stress–strain relations and is affected by the texture, α is
aproportionality factor, ρi is the immobile DD, G is the temperature dependent shear modulus
and b is the Burgers vector.
2.1.2. Flow stress and short-range term contribution. The short-range term σ ∗ in (1) is the
thermally activated flow stress component. It is the stress needed for a dislocation to pass
short-range obstacles and to move it through the lattice. The thermal vibrations of the lattice
assist in overcoming obstacles. The velocity and the density of mobile dislocations are related
to the plastic strain rate according to the Orowan equation
ρm bν̄
p
ε̄˙ =
,
m
(3)
where ν̄ is the average velocity of mobile dislocations having a density ρm . The time taken
by a dislocation to move over a distance between two obstacles consists of a waiting time and
a running time. The moving dislocation has a waiting time in front of an obstacle before it
manages to pass the obstacle (provided it does not become immobile) and then moves to the
next one. The average velocity is related to the waiting time, as the running time is assumed
to be negligible. The waiting time and thereby the average velocity is assumed to depend on
the Gibbs free-energy of activation G for cutting or by-passing of obstacles (see Frost and
Ashby (1982)), and on the temperature T . The average velocity is defined as
G
ν̄ = υa exp −
,
(4)
kT
3
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
where is the mean free path for dislocations between two obstacles, υa is the attempt
frequency which depends on the characteristics of obstacles and k is Boltzmann’s constant.
The previous relations lead to
ρm bνa
G
p
ε̄˙ =
exp −
.
(5)
m
kT
The short-range stress component σ ∗ may have a contribution from the lattice itself as well
as from different kinds of obstacles. In FCC structures the lattice resistance to motion of
dislocations is smaller as compared with the resistance provided by the discrete obstacles.
Large obstacles will make the dislocations immobile and then contribute to the long-range
term. One common relation between the activation energy G and the short-range term
representing a typical barrier encountered by a dislocation is
∗ p q
σ
G = F 1 − ,
(6)
τ
where F is the total free energy required for a dislocation to overcome the lattice resistance
or obstacles without aid from external stress. The quantity τ̂ is the athermal flow strength that
must be exceeded in order to move dislocations across the lattice without the aid of thermal
energy. The exponents 0 < p 1 and 0 < q 2 are related to the shape of energy barriers.
The pre-exponential term in (5), which is approximated following Frost and Ashby (1982) to
be constant, is expressed as
ρm bνa
.
ε̄˙ ref =
m
Combining (5) with (6) and (7) yields
∗ p q σ
F
p
.
1− ε̄˙ = ε̄˙ ref exp −
kT
τ
(7)
(8)
Rewriting (8) we obtain the short-range stress component as a function of the effective plastic
strain rate, see Nemat-Nasser et al (2001):
˙ 1/q 1/p
kT
ε̄ref
∗
.
(9)
ln
σ = τ̂ 1 −
pl
F
ε̄˙
2.1.3. Evolution of DD. The equation for the flow stress (2) requires evolution equations for
internal state variables which are the DD and the vacancy concentration. The mobile DD is
assumed to be much smaller than the immobile one, Bergström (1983), Estrin (2003). The
immobile DD is expressed in terms of hardening (+) and recovery (−) contributions. The
presented model is tracing only the density of immobile dislocations ρi
(−)
(−)
ρ̇i = ρ̇i(+) − ρ̇i(glide)
− ρ̇i(climb)
.
(10)
Mobile dislocations move over a mean free path before they are immobilized or annihilated.
The immobile DD is assumed to increase proportional to the plastic strain rate, which is related
to the density of mobile dislocations, shown in (3), and inversely to the mean free path
m 1 ˙p
ε̄ .
(11)
b
The mean free path is related to the grain size g and the dislocation subcell diameter s
ρ̇i(+) =
1 1
1
= + .
g s
4
(12)
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
The effect of grain size on flow stress, the Petch–Hall effect, is accounted for this way and
contributes to the hardening. The size of subcells is related to the immobile DD by the
parameter Kc
1
s = Kc √ .
ρi
(13)
The DD may be reduced by different processes. Recovery, remobilization and/or annihilation
of dislocations are proportional to the current DD and controlled by dislocation climb and
glide. Recovery by dislocation glide is described by
p
(−)
ρ̇i(glide)
= ρi ε̄˙ ,
(14)
where is a recovery function that depends on temperature and strain rate. This model
accommodates only the dynamic recovery due to the strain rate. Static recovery controlled by
diffusion climb is assumed to have the form
cv Gb3 2
(−)
2
ρ̇i(climb)
= 2cγ Dν eq
.
(15)
ρi − ρeq
cν kT
eq
In (15), cv is the fraction of vacancies, cv is the thermal equilibrium vacancy concentration,
Dv is the self-diffusion coefficient and cγ is a material parameter related to the stacking-fault
energy. The DD decreases towards an equilibrium value ρeq .
2.1.4. Vacancy generation and migration. The calculation of the vacancy concentration is
required for the solution of (15). The generation and motion of vacancies are coupled with the
recovery of dislocations and diffuse solute atoms. The model presented here is only concerned
with mono-vacancies. When a crystal is retained a sufficient time at a given temperature, an
equilibrium level of vacancies is reached. Deforming the material or changing the temperature
generates the excess vacancies. The effect of excess vacancies on the diffusion is accounted
for via (16) as
σy b
cj
0 ˙ p
1
1 ex
eq
ċv = ċv − ċv = χ
+ς 2
(16)
ε̄ − Dvm 2 + 2 cv − cveq .
Qvf
4b
b
s
g
p
The stress in (16) is equal to the flow stress during a plastic deformation, the factor χ σy ε̄˙ is the
fraction of the mechanical work needed for the vacancy formation, Qvf is the activation energy
for forming a vacancy, 0 is the atomic volume and cj is the concentration of thermal jogs. The
eq
parameter ς describes the neutralization effect by vacancy emitting and absorbing jogs, cv is
the equilibrium concentration of vacancies at a given temperature, cv is the non-equilibrium
vacancy concentration and Dvm is the vacancy migration. The equilibrium concentration at a
given temperature is defined in (17) as
Svf
Qvf
cveq = exp
exp −
.
(17)
k
kT
In (17) Svf is the entropy increase when creating the vacancy. The rate of change in the
vacancy equilibrium concentration is related only to the temperature change
eq
eq Qvf
ċv = cv
Ṫ .
(18)
T2
Details of the model for vacancies and diffusion are to be found in Militzer et al (1994) and
Lindgren et al (2008).
5
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
2.1.5. Stress update algorithm. The radial-return operator for the integration of constitutive
equations is used for updating the flow stress. Most plasticity models use the plastic strain as
an internal variable. The DD model uses the immobile DD ρi and the vacancy concentration
cv as internal state variables. The computation of the increment of effective plastic strain,
which fulfils the consistency condition, requires calculation of the yield stress and hardening
modulus for the current iteration of the plastic strain and internal state variables. The evolution
of the internal state variables is governed by the coupled differential equations. The rate of
change in immobile DD given by (19) comes from (11)–(15). The rate of change in vacancy
concentration, (20), is given by (16), (17) and (18).
cv Gb3 2
m1
p
2
ρ̇i =
),
(19)
− ρi ε̄˙ − 2cγ Dv eq
(ρ − ρeq
b
cv kT i
cj 0 ˙ p
1
0
1
∗
eq
eq Qvf
ċv = χ
Ṫ . (20)
(σG + σ ) + ς 2
ε̄ − Dvm 2 + 2 (cv − cv ) + cv
Qvf
4b b
s
g
T2
These equations are solved incrementally by an implicit Newton–Raphson method. Once the
DD and vacancy concentration are known, the hardening modulus and the flow stress can be
evaluated. During the increment iteration the plastic strain rate is assumed constant.
The hardening modulus in the incremental form is given by
dσy
∂σ ∗
∂σG ∂ρi
∂ρi ∂cv
H = pl =
+
+
.
(21)
dε̄
∂ρi ∂ ε̄ pl ∂cv ∂ ε̄ pl
∂ ε̄ pl
Details of the stress update algorithm are shown in Lindgren et al (2008).
2.2. The JC plasticity model
The phenomenological Johnson and Cook (1983) plasticity model is commonly used in FE
simulations of metal cutting. In the study by Shatla et al (2001), Özel and Zeren (2004)
a modification of the JC model was applied to develop the methodology of flow stress
determination for metal cutting simulations. Mabrouki and Rigal (2006) investigated thermomechanical effects during chip formation using this model. Nasr et al (2007) used the JC model
to simulate the effects of tool-edge radius on residual stresses. Calamaz et al (2008) introduced
a strain softening effect in a modified JC model to predict the saw-tooth chip formation.
The JC model assumes an independent effect of the strains, strain rates and the
temperatures. The model provides numerical stability and shows good fit to measured stress–
strain data within a limited experimental range of strains, strain rates and temperatures. The
model is lacking the capability to capture the complex effects that are common in machining,
see Guo et al (2006).
In the classical form, the flow stress σy is expressed as
p T − Troom m
ε̄˙
pn
σy = A + B ε̄
1−
.
(22)
1 + C ln
Tmelt − Troom
ε̄˙ ref
The parameter A is related to the initial yield strength at room temperature, ε̄p is the effective
p
plastic strain, ε̄˙ is the effective plastic strain rate, ε̄˙ ref represents the highest strain rate for
which no rate adjustment to the flow stress is needed, Tmel t is the melt temperature.
Parameters A, B, C, n, m are user-defined material constants. The expression in the first
bracket of (22) gives a nonlinear hardening law; the expressions in the second and third bracket
represent the effects of strain rate and temperature, respectively.
6
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
Table 1. Chemical composition of SANMAC 316L wt (%).
C
Si
Mn
P
S
Cr
Ni
Mo
V
N
0.009
0.31
1.71
0.031
0.023
16.86
10.25
2.04
0.048
0.040
Table 2. Cutting data in experiments and simulations.
Test no
Cutting speed vc (m min−1 )
Feed fn (mm rev−1 )
Cutting depth ap (mm)
1
2
3
4
180
180
240
240
0.05
0.15
0.05
0.15
3.0
3.0
3.0
3.0
3. Experimental procedures
An austenitic stainless steel SANMAC 316L which is a machinability improved steel grade of
AISI 316L was used in the experimental study. The chemical composition is given in table 1.
The material parameters at room temperature are presented in section A.1.
The hot worked steel bar was annealed, straightened, peel turned and polished to final
delivery dimensions with a diameter of 190 mm. Several orthogonal turning and quick-stop
tests were performed in order to measure cutting forces and to capture morphology of the chip.
All cutting tests were performed on specimens with an outer diameter of 170 mm after the
hardened surface was removed.
In order to characterize the mechanical properties at the position where the forces were
measured, quasi-static uni-axial tension tests were performed in the longitudinal and the
transversal directions of the test bar. The average values of the yield strength and the tensile
strength were found to be 276 MPa and 541 MPa, respectively, in the longitudinal direction
and 280 MPa and 538 MPa, respectively, in the transversal direction.
3.1. Measurements of cutting forces and quick-stop tests
All tests were performed using the triangular CVD TiCN–Al2 O3 –TiN coated insert TNMG
160408-QF in grade 4015. This coating and substrate is preferable for continuous as well as
intermittent turning operations. The insert having a wedge angle of 90◦ was mounted on a
tool holder PTGNL 3225P in order to achieve a clearance angle of 6◦ . The insert used had
a grooved geometry made up by a radius, a primary land of 0.15 mm and a rake face with
an active rake face angle of −6◦ due to the tool holder. There are always some dimensional
deviations associated with manufacturing. In order to characterize the actual edge geometry,
the radius was measured at three positions along the cutting edge. Inserts with the edge radius
of 60 µm were subsequently chosen for the cutting force measurements and the quick-stop
tests.
Orthogonal cutting operations were performed on a George Fischer CNC turning lathe on
prefabricated tubes with one closed end having an outer diameter of 170 mm, a thickness of
3 mm and a length of 115 mm. A three component Kistler dynamometer of type 9263 was used
together with a 300 Hz low pass filter for measuring the cutting and feed force components.
Cutting data used for both experiments and simulations are presented in table 2.
The measured history of the cutting forces can be divided into four regions: the ramp-in; the
dynamic response; the steady state and the ramp-out, see Ivester et al (2007). All evaluations
were done by averaging the signal in the steady state region where the force oscillations
7
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
Figure 1. Quick-stop test—case no. 2: (a) chip geometry, (b) evaluation of the shear plane angle
from the quick-stop specimen.
converged towards a plateau of stable oscillations, independent of time. During the dynamic
response section, which follows after the short ramp-in, the forces are reduced down to steady
state conditions. The time histories of measured and simulated cutting forces are presented in
figure 6.
A number of quick-stop tests were performed in order to study the chip morphology.
One important prerequisite and challenge of the test is to freeze the chip formation process
during cutting without affecting the cutting action and distorting the chip. This can be done by
instantaneous reduction of the relative velocity between the insert and the work material to zero.
In this study, in-house developed quick-stop equipment based on the ‘shear pin’ design was
used; see Childs et al (2000) and Satheesha et al (1990). The tests were performed on the same
specimen geometry as used for the measurements of the cutting forces. When the steady state
conditions in the orthogonal turning operation were reached the process was instantaneously
stopped and the chip geometry evaluated. The chip was cut off and then embedded, grounded,
polished and etched for further examinations in a microscope for estimations of the thickness
and the shear plane angle. The chip geometry from quick-stop tests showed some serrated
form. Hence, the average chip thickness was obtained from the estimated chip area and the
chip length which were calculated from the microscope image.
The chip geometry for test case no. 2 is shown in figure 1(a). The evaluation of the shear
plane angle was done by drawing a tangent line from the cutting edge to the intersection of the
un-deformed surface of the work material and the chip as shown in figure 1(b).
However, deformations along the shear plane do not generally follow straight lines which
influenced the accuracy of the evaluation.
3.2. High strain rate testing
The Split-Hopkinson bar is commonly used when determining the inelastic response of metals
at high strain rates. This technique has the capability to measure the response of metallic
materials at strain rates up to 104 s−1 , see Nemat-Nasser (2000). In this work the high strain
rate testing was performed using the Split-Hopkinson pressure bar (SHPB) over a broad range
of temperatures and strain rates. The temperatures were ranging from room temperature up to
950 ◦ C and the strain rates varied up to maximum 9000 s−1 .
Characteristics of Split-Hopkinson technique are that the strain and the strain rate both
are affected by the impact velocity of the striker and the initial length of specimen. Different
8
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
impact velocities of the striker, initial length of specimens and length of the striker were used
to obtain different strains and strain rates. The impact velocity of the striker could in this case
be adjusted by different settings of pressure. A more detailed description of the test equipment
used in this work is given in Apostol et al (2003).
The tests at elevated initial temperature were performed by pushing an arm, loaded with
the specimen, into a stand alone furnace. When the desired temperature was reached, the
specimen was pulled out and placed between the incident and the transmitted bars but still
without physical contact. The transmitted bar was then pushed towards the specimen with a
pneumatic manipulator establishing contact just before the stress pulse reached the transition
between the incident bar and the specimen, all in accordance to prevent thermal gradients in
the specimens.
The signal was measured using strain gauges in front of and beyond the specimen in the
incident respectively the transmitted bar. The recorded signal consisted of several frequency
components. To make up for the dispersion, a fast Fourier transform (FFT) was used in the
evaluation process where the frequency components were shifted, see Gorman (1983). Passive
damping was used to reduce the elastic oscillations.
3.3. Calibration of material models
Very high strain rates exist in the primary and secondary shear zones while the rest of the
workpiece deforms at low strain rates. The primary shear zone is where the major shearing
of workpiece material occurs, in a secondary shear zone adjacent to the tool–chip interface
shearing due to contact conditions takes place. Thus, tests with both low and high strain
rates for varying temperatures are needed. The calibration of the DD model in Lindgren et al
(2008) was based on uni-axial compression tests conducted at strain rates ranging from 0.001
to 10 s−1 and temperature in the range 20–1300 ◦ C. Therefore, recalibration of the model using
additional high strain rate tests was needed.
Deformation behaviour of AISI 316L steel during machining and optimization of material
parameters for the JC model was addressed in several papers. In the study conducted by
Umbrello et al (2007), the influence of different sets of material constants for the JC model
on machining simulation of AISI316L is discussed. Determination of parameters for the JC
model based on orthogonal cutting tests and analytical modelling is presented by Özel and
Zeren (2004).
The set of stress–strain curves used in this study for the calibration of the JC model and
the response of the model to varying load paths is presented in figure 2.
The JC model cannot give a good accuracy when subjected to the entire test set. Therefore,
only data from the high strain rate tests were used to calibrate the model as shown in figure 2.
Furthermore, the JC model response was never extrapolated outside the experimental range.
It was found that the hardening due to large strains but still low strain rates at the start
of the simulation was too large. This prevented the initiation of chip formation. The chip
formation did not proceed in a similar way as shown in figure 4(b) for the test case no. 2 and
the DD model. Therefore, in the implementation of the JC model into MSC.Marc software,
strain rates larger than 9000 s−1 were set to this cut-off value. The same was done with plastic
strain. This cut-off value was 60%. The use of cut-off means that when evaluating the flow
stress according to (22), the effective plastic strain was never larger than max(ε̄ p , 0.6) and its
p
rate never higher than max(ε̄˙ , 9000). Thus, the actual strain and strain rate in the simulation
were overruled by the cut-off limits when evaluating (22).
The obtained parameters for the JC model that were used for simulations are given in
table 3.
9
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
1400
1200
RT 8345 s–1
RT 2850 s–1
1000
True stress [MPa]
200˚C 2930 s–1
400˚C 2965 s–1
–1
800
500˚C 8670 s
500˚C 1020 s–1
600
–1
700˚C 2965 s
400
950˚C 8810 s–1
950˚C 1075 s–1
200
0
0
0.1
0.2
0.3
0.4
True strain [–]
0.5
0.6
0.7
Figure 2. Flow stress curves used for calibration of the JC model. The lines are computed values.
Table 3. Material parameters for the JC model.
A (MPa)
B (MPa)
n
C
ε̄˙ ref (s−1 )
m
Tmelt (◦ C)
111
358
0.49
0.257
1
0.805
1358
For the calibration of the DD model, the entire test set (both low and high strain rate data)
was used. The response of the DD model to varying load paths over the entire test range
shows acceptable accuracy, see figure 3. The figure shows only half of the calibration set in
order to make the figure readable. The excluded curves show the same quantitative agreement
between model and measurements as those included in the figure. No cut-off was used in the
implementation of the DD model.
The parameters obtained from low strain rate tests, including known or assumed parameters
from the literature, are presented in Lindgren et al (2008). Temperature-dependent parameters
determined by the calibration of the model with SHPB tests are presented in sections A.2
and A.3.
4. Finite element formulation
The implicit FE software MSC.Marc was used together with an updated Lagrangian
formulation accounting for large deformation and strains. The ability of Marc software to
adaptively remesh is crucial in the used approach. Then, it is possible to maintain a reasonable
shape of elements and also capture gradients of strain, strain rate and temperature. This
remeshing strategy does not require a criterion for modelling of the chip separation from the
workpiece.
A staggered step approach for coupled transient mechanical and heat transfer analysis
was utilized. The full Newton–Raphson method was applied to solve the global system
10
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
1400
RT 8340 s–1
–1
RT 2845 s
RT 1365 s–1
RT 0.01 s–1
1200
True stress [MPa]
1000
500˚C 8225 s–1
500˚C 3010 s–1
800
700˚C 2965 s–1
700˚C 1.0 s–1
600
400
–1
900˚C 10 s
900˚C 1 s–1
–1
950˚C 8825 s
950˚C 1075 s–1
200
0
0
0.1
0.2
0.3
0.4
True strain [–]
1100˚C 10 s–1
–1
1300˚C 1s
0.5
0.6
0.7
Figure 3. Flow stress curves used for calibration of the DD model. The lines are computed values.
of equilibrium equations and the nonlinear system of equations for state variables,
respectively.
The friction condition is an important factor that influences chip formation. Friction on
the tool–chip interface is not constant and is a function of normal and shear stress distribution.
Normal stresses are largest in the sticking contact region near the tool tip. The stress in the
sliding zone along the contact interface from the tool tip to the point where the chip separates
from the tool rake face is controlled by frictional shearing stress. A variety of complex friction
models exists; however, the lack of input data to these models is a limiting factor. The model
for tool–chip interface employed in this study is a generalized Coulomb friction model. A
value for the friction coefficient was estimated from the measured force components at the
tool rake angle using Merchant’s theory of chip formation; see Childs et al (2000). The
friction coefficient µ = 0.5 was taken as an average value from the whole range of cutting
experiments.
The heat generated in metal cutting has a significant effect on the chip formation. The
heat generation mechanisms are the plastic work done in the primary and secondary shear
zones and the sliding friction in the tool–chip contact interface. Generated heat does not
have sufficient time to diffuse away and temperature rise in the work material is mainly
due to localized adiabatic conditions. Experimental and modelling efforts were devoted to
determine the amount of dissipated energy during severe plastic deformation since Taylor and
Quinney (1937) presented their work. They found that the fraction of dissipation was typically
in the range [0.8, 1.0] for the materials they studied. Vivier et al (2009) investigated the
amount of dissipated and stored energy in structures containing frictional cracks and elastoplastic zones and stated that the Taylor–Quinney coefficient was not a constant. Stainier
and Ortiz (2010) presented the variational theory of thermo-mechanical coupling that results
in precise predictions of the rate of heating due to dissipation and its dependence on the
strain and the strain rate. A standard practice in FEM simulations of mechanical cutting is to
11
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
Figure 4. 2D plane strain FE model of orthogonal cutting: (a) initial mesh, (b) initiation of the chip.
assume the fraction of plastic work transformed into heat equal to 0.9, e.g. Shi et al (2002),
Özel (2006).
In this study the fraction of plastic work converted into heat was assumed to be constant
and equal to 90%. The heat generated by friction was also calculated and applied as a surface
flux; see also Grzesik and Nieslony (2004), Abukhshim et al (2006) and Dogu et al (2006).
The heat generation conversion factor between energy due to friction and heat generated in a
coupled contact analysis was assumed to be equal to 0.9, default value in MSC.Marc software.
Chip formation is also influenced by the contact heat transfer coefficient. Different values
were reported, such as 1×105 W m−2 ◦ C−1 in Özel (2006) or 50×106 W m−2 ◦ C−1 in Mabrouki
et al (2006). A number of numerical trials were run with varying values of the heat transfer
coefficient to obtain a homogeneous temperature field across the nearest region of the contact.
The value 5 × 106 W m−2 ◦ C−1 was used in the presented work.
An orthogonal cutting operation was employed to mimic 2D plane strain conditions. The
depth of cut, used for all test cases, was equal to 3 mm. Coupled thermo-mechanical plane strain
elements were used for the discretization of the workpiece and the tool, respectively. Fully
integrated bilinear quadrilateral thermo-mechanical elements were used for the workpiece. The
mean dilatational formulation and constant temperature distribution over the element were used
to avoid a problem with incompressibility. The initial mesh was defined by 900 elements, see
figure 4(a). The tool geometry was discretized by 2298 tree-node thermo-mechanical elements.
The dimension of the workpiece in the FEM model was 8 × 1.6 mm. A horizontal velocity
corresponding to the cutting speed was applied to the nodes at the bottom of the workpiece.
The nodes along the horizontal and vertical outer edges of the tool were fixed.
Due to adaptive remeshing, the average number of elements increased up to 14 000. The
effect of remeshing near the tool tip is illustrated in figures 4(b), 5(a) and (b) for the test case no.
2 with the DD model. The mesh adaptivity was controlled by the global remeshing parameters.
Elements having the effective size larger then 0.1 mm were automatically remeshed. The lower
bound for element size was 2.5 µm. This element size is approximately equal to the shortest
tool element edge in the tool–chip contact. Furthermore, the mesh was updated every fourth
increment to prevent the severe element distortion.
The adaptive time stepping procedure with multi-criteria such as automatic
reducing/increasing time step size was used to control the analysis. Initial time step size
equal to 10−10 s was increased to 10−6 s during the steady state part of the simulation.
Material properties for the workpiece material are shown in section A.1. Material
properties of the tool were assumed thermo-elastic.
12
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
Figure 5. Intermediate stages of the chip formation: (a) time 9.251×10−5 s, (b) time 2.891×10−4 s.
Figure 6. Cutting force Fc and feed force Ff —test case no. 2; (a) measured, (b) simulated.
5. Results
Finite element simulations were carried out for the cutting data that were used for the
experiments. Observations on the formation of chips and forces needed to create them are
compared with simulation results in the next section.
5.1. Cutting forces
The measured histories of cutting force Fc and feed force Ff for the test case no. 2 are shown
in figure 6(a).The force values presented in table 4 were evaluated by averaging the measured
signal in the steady state region. The steady state in the experiment was reached approximately
during 3 s after the start of measurements.
The loading histories of simulated forces, see figure 6(b), were evaluated at the chip tool
interface. Average values of the computed forces in the steady state region are compared with
the experimental results in table 4. The error used for the evaluation of the computed results
13
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
Table 4. Measured and simulated cutting forces.
Measured
Simulated JC
Test
no
Fc (N)
Ff (N)
Fc (N)
1
2
3
4
446
960
442
929
437
592
439
560
425
1065
440
1050
e (%)
−4.7
10.9
−0.5
13.0
Simulated DD
Ff (N)
e (%)
Fc (N)
e (%)
Ff (N)
e (%)
350
515
355
510
−19.9
−13.0
−19.4
−8.9
485
1100
480
1075
8.7
14.6
8.6
15.7
400
545
390
530
−8.5
−7.9
−11.2
−5.4
Figure 7. Effective plastic strain rate; (a) the DD model, (b) the JC model.
is defined as
e=
Computed − Measured
× 100%.
Measured
(23)
5.2. Chip geometry and shear plane angle
The formation of the chip involves a shearing of the workpiece material at the shear plane, the
region of progressive plastic deformations extending from the tool edge to the position where
the upper surface of the chip leaves the work surface. The highest strain rates are localized
into a narrow band starting in front of the tool edge radius. They are significantly larger in the
primary shear zone than in the remaining material, see figure 7.
The chip formation was simulated as a continuous process without taking into account any
damage criteria. The chip geometry from the quick-stop test was compared with simulated
results achieved by the DD model and by the JC model. The simulated chip thickness was
obtained by averaging values evaluated at three subsequent sections along the chip edge.
Measured and simulated chip thicknesses are presented in table 5.
The shear plane angle was estimated using strain rate plots according to figure 7. Measured
and simulated values of the shear plane angle are also presented in table 5. The error in the
computed chip thickness and shear plane angle was evaluated using (23).
5.3. Material response
All figures presented in this section correspond to the steady state conditions. The results
shown are for the cutting velocity of 180 m min−1 and feed of 0.15 mm. Results from the DD
model and the JC model are presented side by side. Figure 7 illustrates distribution of effective
14
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
Table 5. Measured and simulated chip thickness t and shear plane angle φ espectively.
Measured
Simulated JC
Test
no
t (mm)
φ
1
2
3
4
0.129
0.198
0.119
0.203
20
23
27
28
(◦ )
(◦ )
t (mm)
e (%)
φ
0.095
0.294
0.111
0.293
−26.4
48.5
−6.7
44.3
34
33
26
32
Simulated DD
e (%)
70.0
43.5
−3.7
14.3
t (mm)
e (%)
φ (◦ )
e (%)
0.102
0.241
0.097
0.240
−20.9
21.7
−18.5
18.2
23
25
31
32
15.0
8.7
14.8
14.3
Figure 8. Temperature distribution; (a) the DD model, (b) the JC model.
plastic strain rates in the primary and secondary shear zones. Figure 7(a) presents a maximum
plastic strain rate value of 68 000 s−1 calculated by the DD model. The distribution of plastic
strain rates having a maximum value of 70 000 s−1 calculated using the JC model is shown in
figure 7(b).
Temperature fields are presented in figure 8. Maximum temperature was generated
in the contact between the chip and rake face of the tool.
Figure 8(a) shows
the temperature field generated by the DD model with a maximum value of 695 ◦ C,
temperature distribution related to the JC model having maximum of 673 ◦ C is shown in
figure 8(b).
The material response simulated by the DD model during the dynamic loading history is
controlled by the state variables. The vacancy concentration and the DD are shown in figure 9.
In the domains of highest temperature concentrated close to the tool rake face, see figure 8(a),
the significant generation of vacancies coupled with the dislocation recovery is present, see
figure 9(a). In the area close to the outer surface of the formed chip with lower temperature
level, the increased DD controls the hardening, figure 9(b). Initial value of DD 1012 m mm−3
was increased up to 1015 m mm−3 .
6. Discussions and conclusions
The aim of the presented work was to compare the predictive capability of the DD model with
the JC model. The scope of the evaluated models was to predict chip formation in terms of
chip thickness and shear plane angle, and cutting forces.
The numerical problems in modelling the cutting process are numerous, particularly the
need for reliable remeshing and contact algorithms. The solution for cutting forces and chip
formation are strongly influenced by criteria for the element size as well as time stepping.
15
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
Figure 9. State variables for the DD model; (a) vacancy concentration, (b) DD.
A number of numerical experiments were run to determine the parameters for a stable solution.
Remeshing of the work material was controlled by a set of control parameters. Tool penetration
into the work material was suppressed. When it occurred, the work material was remeshed and
a new boundary in the tool–chip contact was determined. A minimum element edge length
2.5 × 10−3 µm was used for the remeshing of the workpiece.
Other parameters influencing the chip form and the level of simulated forces are
convergence criteria. Default values in Marc software, residual checking (relative force
tolerance 0.1) together with displacement checking (relative displacements tolerance 0.1)
resulted in cutting force Fc equal to 1155 N and feed force Ff equal to 555 N for the test
case no. 2. The modification of these criteria, relative force tolerance 0.01 and relative
displacements tolerance 0.01, provided 5% change in the cutting force and 2% change in the
feed force. The latter values of convergence criteria proved to be acceptable for sufficient
accuracy.
It can be concluded that using mesh and convergence criteria that capture local effects in
shear zones are important factors for a sufficient valid and accurate model.
Forces were slightly better predicted by the DD model. The feed force was underestimated
by the DD model by about 10%, largest deviation of the feed force predicted by the JC model
was 20%. The largest error in the cutting force predicted by the DD model was 15.7%. The
JC model both overestimated and underestimated the cutting force with a largest deviation
of 13%.
The DD model shows also good prediction capability of the chip morphology in
comparison with the JC model. Largest error of the chip thickness was predicted by the
JC model for the cutting speed 180 m min−1 and feed 0.15 mm. The shear plane angle was
overestimated the most by the JC model by about 70% for cutting speed 180 m min−1 and
feed 0.05 mm. The discrepancy in the simulated results can be related to the fact that the
chip formation was simulated as a continuous process and no strain softening or damage was
included.
The question whether the models are useful for evaluating cutting forces and chip geometry
for a given tool design depends naturally on the context. Sometimes it is sufficient to have
a model that gives a qualitative correct answer. However, during later design stages, then
more quantitative correct results are required. The errors in tables 4 and 5 must be related to
the context where they will be used, namely the cutting tool manufacturing industry. Results
in Lindgren et al (2009) indicate that variations in tool radius in the industrial production of
16
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
nominally identical cutting tool as well as variations in material properties of nominally the
same workpieces can cause variations around 10% in forces. We then estimated that errors
below 10% indicate a model fit for the purpose. The DD model can predict forces near the
wanted precision.
Both material models were calibrated with SHPB test for strain rates up to 9000 s−1 .
It is obvious that the calculated values of the plastic strain rate gained by both material
models are significantly larger than the experimental range. Physically based plasticity
models are supposed to have larger domain of validity and can be extrapolated outside their
range of calibration. But, this requires that the mechanisms controlling the deformation
in the calibration set are still dominating in the extrapolated range. The model presented
in Lindgren et al (2008) was used without any particular mechanisms to account for
very high strain rates. Future work will address implementation of those particular
mechanisms which are significant for high strain rate deformations such as phonon or
electron drags as discussed in Frost and Ashby (1982), Regazzoni et al (1987), Guo and
Nemat-Nasser (2006).
Acknowledgments
The presented work is a part of the project financed by the Swedish Research Council and
AB Sandvik Coromant. The authors would also like to thank Veli-Tapani Kuokkala and Jari
Rämö, Tampere University of Technology Finland, for their technical assistance with the high
strain rate testing.
Appendix A.
Appendix A.1. SANMAC 316L—material parameters
Young’s
modulus E
(GPa)
Poisson’s
ratio ν
Heat
capacity
Cp (J kg−1 ◦ C−1 )
Thermal
expansion
α (1/◦ C)
Thermal
conductivity
λ (W ◦ C−1 m−1 )
Density
δ (kg m−3 )
201
0.3
445
15.5e-6
14
7900
Appendix A.2. Parameters determined by the calibration
Notation
Equation no
Value
Dimension
α
ρi0
2
Initial
immobile DD
13
15
7
7
14
7
7
0.4
1×1012
m mm−3
Kc
cg
τ̂
F
P
q
0.333
1.902
17
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
Appendix A.3. Parameters determined by calibration that are temperature dependent
Temperature (◦ C)
Kc
cg
τ̂
F
20
200
400
600
700
800
900
1100
1300
59.73
32.64
56.09
48.15
45.69
46.69
55.13
121.1
200.0
0.0275
0.0275
0.0275
0.0275
0.0275
0.0275
0.0275
0.0275
0.0275
0.010 390
0.010 000
0.007 392
0.005 643
0.006 256
0.007 939
0.005 409
0.032 410
0.050 000
0.5056
0.8609
0.9117
0.9913
0.9297
0.7789
0.6795
0.6242
0.6615
4.499
7.329
10.58
10.29
13.68
13.54
13.05
14.73
28.62
References
Abukhshim N A, Mativenga P T and Sheikh M A 2006 Heat generation and temperature prediction in metal cutting:
A review and implications for high speed machining Int. J. Mach. Tools Manuf. 46 782–800
Apostol M, Vuoristo T and Kuokkala V-T 2003 High temperature high strain rate testing with a compressive SHPB
J. Physique IV 110 459–64
Bammann D J, Chiesa M L and Johnson G C 1996 Modeling large deformation and failure in manufacturing processes
Theoretical and Applied Mechanics ed T Tatsumi et al (New York: Elsevier) pp 359–76
Bergström Y 1970 Dislocation model for the stress–strain behaviour of polycrystalline alpha-iron with special emphasis
on the variation of the densities of immobile and immobile dislocations Mater. Sci. Eng. 5 193–200
Bergström Y 1983 The plastic deformation of metals—a dislocation model and its applicability Rev. Powder Metall.
Phys. Ceram. 2/3 79–265
Calamaz M, Coupard D and Girot F 2008 A new material model for 2D numerical simulation of serrated chip formation
when machining titanium alloy Ti–6Al–4V Int. J. Mach. Tools Manuf. 48 275–88
Childs T H C, Maekawa K, Obikawa T and Yamane Y 2000 Metal Machining—Theory and Applications (London:
Butterworth-Heinemann)
Dogu Y, Aslan E and Camuscu N 2006 A numerical model to determine temperature distribution in orthogonal metal
cutting J. Mater. Process. Technol. 171 1–9
Estrin Y 1998 Dislocation theory based constitutive modelling: foundations and applications J. Mater. Process.
Technol. 80–81 33–9
Estrin Y 2003 High temperature plasticity of metallic materials Thermodynamics, Microstructure and Plasticity (NATO
Sciences Series vol 108) ed A Finel et al (Dordrecht: Kluwer Academic) pp 217–38
Follansbee P S and Kocks U F 1988 A constitutive description of the deformation of copper based on the use of the
mechanical threshold stress as an internal state variable Acta Metall. 36 81–93
Frost H J and Ashby M F 1982 Deformation Mechanism Maps—The Plasticity and Creep of Metals and Ceramics
(Oxford: Pergamon)
Gorman D 1983 A numerical method for the correction of dispersion in pressure bar signals J. Phys. E: Sci. Instrum.
16 477–9
Grzesik W and Nieslony P 2004 Physics based modelling of interface temperatures in machining with multilayer
coated tools at moderate cutting speed Int. J. Mach. Tools Manuf. 44 889–901
Guo W-G and Nemat-Nasser S 2006 Flow stress of Nitronic-50 stainless steel over a wide range of strain rates and
temperatures Mech. Mater. 38 1090–103
Guo Y B, Wen Q and Woodbury K A 2006 Dynamic material behaviour modelling using internal state variable
plasticity and its application in hard machining simulations J. Manuf. Sci. Eng. ASME 128 749–59
Ivester R, Whitenton E, Heigel J, Marusich T and Arthur C 2007 Measuring chip segmentation by high-speed
microvideography and comparison to finite element modelling simulations Proc. 10th CIRP Int. Workshop
on Modelling of Machining Operations (Reggio Calabria, Italy)
Jaspers S P F C and Dautzenberg J H 2002 Material behaviour in conditions similar to metal cutting: flow stress in
primary shear zone J. Mater. Process. Technol. 122 322–30
18
Modelling Simul. Mater. Sci. Eng. 18 (2010) 075005
A Svoboda et al
Johnson G R and Cook W H 1983 A constitutive model and data for metals subjected to large strains, high strain rates,
and high temperatures Proc. 7th Int. Symp. on Ballistics (Hague, The Netherlands) pp 541–7
Lindgren L-E, Domkin K and Hanson S 2008 Dislocations, vacancies and solute diffusion in a physically based
plasticity model for AISI 316L Mech. Mater. 40 907–19
Lindgren L-E, Wedberg D and Svoboda A 2009 Verification and validation of machining simulations for sufficient
accuracy Proc. Int. Conf. on Computational Plasticity COMPLAS X 2009 (Barcelona, Spain)
Mabrouki T and Rigal J-F 2006 A contribution to a qualitative understanding of thermo-mechanical effects during
chip formation in hard turning J. Mater. Process. Technol. 176 214–21
Meyers M A, Benson D J, Vöhringer O, Kad B K, Xue Q and Fu H-H 2002 Constitutive description of dynamic
deformation: physically-based mechanisms Mater. Sci. Eng. A 322 194–216
Militzer M, Sun W and Jonas J 1994 Modelling the effect of deformation-induced vacancies on segregation and
precipitation Acta Metall. Mater. 42 133–44
Nasr M N A, Ng E-G and Elbestawi M A 2007 Modelling the effects of tool-edge radius on residual stresses when
orthogonal cutting AISI 316L Int. J. Mach. Tools Manuf. 47 401–11
Nemat-Nasser S 2000 Introduction to High strain rate testing ASM Handbook Volume 8 Mechanical Testing and
Evaluation (OH, Materials Park: ASM International)
Nemat-Nasser S, Guo W-G and Kihl D P 2001 Thermo-mechanical response of AL-6XN stainless steel over a wide
range of strain rates and temperatures J. Mech. Phys. Solids 49 1823–46
Özel T 2006 The influence of friction models on finite element simulations of machining Int. J. Mach. Tools Manuf.
46 518–30
Özel T and Zeren E 2004 Determination of work material flow stress and friction for FEA of machining using
orthogonal cutting tests J. Mater. Process. Technol. 153–154 1019–25
Regazzoni G, Kocks U F and Follansbee P S 1987 Dislocation kinetics at high strain rate Acta Metall. 35 2865–75
Satheesha M, Jain V K and Kumar P 1990 Design and development of a quick-stop device (QSD) Precis. Eng.
12 205–12
Shatla M, Kerk Ch and Altan T 2001 Process modelling in machining. Part I: determination of flow stress data Int. J.
Mach. Tools Manuf. 41 1511–34
Shi G, Deng X and Shet Ch 2002 A finite element study of the effect of friction in orthogonal metal cutting Finite
Elem. Anal. Des. 38 863–83
Stainier L and Ortiz M 2010 Study and validation of a variational theory of thermo-mechanical coupling in finite
visco-plasticity Int. J. Solids Struct. 47 705–15
Taylor G I and Quinney H 1937 The latent heat remaining in a metal after cold working Proc. R. Soc. Lond. Ser. A
163 157–81
Umbrello D, M’Saoubi R and Outeiro J C 2007 The influence of Johnson–Cook material constant on finite element
simulation of machining of AISI 316L steel Int. J. Machine Tools Manuf. 47 462–70
Vaz M Jr, Owen D R J, Kalhori V, Lundblad M and Lindgren L-E 2007 Modelling and simulation of machining
processes Arch. Comput. Methods Eng. 14 173–204
Vivier G, Trumel H and Hild F 2009 On the stored and dissipated energies in heterogeneous rate-independent systems:
theory and simple examples Contin. Mech. Thermodyn. 20 411–27
Zerilli F J and Armstrong R W 1987 Dislocation-mechanics-based constitutive relations for material dynamics
calculations J. Appl. Phys. 61 1816–25
19
Fly UP