CAD of ECM tools

CAD of ECM tools V K Jain and P C Pandey* The last two decades have witnessed a rapid development in the appfication of electrochemical machining (ECM) techniques in production industries. However, there are tool design problems which hinder its application. Available analytical methods for the design of ECM tools do not yield satisfactory results. In this paper, the finite-element technique is used which has proved superior to all other techniques. Also, details are given of the method employed for the computation of metal removal and gap in ECM. Relevant computer programs are included. Due to rapid development and continuous research, electrochemical machining (ECM) is finding applications z-3 in sophisticated industries such as aerospace, nuclear and in mass production jobs like automobile manufacture. The popularity of ECM is mainly due to its ease of adaptability to different types of machining operations and production of stress-free surfaces with superior finish. One of the problems which hinders free application of ECM is obtaining an accurate tool design. A number of models 4-1° for analytical design of ECM tools have been developed, but these have had little success. This is because the complex problem 1] of ECM has been attempted using over-simplified assumptions. Tools for use in ECM are often designed by 'trial and error'. Finite-element techniques have been used to analyse metal removal and generation of interelectrode gap (I EG) in ECM 12,13. This demonstrated the superiority of finite elements over the finite difference technique. This paper presents the detailed computational aspects of metal removal and gap generation in ECM, using finiteelements and discusses the development of relevant computer programs.

electric field potential density machining efficiency

p ~/

Subscripts

e i, j . . . m o t

electrolyte, equilibrium node numbers workpiece material initial condition condition at machining time 't'

METAL

REMOVAL

IN ECM

Figure 1 shows a schematic diagram of the EC cutting operation. The amount of material removed in time t, can be obtained from the Faraday's law of electrolysis

m = rl EIt/F

(1)

This is valid under idealized conditions only. In practice, it yields inaccurate results because of joule's heating and consequent changes in the electrolyte conductivity. Further, equation (1) does not account for non-uniform heat generation and temperature distribution within the gap, variations in current density whithin lEG, liberation of gas bubbles and the presence of the reaction products; this makes analytical computation difficult. The electrolyte conductivity in ECM, for a given temperature and void fraction is given by:

K = K o (1 + ~ A T)(1-~v)n'

(2)

whereas, the temperature rise of the electrolyte can be computed from:

NOTATION

a,b,c E Ev F FF I J K m n n' t T V Y ~v

constants electrochemical equivalent of work material effective voltage Faraday's constant feed rate function current density electrolyte conductivity metal removed normal to the work surface exponent machining time temperature electrolyte flow velocity inter-electrode gap (lEG) temperature coefficient of electrolyte conductivity void fraction

Mechanical EngineeringDepartment, Motilal Nehru Regional EngineeringCollege, Allahbad-211004, India *Mechanical and Industrial Engineering,University of Roorkee, Roorkee-247672, India

volume 12 number 6 november 1980

dT/dx = j 2 I (K V p e %)

(3)

The potential distribution within the working gap, for the case of two dimensional machining problems, can be obtained from the Laplace equation 13,14: a2~ - -

a x2

a2~ + -= 0

(4)

ay 2

and using the boundary conditions, ~ = 0 at the cathode and ~ = E v at the anode. The current density at the work surface, within the lEG can be obtained from: J = K (a~lan)

(s)

The field vector ~ should be determined 13,1~-19 such that it minimizes the function I(~) and satisfies the boundary conditions above:

I(~) = ~

[(aelax) 2 +

(a~/ay) 2 ] dx dy

0010-4485/80/060309-07 $02.00 © 1980 I PC Business Press

(6) 309

For computation, the lEG has been approximated by triangular elements (Figure 1) which uniquely and continously represent the field variable ~b(x,y) within the solution domain and varies linearly within the elements as follows:

¢(x,y) : N,.,; + Nj /+

: [Nol {¢ o}

I

Readand print inputdata I t J Computeno of computation cycles NC ! [ GenerateK,Vand Yfor different nodes J t JCompute feed rate at every nodeI t J Determine band width MB j I [ Initializeparameters J =i Generateand assemble element matrices !

(7)

and N e is defined as,

N B =a B +bBX+CBY (B=i,j,k. . . . )

Substituting equation (7)in (6) and differentiation with respect to the nodal points, yields a set of equations which can be written as2°-22:

I I

Form ,~obo, mat..

I

i Applyboundaryconditions ! Solve set of eqns to determine _9

-I I Compute~

I

I

I Compute T

Coefficients of the stiffness matrix (Kij) etc., could be evaluated from:

K/i = (b/bj

+

c/cj)/(4 ,~e)

I

Compute heat transfer to air t J Compute decrease in electrolyte temperature

(9)

Assembly of element equations [ of the form (8) ] would give the global matrix. The set of sparse, banded and symmetric matrices thus obtained has been solved by Gaussian elimination to obtain the potential distribution within the lEG. Knowing the distribution of ¢, the values of J, K, Tetc at each node located on the workpiece are evaluated as discussed above. The value of the lEG is computed for the conditions of finite feed and zero feed respectively from (10) and (11 ). (10)

Vt : (2Ct + Y2o)½ where,

t

I Compute K,Y, MR/,MR v and MRs

q

Reass=g.n x - y

outward flow

Reassign x -y side flow Yes

coordinates inward

-! J

I El(~'Wolyte'--~---

T

Tool

Node no =/- 21 Element no =(]~-~

b

Z

0.7

io

L3

1.2 1.7 Distance, mm

~ 19

2.2 2.4

Figure l.(a) Principle of electrochemical machining process (b) two dimensional discretization of inter-electrode gap

310

ComputeV for new value of Y

J

No

I I

L

Work

t 4_

]

Yes

(11)

For computation is of Yt, equation (12) is less time consuming and is more accurate than equations (10) and (11 ).

0.0 0.2

]

JCompute effective electrolyte temperatureI

re = f l E E r K/(F PmFF) = ClEF

I

J

Yes

1 Yo- Ye t = FT [ Yo-Yt+Ye In ( ~ ) ]

a

1

t

Print resulte

I

Figure 2. Scheme of computation employed in mode/ SGFET-22 Table 1 shows that most of the experimental temperature values1~ are closer from (12) than from (10). However, other ECM parameter values are about the same from both (10) and (12).

C '= ~?J E (FPm) where,

(12)

C'= ~J E (FPm) Linear MRR at each node and the volumetric and specific metal removal rate within each element can be computed from:

MRI(i) = At F F + Y(i)-FPY(i); If KIN=I,FPY(i)=YI (13)

computer-aided design

MRv(i ) = W[(MR/(i) + MR/(i+I)) MR,(;) : [MR (i) /

+

x(i+ l ) - x(i) (14)

2

J(;t ]

The ratio of the actual I EG to the corresponding equilibrium gap at each node is given by:

CYte = Yt / Ye

(Is)

(16)

Table 1. Comparison of equations (10) and (12) with data for anode shape (t = 300s) x mm

Y(mm) x 10-1

Temperature rise (°C) Eq(10) £q(12) yo=0.8 Yo=0.8

MR/(mm)

K(~,mm) -1 x 10"1

J(A/mm 2 ) (10)

(12)

(io)

(12)

(10)

(12)

5.0 11.625 15.0"

0.3709 0.3706 0.3711

0.3712 0.3704 0.3708

0.1560 0.1799 --

0.1561 0.1799 0.1933

0.4408 0.4415 0.4413

0.4410 0.4410 0.4410

5.3100 15.0073 20.4389

15.8687 22.45*

0.5034 0.5023 0.5024

0.4854 0.5060 0.5039

0.1573 0.1834 0.1982

0.1567 0.1840 0.1983

0.6000 0.6000 0.6000

0.6050 0.6058 0.6030

6.3067 16.6175 22.5336

10.0875 23.6*

0.4687 0.4676 0.4679

0.4687 0.4675 0.4680

0.1590 0.1881 0.2048

0.1590 0.1881 0.2048

0.5640 0.5640 0.5640

0.5640 0.5639 0.5641

-

9.9375 22.5 26.0*

0.4994 0.4978 0.4985

0.4995 0.4977 0.4988

0.1622 0.1974 0.2181

0.1622 0.1974 0.2182

0.6089 0.6090 0.6090

0.6090 0.6089 0.6095

--

10.2375 23.0 28.4*

0.5015 0.5001

0.5016 0.4997 0.5008

0.1612 0.1943 0.2136

0.1612 0.1943 0.2137

0.6089 0.6090 0.6090

0.6099 0.6084 0.6093

(10)

(12)

7.0 17.0 22.0

7.8500 9.0419 9.7153

7.8477 9.0452 9.7192

13.7024 18.4565

5.2623 13.6954 18.4422

7.0 17.0 22.0

3.4808 4.0584 4.3841

3.3670 4.0890 4.3980

5.7189 14.9660 20.1821

7.0 17.0 22.0

5.0027 5.9176 6.4426

5.0018 5.9177 6.4434

6.3056

16.6205 22.5309

7.0 17.0 22.0

5.9083 7.1881 7.9424

5.9097 7.1877 7.9452

7.4490 19.9007 27.2433

7.0 17.0 22.0

4.6961 5.6616 6.2235

4.6955 6.6587 6.2265

7.0750 18.8207 25.6587

5.2595

Experimental

0.5005

* x = 21.5 mm

Table 2. Machining conditions 1

FF(mm/s )

2 0.011

4 0.011

5 0.011

7 0.015

8 0.019

12 0.019

13

0.02

15

0.011

17

0.011

18

0.015

19

0.015

0.019

Ev(V )

16.0

12.0

20.0

16.0

12.0

16.0

20.0

12.0

20.0

16.0

20.0

12.0

Q(mm3/s) x 103

31.67

25.0

35.0

25.67

21.5

24.17

26.33

44.5

56.17

49.67

56.83

61.17

C = 4.18 W/g°C F = 96500.0 coulombs

Ko = 0.01412 (Qmm)-I To = 30°C

A t = 30s t = 300 s

Yo = 0.8 mm

/, I I

I I I

/

/

/

E

/

/

/J

I I

:="

/

K)

®

2O

"0

I0

20

'ZJ I0

I

.,,,..~

5

I

CIlia,"-

I I0

0

I 20

I

c 0

IO Distance, mm

Distance, mm

. . . . V, m/s • ?..65 o 4.49

FF, mm/s 0 . 0 0 2 0 x 10 0.0011 x I0

Experiment ~

V, m/s O 5.66

• 5.73

J

0 v 0

20

"~

I IO Distance, mm

I 20

SGFET-22

FF, rnm/s 0.0011 x I0 O.O015x I0

V, rn/s x 5.01

FF, mm/s 0.0OI5 x I0

o 6.17

0.0019 x I0

16

Figure 3. Comparison of theoreticel (model SGFET-22) end experimental results volume 12 number 6 november 1980

311

COMPUTE R APPLICATIONS

This validates the analysis and the computer program. Figures 4 and 5 show the computed results pertaining to lEG, current density, electrolyte conductivity, linear and volumetric metal removal rates and the ratio between lEG and equilibrium gap. Current density is high at the corners of the tool as compared to that at the centre (Figure 4). Due to this, a larger lEG is obtained at the exit end. However, at inlet the effect is not so pronounced because of the nature of the boundary conditions assumed. However, for machining time in excess of 240s, the lEG tends to become constant and virtually equilibrium conditions are reached. A similar trend can be seen in Figure 5.

Computation of ¢, J, K, Yt etc, is discussed in Figure 2 through flow chart and computer programs (Appendix 1). Subroutines related to READ, WRITE or/and PRINT, DIMENSION, COMMON etc., have been omitted.

RESULTS A N D DISCUSSION Figure 3 shows a comparison between the experimental and computed (model SGFET-22) values assuming an anodic dissolution efficiency of 100 per cent (other mach in ing conditions are given in Table 2). However, in actual practice machining efficiency varies between 75 and 100 per cent. In a few cases (Figure 3) some discrepancy between the computed and experimental temperature values can also be noticed. This is due to defects in the temperature measurement, ineffective filtration of the electrolyte before recirculation, and the assumption for ECM efficiency. However, there is a good correlation between the computed and analytical results (Figure 3).

o.ozz a

CONCLUSIONS The computer method for the analysis of the ECM process discussed yields results which are close to the experimental data. However, the accuracy of the results is decided by the correctness of the assumed values for the anodic dissolution efficiency and other parameters for which little data is

available.

°~'b

7-. 0.020

o.851-

/~

\

y.

oo,8

=

y

0.016

O.

/

-

~E

./

~ 0.75

036-

" "--

••~

""

i

.

/ /

0.65

0.34

/ / /

-

oo,,~ • D

I0

5

Distance,

, 15

, 20

G 25

05

~o

mm

5

032n

"" ~

_= ~

5

Distance, mm - ......

30s 60s

90s ........ 120s - - o - - - o - - - 150s

= x

. . . . .

I

,;

,5

~o

G

25

Distance, mm = x

210s 240-300s

Figure 4. Variation of (a) electrolyte conductivity (b) inter-electrode gap and (c) current density along electrolyte flow direction a

C

b

1.2

.

\ /

%(3

4.S

x=0.0

,=,2o

'?, o 4.~=

I.I

5.0 m

1.0--

~

0.9

4.1

~-.0

3.7

1.0

® O.voR~

I 60

t I 120 180 Machining time, s

I

240

~

300

3.$

I 5

I K)

I 15

I 20

r ~" 25

I 60

Distance, mm - .......

30s 60s

. . . . . .........

90s 120s

I

i

120

180

,

240

F."

M a c h i n i n g time, s

~ ~

x

210s 240-~K)Os

Figure 5. (a) Relationship between machining time and ratio Yt/Ye (b) variation in penetration along the electrolyte flow direction (c) relationship between volumetric metal removal and machining time

312

computer-aided design

Thus, the computer method of ECM tool design would be successful only when a large amount of experimental data concerning the process is available.

REFERENCES 1 Wilson, J F Practice and thory of electrochemical machining Wiley I nterscience, NY, USA (1971) 2 0 p i t z , H 'Electrical machining processes'proc. Int. Production Eng. Res. conf. (1963) pp 225-231

reproduction in ECM' Proc. 19th Int. MTDR conf. (1978) pp 533-540 11 Jain, V K and Pandey P C 'On some aspects of tooling design in ECM'Mech. Eng. Bull. Vol 8 No 3 (1977) pp 66-72 12 Jain, V K and Pandey P C 'Design and analysis of ECM tooling' Precision Eng. Vol 1 (October 1979) pp 1 9 9 206 13 Jain, V K and Pandey P C 'Finite element approach to the two dimensional analysis of electrochemical machining process' Precision Eng, Vol 2 No 1 (January 1980) pp 23-28

3 Krabacher, E Jet al 'Electrical methods of machining' ibid pp 232-241

14 Tipton, H 'The determination of shape of tools for use in electrochemical machining' Report No 40, MTI RA UK (1971)

4 Ki~nig, W and Pahl D 'Accuracy and optimal working conditions in ECM' Annals CIRP Vol 18 (1970) pp 223-230

] 5 Jain, V K and Pandey P C 'Anode profile in ECD' (to be published)

5 Hopenfield, J and Cole R R 'Prediction of the one dimensional equilibrium cutting gap in electrochemical machining' Trans. ASME J. Eng. Ind. (1969) pp 755-765 6 Kawafune, K et al 'Effect of the working parameters on the working gap during electrochemical machining' Annals CIRP Vol 18 (1970) pp 305-317 7 Tipton, H 'The dynamics of electrochemical machining' proc. 5th Int. MTDR conf. (1964) pp 509-522

16 Bhatia, S M 'Effects of electrolyte conductivity on ECM process' M.Tech. Dissertation, I IT Bombay, India (1971) 17 Zienkiewicz, O C The finite element method McGraw Hill, UK (1977) 18 Desai C S and Abel, J F Introduction to the finite element method Affiliated East West Press, New Delhi, India (1977) 19 Bathe, K J and Wilson E L Numerical methods in finite element analysis Prentice Hall, New Delhi, India (1978)

8 McGeough, J A and Rasmuseen A Theoretical analysis of electrochemical machining process' ASME paper No 69-WA/Prod-10

20 Segerlind, L J Applied finite element analysis John Wiley and Sons, NY, USA (1976)

9 Collet, D E e t al 'A complex variable approach to electrochemical machining problems' J. Eng. Math. Vol 4 No 1

21 Martin, H C et al Introduction to finite element analysis Tara McGraw Hill, New Delhi, India (1975)

(1970i pp 29-37 10 Larsson, C N e t al 'Electrochemical effects on shape

22 Huebener, K H The finite element method for engineers John Wiley & Sons, UK (1975)

APPENDIX I. NOTATION AND PROGRAM

IBC

ABC

Subroutine related to the application of boundary conditions AGCSG Subroutine that automatically generatesthe coordinates of global stiffness matrix A~ Electrolyte's electrical conductivity AKI Initial electrolyte's electrical conductivity AJ Current density ALF Temperature coefficient of specific conductivity of the electrolyte AMRR Linear metal removal rate AMRV Volumetric metal removal rate AV Electrolyte's flow velocity B Electric field potential BANSOL Subroutine that solvessparse,banded and symmetric matrix BWD Subroutine that determines band width of a matrix C Specific heat of electrolyte CDTC Subroutine that evaluatescurrent density, temperature, conductivity etc. DATIN Subroutine that readsand prints out the input data DE Density of electrolyte DELN Length of the normal to a node on the anode surface DELT Area of an element (triangular element in this case) DW Density of the workplace material E Electrochemical equivalent of the work material EPS Precision to which interelectrodegap is computed EV Effective voltage I Connectivity matrix (only in subroutine BWD)

volume 12 number 6 november 1980

Matrix that supplies node numbers for which boundary conditions have been specified INELGP Subroutine that calculatesthe interelectrode gap for a node at any instant of time 't' IT Matrix that assignsa number to a node just above the node on the work surface IW Matrix that transfersallocated numbers to the nodes located on the work surface KSTC Computational cycle time MB Bandwidth of the matrix under consideration METR Subroutine that evaluateslinear, volumetric and specific metal removal rate NE Total number of elements NN Total number of nodes NTST Total specified machiningtime NNUS Number of nodes where 'EV' has been specified NNW Total number of nodes located on the work surface PRNTOT Subroutine that prints out the desired results Q Volumetric flow rate of the electrolyte SINT Numerical value by which interelectrodegap is increased in each iteration when solvingan implicit equation (1 0) SMRR Specific metal removal rate T Nodal temperature VIBC Givesvalue of ith boundary conditions W Width of the workpiece X X-coordinate of a node XTOT Total tool movementat an instant t Y Y-coordinate of a node YE Equilibrium interelectrode gap

313

YSUG Z ZER

If this has positive numerical value then interelectrode gap is computed by equation (12) otherwise by (10) Feed rate Subroutine that makes the coefficients of a matrix as zero

DO 40 LL=L,MB J=J+l A(I,J)=A(I,J )-C(L)*A(K,L L) 40 CONTINUE B(1)=B(1)-C(L)*B(K) 30 CONTINUE

PROGRAM WRITTEN IN FORTRAN IV SUBROUTINE ABC KF=MB-I DO 300 K=I,NNUS N5=IBC(K) DO 400 K=I,KF IR=NS-K1 IF(IR.LT.1) GO TO 450 B(IR)=B(I R)-A(I R,(KI+I ))*VIBC(K) A(IR,(KI+I ))=0.0 450

IR=NS+K1 IF(IR.GT.NN) GO TO 400 B(IR)=B(I R)-A(NS,(KI+I ))*VIBC(K) A(NS,(K1+I))=0.0

10 CONTINUE 100 K=K-1 IF(K.EQ.0) GO TO 200 DO 50 J=2,MB L=K+J-I IF(N.LT.L) GO TO 50 B(K)=B(K)-A(K,J)*B(L) 50 CONTINUE GO TO 100 200 CONTINUE RETURN END SUBROUTINE BWD IX=O DO 100 K=1,NE IYI=IABS(I(K,1 )-I(K,2)) IY2=IABS(I(K,2)-I(K,3)) IY3=IABS(I(K,3)-I (K,1)) I F(IY1 .GE.IY2.AND.IY1 .GE.IY3)IY=IY1 I F(IY2.GE.IY3.AND.IY2.GE.IY1) IY=IY2 I F(IY3.GE.IY1 .AND.IY3.GE.IY2)IY=IY3

400 CONTINUE A(NS,1)=l B(NS)=VIBC(K) 300 CONTINUE RETURN END SUBROUTINE AGCSG DO 101 K=I,NE DE LT=(X(I (K,2))*Y(l (K,3))+X (I(K,I))*Y(I (K,2)) +X(I(K,3))*Y(I(K,I )))-(X(I(K,2))*Y(I(K,I )) +X (I (K,3))*Y(I(K,2))+X(I(K,I))*Y(I(K,3))) BI =(Y(I(K,2))-Y(I(K,3)))/DE LT B2=(Y(I(K,3))-Y(I(K,I )))/DELT B3=(Y(I(K,I ))-Y(I(K,2)))/DELT AI =(X (I(K,3))-X (I (K,2)))/DE LT A2=(X(I(K,I ))-X(I(K,3)))/DE LT A3=(X(I (K,2))-X(I(K,I)))/DELT S(I ,I )=(BI *BI+AI *At )* DE LT/2.0 S(1,2)=(B2" B 1+A2*AI )*DE LT/2.0 S(1,3)=(B3*B I +A3*AI )*DE LT/2.0 S(2,2)=(B2*B2+A2*A2)* DE LT/2.0 S(3,3)=(B3*B3+A3*A3)* DE LT/2.0 S(2,1 )=S(1,2)

s(3,1)=s(1,3)

S(3,2)=S(2,3) DO200 IM=1,3 IR=I(K,IM) DO200 IN=1,3 IC=I(K,IN) ID=IC-(IR-I ) IF (ID.LT.I)GOTO 200 A(I R,I D)=A(I R,ID)+S(I M,IN) 200 CONTINUE 101 CONTINUE RETURN END SUBROUTINE BANSOL N=NN DOIO K=I,NN B(K)=B(K)/A(K,1) IF(K.EQ.N) GO TO 100 DO 20 J=2,MB C(J)=A(K,J) A(K,J )=A(K,J )/A(K,I ) 20 CONTINUE DO 30 L=2,MB I=K+L-1 IF(N.LT.I) GO TO 30 l=O

314

100

IF(IY.GT.IX) IX=IY MB=IX+I RETURN END SUBROUTINE CDTC DO10 II=I,NNW IF(II.EO.NNW) GO TO 12 12=11+1 AC=(Y(IT(12))-Y(IT(I I ))) BC=(X(IT(II ))-X(IT(12))) CC=-X(IT(II ))*AC-Y(IT(I 1))*BC DE LN=ABS((AC*X(IW(I 1))+BC*Y(IW(I I ))+CC) /SORT(AC**2+BC**2)) IF(II.EQ.I) AK(IW(II-I ))=AKI

12 AJ (IW(I1))=((B(IW(I1))-B(IT(I1)))/DE LN) * 7 AK(IW(II-1)) 10 CONTINUE DO15 II=I,NNW V=AV0Wfll)) TAA=(AJ (lW(ll))**2)/(DE*V*C*AK(IW(I I ))) EX I=EX P(TAA*AL F*(X(IW(I I ))-X (lW(ll)))) T(IW(II))=(EXI-I.)/ALF TA(IW(I I ))=T(TI+T(IW(I I )) AK(IW(I 1))=A KI*(1.0+A L F*T(IW(I 1))) 15 CONTINUE RETURN END SUBROUTINE INELGP AKSTC=KSTC DO 10011=1,NNW IND=IW(I1)

Z=AZ(,W(H)) IF(Z),1,55,1 55

BT=E*EV*AK(IW(I1 ))/(F*DW) IF(YSUG.GE.I.0) GOTO 186 Y(IW(I1 ))=SQRT(Y (IW(I 1))*'2+2.* BT*AKSTC) GO TO 8888 NPATH=I IF(YSUG.GE.1.0) GO TO 186 AYT=Y(IW(I1))

computer-aided desigr

BINT=SINT 1=1 YE=EV*E*AK(IW(I1 ))/(F*DW*Z) EY(IW(I1))=YE 2 AA=Y(IW(I1))-AYT+(YE*ALOG((YE-Y(IW(I1 ))) 1 /(YE-AYT)))-Z*AKSTC 3 II:(YE-AYT) 91,91,33 91 BYT=AYT-BINT IF(YE-BYT) 65,163,65 65 CHEN=(YE-Y(IW(I1 )))/(YE-BYT) IF(CHEN) 163, 163, 64 163 BINT=BINT/2. GO TO 91 64 CONTINUE

994 FORMAT (10X, 18HWRONG VALUE OF GYT)

19 CONTINUE IF(AA*BB) 900,900,901 900 CC=ABS(AA)+ABS(BB) RA=ABS(AA)/ABS(CC) Y(IW(II ))=AYT+RA*BINT GO TO 9999 901 IF(AA-BB) 903,903,902 902 Y(IW(I1))=BYT GO TO 9999 903 Y(IW(I1))=AYT GO TO 9999 186 Y(IN D)=Y(IND)+(((E*EV*AK(IN D)) I /(Y(IND)*F*DW)-Z)*AKSTC GO TO 8888

GO TO 56 9999 CONTINUE

33 BYT=AYT+BINT IF(YE-BYT) 66,61,66 66 CHEP=(YE-Y(IW(I1 )))/(YE-BYT) IF(CHEP) 61,61, 62 61 BINT=BINT/2. GO TO 33 62 CONTINUE 56 BB=Y(IW(I1))-BYT+(YE*AIOG((YE-Y(IW(I1))) 1 /(YE-BYT)))-Z*AKSTC 667 J=J+l

IF(J-100) 9,9,12 9 IF(AA*BB) 5,5,4 4 AYT=PYT AA=BB GO TO (3,5), NPATH 5 BINT=BINT/2.0 NPATH=2 IF(BINT-EPS) 11,11,3

8888 AV(IW(I I ))=OI(W*Y(IW(I I ))) 100 CONTINUE

RETURN END SUBROUTINE METR DO4 I1=I,NNW Z=AZ(IW(II )) IF(KIN.GT.1) GOTO 5 AM RR(IW(I I))=Z*AKSTC+Y (IW(I 1))-YII(IW(I I )) GO TO 41 5 AM RR(IW(I 1))=Z*AKSTC+Y(IW(I1 ))-FPY(IW(I 1)) 41 CONTINUE

FPY(IW(I1))=Y(IW(I1 )) CONTINUE DO 6 11=1,NELW AM RV(IW(I 1))=W*(AM RR(IW(I 1))+AMRR(IW(II+I )))* ((X(IW(II+I))-X(IW(I1)))/2.) SMRR(IW(I1 ))=AM RV(IW(I 1))/ ((AJ(IW(II+I))+AJ(IW(I1)))•2.)

11 CONTINUE GOTO 19

RETURN

12 PRINT 994

END

volume 12 number 6 november 1980

315