Novel computations of a moving boundary heat conduction problem applied to EDM technology

Computers them. Engng, Vol. 15, No. 8, pp. -18. Printed in Great Britain.All rightsreserved

OO!I8-1354/91 53.00 + 0.00

1991

CopyrightQ 1991 Pcrgamon Rcsa pk

OF A MOVING BOUNDARY NOVEL COMPUTATIONS PROBLEM APPLIED TO HEAT CONDUCTION EDM TECHNOLOGY M. A.

BARR-T,

M. R. PATFS and P. T. EUBANK

Department of Chemical Engineering, Texas A&M (Received

University, College Station, TX 77843, U.S.A.

23 huly 1990; finaf revision received 31 May

1991; received for publication

11 June 1991)

AWract-A computational procedure is presented to solve a model for thermal conduction which properly describes erosion rates for the anodic material in electric-discharge machining (EDM) processes. The shape and size of melt craters obtained by electric sparks in EDM are calculated numerically by use of a model for thermal conduction. The computational procedure is based on the exact solution of the associated problem of 2-D transient conduction in a semi-infinite body subjected to a uniform time-dependent heat flux over a small circular part of its surf&e. The superposition principle in conjunction with adjustable heat sources is used to cope with the time-dependent boundary conditions. Stability is not a consideration in the present method, which requires less computational time than other methods, particularly for long-time solutions. Because discretization of time or space is not used, the present method can be used to calculate the temperature at any location and time. In contrast, finite-element or finite-difference methods must compute all of the nodal temperatures at each time step. These constraints may lead to schemes that become either inaccurate or unstable, or introducespurious

oscillationsin the solution for a domain of a realisticgrid-size. The methodis in generalmore complicatedto formulateand more tediousto program,but is lesscostly in the numericalimplementationand faster comparedto other numericalmethods.

AN

OVERVIEW

OF THE

EDM

PROCESS

Electric-discharge machining (EDM) has steadily gained importance over the years because of its ability to cut and shape a wide variety of solid materials

difficult to machine by conventional technologies for manufacturing. Carbides, nitrides, ceramics, composites and cermets are examples of materials for which this machining technology has been successful. The removal of material in EDM is associated with the erosive effect occurring when electric sparks take place between two electrodes. Sparks of short duration (l-1000 ps) are generated in a liquid dielectric gap separating the tool and the workpiece. During operation, an applied voltage of =200 V across a typical gap of 40 ,um causes the dielectric to break down. The voltage U falls to a constant value of about 23-25 V, and a constant current lis established for a certain on-time. A plasma channel expands during this on-time or pulse tr. The electrical energy (VI&) delivered by the generator is extremely concentrated because the surrounding dense, liquid dielectric restricts the plasma growth. During this pulse, the high-energy plasma radiates energy to both electrodes followed by thermal conduction and melting within the electrode interior. Scanning electron microscopic observations of EDM machined metal surfaces and of debris or chips formed during the erosion process, prove that the dominant nature of the process is melting, with limited electrode vaporization because of the high plasma pressures (Lhiaubet and Meyer, 609

1981; Benzerga et al., 1985). The anode 6rst melts rapidly due to the absorption of fast-moving electrons at the start of the pulse, but then begins to resolidify after a few microseconds. This is thought to be due to the expansion of the plasma radius at the anode interface which causes a decrease in the local heat flux at the anode surface. Melting of the cathode is delayed in time by one or two orders of magnitude beyond that of the anode due to the lower mobility of the positive ions. At the end of this pulse, a period of pause S begins when power is terminated to the machine. During this time the dielectric is cleaned from the debris formed by the process of erosion. A theoretical rate of erosion can be defined as the volume of molten cavity divided by the sum of the pulse and the pause times for each cycle of operation. HEAT

CONDUCI’ION

MODELS

Several mathematical models have been developed to explain the erosive effect of the sparks. In these models the area of contact between the plasma channel and the electrodes is usually considered as a surface heat source with power and a heat flux independent of time. The melting-point isotherm is calculated at the end of the p&se, and the volume of molten metal is supposedly equal to the material removal per pulse. These assumptions are oversimplifications because the area of contact of the plasma channel is typically transient. In addition,

M. A. BAR~UPET et aI.

610

other models must be developed to relate the molten volume per pulse and the metal removal per pulse (Pate1 et al., 1989). Since the discovery of EDM nearly 4Qyr ago, considerable theoretical and experimental research has been performed to identify the basic physical processes involved. Although thermal conduction and melting are understood to be the dominant mechanisms for erosion of the electrode, the physics of the plasma bubble is not well understood. The dissipated electrical power at the anode and cathode may be represented as heat sources at the surface of the respective electrodes. T’he electrical current also produces joule heating within the electrodes, although the importance of this effect is negligible (Van Dijck, 1973; Comstock, 1959). Therefore, the electrical power in the plasma column, the power at the surface of the anode and the power at the surface of the cathode are key inputs for a general model for thermal conduction. This thermal model must predict crater shapes and volumes in agreement with observations and, on the other hand, the computer time must be minimized. Because of electron emission, the plasma radius at the cathode interface is much smaller than the radius at the anode interface; therefore, at the cathode the heat flux can be approximated by a point source. This model of a point source model at r = 0 in 1-D spherical coordinates is adequate for describing the erosion rate of the cathodic material (DiBitonto et al., 1989). However, in this model the influence of the area of the hot spot is completely neglected, which results in discrepancies between the calculated hemispherically shaped craters and the observed shallow saucer-shaped craters for the anodic material. Therefore, for the anode, more sophisticated models were developed, to include the time-dependent character of the heat source. Several efforts have been made to calculate the removal rate per spark in EDM processes. Some existing models, Zingerman (1959) and Zolotykh (1960), assume a constant heat flux applied to a semi-infinite body. The crater diameter is assumed to coincide with the isothermal melt surface at the end of the pulse because the molten metal is evacuated by the pressure drop and electromagnetic forces when the plasma collapses. Mukoyama (1968) assumed a constant heat flux over a circular disk. Van Dijck and Dutr& (1974) presented a 2-D model although their boundary conditions are different from the ones we used herein. Our objective is to develop more realistic models for the EDM process. Herein the hot spot is considered a circular heat source with time-dependent radius and time-dependent power on the surface of the electrode. EXPANDING

CIRCLE

HEAT

SOURCE

MODEL

(ECHSM)

The governing partial differential equation for the determination of the temperature distribution

T(r, z, t) due to unsteady-state heat conduction into a homogeneous isotropic solid (symmetry in the f? direction) without internal heat production is: 1 i?T @T 1aT --=Yj-p+;ar+p. u at

a2T

Here, CL= K=/&,, is the thermal diffusivity. The associated initial and boundary conditions for the above partial differential equation (PDE) are: IC: BC:

t > 0,

t=O,V(z,r)

T=To

(2)

z = 0, --K=g=Q(r)

forO
z=r=co,

(3)

for r > a(t),

= 0 t > 0,


(4)

T=T,

(5)

where a(t) is a time-dependent circular heat source, To is the ambient temperature of the solid and or the thermal conductivity. Because of the time variation in the radius of the heat source, current density and power are also time dependent. To cope with this feature the principle of superposition was used, with the great advantage of much smaller computer time as compared to finiteelement methods. For a time-independent boundary condition we tried a finite-element heat conduction code called TOPAZ (Shapiro, 1984) however, this package did not handle a time-dependent heat flux boundary which is the appropriate for this problem. To solve the PDE, the following simplifying assumptions were made: 1. The heat flux radius at the anode surface (a) expands at a known rate of growth. 2. Average thermophysical properties of the anodic material apply over the temperature range from solid to liquid melt. 3. The heat of fusion of the anodic material can be neglected for the determination of melt fronts (DiBitonto et al., 1989). Solution to the PDE for a jixed circular heat source As a basis the temperature distribution for a 6xed circular disk source, as obtained from the analytical solution given by Carslaw and Jaeger (1959), was used:

with 1 as a dummy integration variable and X=e-&erf

yz ( --aJolt -+erfc

> (z+A.Jar q/G

>

,

(7)

Jo and 3, are Bessel functions of zero and first-order, respectively, PO = @,,a2 is the effective power of the

Moving boundaryheat conductionproblem energy source, e&(x) is the complementary error function and Q is the radius of the heat source. The axial temperature is provided by a much simpler expression derived as a special solution of equation (6): 2P, J&Y T(0, 2, t, a) = xa2KT x {ierfc(&)-ierfc(s)l-

@)

where ierfc is the integral complementary error function defined as: m ierfic(x)= e&(x) dx = --L e-x’ - xerfc(x). (9) sX fi This is the analytical solution for an infinite cylindrical solid. An order of magnitude analysis of the size of the molten cavities and the size of the anodes and cathodes in EDM operations shows that the latter are several orders of magnitude larger, therefore the infinite cylinder assumption is justified. Traditionally, solutions of the heat conduction equation indicate the time dependence of temperature at a fixed location in space. An alternative solution of the problem, particularly for moving boundary problems, is to calculate how specified temperatures move through the heat conducting medium at specified times, using an “isotherm migration method” (Crank, 198 1). Temperature determinations for a fixed circular heat source Equation (6) is the exact solution for the temperature distribution of the heat conduction equation with a tixed circular heat source as a boundary. At first glance it appears to involve a rather complex

611

numerical integration, which has to be executed in order to evaluate temperatures. In practice, it would be prohibitive if not impossible in terms of computer time and/or round-off error. A series of approximations are proposed, valid for most engineering applications, that allow temperature evaluation while minimixing the computer time by drastically reducing the number of functional evaluations. As part of the overall practical strategy, the Bessel functions J,(x) and J,(x) were computed in terms of rational functions of polynomials in the argument (x). The coefficients of the various rational functions and polynomials are given by Hart (1968) for various levels of desired accuracy. An accuracy of 10m6 was used herein. The wmplementary error function appearing in equation (7) was evaluated with a subroutine from Press et al. (1987), which uses a Chebyshev fitting polynomial and provides erfc(x) with a fractional error everywhere less than 1.2 x 10 - ‘. With subroutines for the Bessel functions and the complementary error function, the integrand of equation (6) can be easily evaluated for auy values of time (if), radius (r), flux radius (a) and depth (2). To evaluate a temperature, integration is required between 0 and OZJ.For practical applications this integral can be truncated by using an acceptable upper limit. The behavior of the integrand was analyzed for several values of the variables (r, u, t, z) and the function was observed to have a decaying oscillatory behavior. Figure 1 illustrates the oscillatory behavior of the function in the least favorable case in terms of truncation errors, that is at the surface (z = 0). Since J, and J, are Bessel functions of zero and first-order, respectively, PO= &,,a2 is the effective power of the

a Fig. 1. Example of the otillatory behaviorof the intetpandof equation(6) at the surface.

M. A. BARRUP~ et al.

612

the computed temperature is proportional to the area under the curve, one may assume that after the second or third zero of J,(X) = 0 or J,(x) = 0 (whichever is first) the areas above and below would compensate. In the program control statements were included to adjust the upper limit of the integral in equation (6) according to an order of magnitude analysis that takes into account: (a) the value of the integrand itself at the upper limit; and (b) the zeros of the Bessel functions. The upper limit in the integral (;I = 2”) was calculated as: a

m_

min

-

zero”(Jo 1 zero”(Jl 1 r ’ a 1’ [

(10)

where n is the nth zero or root of the Bessel function. After a series of calculations, the variation in the molten craters due to extension of the interval of integration beyond the second or third zero of the Bessel functions, Jo and J, , respectively, was found to be less than 1% . Therefore, truncation after the second zero was sufficient for most engineering purposes. The integration was carried out using GaussLegendre quadrature. Control statements are included in the master program to decide on the number of quadrature points to be used, and a subroutine provides the weights and the abscissas accordingly. The evaluation of the axial temperature is straight forward. The integral of the complementary error function can be evaluated directly from the complementary error function using equation (9). So far this is only the evaluation of a temperature for a given set of space coordinates, spark duration t and for a fixed circular heat source. For a timedependent heat flux radius, the principle of superposition, based upon the individual solutions of tixed heat sources of different sizes, can be used (Carslaw and Jaeger, 1959, p. 266). The superposition source

principle

for

an expanding

heat 6E;5;6E;

For discs of radii that change with time, the temperatures can be determined by superposition (Van Dijck, 1973). Here the temperature contributions are summed over discrete disc-sizes a, (a at time 2,) in the corresponding time interval t,to t,, , as:

T[r, z, tf, a(r)1 =

5 (Tb-, z, tf -

h, athll

i-L

The “on-time’* t, is usually divided into n equals intervals (tf = n x At) and a disk size is calculated for each one of these intervals according to its growth formula. One must have a better physical understanding of the plasma dynamics and the plasma geometry during plasma formation and collapse in order to supply an accurate boundary condition to this model. Altema-

I

I

I

I

I

Ej8 8888 I +

++++

8888888686E ++++++++++i

Moving boundaryheat conductionproblem

613

first differences were found to be very small, but then start to increase considerably as the 6nal time is achieved. Because of that, the order of summation was reversed to have control over these incremental differences by increasing the time interval. The last incremental difference in the previous scheme, the disk radius g calculated at ff, be-comes the first in the reversed summation, and this is set as the allowable minimum difference to have along the summation. Therefore, equation (11) can be rewritten as:

tively, these equations may be decoupled by proposing physically realistic radius growth equations and wmparing the simulations to available experimental data. This procedure was used successfully and the calculated craters and erosion rates compared well with experimental data (Pate1 et al., 1989). The time dependence of the heat source radius a, = 0.788 x t’14 was chosen based upon high-speed photography experiments which correspond to EDM conditions (Robinson, here li=i x At. 1967), Equations (6) or (8) were used to evaluate the temperature for this heat source radius a, at two different times: r, - t,andtft,,,. The re.sul&were substracted as indicated in the above equation. If this process is repeated for -n incremental differences, (2 x n - 1) functional evaluations of either equation (6) or equation (8) are required. The above computations require considerable computer time, particularly when the time intervals are small or the number of discrete disc sizes used is large. The behavior of equations (6) and (8) was analyzed for fixed time increments (At = q/n = 0.01 ps) and the

Ttr, 2,

tf.

a(t)1 = 2 {T,[r, 2, t,, ah-i-l

-

?+,[r,z,

tr+l.

a,-

till m.

Table 1 shows how the temperature evolves from this summation using the adjustable interval (At”) scheme. This time interval is adjusted by a prespecified factor in order to keep a uniform rate of growth of T. In Table 1, tatis the time used to calculate the heat-source radius a(r), which starts at tf and it is successively decreased by At, until it reaches the value

Read Parameters and variables

Initialize

dif = big tsup L 0. tleft = tfinal

AtV = Aii I

Do whlk (Atv.lt.tfinal) Do while(dlf.gt.small.and.tleft.gt.Atv) = 0.788’tleft”(3/4) ;i + 1) = tfmal - theft+ AW = tstore tleft = tleft-Atv

( dif too small reset to iam

increase At@

dif = big

AW =

times*

enddo

Add contribution

&v

(rime.9

value and

is a prespecified mUftip#8r)

for last time increment

a = 0.788’Ati”(3/4) Calculate Temperature for this radiususing either Equations (6) or ( 8) at the following times:

t(i+l) = tfinal - wore t0 ) Atv = tfinal - t(i)

Fig. 2. Flow diagramwith the major

(12)

$Z$+j

this is the last value of t(i)

computingstepsusingthe principleof superpositionwithadjustable intervalsof time.

614

M.

I

I

A. BARR-

et al.

I

.Ol

.1

Fig. 3. Variable

Time

distribution

10

‘@s)

of interval and associated

‘5 frequencies.

1.0

0.8 8 _ S m

t 0.6

-

2 : c

0.4 -

i .s t

0.2 (a) 0.0

’

0

100

200

Pulse Fig.

4. Temperature

differences

calculated

Time

300

400

(p8)

from: (a) an adjustable; superposition.

and

(b)

a fIxed interval

in the

Moving

boundary heat conduction problem

of O.Ol-the &ted At used for comparisons. NC.is the number of functional evaluations required. Figure 2 is a flow diagram indicating the major programming steps designed to accomplish these calculations. Figure 3 shows how the pulse time (15 ps) is distributed and it is a complement of Table 1. This on-time or pulse is distributed as: tr=

5

i-l

Ati XL,

(13)

whereJ; is the frequency of interval of size Ati and N is the number of different sixes used. Figure 4 indicates the absolute difference of the temperatures calculated using adjustable time intervals and tixed time intervals. The powers used in this figure are typical of EDM conditions. This scheme would be expected to underestimate the temperatures (the largest difference found was 1 K); however, one should not take the fixed-interval solution as more exact that the solution for an adjustable interval. Since computers carry numbers with a tlnite number of digits, the exact solution is not known and substantial round-off error may accumulate. Figure 5 shows the dramatic reduction in the number of functional evaluations from using the adjustable interval procedure over the fixed interval scheme. The number of functional evaluations [Nreis calculated as the number of times that equation (6) or equation (8) have to be solved and it is not related to the degree of difficulty in solving these equations]. For instance, N, is higher for the axial, r = 0, and for the radial direction (z = 0), than for a nonzero value of z and r. However, the axial temperature (r = 0) is simpler to calculate because it does not require a numerical integration unlike equation (6), and z = 0 (the surface) also simplifies somewhat the evaluation of equation (6).

10

15 Pulse

Fig. 5. Comparison

20

615

RESULTS OF MODEL

CALCULATIONS

The procedure to evaluate temperatures for a time-dependent heat flux has been outlined above. This method was applied in the evaluation of meltisothermal locii for different input parameteM.g. thermal properties, pulse times and powers characteristic of spark-machining processes. To find the coordinates (r,,, , zm) at which T = T,, a bisectional method was used for the axial direction [equation (8)]. To reduce the number of temperature determinations two key points of the melt-front are first calculated: one along the z-axis using equation (8) and the other along the radial coordinate [T(r, z = 0, t)] using equation (6). These two locations can be joined with a straight line. The equation for this line: z=z,----XT,=a? r,

(14)

provides a first estimate for the melt front in the (r, z) plane along uniformly spaced vectors in the radial direction. A number of points chosen along these vectors where the temperature is calculated. The desired isothermal locus, the melting locus in this example (z, , rm), can be determined by interpolation, least squares or bisectional routines. Interpolating polynomials were used which provided results very close (< 1% differences) to those of rigorous bisection at reduced computer CPU cost. Figure 6 indicates the interpolated Z~ location given a set of temperatures in the neighborhood of T = T,,, . The variable heat flux solution covered ranges of power input ranging from 1 to 500 W, with no limitations of convergence or stability for any pulse time used except for unrealistic “on-times” that would correspond to a completely resolidified

0

Flxed Step

0

(r=O,z=lO)pm

3

(r=lO, 240) pm

g

(r=lO,z=O)pm

25

Time (ps)

of the number of functional evaluations using a fixed or a variable interval.

M.

616

A. BARRWET et al.

8

8

T+mperature(K)

Fig. 6. Interpolated

50

locations

of the melt front for copper.

Mett Fronts for Copper (Advancing Power = 130 W

1

,-

1199

1CKKJ

6w

Front)

-

1-r riny l--w

-

IOJ

1I

/ .

(8)

O-' 0

40

20

(urn)

Radius

Melt Fronts for Copper (Receding Power I 130 W

Fronts) .......... .,_._.I._.

tr12op8 tr2oojl8 t _ 240 w t = 300 cu

-

t=560&4S

.. .. 3

-

t=8w

30

N 20

10

,.,

_-

--

40

PO

Radius

Fig_ 7. predicted advancing

‘

(pm)

(a) and receding (b) melt-fronts for copper as a mnctlon a typical power of EDM proccsscs.

OI pulse time for

Moving boundaryheat conductionproblem

617

Predicted Eroslon Rates for Copper

-

A

IncreasIng power 10

0

20

30

40

Pz4oW

P=MW

. . . . *. . ..I _._._._.I.

p_ggw p= ,,gw

-

P=136W 50

60

70

nme (Pa) Fig. 8. Theoreticalrates of erosion for copper materialfor differentpower inputs. material. Figure 7 shows the advancing and receding calculated melt-fronts for copper as a function of time for a typical power input in EDM operation. The melted volume first increases with time as well as the radius of the source. After a “critical value” of this radius is allowed, the melted volume starts to decrease or resolidify. For constant dissipation of total power, the heat flux decreases as the radius of the plasma channel grows. At some instant, the heat loss from the liquid-metal zone by conduction to the surrounding solid material exceeds the heat supplied to the liquid. Consequently, the liquid metal in the deepest part of the crater will solidify and the crater shape will become shallower. If different materials are compared under the same power inputs, those with better thermal conductivity of the solid material will produce this phenomenon sooner. This explains the low relative wear of copper electrodes when EDM operations are performed on steel. The only model able to mimic this behavior is one with an expanding heat-source boundary. The area under these melt fronts is calculated numerically and further is rotated to calculate the volume of the molten cavity. Figure 8 shows the theoretical rates of erosion for copper under different applied powers typical of die-sinking EDM operation. These rates of erosion were calcualted as: t’,=

v, [ t+a

1 ’

05)

where V, is the volume of melt cavity, t is the pulse time and 6 is the pause time, or the time allowed between sparks. It can be calculated from the frequency at which the machine operates. In the operation of a wire EDM machine the work piece is the anode and one should operate at pulse times closer to the optimum to obtain maximum material removal per unit time. For die-sinking operation the polarity is reversed and the tool electrode is the anode, therefore one is interested in minimizing the wear. This is accomplished by allowing the spark

to stay longer which would favor the resolidification of the material. CONCLUSIONS

A series of numerical approximations was proposed, to evaluate isothermal locii from a heat conduction model with a time-dependent, heat-flux boundary condition. This strategy, valid for most engineering applications, minimizes the computer time by reducing considerably the number of functional evaluations. This procedure can be used to calculate temperature profiles and stock removal rates for different applications of spark erosion as in EDM technology and laser welding applications. The present expanding circle heat source model (ECHSM) provides qualitatively correct curves of erosion rate and explains the low rates of erosion that result due to resolidification of the anodic piece for die-sinking EDM during long on-times. At the same time, it also quantifies the rationale behind operating the wire machines with very short on-times. The model is capable of describing all the qualitative features associated with erosion of the electrode, e.g. resolidification of the melt cavities which eventually results in relatively minimal erosion for the long pulse times associated with die-sinking machines_ For short pulse times associated with wire machines, the model is unique in terms of the correct qualitative prediction of the maxima known to exist in anode-erosion curves. NOMENCLATURE n, = C = ECHSd = EDM = I= J, = J, = PDE = P,, =

Radius of the plasma at the anodic surface Heat capacity Expandingcircleheat source.model Electric-dischsrge machining current Besselfunctionor order zero Resselfunctionof order one partial Uferential equation F%lbctive Powerdelinedin equation(5)

M. A. BARR-

618 4 = Heat flux r = Radial coordinate t=Time (,us) T= Temperature U = Voltage V, = Volume of eroded cavity 3, = Rate of erosion (mm3 min - 3 X= Defined in equation (8) z = Depth or axial coordinate Greek letters 01= d= A = JC~= 1= c = p =

Thermal diffusivity Pause time, or time between sparks Interval Thermal conductivity Dummy variable defined in equation (7) Dummy variable defined in equation (5) Density of solid material

Subscripts A, a = C = i, i + 1 = v= m = 0 =

Relative to the anode Relative to the cathode Time counter variable Melting conditions Initial time or at ambient conditions

Acknowlea’gements-We acknowledge the financial support of the State of Texas Advanced Technology Program through the Grant No. 4437 and the technical advice received from a number of physicists and engineers. In particular, we thank Dr Ivan0 Beltrami of AGIE and Professors A. M. Gadalla, J. C. Holste and K. R. Hall of Texas A&M University. REFERENCES

Benzerga L., C. Lhiaubet and R. M. Meyer, J. appl. Phys. !b3, 604 (1985).

et al.

Carslaw H. S. and J. C. Jaeger. Cot&&on of Heat in Solti, 2nd Edn, Chaps 12, 13. pp. 214, 264, eqns (4) and (5). Clarendoa Press, Oxford (1959). Comstock W. D.. Current distribution in the cylindrical source p~e-el&rode configuration. Commun. El’eetron. Jul, (1959). Crank J., How tc deal with moving boundaries in thermal problems. Numerjcal Methodp in Heat Transfer, Chap. 9. Wiley, New York (1981). DiBitonto D. D., P. T. Eubank, M. R. Pate1 and M. A. Barrufet, Theoretical models of the electrical discharge machining process. Part I: a simple cathode erosion model. J. a&. Phys. 66, 40954103 (1989). Hart J. F.. Computer Approximations. Chap. 6.8. p. 141. Wiley, New York {1%8). Lhiaubet C. and R. M. Meyer, J. appl. Phys. 52, 3929 (1981). Mukoyama Y., The mechanisms of EDM. Bull. Jap. Sot. Precir. Engng 2, 288 (1968). Pate1 M. R., P. T. Eubaak, D. D. DiBitoatc and M. A. Barrufet, Theoretical models of the elect&al discharge machining process. Part II: the anode erosion model. J. appl. Phys. 66, 4104-4111 (1989). Press W. H., B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, Numericat Recipes, The Art of Scienti#c Computing, Chap. 6, p_ 164. Cambridge University Press (1987). Robinson J. W., J. appl. Phys. 38, 210 (1967). Shapiro A. B.. A finite element heat conduction code for analyzing 2-D solids. Report UCID-20045, Lawrence Livermore National Laboratory (1984). Van Dijck F., Ph.D. Dissertation, Catholic University of Leuven, Belgium (1973). Vank Dijck F. and W. L. Dutr& J. Phys. D: Appl. Phys. 7, 899 (1974). Zingerman A. S., Regarding the problem of the volume of molten metal during electrical erosion. Sov. Phys.-Sol. 1, 225 (f959). Zolotykh S. N., Sov. Phys.-Tech. Phys. 4, 2 (1960).