Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Plane-strain visioplasticity for dynamic and quasi-static deformation processes Dwivedi, Surendra Nath 1978

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1979_A7 D95.pdf [ 6.33MB ]
Metadata
JSON: 831-1.0100192.json
JSON-LD: 831-1.0100192-ld.json
RDF/XML (Pretty): 831-1.0100192-rdf.xml
RDF/JSON: 831-1.0100192-rdf.json
Turtle: 831-1.0100192-turtle.txt
N-Triples: 831-1.0100192-rdf-ntriples.txt
Original Record: 831-1.0100192-source.json
Full Text
831-1.0100192-fulltext.txt
Citation
831-1.0100192.ris

Full Text

PLANE-STRAIN  VISIOPLASTICITY FOR DYNAMIC AND  QUASI-STATIC DEFORMATION PROCESSES by SURENDRA NATH DWIVEDI M.E., U n i v e r s i t y o f Roork.ee;. 1971  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE  in THE FACULTY OF GRADUATE STUDIES (Department  o f Mechanical Engineering)  We accept t h i s t h e s i s as conforming to the r e q u i r e d standard  THE UNIVERSITY OF BRITISH COLUMBIA November, 1978  (c) Surendra Nath Dwivedi, 19 7 8  In p r e s e n t i n g  t h i s t h e s i s i n p a r t i a l f u l f i l l m e n t of  requirements f o r an advanced degree at the U n i v e r s i t y  the  of  B r i t i s h Columbia, I agree that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r reference  and  study.  I f u r t h e r agree  that p e r m i s s i o n f o r e x t e n s i v e copying of t h i s t h e s i s s c h o l a r l y purposes may  be granted by  ment o r by h i s r e p r e s e n t a t i v e s .  the Head of my  Depart-  I t i s understood t h a t  i n g or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain not be allowed without my  for  written  shall  permission.  (Surendra Nath Dwivedi)  Department of Mechanical E n g i n e e r i n g The U n i v e r s i t y of B r i t i s h Columbia 20 75 Wesbrook P l a c e Vancouver, B.C., Canada V6T 1W5  D  a  t  e  •• 3 o [ f 11 7  g  copy-  ABSTRACT  The v i s i o p l a s t i c i t y approach complete  i s developed t o enable the  s t r e s s h i s t o r y of any steady or non-steady,  quasi-  s t a t i c o r impact, plane s t r a i n p l a s t i c deformation process t o be determined from a r e c o r d o f the deformation p a t t e r n .  The  v e l o c i t y f i e l d i s determined e x p e r i m e n t a l l y and f o r dynamic c o n d i t i o n s h i g h speed photographs  are taken of a g r i d p a t t e r n  marked on the end s u r f a c e of the specimen.  D i g i t i z a t i o n of  the instantaneous g r i d node p o s i t i o n s allows the v e l o c i t y f i e l d s to be o b t a i n e d a t predetermined time i n t e r v a l s out the t r a n s i e n t deformation p e r i o d .  through-  Hence, the s t r a i n - r a t e ,  e q u i v a l e n t s t r a i n r a t e , e q u i v a l e n t s t r a i n and f i n a l l y  stress  f i e l d s can a l l be o b t a i n e d . A three dimensional s u r f a c e f i t t i n g procedure, u s i n g f o u r t h o r d e r p o l y n o m i a l s , i s used t o smooth the s c a l a r component of the e x p e r i m e n t a l l y determined v e l o c i t y f i e l d . c o n d i t i o n of c o n t i n u i t y  The  ( -^=- ^ ^ • f o r plane strain), i s e  e  r  imposed on both s u r f a c e s thereby r e d u c i n g the number o f i n dependent parameters  from 30 to 10.  Besides smoothing  experimental p o i n t s t h i s procedure has the d i s t i n c t  the  advantage  t h a t the polynomials can be r e a d i l y d i f f e r e n t i a t e d f o r det e r m i n i n g s t r a i n - r a t e s and t h a t deformation can be  referred  to a master r e f e r e n c e g r i d t h a t i s f i x e d w i t h r e s p e c t to time. P l a n e - s t r a i n u p s e t t i n g t e s t s , conducted a t a speed of 0*02 f t / m i n g i v e r e s u l t s t h a t agree c l o s e l y w i t h the w e l l  docu-  merited ' f r i c t i o n quasi-static  hill  1  type of normal s t r e s s  r a t e s o f s t r a i n . However.,with  at a speed o f 15;7 f t / s e c the normal s t r e s s radically different exhibiting  distribution for  the specimen deformed distribution i s  a saddle type d i s t r i b u t i o n .  The e f f e c t o f s t r a i n r a t e on the i n t e r f a c e  and body  stresses  w i l l have s i g n i f i c a n t b e a r i n g on a number o f metal forming operations.  iv'  ACKNOWLEDGEMENTS  The G.W.  author s i n c e r e l y thanks the t h e s i s s u p e r v i s o r , P r o f e s s o r  V i c k e r s f o r h i s unending words o f encouragement and whose  amiable a t t i t u d e made t h i s work a p l e a s a n t  experience.  S i n c e r e thanks must a l s o go t o the other members o f the t h e s i s committee i . e . , P r o f e s s o r H. Ramsey, I.S. Gartshore, N.G. E l y , R.E. McKechnie and Norman Franz f o r t h e i r time t o time f r u i t f u l d i s c u s s i o n s and expert The  opinions.  author i s most g r a t e f u l t o a l l the other members o f  the department f o r t h e i r i n s p i r i n g i n t e r e s t and a s s i s t a n c e . Thanks are a l s o due t o my wife Shashi help and f o r her moral Above a l l ,  support.  I thank God f o r g i v i n g me the o p p o r t u n i t y to  b e t t e r myself academically the  research.  f o r some t y p i n g  and s p i r t u a l l y d u r i n g the course o f  V  TABLE OF CONTENTS  Page ABSTRACT. . . .  i i  ACKNOWLEDGEMENTS  i v  TABLE OF CONTENTS  V  LIST OF FIGURES  V,ijLi  NOMENCLATURE  xi'.  CHAPTER ONE - INTRODUCTION AND LITERATURE SURVEY 1.1  1.2  1  PLASTICITY IN METAL WORKING  3  1.1.1  S l a b Method  3  1.1.2  Uniform Deformation Energy Method  4  1.1.3  S l i p Line Method  4  1.1.4  Limit Analysis  5  1.1.5  Finite  Element Method  6  1.1.6  Finite  D i f f e r e n c e Method  . ..  VISIOPLASTICITY  8  CHAPTER TWO - NON-STEADY, PLANE-STRAIN, DYNAMIC AND QUASI-STATIC VISIOPLASTICITY 2.1  EQUATIONS FOR QUASI-STATIC VISIOPLASTICITY 2.1.1  2.1.2  2.1.3  2.2  9 10  Equations f o r Three Dimensional Nonsteady S t a t e , Q u a s i - s t a t i c , v i s i o plas t i c i t y  10  Equations f o r P l a n e - S t r a i n , Non-steady State, Quasi-static V i s i o p l a s t i c i t y  13  Determination o f the S t r e s s F i e l d the S t r a i n F i e l d  14  from  EQUATIONS FOR DYNAMIC VISIOPLASTICITY 2.2.1  7  Equations f o r Three Dimensional, Nonsteady S t a t e , Dynamic, V i s i o p l a s t i c i t y . . .  16  16  vi  Page 2.2.2  2.2.3 2.3 2.4  2.5  Equations f o r Plane S t r a i n Non-Steady S t a t e Dynamic V i s i o p l a s t i c i t y  18  Determination o f the S t r e s s F i e l d from the S t r a i n F i e l d  18  GENERAL PROCEDURE FOR THE STUDY OF DEFORMATION USING VISIOPLASTICITY  19  GENERALIZED EQUATIONS FOR PLANE-STRAIN VISIOPLASTICITY FOR COMPUTATIONAL PURPOSES  23  THE COMPUTER PROGRAM  29  CHAPTER THREE - VARIABLE SPEED, CONTROLLED VELOCITY PROFILE, SINGLE CYCLE IMPACTING PRESS.. 3.1  3.2  DESCRIPTION OF EQUIPMENT  37  3.1.1  Background Information  37  3.1.2  D e s c r i p t i o n o f the Press  39  3.1.3  Operation o f the Press  39  KINEMATIC ANALYSIS OF THE MECHANISM 3.2.1 D e r i v a t i o n o f E x p r e s s i o n s f o r V e l o c i t y and A c c e l e r a t i o n  40  3.2.2  3.3  36  Maximum V e l o c i t y and A c c e l e r a t i o n of the Mechanism  CONTROL  40  43 47  3.3.1  C o n t r o l f o r S i n g l e Cycle Operation  3.3.2  C o n t r o l f o r S y n c h r o n i z i n g the High Speed Camera  CHAPTER FOUR - EXPERIMENTAL PROCEDURE AND DISCUSSION.. . .  47  50 51  4.1  EXPERIMENTAL PROCEDURE  51  4.2 4.3  DISCUSSION SOURCES OF ERRORS  53 84  vii  Page  CONCLUSIONS  86  SUGGESTIONS FOR FURTHER, WORK  8  7  REFFERENCES-,  88  APPENDIX  95  viii  LIST OF FIGURES  Page  1.  E l e m e n t a l cube f o r d e r i v a t i o n of the equation of motion  2'.  17  P o s i t i o n of p o i n t P a t " i n s t a n t t ' and't +i n  3 .-- G r i d - l i n e s -on the specimen. .4.  Flow c h a r t of computer program  5.  V a r i a b l e speed, c o n t r o l l e d v e l o c i t y  n  . ...  : .V  2 0  31 32  profile,  s i n g l e c y c l e impacting p r e s s  38  6.  D r i v e mechanism of impacting press  41  7.  Motion of the c y c l o i d a l cam  44  8. 9.  Details of e l e c t r i c a l c i r c u i t for control Plow diagram f o r r e l a y sequence f o r s i n g l e control  10. 11.  12.  13.  14.  15.  16.  46 cycle 48  C i r c u i t o p e r a t i o n b l o c k diagram f o r high speed camera s y n c h r o n i z a t i o n  49  D i s t o r t i o n of g r i d l i n e s d u r i n g deformation f o r 0.02 f t . / m i n . deformation speed (a) F i r s t s t e p ; and (b) second step  55  Distortion of g r i d l i n e s d u r i n g deformation f o r 0.02 f t . / m i n . deformation speed (c) T h i r d step; and (d) f o u r t h s t e p . . . .  56  D i s t o r t i o n of g r i d l i n e s d u r i n g deformation f o r 15.7 f t . / s e c . deformation speed. ( F i r s t step)  57  D i s t o r t i o n of g r i d l i n e s d u r i n g deformation f o r 15.7 f t . / s e c . deformation speed. (second s t e p ) . . . .  58  D i s t o r t i o n of g r i d l i n e s d u r i n g deformation f o r 15.7 f t . / s e c . deformation speed. ( t h i r d step)  59  D i s t o r t i o n o f g r i d l i n e s d u r i n g deformation f o r 15.7 f t . / s e c . deformation speed. (fourth step)....  60  i x Page 17.  Distortion  o f g r i d l i n e s d u r i n g deformation f o r  1 5 . 7 f t . / s e c . deformation speed. 18. 19. 20. 21. 22.  23.  24.  25.  26.  27.  28.  29.  ( f i f t h step)....  61  D i s t o r t i o n o f g r i d l i n e s d u r i n g deformation f o r 1 5 . 7 f t . / s e c . deformation speed. (sixth step)....  62  D i s t o r t i o n o f g r i d l i n e s d u r i n g deformation f o r 1 5 . 7 f t . / s e c . deformation speed. (seventh s t e p ) . .  63  G r i d node p o i n t movements d u r i n g deformation f o r 0.02 f t . / m i n . deformation speed  64  G r i d node p o i n t movements d u r i n g deformation f o r 1 5 . 7 f t . / s e c . deformation speed.  65  Three dimensional p l o t o f h o r i z o n t a l v e l o c i t y (u) as a f u n c t i o n o f x and y f o r 0.02 f t . / m i n . deformation speed  66  Three dimensional p l o t o f v e r t i c a l v e l o c i t y as a f u n c t i o n o f x and y f o r 0 . 0 2 f t . / m i n . deformation speed  67  ( v)  Three dimensional p l o t o f h o r i z o n t a l v e l o c i t y (u) as a f u n c t i o n o f x and y f o r 1 5 . 7 f t . / s e c . deformation speed Three dimensional p l o t o f v e r t i c a l v e l o c i t y as a f u n c t i o n o f x and y f o r 1 5 . 7 f t . / s e c . deformation speed  68  (y0 69  Three dimensional p l o t o f e f f e c t i v e s t r a i n - r a t e as a f u n c t i o n o f x and y f o r 0.02 f t . / m i n deformation speed  70  Three dimensional p l o t o f e f f e c t i v e s t r a i n - r a t e as a f u n c t i o n o f x and y f o r 1 5 . 7 f t . / s e c deformation speed  71  Three dimensional p l o t o f t o t a l e f f e c t i v e s t r a i n as a f u n c t i o n o f x and y f o r 0.02 f t . / m i n deformation speed  72  Three dimensional p l o t o f t o t a l e f f e c t i v e as a f u n c t i o n o f x and y f o r 1 5 . 7 f t . / s e c  73  strain  X  Page  30.  31.  32.  33.  34.  35.  36. 37.  38.  39.  a  Three dimensional p l o t of normal s t r e s s ( ) as a f u n c t i o n of x and y f o r 0.02 f t . / m i n . deformation speed v  7 4  a  Three dimensional p l o t of normal s t r e s s ( ) as a f u n c t i o n of x and y f o r 0.02 f t . / m i n . deformation speed  75  Three dimensional p l o t o f shear s t r e s s ( T y ) as a f u n c t i o n of x and y f o r 0.02 f t . / m i n . deformation speed  76  Three dimensional p l o t of normal s t r e s s ( Oy) as a f u n c t i o n o f x and y f o r 15.7 f t . / s e c . deformation speed  77  x  X  a  Three dimensional p l o t of normal s t r e s s ( ) as a f u n c t i o n of x and y f o r 15.7 f t . / s e c . deformation speed. x  Three dimensional p l o t o f shear s t r e s s ( x ) as a f u n c t i o n of x and y for, 15.7 f t . / s e c . deformation speed x  7 8  v  Two dimensional p l o t of o. ' as a f u n c t i o n of x f o r 0.02 f t . / m i n . deformation speed...  7  ^  1  y  8 0  Two dimensional p l o t of ' T y ' as a f u n c t i o n of x f o r 0.02 f t . / m i n . deformation speed X  81  a  Two dimensional p l o t of ' - y' as a f u n c t i o n of x f o r 15.7 f t . / m i n . deformation speed Two dimensional p l o t of x y as a f u n c t i o n of x f o r 15.7 f t . / m i n . deformation speed  8  ^  X  83  NOMENCLATURE yt  t h e component o f v e l o c i t y i n x - d i r e c t i o n  v  t h e component o f v e l o c i t y i n y - d i r e c t i o n  W  t h e component o f v e l o c i t y i n z - d i r e c t i o n  £  x  s t r a i n rate alone x - d i r e c t i o n  *  £y ;  V ixy  . s t r a i n rate along y - d i r e c t i o n shear s t r a i n - r a t e  £  effective strain-rate  £  total effective strain  0  effective stress  0"  normal s t r e s s a l o n g x - d i r e c t i o n  a  normal s t r e s s along y - d i r e c t i o n y  0"  normal s t r e s s a l o n g  T  xy  z-direction  shear s t r e s s  p  the density  V  t h e volume o f t h e p l a s t i c a l l y d e f o r m i n g body  £  volumetric strain-rate  v  F  o f the m a t e r i a l  t h e t r a c t i o n on t h e p a r t o f s u r f a c e S r  0 A  t h e v e l o c i t y p r e s e n t e d on t h e remainder s u r f a c e S. Lagrange m u l t i p l i e r  a  the penalty  function  0  t h e a n g l e made fcyrthe d r i v i n g arm w i t h  vertical  C H A P T E R  ONE  INTRODUCTION & LITERATURE SURVEY  1  2 CHAPTER ONE  INTRODUCTION  The mechanism o f p l a s t i c deformation p l a y s a v i t a l r o l e i n many i n d u s t r i a l metal-working  processes.  However, i t has not  proved p o s s i b l e to analyse completely many o f these processes u s i n g the g e n e r a l b a s i c equations "derived from the - theory - of plasticity.  T h i s i s p r i m a r i l y due to u n c l e a r l y d e f i n e d boun-  dary c o n d i t i o n s ; f o r example, t h e a c t u a l f r i c t i o n a l c o n d i t i o n s p r e s e n t a t the metal-die i n t e r f a c e are f r e q u e n t l y unknown. Many s i m p l i f i e d a l t e r n a t i v e methods have been developed and used to study c e r t a i n o f the metal forming p r o c e s s . these analyses c e r t a i n assumptions  In  and s i m p l i f i c a t i o n s a r e  made r e g a r d i n g the processes and the behaviour o f the m a t e r i a l s d u r i n g deformation.  However, i n s p i t e o f these  the s o l u t i o n s o f t e n l a c k uniqueness One  approach  idealizations  and completeness.  c a l l e d v i s i o p l a s t i c i t y has been used with  some success t o determine  the complete  stress picture  through-  out the deformation zone i n c e r t a i n s t e a d y - s t a t e e x t r u s i o n and f o r g i n g o p e r a t i o n s .  I t r e q u i r e s t h a t the v e l o c i t y  field  be determined e x p e r i m e n t a l l y and hence the s t r a i n - r a t e s ] .;. A .v and f i n a l l y s t r e s s f i e l d s can a l l be obtained.  This  method has been shown t o give r e a l i s t i c s o l u t i o n s and i t s a p p l i c a t i o n has been extended  d u r i n g the l a s t decade.  In t h i s work the v i s i o p l a s t i c i t y approach  has been used  to study m a t e r i a l deformation i n dynamic and non-steady  condi-  3 tions.  The r e l e v a n t equations and procedure have been embodied  i n t o a s p e c i a l l y developed computer program,. so t h a t the comp l e t e s t r e s s h i s t o r y o f any steady o r non-steady,  quasi-static  or impact plane s t r a i n , deformation can be determined r e c o r d o f the deformation p a t t e r n .  from a  S p e c i a l a t t e n t i o n has been  given i n t h i s work to smoothing the e x p e r i m e n t a l l y determined v e l o c i t y f i e l d s , a p o i n t which has caused some d i f f i c u l t y i n the p a s t .  R e s u l t s from t h i s work have been s u i t a b l y compared  with p r e v i o u s s t e a d y - s t a t e r e s u l t s f o r v e r i f i c a t i o n 1.1.  purposes.  PLASTICITY IN METAL-WORKING While the analyses of metal-working  processes has been r e -  s t r i c t e d by the c o m p l e x i t i e s i n v o l v e d , some approaches been made.  have  A number o f these i n common use are the s l a b (or  e q u i l i b r i u m ) method, uniform deformation energy method, s l i p l i n e s o l u t i o n s , upper and lower bound s o l u t i o n s , d i f f e r e n c e and f i n i t e element methods.  F o r completeness  b r i e f d e s c r i p t i o n o f these common approaches 1.1.1.  finite a  i s given below.  Slab o r E q u i l i b r i u m Method.  The method i n t r o d u c e d by Sachs (1) i n 19 31, c o n s i s t s o f i s o l a t i n g a s m a l l elemental volume o f the m a t e r i a l under goihrg deformation and o b s e r v i n g the behaviour o f t h i s element as i t moves through the working  zone.  Since t h i s element i s an i n -  t e g r a l p a r t o f the m a t e r i a l , i t should always be i n a s t a t e of e q u i l i b r i u m .  The assumption  i s made t h a t s t r e s s e s on a  plane s u r f a c e p e r p e n d i c u l a r to the d i r e c t i o n o f the flow are p r i n c i p a l s t r e s s e s and t h a t these do n o t vary on t h i s plane.  4 A n a l y s i s of the e q u i l i b r i u m c o n d i t i o n r e s u l t s i n one  or  more d i f f e r e n t i a l equations which together with the necessary boundary c o n d i t i o n s , give the deformation  stresses.  Since the e f f e c t o f redundancy, f r i c t i o n and p a t t e r n of flow are not c o n s i d e r e d , t h i s method g i v e s an of the deformation s t r e s s e s .  underestimate  However, the a n a l y s i s i s s t r a i g h t -  forward and i t has been w i d e l y used i n wire and tube  drawing  problems as w e l l as hot and c o l d r o l l i n g of s t r i p and sheet (1). 1.1.2.  Uniform Deformation  Siebel  Energy  Method:  (2) proposed t h i s approach  i n 1932  i n which the  amount o f deformation i s determined by c o n s i d e r i n g the shape of an element o f m a t e r i a l b e f o r e and a f t e r deformation. hence g i v e s o n l y the average s p e c i f i c i n t e r n a l energy s t a t e metal-working 1.1.3.  It  forming p r e s s u r e as a f u n c t i o n of  and i s g e n e r a l l y used f o r steady-  processes.  S l i p L i n e Method:  Hencky(3) i n t r o d u c e d the s l i p l i n e theory i n 1923. can be used f o r determining the l o c a l s t r e s s and  It  velocity  d i s t r i b u t i o n d u r i n g deformation, although i t i s r e s t r i c t e d t o plane s t r a i n c o n d i t i o n s and r e q u i r e s a predetermined p a t t e r n of flow. The  s l i p l i n e s o l u t i o n c o n s i s t s of f a m i l i e s o f c u r v i l i n e a r  o r s t r a i g h t l i n e s , which are p e r p e n d i c u l a r t o each o t h e r and correspond to the d i r e c t i o n s of maximum and minimum constant shear s t r e s s .  These l i n e s s a t i s f y the s t a t i c e q u i l i b r i u m  con-  d i t i o n , y i e l d c o n d i t i o n and the p a t t e r n of flow everywhere i n the p l a s t i c zone of the m a t e r i a l .  These shear or s l i p  lines  5  are known as c h a r a c t e r i s t i c s o f the d i f f e r e n t i a l equations of equilibrium.  In the s l i p l i n e method the forming t o o l and the  m a t e r i a l o u t s i d e the s l i p l i n e are c o n s i d e r e d as r i g i d ( i . e . the metal ahead and behind the p l a s t i c zones and the t o o l m a t e r i a l have an i n f i n i t e modulus of e l a s t i c i t y ) .  The  slip  l i n e s o l u t i o n i s not o p t i m a l or unique and a l s o g i v e s values higher than the t r u e s o l u t i o n . T h i s method has been w i d e l y used f o r the study o f many metal deformation processes (4-14), some of the l a t e s t work has i n v o l v e d s l i p l i n e s o l u t i o n s f o r a n i s o t r o p i c m a t e r i a l s (15, 16) and has taken account o f f r i c t i o n on the die-workpiece interface  (17, 18).  have produced  A l s o Ewing  (19) and l a t e r C o l l i n s  s l i p l i n e s o l u t i o n using numerical  (20)  computation  by power s e r i e s and by matrix o p e r a t i o n a l methods. 1.1.4.  Limit Analysis:  The mathematical  model o f l i m i t a n a l y s i s p l a c e s upper  lower estimates on the l o a d r e q u i r e d f o r deformation.  and  This  l i m i t a n a l y s i s i s based on two extremum theorems put forward by Prager and Hodge (8), and Drucher and Prager gave the mathematical on the assumption plastic. a)  (21).  Hill  (7)  p r o o f o f these theorems, which are based  t h a t the m a t e r i a l i s r i g i d and  perfectly  They can be s t a t e d as:  Upper Bound Theorem.  I f a kinematically admissible v e l o c i t y  f i e l d e x i s t s , the loads r e q u i r e d to be a p p l i e d to cause the v e l o c i t y f i e l d to operate c o n s t i t u t e an upper bound s o l u t i o n . b)  Lower Bound Theorem.  I f a s t a t i c a l l y admissible stress  f i e l d e x i s t s such t h a t the s t r e s s e s are everywhere j u s t below  6 those necessary to cause y i e l d i n g , then the loads a s s o c i a t e d w i t h t h a t f i e l d c o n s t i t u t e a lower bound s o l u t i o n . These techniques have been used e x t e n s i v e l y study metal-working p r o c e s s e s , such as f o r g i n g ,  (22-39) to extrusion,  wire drawing and tube drawing. 1.1.5.  F i n i t e Element- Method:  The f i n i t e element method i s one o f the most powerful t e c h niques f o r s o l v i n g two dimensional problems i n metal-working but  at p r e s e n t has a l i m i t e d p o t e n t i a l f o r complex  to economic  constraints.  (40) i n 19 54.  This method was  problems  due  i n t r o d u c e d by A r g y r i s  In t h i s approach, the deforming area o r continuum  i s s u b d i v i d e d i n t o an e q u i v a l e n t system o f elements, known as f i n i t e elements. of  The f i n i t e elements may  be t r i a n g l e s ,  group  t r i a n g l e s , q u a d r i l a t e r a l e t c . f o r two dimensional s t u d i e s  and t e t r a h e d r a , r e c t a n g u l a r prisms o r hexahedra e t c . f o r t h r e e dimensional s t u d i e s .  Once a displacement model i s s e l e c t e d ,  an element s t i f f n e s s m a t r i x i s d e r i v e d u s i n g v a r i a t i o n a l ciples.  prin-  The a l g e b r a i c equations f o r the whole continuum are  then assembled and s o l u t i o n s f o r unknown displacements a t the nodal p o i n t s can be o b t a i n e d .  By use o f the computed d i s p l a c e -  ments and the s t r e s s and s t r a i n r e l a t i o n s , the s t r e s s e s a t the nodal p o i n t s may  be determined.  The s o l u t i o n i s based on ex-  tremum p r i n c i p l e a c c o r d i n g to which the a c t u a l s o l u t i o n m i n i mizes the f u n c t i o n a l <f>= /  a  $  e d V r  with the c o n s t r a i n t t h a t & = 0 v  , where /  I  .. y dS  7  where  v  = volume o f the p l a s t i c a l l y deforming body  cP = e f f e c t i v e s t r e s s e  = effective strain  i  = volumetric  v  F and  rate  strain  rate  = the t r a c t i o n o f the p a r t o f s u r f a c e  ii = the v e l o c i t y p r e s c r i b e d on A modified  (41)  f u n c t i o n a l has  1  where wicz  a penalty  f u n c t i o n , a , as  cb = / a e d V + /  ^- (e )  By  = very  that  While Godbole and  Zienkie-  f u n c t i o n a l be m o d i f i e d  using  follows:  dv - / F • u dS S  where a  & Kobayashi  F • udS  have suggested t h a t the  2  surface  Op  A = the Lagrange m u l t i p l i e r . (42)  p  the remainder  been given by Lee  using the Lagrange m u l t i p l i e r so  cb = / CT e dV + 7 A e dV v V  S  F  l a r g e number.  use of above f u n c t i o n a l s , many metal-working processes such  as u p s e t t i n g drawing, p i e r c i n g e t c . 1.1.6.  (40-57) have been  studied.  F i n i t e D i f f e r e n c e Method:  This method i s one  of the most r e c e n t techniques to  used f o r the study of metal-working processes.  be  I t requires  t h a t the continuum be d i v i d e d i n t o a number o f g r i d s and 'difference  1  ( i . e . f i n i t e ) q u a n t i t i e s are s u b s t i t u t e d f o r  d i f f e r e n t i a l q u a n t i t i e s across d i f f e r e n t i a l equation with taneous equations can be numerically  that  the g r i d s .  Thus f o r a  given  boundary c o n d i t i o n s a s e t o f  s u b s t i t u t e d , which can be  u s i n g a computer.  The  s i z e o f the g r i d  determines the accuracy of the s o l u t i o n .  The  simul-  solved spacing  f i n e r the  grid,  8  the b e t t e r i s the accuracy o b t a i n e d , but t h i s a t the expense of computer c o s t . S t u d i e s which d e s c r i b e the a p p l i c a t i o n o f t h i s to  technique  f o r g i n g , e x t r u s i o n and sheet-metal processes are given i n  references 1.2.  (58-64).  VISIOPLASTICITY The v i s i o p l a s t i c i t y was  i n t r o d u c e d by Thomsen (56, 66, 6 7)  and l a t e r developed and extended by Shabaik, Kobayashi e t a l . (68, 69, 70, 71, 72). are  photographed  In t h i s method, the g r i d l i n e p a t t e r n s  f o r each i n c r e m e n t a l step of deformation and  thus the movement o f g r i d p o i n t s can be determined. l a r g e d photographs  From en-  o f c o n s e c u t i v e g r i d p a t t e r n s the i n s t a n t a n -  eous v e l o c i t i e s o f a l l g r i d node a c r o s s the s u r f a c e s can be found.  The s t r a i n s ,  thus be determined  strain rates, t o t a l effective  strain  can  f o r a l l p o i n t s i n the deformation r e g i o n  and f i n a l l y the s t r e s s f i e l d and forming loads may  be  found.  In t h i s method the instantaneous flow f i e l d i s an a c t u a l one and g i v e s i n f o r m a t i o n o f a l l s t r a i n s e n t i r e deformation r e g i o n .  I t may  and s t r e s s e s over the  be used f o r both work-  hardening and non-workhardening m a t e r i a l s . D e t a i l s o f the b a s i c equations used i n v i s i o p l a s t i c i t y are given i n the next chapter. The v i s i o p l a s t i c i t y method has been a p p l i e d to f o r g i n g and axisymmetric cesses.  and plane s t r a i n e x t r u s i o n and r o l l i n g pro-  (Reference 68 to 76).  Recently i t has been used f o r  i n v e s t i g a t i n g the r e l a t i o n s h i p between s t r a i n and (77),  microhardness  crack propagation and f o r the d e r i v a t i o n of c r i t e r i a f o r  ductile  r u p t u r e o f f u l l y p l a s t i c notched bars (78).  9  C H A P T E R  TWO  NON-STEADY PLANE-STRAIN DYNAMIC AND QUASI-STATIC VISIOPLASTICITY  10  CHAPTER  TWO  NON-STEADY PLANE-STRAIN DYNAMIC AND QUASI-STATIC VISIOPLASTICITY  2.1.  EQUATIONS FOR  2.1.1  QUASI-STATIC VISIOPLASTICITY  Equations f o r Three Dimensional Non static Visioplasticity  Steady S t a t e Quasi-  The f o l l o w i n g equations i n three dimensions  are used t o  d e s c r i b e the mechanism o f p l a s t i c deformation of an  isotropic  solid. •  •  •  The s t r a i n r a t e £ , e . x  ,y  can be given terms of  y  xy  v e l o c i t y components as f o l l o w s :  •  _  x  | 2 y  3u 3  x,  Y '/ _ 3 u  3v  +  xy ~ 3 y  7Y.  +  E  3 V  =  z  3y,  Y  ax',  _3_u 3z  =  zx  ;  _  3w  yz  ^y  +  ,  9 W  3z  3v 3z  ; and  3w 3x  (2.1)  where u, v, w are the components of v e l o c i t y i n the x, y and z directions respectively.  The equations of s t a t i c  n e g l e c t i n g a l l body f o r c e s are: 3a 3x  8 T  xy 3x  3T  3i: 3y  +  i i x gy  9z  +  5 T  y  z  = o  3z  equilibrium,  11 da  3T 3  T  x  +  z  IJl  3 x  rule)  a  stress  and s t r a i n  •  x  o  9 z  rate  relationships (or  i s g i v e n by  • E  =  2  3y  The L e v y - V o n M i s e s flow  +  x +p  £  _ ""'  ^  a  y  • £  y _ +p  a  •  z z  Y  _ +P  ~  *  xy _ 2  ^  xy  •  ^yz _ ~  2T  yz  ^zx _ \ 21 zx  ~  ! where p = -  (2.3)  j (  CT v  + x  CT  „ +  ° ) / rr  2 <-a e=  effective strain  rate,  and  o = effective stress. The Von M i s e s y i e l d (  a x  a ) y  2z  condition i s  2  + ( a - C T ) ^ + ( Y z  a z  2  2  2  2  a ) +6 ( r +T +T ) = x xy zy zx  -2 20 —  (2.4) The  above e q u a t i o n  stresses  ( a  where  ^,  a  c  a )  V  and  2  2  2  +  (  a  2  =  a  a  - a ) 3  2 + (  a^  CT  a  =  £  ()  ( e  strain  From e q u a t i o n  f  f  and  2  o ^ ) = 2a  2  i s constant f o r perfectly  materials  crCe) ^ f  =  of principal  as  = e f f e c t i v e s t r e s s , which plastic  rate,  ( 2 . 4 ) can be e x p r e s s e d i n terms  o  r  non-workhardening  plastic  f o r workhardening  plastic  e , T ) , f o r material  which  temperature.  (2.3),  we  may  write  material  material  i s a f f e c t e d by  strain  12  •  3— e X  x  25  •  ^  =  e  +p)  X(a  3  G  Y  >a  *  =  +p)  y  = .  A (a  A(a e  +p)  y  +p)  z  2a  z r  Y  3~ t  _  xy  T  2*T  ^  xy  2^  , xy  • •  Y  si  =  2Ax  y z  1 A T  zx  a  . zx  Now  =  j •=a  z"  zx  (2.5)  s u b t r a c t i n g the second equation  (2.5)  (2.5) from f i r s t  equation  gives  _  e  e  and s i m i l a r l y  e  =Ar(a  Y  x  -  7  f o r the o t h e r equations  _  = 3~  Y  a )  x  2a.  ( a  —  z?  Y  (2.5)  ~ a .) z  2a  = Z  3  e  ( a  X  ~ a Z  ) X  2a  3Y T  Y  2a  3y  3  Y  T  x ^  ^ X  x  2f  Z =  A zx T  y  r z 2a  = V-  x  2  t z  x  z  ^  7  (2-6)  2a  Squaring the l e f t hand s i d e and equations (e  "  x  (2.6) and adding, we 2  e ) +(  y  r i g h t hand s i d e o f a l l the  e -  y  2  e ) +(  z  e -  z  get ) +  e  2  x  \(  y  2  xy  +  Y  7  2 Z  +  Y  ,  2 z x  7  13  /3 e'  I(a  L  2  2  - a ) + ( a  x  y .  - a ) + ( a  y  z  ) +6T 2  -a  z  X  2  x  +6T  *  7  2  +6T  ]  2  *J  z  z  (2.7) Combining equation .  ( £  x  -  •  , 9 £  y  )  2  (2.4) and (2.7) gives  • v 2  3  , • • >2 , '2 (e - e ) ( Y + xy  2  (e - e ) . v  +  z  J  +  z  y  f -2 |x 2a — 0  9 3  2  2  x  T  e  ! '< (x e —e e) + ( IF 1 y  ;  +  s t r a i n r a t e may be given as e  "  2  e  y.  )+(  e  z  2  e ) ^  -  2  2  "z\ y  x  e z  e ) x  e) +( I  2  e) +( e y -  2  -  o  'e)  z  y  ^  y  +  z  Y  .  l\[  zxif  2  +  e -  +  x  z  e ) +( e - e ) +( y y z  ( y  3  I  (  -2, Y ) zx  -2  So t h a t the e q u i v a l e n t  >  x  '2 Y + yz  TJ-  .  2 - 2  •  ( (yY 'xy  .2  + Yy z  +  Y z x ),  1 2  2  3(Y '4* Y X  2 +  ' + Y Y  2 +  + z  * Y z  2  :  x  (2.8) 2.1.2.  Equations f o r Plane S t r a i n Non-Steady State Visioplasticity  For the plane s t r a i n c o n d i t i o n ,  e  , Z  Y  ,  ' yz  y  ,  zx  x  yz  p l a s t i c material). given  and  e  , 2  Quasi-Static  Y„„# JL  y  „„r ZX  ~ ', a l s o - C L = i ( f o r a r i g i d p e r f e c t l y zx  Hence the b a s i c equations f o r plane s t r a i n are  as f o l l o w s : The equation f o r s t a t i c e q u i l i b r i u m from equation da  • x  9T  . .xy = 0  +  ' 3x  ",9y  9T  3a +  3x  (2.2) i s  Y _  •  9y  (2.10)  14  The  flow r u l e from equation (2.3) now becomes  8  a  where  x  x  E  _  ,y  +p  A  Y  _  a +p c  y  3  =  xy  _  •  (2.11)  2T  xy  | e a  Von Mises y i e l d c r i t e r i a  2  x -:,ay ) +  ( a  +6(  2  x  or a'-=  xy  +  1  from e q u a t i o n (2.4) becomes  a ] y -2i ( ax + ay )l + [2i ( xa +y a )~ x 2  \a  J  2  T  4  u  u  2  J  2  0 + 0) = 2 a  [ -nfa - a ) + 3 * x y  L  2  xy  ] '  5  J  (2.12)  The e f f e c t i v e s t r a i n r a t e can be w r i t t e n from e q u a t i o n (2.8) as 2 -  L E  T  =  2.1.3.  [  3  2  " 2 3 * e  + T  x  Y  ^ y  J ;  ,^-  :  (  Determination o f the s t r e s s f i e l d  2.i3)  from the s t r a i n  field  Once a l l the normal and shear s t r a i n s are known throughout the  deforming r e g i o n from e q u a t i o n (2.1) the s t r e s s e s may be  calculated. the  To determine the s t r e s s f i e l d  from the s t r a i n  f o l l o w i n g steps are r e q u i r e d . From equation (2.3)  e  x  e  y  _  3  P  - -52.a  ( a +p^ a -p)  x ^  y  r  3  e  = -2a  ( a ~ a ) u  x  °y  field,  15 2.0a  Differentiate •3a,  (e  x~  3a  x 3x  3  From the e q u i l i b r i u m  w i t h r e s p e c t to x g i v e s  e —e x> y L A _J  (2.15)  (2.10)  equation  i n equation 3T  1  (2.15)  3  xy 3y  3x  E  Tx  we g e t  -e x .y  (2.16)  (2.10)  From equation 3a,  9T.  (2.17)  xy 3x  3y  (2.16)  Using equation  and ( 2 . 1 7 )  a t x = 0 and y = bt i . e . 3  X-  T  (0,a) - • •*( - xy dy - f  a Y  y  w i t h the known value o f  0"(x,y)  a(o,a) we g e t ,  V  (x,y) =  a  (2.14)  £x-£y x  xy 3y  Substituting  Y  X  a  3T  x  3X  3a  - e ) = y  (2.14)  equation  . 3x  3a  x  a'  av  J  '3T  ,£  -  £  xy+ 3 \ x y \ 3y 3x\ ^ /—1 y=a  dx  (2.18) (2.14)  From equation  (2.19)  x  f  y+  and from equation  xy 2x  xy  ^  F  ^  (2.11)  (2.20)  where A- =  T  c  Z  (2.21)  2.2  EQUATIONS FOR  2.2.1  DYNAMIC VISIOPLASTICITY  E q u a t i o n s f o r Three Dimensional Dynamic V i s i o p l a s t i c i t y  Non-Steady S t a t e  For the case of dynamic d e f o r m a t i o n s t r e s s and  s t r a i n r a t e r e l a t i o n s h i p and Von M i s e s y i e l d  are s i m i l a r t o the t h r e e d i m e n s i o n a l case. by  the Levy-Von M i s e s  non-steady  criteria  quasi-static  However, the s t a t i c e q u i l i b r i u m e q u a t i o n i s r e p l a c e d  the e q u a t i o n o f motion. The  e q u a t i o n o f motion i s o b t a i n e d by c o n s i d e r i n g a g e n e r i c  elemental shear  cube s u b j e c t to three normal and t h r e e i n d e p e n d e n t  s t r e s s e s as shown i n F i g . ( 2 . 1 ) .  coordinates of a p a r t i c l e  x  =  y where x /  y  0  x  =•  (  y  x  ( x 0  Y  '  o'  then  h )  c o o r d i n a t e s a t time t = 0 .  are the i n i t i a l  Q  a r e the c u r r e n t  t }  V  o '  I f x,y  The  components o f v e l o c i t y along the x, y, z axes a r e g i v e n by  u,  v, 'Wrespectively. W r i t i n g . t h e e q u a t i o n of motion a l o n g x a x i s g i v e s (c  + x  x dx) ~Ti  x  -T.,_ yx  dz  dx  +  dy dz -  a  dy dz + (  x  (T  z  x  (dx dy d l - ^ P {  y  + -y^-  )  dx  dy  - T  3T  3o  xx j. 3~>:  by dx dy  3T  yx + 3y  z x  dx  dx  dy  (2.22)  u }  2  D i v i d i n g throughout  °- yx dy)dz ay  x  . zx 3z  =  dz  . • ,. £_ * p<•  / »> o o \  ( p u )  S i m i l a r l y e q u a t i o n of motion a l o n g y and  z d i r e c t i o n s we  get  17  FIG. 2.1  E L E M E N T A L C U B E F O R DERIVATION O F EQUATION OF M O T I O N  THE  18  Hy*  Jlxi  +  • 3x and - ZX 9x  2lx*  +  3y 3  3 T  7  3 a  zy  Considering  xx 3x  yjs  zx  g:u 9-t  p —  Z  9-  ,  3x  3T  .. yz 9z  +  p• 9t  =  9a zy , .. zz  ?y  '2.2.2.  (2.25)  P  3T  9a . yy  +  J j (v.-) ~ 9 i  =  yx . zx = — + ->— 9y -,9  •  9T  zz 9z  p as c o n s t a n t , we g e t  9x  ,  C  9 z  - 3y  9a  (2.24)  = j| pv)  (2.26)  3.w  9z  e q u a t i o n s f o r Plane Visioplasticity  F o r the p l a n e ^  strain  S t r a i n Non-Steady S t a t e  condition (2.26)  vso t h a t the e q u a t i o n o f m o t i o n 9c . xx 3X  +  —  . xy 9y K~TT~  •  y  =  3a Vx , . yy - + — = 9x 3y J  Determination  Proceeding  Y =  = T . zy  = - — — 0. 9t  can be w r i t t e n a s  9u 91 .  (2.27)  P  „ 9v p — 91 .  .  of the Stress F i e l d  from  the S t r a i n  Field  i n t h e same way as i n t h e case o f q u a s i - s t a t i c  v i s i o p l a s t i c i t y method, 3c  zx  9T  .  3T  2.2.3.  T  Dynamic  9a X _ _2  the e q u a t i o n ' c  (2.15) can be w r i t t e n as  I  ,x- •- y  •  (2.15)  19  Now from the equation o f motion (2.2 7) we get 3a  3.T  x  3X  xy  (2.28)  +  dY  p  3t  S u b s t i t u t i n g i n equation 3a.  .3x  3x  (2.15) g i v e s  £  3 u _ _ j _ [ x3 t 3x (_ .  3y  £  y  J  (2.29)  S i m i l a r l y from the e q u a t i o n o f motion (2.2 7) ,8 a  3  3y  3x  Using equation o(x,y) G  T  , + p  (2.30)  .3V  Tt  (2.29) and (2.30) w i t h the known v a l u e s o f  a t x=o and y=a, i . e .  y(x,y)  X  "  =  a  y(o,a)  3T  -  /  I  xy.  3 x  3y  a(o,a) we get y  (  8 x  3 v>  xy  x- y  3u  +P 3  From e q u a t i o n a  x  t  _ j  dy  dx  y  (2.31)  (2.14)  = a  y  +  and from equation  (2.32)  £ x-ey (2.li)  (2.33)  where  2.3.  A =  3s  (2.34)  —  GENERAL PROCEDURE FOR THE STUDY OF DEFORMATION USING VISIOPLASTICITY  The instantaneous g r i d v e l o c i t i e s  are determined  p e r i m e n t a l data and thus the s t r a i n rates,*.'equivalent r a t e s and f i n a l l y  s t r e s s e s can be determined.  from exstrain  20  F I G . 2.2 P O S I T I O N S O F P O I N T P A T I N S T A N C E S T „ A N D  T  n  +  1  21  G r i d l i n e s are marked on an end face of the specimen, which i s deformed a t a predetermined  speed.  The  deforming  p a t t e r n i s photographed u s i n g a high speed camera. p o i n t s at a l l stages of deformation  grid  are d i g i t i z e d from e n l a r g e d  photographs and the d i g i t a l p o s i t i o n a l determine the instantaneous  The  grid  data used as i n p u t to  g r i d nodes v e l o c i t i e s .  The  procedure  i s i l l u s t r a t e d i n F i g . (2.2), where a g r i d p o i n t formed by I row  and J column has  of time t t  , ,. n+1  and g r i d - c o o r d i n a t e s of x , , , y , , at the i n s t a n t ^ n+1 n+1 J  n  The  velocity  c o o r d i n a t e s x , y ..at a p a r t i c u l a r i n s t a n t  instantaneous  horizontal velocity J  u and the v e r t i c a l  v i s then g i v e n by x  u. .  , , -x n+1 n -t  :  AX A  t  v. (2.35) Since the above components o f v e l o c i t y  are o b t a i n e d  from  the d i g i t i z e d c o o r d i n a t e s o f experimental  g r i d p o i n t s an  ent smoothing procedure i s r e q u i r e d .  smoothing procedure  mentioned by Shabaik the p o i n t s .  (71)  The  i s based on a simple averaging  This procedure has  r e s p e c t to time ( c a l l e d a master g r i d ) .  with time.  . to  F u r t h e r f o r non-steady s t a t e  c o n d i t i o n s a r e f e r e n c e g r i d i s needed, which can be f i x e d  technique  of  caused d i f f i c u l t i e s i n the  past i n t r e a t i n g data which i s i l l - d e f i n e d and a l s o tends be time consuming i n o p e r a t i o n .  effici-  The  simple  with  averaging  g i v e s g r i d node p o s i t i o n s t h a t c o n t i n u a l l y change In order to surmount these d i f f i c u l t i e s a number  of a l t e r n a t e methods were c o n s i d e r e d and  f i n a l l y a three dimen-  22  s i o n a l s u r f a c e smoothing procedure was  adopted, which t r e a t s  x, y and u and a l s o x, y and v as separate s u r f a c e s and a complete  fits  f o u r t h o r d e r polynomial through the experimental  p o i n t s i . e . smoothing i s done i n three dimensions  to a s u r f a c e  formed from the s c a l e r components o f the v e c t o r f i e l d . (i.e. I  c o n d i t i o n of c o n t i n u i t y  * ^ X  *- or ¥  -  3u 3x  The  3 v for 3 y  plane s t r a i n ) can a l s o be imposed w i t h i n the s u r f a c e f i t t i n g procedure,  thereby r e d u c i n g the number of independent  para-  meters from 15 f o r each s u r f a c e (a t o t a l of 30) to 10 f o r both surfaces.  The  smoothing procedure mentioned by Shabaik  does not take account o f c o n t i n u i t y and merely i f the e r r o r i s l e s s than Besides f i t t i n g  checks  to see  15%.  a smooth s u r f a c e , the polynomials have  the d i s t i n c t advantage t h a t they can be r e a d i l y for  (71)  differentiated  determining s t r a i n s rates, and t h a t the deformation can be  a u t o m a t i c a l l y r e f e r r e d to a master r e f e r e n c e g r i d . t h a t s t r a i n s and s t r e s s e s can be determined w i t h i n the non-steady  deformation zone.  T h i s means  for fixed points  A l s o s t r e s s e s can be  e v a l u a t e d d i r e c t l y a t any p o s i t i o n o f x and y p u r e l y by subs t i t u t i n g the c o o r d i n a t e s o f a p o i n t (not n e c e s s a r i l y a g r i d point) r e q u i r e d , whereas the simple averaging technique r e q u i r e s an i n c r e m e n t a l e v a l u a t i o n o f s t r e s s e s over contiguous  grid  p o i n t s u n t i l the r e q u i r e d g r i d p o i n t i s reached. After calculating  e x  and  e  ^  u s i n g equation  (2.1) , the  e f f e c t i v e s t r a i n r a t e s at a l l g r i d p o i n t s f o r a l l i n s t a n c e s o f deformation can c a l c u l a t e d from the equation  (2.13).  23  T  In order to c a l c u l a t e 2.31)  and  x  ^.  (equation 2.32)  (equation 2.33),  (equation  the value o f A . i s r e q u i r e d .  For  non-workhardening m a t e r i a l s the value of A A i s p u r e l y a f u n c t i o n of e f f e c t i v e s t r a i n rate ( constant.  e  ) as e f f e c t i v e s t r e s s  For a workhardening m a t e r i a l , the value of  be o b t a i n e d from experimental  a vs  e  the r e l e v a n t s t r a i n r a t e c o n d i t i o n s . curve such as Thus i f 7  a = ^  e  "5 must  m a t e r i a l data taken a t I t i s normal to f i t a  where c and n are m a t e r i a l c o n s t a n t s .  i s known at a l l i n s t a n c e s of time, the value of  equivalent strain may  ( a ) is  e  be determined  at any deformation  time t and hence a  i n c r e m e n t a l l y (assuming small i n t e r v a l s ©'f  time) from the e x p r e s s i o n  F .. i  t  = '  f  where 7 ^  7  o  D  i s the  dt  ± j  effective  s t r a i n rate of a p a r t i c u l a r  p o i n t as i t moves along i t s deformation  2.4.  (2.36)  e  grid  path.  GENERALIZED EQUATIONS FOR PLANE-STRAIN VISIOPLASTICITY FOR COMPUTATIONAL PURPOSES  The equations  i n a form s u i t a b l e f o r the development o f  the computer program are given below. The  c a l c u l a t i o n o f u and v i s done u s i n g equation  X  u. . _ n + 1 13 -  v. .  =  J  W  1  - X  _  n  ^  AX A t  y -Y ^n+1 t ,.-t n+1 n  n  =  Ay_ At  (2.35)  24  Curve f i t t i n g o f the v e l o c i t i e s i s done by a l i b r a r y subr o u t i n e c a l l e d DLSQHS.  This f i t s a complete f o u r t h o r d e r  polynomial i n x and y and determines the 15 parameters (constants) f o r u and v f o r the equations given below."  The equa-  t i o n s used f o r computational purposes can be g i v e n as 2 3 4 u(x,y) = a^+a x+a x +a^x +a^x 2  + a  3  Y + a  6  Y  7  + a  8  Y  + a  9  Y  + a  10  X y  2 2 +a xy +a xy +a x y+a x y 11  12  1 3  14  +a x y  (2.37)  1 5  Similarly: 2 v(x,y) = b + b x + b x b 1  2  3  x 4  3 4 +b x 5  2 3 4 +b y+b y +b y +b y 6  y  g  g  2  +b xy+b xy +b xy 1 0  1 1  3  1 2  2 2 2 3 + b x y + b x y +t> x y 13  14  (2.38)  15  Thus the v e l o c i t y component u anv v can be expressed •  •  s e p a r a t e l y and the c a l c u l a t i o n and  f o r s t r a i n rates  •  Y  e  e  x,  y,  ' xy  can be done as f o l l o w s : "K  =  If  =  2  3  a 2a x 3a x 4a x a y +a y +a y +2a xy+2a xy 2 +  3  +  4  2  5  +  1 0  3  1 1  1 2  2  +3a x y 1 5  +  1 3  2  1 4  (2.39)  25  '£y  If  =  =  V  2  b  7  y  +  3  b  Y  8  2 +  4 b  9  y 3 + b  10  X  2 2 + 2 b x y + 3 b x y +k> x 11  12  13  2  +2b x y+b x 1 4  and  ^  =  (2.40) 2  + |Z  = a +2a_y+3a y +4a y 3x . 6 7 S" 9^ c  Q  J  3y  'xy .  3  1 5  3  Q  7  +a x+2a xy+3a xy 1 Q  1 1  2 + a  13  x  +  2  a  i4  2  1 2  2 y  x  3 + a  l5  x  2  +b +2b x+3b .x +4b x 0  + b  0  10  y + b  /  ll  y 2 + b  3  n  12  y 3  2  2  +2b xy+2b xy +3b x y 1 3  The c o n d i t i o n o f c o n t i n u i t y J  1 4  (  e  *  = X  (2.41)  1 5  e  *  or T — y  3u  3v  =--—  3x  )  is  3y ;  imposed on the curve f i t t i n g by r e q u i r i n g the c o e f f i c i e n t s . t o be r e l a t e d as f o l l o w s . a  2  =  b  " 6  2a = - b  1 Q  3a = - b  1 3  4a = - b  1 5  3  4  5  a±1= 3a  I 5  3b  =-2b  8  14  •  (2.42)  26  a  12  =  "  4 b  9  a  10  =  -  2 b  7  2 a  13  =  "  2  b  l l  2a = - 3 b 14  1 2  The p a r t i a l d e r i v a t i v e of -y w i t h r e s p e c t t o x and y i s given by xy 3* xy 3x Y  = a +2a y+3a y 1 0  1 1  2  1 2  +2a x+4a xy 13  14  2 2 + 3a, x +2b-+6b x+12b -x 15 3 4 5 c  /1  1  2  +2b y+2b y +6b xy 1 3  1 4  (2.43)  1 5  and hence 8Y  xy T'iT * a  =  ( a  10  + 2 b  3  ) y  +  ( a  ll  + b  13  ) Y  +  (  3 a  12  + 2 b  + (2a +6b )xy. + ( 4 a + 6 b ) x y 13  4  1 4  2  + (3a +12b )x y 1 5  )  14 y  2  1 5  (2.44)  5  Similarly 3Y 2 ' xy = 2a +6a y+12a y + 2 a x + 6 a x y 7  g  g  11  2  +2a x +b +2b 1 4  1 0  i ; L  12  y+3b y  2  1 2  2 +2b x+4b xy+3b x 13  14  15  (2.45)  27  and  /  _i]|y  d  x  ^a^b^)  =  x + (6a 2b )xy 8 +  ( 2 a 2 +  ( a  F u r t h e r equation .. »• ( e i+ l -  £  i +  ll  + b  13  }  x  +  14 3  +  (12a +3b ) y x  +  (3a +2b )x-y  g  1 1  +  3 b  15  }  x  3  12  1 2  (2.46)  1 4  (2.36) can be expressed as £  A  i+1)  t -4- £. (2.47)  K..  2  The normal s t r e s s -<?y can be c a l c u l a t e d u s i n g equation :  (2.31) ,  i.e. y  3 T  y ( o , a ) - / •. x y + a ax  Y -  v  p  X  J_v at  ; y  / o  xy. 3y  < 3 3x  r  x-"• y -, ..' 3 u dx ^ ^3 t  y  = a  (2.31)  The terms f o r t h i s equation may be c a l c u l a t e d s e p a r a t e l y .  The  f i r s t termo" y (o,a) i s determined e x p e r i m e n t a l l y at each i n t e r v a l of time.  The second term can be c a l c u l a t e d from the equa-  tion. Y T  ,;  xy  so  xy  (2.33)  =  2A  that 3T  xy  (  =  7 ^  ( _  Y  x y  2T  Y xy  ^  3x  )  3  £.  a  = 1 3 3 3x  -Y  xy  c  —n e  28  Y  = £ —a_  _n  x  3  9X  c  Y  y  fe  r  Y  n  c 3  xy  e  e _ Y xy  xy 3>^1  I  3  S, ax £  . 2 e  '  gX^  3  3Y  - n-1 „ •3 3  7  X  n 8  Y„„  T  +  MHY  (2.48)  .3x  3x  Similarly 3T  Y  _ c " 3  XV  'V  — n-1  xy .  —  3  F  n-  e  _ n  3Y /—<n Y ^xy _(; ) xy  + _s_  £  sy  (i) X y  Using equation  (2.1), the t h i r d term o f equation  2  /— \ ^ ^  3 3  e  Y  .  (2.49)  (2.31) can be  computed as £  ( x^  3  <y)  £  £  X  y  3 A  A  "^X  (  £  3X  ,3 • (/ 7 ) . 1 3x - n TA  £  x  + -r A  J  y) .2  3  E  +  A  L 3x  2Ce  But  n  [7(7) 3X  2c  2c n  e ( e)  L t  J?  (2.50)  - *  ]  3s  -  +  3e 3x  n  +  ( e ) _3_ ( e ) ]  3x 3  + 2c ( i ) n ^  3x  (7 )  3x  3s_  _A  -  3X  }  3E  1  (f^T  n A  ax  E;  IL_ e  differentiating  37 3X  1 +  -  3s 3x  (2.51) :J .  s  equation  (2.13) w i t h r e s p e c t t o x, we get l  3s 3x  e )TI  -x  J  -y -  gx -n-1  •e A[_  "  t  E  = _ ! r t~ 2c 3 n " 2c"  V b  3x • x  3_ 3X  r2 (3 s x  2  2  4  x y  = 1 3  [ (3 e x  2  2  i 3 ; ) 4 x y  2  (6 e _ _ _ x 3 x 2y ' e 8x 4 e  5 y  xy) ' 9x  v  x  y  -h  = (3 e  V  2  + 3i  x  4  2  x y  )  H e  /  |_  x  ^x + 1 i  _!_xy7T  2 . 5 -TirJ xy  —  (2.52) So t h a t the t h i r d term o f equation (2.50-2.52) 2  £  ~  (2.31) , u s i n g  equations  can be computed as  £  ~K  (e„ - e..){ —  -g^ + - — £  > + — ^ d  -X  y  d  (2.53) Thus the normal s t r e s s e s equations  and  c? can be computed u s i n g x  (2.31, 2.32, 2.48, 2.49, 2.50, 2.51, 2.52 and 2.53)  w i t h the known v a l u e s o f  "2.5  a  a(o,a).  COMPUTER PROGRAM  The flow c h a r t o f the computer program developed  f o r plane-  s t r a i n dynamic and q u a s i - s t a t i c v i s i o p l a s t i c i t y i s given i n F i g . (2.4)  and a l i s t i n g o f the program i s g i v e n i n Appendix I. The running i n s t r u c t i o n s and the main steps i n t h e program  are d e s c r i b e d below: (1)  Input a l l data r e q u i r e d f o r the c a l c u l a t i o n s (a)  Read  IX  = No. o f experimental g r i d l i n e s p a r a l l e l t o y a x i s  IY  = No. o f experimental g r i d l i n e s p a r a l l e l t o x a x i x The above are given i n FORMAT (312)  30 IT  = No. o f the time  steps  DT  = Time i n t e r v a l between two c o n s e c u t i v e photographs ( t h i s need not be constant time  interval)  FORMAT (8F10.0) NIX  = No. o f g r i d l i n e s p a r a l l e l t o y a x i s i n master g r i d  NIY  = No. o f g r i d l i n e s p a r a l l e l t o x a x i s i n master g r i d Both i n FORMAT (312)  CA  = Constant  a used f o r a (o,a)  CC  = Constant  used i n t r u e s t r e s s and t r u e s t r a i n  CN  = Constant  (or index)  SIGOA =  relation  f o r workhardening  a (o, a ) , known value o f y a t x = o and y = a  A l l above i n FORMAT (8F10.0) (b)  Read, the instantaneous  c o o r d i n a t e s o f the g r i d p o i n t s  of c o n s e c u t i v e photograph.  The format  i s g i v e n as  FORMAT (5X, 2F6.3, IX, 2F6.3, IX, 2F6.3, IX, 2F6.3, IX, 2F6.3) )  C a l c u l a t e the values o f h o r i z o n t a l v e l o c i t y u and v e r t i c a l v e l o c i t y v a t each g r i d p o i n t u s i n g equation  )  F i t a 4th degree polynomial through  the three  (2.35) dimensional  curves f o r u and v as a f u n c t i o n o f x and y u s i n g the computer l i b r a r y r o u t i n e c a l l e d  'DLSQHS'.  This program p r o v i d e s a l e a s t square t i o n o f M parameters  f i t to a l i n e a r  func-  ( i . e . M independent v a r i a b l e s ) and  N data p o i n t s by the Householder t r a n s f o r m a t i o n  techniques.  r  22 21 20 19 18 17 16 15 14 13 12 11 10 9 8  L .  0.1  2  3  4  5  6  7  IX  FIG. 2.3.  GRID . L I N E S ON  THE  SPECIMEN  8  32  FLOW CHART OF COMPUTER PROGRAM  Read IX, IY, IT, DT, NIX, NIY, CA, CC, CN and (o,a)  ( i ) Set Master G r i d f o r f i n a l p l o t and f u r t h e r calculations ( i i ) F i l l zer.o f o r the p o i n t s , o u t s i d e t h e boundary o f the deforming zone by . use o f s u b r o u t i n e F I L L  T Read ( x , y ) , t h e g r i d p o i n t s a t each time i n t e r v a l . P l o t them  X C a l c u l a t e component v e l o c i t i e s u and v f o r each time increment using the equations X -x n+1 n t -t i j n+1 n y -y (equ. 2 . 3 5 ) n+1 n -t .V « = t n+1 n U  Ploti u .v  £, 6  3  00  F i t t h e c u r v e s t o t h e u and v v a l u e s w i t h complete f o u r t h degree p o l y n o m i a l i n x and y ( S u b r o u t i n e DLSQHS)  = f(x,y),  i n two d i m e n s i o n a l form u s i n g s u b r o u t i n e PLOT, PLOTAND & SYMBOL i n three dimensional form u s i n g s u b r o u t i n e PERS  Plot:• ( i ) <7~ as a f u n c t i o n o f x and y ( i i ) (J"^ as a f u n c t i o n o f x and  (equ.  2.1, 2.13,2.38,  y  (iii)T^*  y ay -v-tr  = f(x,y),  1  2.39,2.40)  1  y  = f(x,y),  C a l c u l a t e normal & shear stress using e q u a t i o n ' 2 . 1 , 2 . 1 3 , 2 . 3 8 , 2 . 3 9 , 2.40  C a l c u l a t e and s t o r e s t r a i n rates f o r a l l points, using s u b r o u t i n e DERIV  3U Y xy J~ ~sy  x  t (a)  E  x~ 3X  = f(x,y), = f(x,y),  as a function of x xy  and y (a) i n two dimensional form  o»V  +ax  using subroutine PLOT, PLOTAND & SYMBOL  (b) i n three dimensional form  xy  using subroutine PERS  E C a l c u l a t e and s t o r e t h e t o t a l effective strain at a l l s p e c i f i e d node p o i n t s  STOP  FIG.  2.4  33 That i s i t minimises  N  2  M (v. 1  E i=l  a.x..) E 1 j=l  DLSQHS transforms the matrix X to an upper t r i a n g u l a r v i a Householder t r a n s f o r m a t i o n s , and then s o l v e s by backward s u b s t i t u t i o n . by  'TRUE':  the system  I f the command REFINE i s d e f i n e d  a c o r r e c t i o n v e c t o r i s computed from the r e s i d u a l  e r r o r s between the dependent v a r i a b l e s Correction  form  v e c t o r s are then a p p l i e d  and the f i t t e d  values.  t o the s o l u t i o n and r e -  computed i t e r a t i v e l y u n t i l convergence i s o b t a i n e d .  DLSQHS  i s most e f f e c t i v e f o r problems, where the c o r r e c t i o n o f the m a t r i c e s i s unknown and the s c a l e o f d i f f e r e n t varies (4)  (5)  widely.  impose t h e " c o n d i t i o n  i.e.  au^  -  of continuity • • --  _§v  Calculate  the s t r a i n - r a t e s £ "/ "£  effective strain-rate and  variables  (2.73)  :  1  e = - e x y "-  or - '  , shear s t r a i n - r a t e  u s i n g equations  ;  xy and  (2 . 3 9 ) , ( 2 . 4 0 ) , ( 2 . 4 1 )  1  ,are combined i n t o a DERIV' ( d e t a i l s are g i v e n  i n Appendix) . (6)  Calculate  the t o t a l e f f e c t i v e s t r a i n by i n t e g r a t i n g  along  the path o f p a r t i c u l a r p a r t i c l e o r g r i d node ( r e f e r equat i o n s 2 . 3 6 and 2 . 4 7 ) . (7)  Calculate 2.41).  shear s t r e s s  T x y u s i n g equations  ( 2 . 3 3 , 2 . 3 4 and  34  (8)  C a l c u l a t e independently the terras from equations (2.48),  (2.49),  (2.50),  (2.51),  (2.31),  (2.52), and (2.53) and  thereby c a l c u l a t e  a  routine  T h i s subroutine i n t e g r a t e s a f u n c t i o n  'DQUANK'.  y u s i n g the computer l i b r a r y  sub-  "3D !  f(x)  when the l i m i t s a and b are given.  dx.  I t i s b a s i c a l l y based on Simpon's three p o i n t s i n -  tegration fifth The  &  i . e . I = j f(x)  improved by u s i n g an adjustment  term o f  degree i n p l a c e o f the three degree term.  absolute e r r o r can be l i m i t e d to any a r b i t r a r i l y  s p e c i f i e d value. C a l c u l a t e s t r e s s e s ° x u s i n g equations (9)  P l o t the d i f f e r e n t q u a n t i t i e s , u, v, 7f  CT  y , ° x,  1  e  (2.31) and (2.31). x,  e  y,  Y  xy, F ,  xy, e t c . f o r d i f f e r e n t x, y v a l u e s .  The  p l o t t i n g v e c t o r s are taken a t d e f i n e d master g r i d node points.  Any p o i n t w i t h i n the master g r i d b u t o u t s i d e the  deformation  zone are g i v e n zero v a l u e s by the subroutine  'FILL'. (a)  These v a l u e s can be p l o t t e d i n two dimensions  as a f u n c t i o n o f x o r a f u n c t i o n o f y. the subprogram 'PLOT' can be used. the b a s i c p l o t s u b r o u t i n e .  either  F o r t h i s purpose  Thus subprogram i s  I t generates  the pen movement  r e q u i r e d t o move the pen i n a s t r a i g h t l i n e from i t s p r e sent p o s i t i o n t o the p o s i t i o n i n d i c a t e d i n the c a l l . i s a l s o used t o r e l o c a t e the o r i g i n of the p l o t t e r dinate, system i n the X d i r e c t i o n . i s complete,  a second  To ensure  It  coor-  that p l o t t i n g  subprogram 'PLOTND' i s used.  For  35 c l e a r d i s t i n c t i o n of d i f f e r e n t l i n e s i n a p l o t , a t h i r d subprogram 'SYMBOL' can be used. c h a r a c t e r s and s p e c i a l  This plots  alphanumeric  symbols.  (10) Three dimensional o r t h o g r a p h i c d i s p l a y s can be made u s i n g the subroutine  'PERS'.  The above values are taken as 'Z'  values and are p l o t t e d i n three dimensional f u n c t i o n o f x, and y.  form as a  C H A P T E R  T H R E E  VARIABLE SPEED, CONTROLLED VELOCITY PROFILE, SINGLE CYCLE, IMPACTING  PRESS  37  CHAPTER THREE  VARIABLE SPEED, CONTROLLED VELOCITY PROFILE, SINGLE CYCLE IMPACTING PRESS  3.1  DESCRIPTION OF EQUIPMENT  3.1.1.  Background  Information  For experimental as heading  i n v e s t i g a t i o n of forging operations  such  and u p s e t t i n g , where s t r a i n - r a t e , forming-speed  and  forming l o a d are important,  a device w i t h s p e c i a l  requirements  i s needed. Obviously as wide a range o f impact  speed as p o s s i b l e i s  r e q u i r e d t o g e t h e r With a v e l o c i t y p r o f i l e which i s s e n s i b l y i n dependent of forming l o a d and which may  be a d j u s t e d to s u i t  circumstances. The commercial a l t e r n a t i v e s t h a t are a v a i l a b l e have been developed  f o r s p e c i f i c a p p l i c a t i o n s and o f n e c e s s i t y have a  limited f l e x i b i l i t y .  For example, the crank p r e s s , used i n  many f o r g i n g o p e r a t i o n s has a v a r i a b l e s t r o k e and maximum v e l o c i t y and has a w e l l c o n t r o l l e d v e l o c i t y p r o f i l e , which i s s e n s i b l y independent o f forming  load.  v e l o c i t y i s l i m i t e d to approximately  However, the maximum 15 to 20 f t . / s e c .  With  the drop hammer and h i g h v e l o c i t y forming machines (such as the p e t r o f o r g e ) a h i g h e r v e l o c i t y can be reached for  (15-30 f t . / s e c .  the drophammer, 90-100 f t . / s e c . f o r the p e t r o f o r g e ) but  stroke and v e l o c i t y p r o f i l e d u r i n g deformation p u r e l y by the r e s i s t a n c e of the workpiece.  are  the  determined  That i s the h i g h  SOLENOID.  FIG 3.1  VARIABLE S P E E D , C O N T R O L L E D VELOCITY PROFILE, SINGLE C Y C L E IMPACTING P R E S S  39  v e l o c i t y r e s u l t s i n h i g h energy very l i t t l e o f the energy  and with low s t r e n g t h t a r g e t s  i s absorbed  i n p l a s t i c work.  Further  the h i g h v e l o c i t y d e v i c e s cannot be used f o r low v e l o c i t y work. With these design c r i t e r i a in:mind a v a r i a b l e speed, controlled velocity profile,  s i n g l e c y c l e , impacting press was  designed and b u i l t w i t h i n the department. 3.1.2  D e s c r i p t i o n o f the Press:  A diagram o f the impacting press i s shown i n F i g . 3.1. The d r i v e comprises  a m o d i f i e d Whitworth q u i c k - r e t u r n mechanism  c o n s i s t i n g o f a crank and a d r i v e arm together w i t h a v a r i a b l e speed D.C. motor, a f l y w h e e l , bearings e t c .  The end o f the  d r i v e arm i s attached by a connecting r o d t o a c y c l o i d a l cam. In s i n g l e c y c l e o p e r a t i o n , the cam i s made t o engage with an upper p l a t e n (or ram) which impacts  the workpiece.  The upper  p l a t e n and cam are both mounted on m u l t i r o d supports w i t h l i n e a r b a l l bushings. f o r emergency 3.1.3  A brake  i s p r o v i d e d on the f l y w h e e l  purposes.  Operation o f the P r e s s : The d r i v e wheel i s r o t a t e d a t a p a r t i c u l a r speed by ad-  justment  o f the D.C. motor c o n t r o l l e r causing the d r i v e arm  to o s c i l l a t e about i t s lower p i v o t . mechanism connects  The s i n g l e c y c l e tripping.?  the d r i v e arm w i t h the cam and the cam en-  gages w i t h the cam f o l l o w e r s on t h e upper-platen.  The upper-  p l a t e n i s thus f o r c e d down towards the workpiece.  The p l a t e n  achieves a maximum v e l o c i t y when the d r i v e arm i s i n a c e n t r a l p o s i t i o n a f t e r which time the p l a t e n i s brought  to rest.  the r e t u r n s t r o k e o f the cam the p l a t e n i s r e t u r n e d t o i t s  On  40 i n i t i a l position. connecting  The t r i p p i n g  mechanism then disengages the  r o d from the cam and the d r i v e arm continues t o  o s c i l l a t e f r e e l y about i t s lower p i v o t . The  s t r o k e , v e l o c i t y and a c c e l e r a t i o n p r o f i l e o f the upper-  p l a t e n are determined s o l e l y by the cam contour s e t t i n g o f the D.C. motor.  and the speed  A c y c l o i d a l cam i s used f o r h i g h  v e l o c i t y work t o minimise e x c e s s i v e wear, through shock and vibration. The  lower p l a t e n h e i g h t can be a d j u s t e d r e l a t i v e t o t h e  upper p l a t e n thereby  determining  the working p o r t i o n o f the  stroke.  3.2  KINEMATIC ANALYSIS OF THE MECHANISM  3.2.1. The The  Derivation of Expressions  f o r V e l o c i t y and A c c e l e r a t i o n :  d r i v e mechanism i s shown s c h e m a t i c a l l y i n F i g . (3.2.).  centre l i n e o f the d r i v e arm i s i n d i c a t e d by l i n e AB, where  p o i n t A denotes the f i x e d lower p i v o t . d i c a t e d by the d o t t e d l i n e with B points.  1  The path o f B i s i n -  and B", showing the extreme  The p o i n t C i s the c e n t r e o f r o t a t i o n o f t h e drivewheel,  a d i s t a n c e P above the f i x e d p i v o t o f the d r i v e arm.  The eccen-  t r i c i t y o f p o i n t E, the cam f o l l o w e r , about C i s given by a d i s t a n c e Q.  The l i n e CD i s a r e f e r e n c e f o r angle e.  of AB i s L.  The path o f p o i n t F denoting  shown by d o t t e d l i n e .  The l e n g t h  the cam p o s i t i o n , i s  The l e n g t h o f connecting  rod, BF = R and  the d i s t a n c e o f cam p o s i t i o n F from l i n e AC (or AH) be X. angle made by AB with l i n e AC i s 4>  The  42 Then from geometry = arc tan QCosQ P+QSinQ  (3.1) 1  D i f f e r e n t i a t e <\> with r e s p e c t to t , we ,= _ PQSinQ+Q d t  d&  2  2  '"~ " P +2PQSin9+q  2  D i f f e r e n t i a t i n g again, we  (  dt Now  (2.3)  get  "> : (P +2PQSine+Q ) 2  from F i g . (3.2)  2  HF = HG + GF = L S i n  BG )1 2  but GF = [(BP'  ^  )  dt  get  h  =  [R  2  - ' { L (1-Cos cb) }  D i f f e r e n t i a t i n g with r e s p e c t to t , we dx ± dt  L d_j dt  [ R 4 L ( 1-Cos 2  = L dj> dt  [Cos  - W V  cb  where W = L (1-Cos and V = [(R  2  Sin  1  ]^  sin * 2  cb)} ;^  <b  h  get  - L ( l - C o s cb )  Cos*  cb + GF = X 2  O  [R  (3.3)  -{• L ( l - C o s <f> ) } ]  O  Hence, X = L S i n <j) +  2  ^  (3.4)  ]  cj> )  2  - W )]^  r  D i f f e r e n t i a t i n g again and r e a r r a n g i n g the terms, we 2  d x 2 dt  = L d  2  cb  ,2  dt'  (Cos  <f>  +  W) - L ,d<b. v dt c  [(Sin  ;  + 1)  (W V  2  <b  get  + L Sin V  + W Cos V  <b  tb  )  ]  (3.5)  43  3.2.2.  Maximum V e l o c i t y and A c c e l e r a t i o n o f the Mechanism:  The press has f o l l o w i n g dimensions:P = 14 i n c h . Q =  6 inch.  L = 24 i n c h . R = 26 i n c h . Ratio o f d r i v e r t o driven, p u l l e y d i a . = 0.573 For a flywheel  speed o f 250 RPM,  ( i . e . a motor speed of 434.4  rpm) the angular v e l o c i t y o f the f l y w h e e l  2 x 250 = . ,—r = 15 rad/sec. 6 Q  (3.6) The p o s i t i o n of the maximum v e l o c i t y of the ram (hence the upper-platen) occurs when Hence from equation d£ dt  <j> = o o r 0 = - 90'  (3.2) we get 2  =  _•• 'd9; PQ S i n e + Q dt p Qsine+Q 2  2  + 2 P  = (- a)  = 0. 75 a)  (14x6) (-1)+(6x6) (14x14) +(2x14x6 (-1) ) + (6x6) (omega)  From equation (3.4) dx dt  d (j) (Cos (j) - W S i n cf> ) ; here dt V  _ L  _ _ dA " dt L  11  c>j = o  (1 - W x o ) -J d(j) v " dt  It  12 dt  T  L  = 2 x .75 u = 1.5 u ( f t . / s e c . ) ' -  (  3  *  7  44  45 Now the displacement c y c l o i d a l cam h —  ^  o f p l a t e n , y, i n r e l a t i o n to the  ( r e f e r . F i g . 3.3) i s given by TT£_ _ 1 S i n 2ja_r 2  L'  IT  -  L  (3.8)  1  if  where I = d i s t a n c e t r a v e l l e d along x a x i s at a p a r t i c u l a r time L  1  = length o f c y c l o i d a l cam i n x d i r e c t i o n  h = maximum d i s t a n c e t r a v e l l e d by upper-platen o r f o l l o w e r i n the y d i r e c t i o n V e l o c i t y o f f o l l o w e r i . e . p l a t e n i s given as (£4.) v  v  h dt  =  ;  , .,  -jTT-f  =  v  max  T  COS2TT£  1 -  —  -  1  (3.9)  )  dt ,  (3.10)  S i m i l a r l y the equation f o r a c c e l e r a t i o n of the u p p e r - p l a t e n , a, i s given by _ . . a- v = 3  (£-*) 2h TT d t n  y  S  ^ i . . . 2 v.I  i  n  j i  2  , h  + — (11  L  L  L'  2 Cos2 -a d Z ) —xT I 2 L  . . (3.11)  dt  For the present cam, ,d£. h = 4.0 i n c h ,, dx. and L ' = 8", so t h a t 2 x 4 . 0 x d t _ 2x4.0x dt' max 8 8 l  ;  (  v  2x4. 0x1.5 ui 8 = 1.5 a) Using the value o f u = 1 5 v  max = 22.5  rad./sec. from equation (3.6)  ft./sec.  PUSH-BUTTON R E D  "Eft n k f f  68  k  off"  3  P<wer i n  500  i  ORANGE •  9  BRN *ED  JT  2 0u  -/Fuse Power  BRN  9 9  l  ILl l a t c h et  2  IRL, l l  TTT  JLACK  20u  >( B L A C K 6  METAL-CAS  500V  10  BLUE  NO I  £l50  BLUE  C  Reset'  7 BRN  (LACK  t 1 1 7 V , 60A.  BLUE  415/300V ORANGE  4  6  RL3 '  ORANGE  RED  ^C  T 3 DC I I R L  MICROSWITCH  t  + 170V-p  fSolenoid iPower  Power  5K  HJETAL-CHASIS '  I o  1 - -L- +  SYNCHRN.  r-p—Urj—° 1  J-100 T  1  FIG. 3.4 ELECTRICAL  u  -M T  <£k  0 6 6 6  CAMERA  ND  L N  110V AC(OUT) SOLENOID  "  I  CIRCUIT  _o  FOR  SYNCHRONIZATION U  N  I  T  SINGLE - CYCLE  j  I  L_o I 0 9  u 10V T o CAMERA  Camer  .  s  L 0  I, 0  I—  OPERATION  CONTROL  O  From cam. Control  47  3.3  CONTROL  With s i n g l e c y c l e o p e r a t i o n i t i s e s s e n t i a l t h a t engagement o f the d r i v e arm with the cam occur a t the c o r r e c t p o i n t i n t h e c y c l e and have s u f f i c i e n t time t o engage p r o p e r l y .  This i s  achieved by r e q u i r i n g t h a t a sensing device p l a c e d a t the extreme p o i n t o f the d r i v e arm motion be a c t u a t e d i n c o n j u n c t i o n with a push-button s t a r t switch before engagement occurs.  A solenoid  then r e t r a c t s and allows the maximum time f o r the two p a r t s t o engage. S y n c h r o n i z a t i o n o f the h i g h speed camera with the specimen deformation  i s a l s o accomplished  the remote sensing switch.  . u s i n g the  signal  from  A time delay device i s used t o vary  the s t a r t o f f i l m i n g so t h a t d i f f e r e n t impact speeds can be accommodated.  A f u r t h e r sensing device p l a c e d a t the o p p o s i t e  extreme o f the d r i v e arm movement i n d i c a t e s when the event i s completed and t r i g g e r s the magnetic camera brake. D e t a i l s o f the e l e c t r o n i c c i r c u i t r y  f o r engagement and  disengagement o f the d r i v e arm and f o r s y c h r o n i z a t i o n o f the high speed camera are given i n the f o l l o w i n g s e c t i o n s . 3.3.1.  Control f o r Single Cycle  Operation:  The e l e c t r i c a l c i r c u i t used f o r c o n t r o l l i n g the engagement of the d r i v i n g arm o f the cam i s given i n F i g . (3.4) and the flow diagram o f the r e l a y sequence i s given i n F i g . (3.5). sequence o f events (i)  The  i s as f o l l o w s :  S i n g l e c y c l e s t a r t switch  (push button) i s pressed which  s e t s the l a t c h r e l a y (RL2) through r e s e t r e l a y (RL3),  faormally  OPTICAL-SENSORS 1  i  JL  HIGH SPEED AC CAMERA LOGIC  START PUSH BUTTON  SOLENOID  2  SYNC. < CONTOL  AC AC  POWER RELAY L4  HIGHSPEED CAMERA  ^  LATCH RELAY HIGH SPEED CAMERA RELAY  ON RATCHET OFF TRIGGER ^ R PULSE  RL  2  L1  POWER ON RESET  LATCH CONTROL  RESET RELAY RL 3  117 V AC  POWER SUPPLY  POWER ON PULSE FOF.MER  170V DC  PULS E FORTVtER i  MICRO SWITCH ( A )  F I G . 3.5  FLOW DIAGRAM FOR  R E L A Y S E Q U E N C E FOR S I N G L E C Y C L E C O N T R O L  CAMERA RELAY  OPTICAL SENSOR 1  | SINGLE | CYCLE CONTROL UNIT CAMERA RELAY  OPTICAL SENSOR 2  FIG.3.6  LAMP DRIVER  SAMPLE PULSE 1ms  CIRCUIT OPERATION BLOCK  DIAGRAM FOR  LAMP INDICATOR  HIGH SPEED CAMERA SYNC.  50 closed contact p o s i t i o n ) . (iii)  The  r a t c h e t r e l a y steps to 'ON'  the micro switch  as the d r i v e arm  (A) a t the bottom o f the s t r o k e .  a s i g n a l on the s o l e n o i d r e l a y  (RL4)  through  contacts  This puts  the l a t c h  relay  (RL2), and e n e r g i z e s the s o l e n o i d and hence the t r i p p l i n g mechanism. (iv)  When the d r i v e arm  c o n t a c t s the micro  second time the r a t c h e t r e l a y steps to the  switch 'OFF'  (B) f o r the  position.  reset relay  (RL) momentarily e n e r g i z e s and the l a t c h r e l a y  deenergizes  via  3.3.2  The (RL2)  RL3.  C o n t r o l f o r S y n c h r o n i z i n g the High Speed Camera: The b l o c k diagram o f the s y c h r o n i z a t i o n c o n t r o l o f high  speed camera i s g i v e n i n F i g . (3.6). sensor microswitch  I t c o n s i s t s o f an  (A) to i n i t i a t e a f i l m i n g s i g n a l and  i n t e g r a t e d c i r c u i t timers type 556  optical an  to give a v a r i a b l e time  delay to the a c t u a t i o n of the camera r e l a y and hence the camera. Reference  p u l s e s of 1 m.sec d u r a t i o n are compared with  the p u l s e s r e c e i v e d from the o p t i c a l sensor no. (3.6)) and when these correspond, showing c o r r e c t s y n c h r o n i z a t i o n .  2.  (in F i g .  a lamp i n d i c a t o r i s t r i g g e r e d  CHAPTER FOUR  EXPERIMENTAL PROCEDURE AND DISCUSSIONS  4.1.  E XPE RIMEN TAL  P ROCE DURE  A plane s t r a i n pact  s p e e d 15.7 f t . / s e c .  from p l a s t i c i n e l"xl-l/2"x2" cone  u p s e t c o m p r e s s i o n - t e s t was made a t an i m -  using  The p l a n e - s t r a i n  a m e t a l mould.  specimen  The s p e c i m e n  was made measuring  was c u t w i t h a f i n e w i r e and a l u b r i c a n t o f s i l i -  g r e a s e was u s e d t o p r e v e n t t h e s p e c i m e n  from s t i c k i n g i n  t h e mould. The  end s u r f a c e o f t h e specimen  was s p r a y e d w i t h b l a c k  paint  a n d a s q u a r e g r i d p a t t e r n was made by s p r a y i n g  paint  on  inch.).  the black surface  on t h e l o w e s t p l a t e n o f t h e h i g h s p e e d lubricated pexiglas plates  apparatus).  The u p p e r  the specimen  a t t h e maximum v e l o c i t y  and f i l m i n g  appropriately  contrast  c o n d i t i o n s were o b t a i n e d by p l a c i n g t h e  between two p a r a l l e l  that  gave good  (14 mesh,  photography.  Plane-strain  trols  a w i r e mesh g r i d  The u s e o f b l a c k a n d w h i t e p a i n t  f o r h i g h speed  specimen  through  white  the objective  t h a t o f t h e specimen  (called  a Kudo  p l a t e n o f t h e ram was a d j u s t e d t o i m p a c t  synchronization  set.  impacting press  i n the cycle  (as d i s c u s s e d  and a l l c o n -  i n C h a p t e r 3)  The h e i g h t o f t h e camera was k e p t , s u c h l e n s o f t h e camera was on t h e same h e i g h t a s and t h e p l a n e o f specimen  the p l a n e o f t h e l e n s . ( s e e F i g . A )  was p a r a l l e l t o  52  FIG.  (A)  LOCATION  OF THE SPECIMEN IN KUDO APPARATUS  53  4.2.  DISCUSSION The g r i d p o i n t s ,  f o r the dynamic plane s t r a i n upset com-  p r e s s i o n t e s t , were d i g i t i z e d at 0.00133 second time (every 4th frame of the h i g h speed f i l m ) . a quasi-static (71)  plane-strain  compression  (see F i g . 15) a t a speed of 0.02  4.3  to 4.9  t e s t done by  from  Shabaik  f o r the f o u r steps of  deformation i n F i g . 4.1  and 4.2  and i n F i g s .  f o r the seven steps o f dynamic deformation.  movement of is plotted  Also results  f t . / m i n . were d i g i t i z e d .  The d i g i t i z e d g r i d p o i n t s are p l o t t e d quasi-static  intervals  The  . c e r t a i n grid-node p o i n t s d u r i n g deformation i n Figs.  4.10  and  The smoothed h o r i z o n t a l (v) are given as a f u n c t i o n  4.11. velocity  (u) and v e r t i c a l v e l o c i t y  of x and y i n F i g s .  4.12  to  4.15  f o r l a s t time i n t e r v a l i n both s t a t i c and dynamic t e s t s . Plots 4.16  of e f f e c t i v e s t r a i n rate  and 4.17,  w h i l e t o t a l e f f e c t i v e s t r a i n accumulated  m e n t a l l y f o r a l l time i n t e r v a l s 15.7  ( e~ ) are g i v e n i n F i g s .  f o r the 0.02  incre-  f t . / m i n . and f o r  f t . / s e c . deformation speeds are shown i n F i g s .  4.18  and  4.19. The normal  and shear s t r e s s e s  as f u n c t i o n s of x and y f o r the 0.02 in Figs.  4.20,  4.21  t i o n speed i n F i g s . f o r 0.02  0  ( y,  0  T  x,  4.24  and 4.25.  f t . / m i n . deformation speed was  _ and s t r a i n as  The v a l u e s of a (o,a) taken from r e f e r e n c e between  — o 25  a = 31000 e  a (o,a„) f o r the p l a s t i c i n e specimen  speed  f t . / s e c . deforma-  (71) as 31000 p s i a t the o r i g i n w i t h the r e l a t i o n s h i p  effective stress  plotted  f t . / m i n . deformation  and 4.22,and f o r the 15.7 4,23,  xy) are  was  "  .  The v a l u e s o f  taken as 20 p s i a t the  54  o r i g i n from r e f e r e n c e (79).  The s t r e s s e s  ° x and-T'cxy are a l s o  p l o t t e d as a f u n c t i o n o f x f o r the CL 02 f t . / m i n . deformation speed i n F i g s . 4.26 and 4.27 and f o r 15.7 f t . / s e c .  deformation  speed i n F i g s . 4.28 and 4.29. The  r e s u l t s from the specimen deformed a t 0.02 f t . / m i n .  agree c l o s e l y with w e l l documented r e s u l t s f o r q u a s i - s t a t i c deformation showing a ' f r i c t i o n h i l l  1  type o f normal  a  ( y ) inter-  face s t r e s s d i s t r i b u t i o n w i t h maximum s t r e s s o c c u r r i n g a t the c e n t r a l p o r t i o n o f the specimen.  The normal s t r e s s  distribution  f o r the specimen deformed a t 15.7 f t / s e c i s r a d i c a l l y  different,  showing a saddle type d i s t r i b u t i o n o f normal s t r e s s , with t h e maximum s t r e s s o c c u r r i n g near the p e r i p h e r y o f the c o n t a c t zone. The i n t e r f a c e shear s t r e s s d i s t r i b u t i o n s a l s o change form w i t h strain  rate.  The dramatic change o f normal s t r e s s d i s t r i b u t i o n w i t h s t r a i n r a t e i s t o t a l l y a t v a r i a n c e with c u r r e n t l y h e l d views and furthermore i t occurs a t q u i t e moderate v e l o c i t i e s which are c e r t a i n l y w e l l w i t h i n t h e range o f those encountered i n many metal  forming o p e r a t i o n s .  X  FIG.4.1  DISTORTION OF GRID LINES DURING DEFORMATION FOR 0,02 FT./MIN. DEFORMATION SPEED,  o >-oJ  FIG.4.2  DISTORTION OF GRID LINES DURING DEFORMATION FOR 0.02 FT./MIN DEFORMATION SPEED.  c  \  '  ) \  <s' 7k  \^ 7  '  \S  \ 7^  \7 \  -7 <  \' 7 s.  \7 7\  X r ) \  >  V.  in  )  >  k  )V.  fV  ^ 7 \  '  \Y  f  v- ' \ / \• \ —7s. / \/  R /•  \ 7  S  f ^  )  7  f  / k i-  7  \ f 7  0 i.^  /_  (  >  s  >  ( /  7< \  i  \  3«C  )  {  /\  )  <=> N  \  7\  "sf  : — 7 .  ' < f  )^  '  ;  \  NK \{ 7k  /\.  7  \ /• 7 k—7  }  -7 \  \  w  /  S7  \  )  Ln  7\ 7  \ 7v  '  k  >  C '^  7  >7  / 7v  7  f  f  .\ f 7  /  ' ^  -s t /  L  > J  £ -A c  /  N/  >,  )c  >\  \* 7  VV  i  c  \ 7\  /• C— "  ^'  K—~7  >\ 7\  v  > J ^  f '\  7\  ; /  V r  )  7  \f 7 <  >  7\  \' 7**  /  2  7K  (  •>  '  ^ f  ^c. > 7  /^  V  >/  >  f  -)  in  a  —j~~  m  a  0.0  0.5  1.0  1.5  2.0  2.5  3.0  X  FIG4.3 DISTORTION OF GRID LINES DURING DEFOLIATION FOR 1 5 . 7 FT./SEC. DEFORMATION SPEED.  58  FIG.44 DISTORTION OF GRID LINES DURING DEFOMATION FOR 15.7 FT./SEC. DEFORMATION SPEED.  59  m  in o  0.0  0.5  1.0  1.5  2.0 X  2.5  FIG. 4.6 DISTORTION OF GRID LINES DURING DEFOMATION FOR 1 5 . 7 FT./SEC. DEFORMATION SPEED.  3.0  61  62  63  FIG.4.9 DISTORTION OF GRID LINES DURING DEFOLIATION FOR 15.7 FT./SEC. DEFORMATION SPEED.  64  68  prG4.14  THREE DIMENSIONAL PLOT OF HORIZONTAL VELOCITY (u) AS A FUNCTION OF X AND Y FOR J5.7 FT./SEC. DEFORMATION SPLED.  69  FIG.4.15  THREE  DIMENSIONAL PLOT OF VERTICAL VELOCITY (V) AS A FUNCTION  OF X AND Y FOR  15.7 FT./SEC. DEFORMATION SPEED.  FIG.4.16 THREE DIMENSIONAL PLOT OF EFFECTIVE STRAIN-RATE (£) AS FUNCTION OF X AND Y FOR 0 . 0 2 FT./MIN-DEFORMATION SPEED.  71  FIG.4.17 THREE DIMENSIONAL PLOT OF EFFECTIVE STRAIN-RATE (t) AS FUNCTION OF X AND Y FOR 15.7 FT./SEC. DEFORMATION SPEED.  FIG.4.18 THREE DIMENSIONAL PLOT OF TOTAL EFFECTIVE STRAIN (§) AS FUNCTION OF X AND Y FOR DEFORMATION SPEED OF 0.02 FT./MIN.  73  FIG.4.19 THREE DIMENSIONAL PLOT OF TOTAL EFFECTIVE STRAIN . ( . £ ) . AS FUNCTION OF X AND Y FOR 15,7 FT./SEC DEFORMATION SPEED.  FIG.4.20THREE DIMENSIONAL PLOT OF NORMAL STRESS (a ) AS FUNCTION OF X. AND Y FOR 0.02 FT./MIN DEFORMATION SPEED.  76  FIG.422 THREE DIMENSIONAL  PLOT OF SHEAR STRESS (T  ) AS FUNCTION  OF.X AND Y FOR 0.02 FT./MIN DEFORMATION SPEED.  80  81  82  FIG.4.29 TWO DIMENSIONAL PLOT  DEFORMATION SPEED.  OF 7^ AS A FUNCTION OF X FOR 1 5 . 7 FT./SEC. f  4.3.  84  SOURCES OF ERROR  The q u a l i t y  o f the g r i d node p o s i t i o n a l  to the program i s the most important s i n g l e results.  data used as i n p u t item a f f e c t i n g the  Sparse o r p o o r l y d i g i t i z e d data i s l i k e l y t o cause  inconsistencies  i n the output s t r e s s  distributions.  The  surface f i t t i n g routines w i l l  smooth c e r t a i n i r r e g u l a r i t i e s  but t h e r e i s a l i m i t t o t h e i r  capabilities.  The v e l o c i t y o f the upper p l a t e r  varies  d u r i n g the c y c l e ,  and d u r i n g specimen deformation, a c c o r d i n g t o the equations 3.4 and 3.9.  A further  fluctuation  of plateja v e l o c i t y may  occur as the drivewheel speed changes d u r i n g the c y c l e . energy balance c a l c u l a t i o n s  Initial  show this change i s l i k e l y t o be  very s m a l l p a r t i c u l a r l y with low s t r e n g t h p r o j e c t i l e s . p r i o r knowledge o f the a c t u a l  upper-platen v e l o c i t y  A  profile  d u r i n g deformation i s not r e q u i r e d as t h i s i s o b t a i n e d automatically  from the d i g i t i s e d displacement d a t a and a knowledge  of t h e time increment between frames o f the h i g h speed  photo-  graphs. P l a n e - s t r a i n deformation was achieved u s i n g a Kudo apparatus.  While t h i s assured p l a n e - s t r a i n  condition i t did  i n t r o d u c e a f r i c t i o n a l drag on the end faces o f the specimen. The e f f e c t was minimized  u s i n g s i l i c o n grease as a  lubricant  and from examination o f deformed specimen i t was concluded the e f f e c t was not important. In the a n a l y s i s  the m a t e r i a l was assumed t o be s t r a i n - r a t e  i n s e n s i t i v e which i s a common assumption  and not unreasonable  85  f o r many m e t a l - f o r m i n g m a t e r i a l s .  I t i s possible to r e l a t e  e f f e c t i v e s t r e s s , a , t o b o t h e f f e c t i v e s t r a i n e" and e f f e c t i v e *  s t r a i n r a t e e.  With m o d i f i c a t i o n s t o t h e a n a l y s i s s t r a i n - r a t e  s e n s i t i v e m a t e r i a l s c o u l d be accommodated.  86 5.  CONCLUSIONS  The v i s i o p l a s t i c i t y a p p r o a c h dynamic and q u a s i - s t a t i c ,  steady  have been developed f o r o r non-steady  deformation  processes. The tion using  effect  during different t h i s work.  upsetting, with  o f impact  velocity  metal-working  I t i s clear  from  t h a t s t r a i n and s t r e s s  strain-rate.  on t h e mechanism processes the i n i t i a l  o f deforma-  can be s t u d i e d study o f  d i s t r i b u t i o n vary  significantly  6.  87 SUGGESTIONS FOR FURTHER WORK  The method developed enables the s t r e s s d i s t r i b u t i o n s to be determined i n many dynamic metal-forming o p e r a t i o n s .  A  s t a r t i n g p o i n t f o r t h i s i s to determine the change o f s t r e s s d i s t r i b u t i o n with s t r a i n rate  (or impact v e l o c i t y ) ,  material  d e n s i t y and s u r f a c e geometry f o r p l a n e - s t r a i n u p s e t t i n g operations. M o d i f i c a t i o n s can be made to the s u r f a c e f i t t i n g  routines  to accommodate the c o n s t r a i n t s that the v e l o c i t y g r a d i e n t i s zero along  the y a x i s , and t h a t the v e r t i c a l component o f  v e l o c i t y , v, i s equal to the p l a t e n v e l o c i t y f o r p o i n t s on the platen-workpiece i n t e r f a c e .  It i s likely  t h a t a 5th order  polynomial would then be needed f o r s u r f a c e  fitting.  88  R E F E R E N C E S  89  References 1.  Hoffmann, 0 . , and Sachs, G . , " I n t r o d u c t i o n t o the Theory o f P l a s t i c i t y f o r E n g i n e e r s , " McGraw H i l l , New York, 1 9 5 3 .  2.  S i e b e l , E., "Die Formgebung i n b i l d s m e n Zustande", V e r l a n g S t a h l e i s s e n , Dusseldorf, 1932.  3.  Henky, H. , "Tiber e i n i g e s t a t i s c h bestimmte F a l l e des G l e i c h g e w i c h t s i n p l a s t i s c h e n K o r p e r n " , Z. Angew. Math. Mech. 3, 241-251, 1 9 2 3 .  4.  P r a n d t l , L., "Anwendungs b e i s p i e l e zu einem Henckyschen S a t z uber das G l e i c h g e w i c h t " , A. Angew. Math. Mech., _3_, 401-406,  5.  1923.  C a r a t h e o d o r y , C , and Schmidt, E., "Uber d i e HenckyP r a n d t l s c h e n Kurven", Z. Angew. Math. Mech., 468475,  1923.  6.  G e i r i n g e r , H., " B e i t r a n g zum v o l l s t a n d i n g e n ebenen P l a s t i z i t a t s - Problem," P r o c . T h i r d I n t e r n a t i o n a l Congress o f A p p l i e d Mechanics, 1 . , 185-190, 1 9 3 0 .  7.  H i l l , R., The M a t h e m a t i c a l Theory o f P l a s t i c i t y Clarendon Press, Oxford, 1950.  8.  P r a g e r , W. and Hodge, P.G. J r . , Theory o f P e r f e c t l y P l a s t i c S o l i d s . W i l e y , New Y o r k , 1 9 5 1 .  9.  L e e , E.H., "The T h e o r e t i c a l A n a l y s i s o f M e t a l Forming Problems i n P l a n e S t r a i n " , J . A p p l i e d Mechanics, 1 9 , 97-103,  1952.  10.  Johnson, W. and M i l l e r , P.B., P l a s t i c i t y f o r M e c h a n i c a l E n g i n e e r s , D. Van N u s t r a n d Co., L t d . , 1 9 6 2 .  11.  B i s h o p , J . F . , "On t h e E f f e c t o f F r i c t i o n on Compression and I n d e n t a t i o n Between F l a t D i e s " , J . Mech. Phys. S o l i d s , 6, 132-144,  1958.  12.  Green, A.P., "On Unsymmetrical E x t r u s i o n i n P l a n e S t r a i n " , J . Mech. Phys. S o l i d s , 3, 189-196, 1955-  13.  Oxley, P.L.B. and Farmer, L.E., " S l i p - L i n e F i e l d f o r Plane-Strain Extrusion of a Strain-Hardening M a t e r i a l " , J . Mech. Phys. S o l i d s , V 1 9 , N6 Nov. 1971, 3 6 9 - 3 8 8 .  14.  C h i t k a r a , N.F. and C o l l i n s , I . F . , "A G r a p h i c a l Technique for Constructing Anisotropic S l i p - L i n e F i e l d " , I n t . J . Mech. S c . , 1 9 7 4 , pp. 241-248  90  15.  R i c e , J.R.,-"Plane S t r a i n S l i p - L i n e Theory f o r A n i s o t r o p i c R i g i d / P l a s t i c M a t e r i a l s " , J . Mech. Phys. S o l i d s , V21 N2, March, 1973, 6 3 - 7 4 .  16.  Booker, J.R. and Davis, E . I . , "A General Treatment o f P l a s t i c A n i s o t r o p y under c o n d i t i o n s of Plane S t r a i n " , J . Mech. Phys. S o l i d s , V o l . 2 0 , 1972, pp. 239-250.  17.  Shabaik, A., " E f f e c t o f F r i c t i o n and Degree of Deformation on Buldge P r o f i l e During Compression", Proc. North American Metal-Working Conference, McMaster U n i v e r s i t y , Canada,  1973,  18.  PP. 221-238.  Wilson, W.R.D., " S l i p - L i n e S o l u t i o n s f o r S t r i p Drawing with A r b i t r a r y F r i c t i o n C o n d i t i o n s " , Proc. 5 t h . NAMRC,  SME,  1977, PP. 8 0 - 8 6 .  19.  Ewing, D.J.F., "A S e r i e s Method f o r C o n s t r u c t i n g P l a s t i c F i e l d s , " J . Mech. Phys. S o l i d s , V o l . 15, 196?, pp. 105-114.  20.  C o l l i n s , I.F., "The Algebraic-Geometry o f S l i p - L i n e F i e l d s w i t h A p p l i c a t i o n s to Boundary Value Problems", Proc. Roy.  S o c , Ser. A, V o l . 3 0 3 , 1968, pp. 317-338.  21.  Drucker, D . C , Prager, W. and Greenberg, H.J., "Extended L i m i t Design Theorems f o r Continuous Media", Quart. Appl. Math., 2, 1952, pp. 381-389.  22.  A v i t z e r , B., Metal Forming; Process and A n a l y s i s , H i l l , New York, 1968.  23.  Johnson, W., "Estimate o f Upper Bound Loads f o r E x t r u s i o n and C o i n i n g Operations", P r o c I n s t . Mech. Engrs. (London),  173,  McGraw  pp. 6 1 - 7 2 , 1957.  24.  Kudo, H. , "An Upper Bound Approach to Plane S t r a i n F o r g i n g and E x t r u s i o n - I " , I n t . J . Mech. Science, 229-252, i 9 6 0 .  25.  Kudo, H., "An Upper Bound Approach to Plane S t r a i n F o r g i n g and E x t r u s i o n - I I " , I n t . J . Mech. S c i . , 1:  229-252,  i960.  26.  Kudo, H., "An Upper Bound Approach to Plane S t r a i n and E x t r u s i o n - I I I " , I n t . J . Mech. S c i . , 1: 366-368, i 9 6 0 .  27.  Kobayashi, S., "Upper Bound S o l u t i o n of Axisymmetric Forming Problems - I " , Presented a t the P r o d u c t i o n E n g i n e e r i n g Conference o f ASMS, May, 1963.  28.  Kobayashi, S., "Upper Bound S o l u t i o n s of Axisymmetric Forming Problems I I " , Trans. ASME, S e r i e s B, J . Eng. Ind., 86 No. 4 , 196^.  91  29.  Nagpal, V., Lahoti, G.D. , Altan, T., "A Numerical Method for Simultaneous Production of Metal Flow and Temperatures i n Upset Forging of Rings", ASME Paper No. 77-WA/Prod.-35.  30.  Lahoti, G.D. and Altan, T., "Prediction of Temperature Distributions i n Tube Extrusion using a Velocity F i e l d without Discontinuities", Proc. 2nd. NAMRC, University of Wisconsin, Madison, 1974, pp. 207-224.  31.  Lahoti, G.D. and Altan, T., "Computer-Aided Analysis of Metal Flow and Temperatures i n Radial Forging of Tubes", Proc. Int. Prod. Eng. Res. Conference, New Delhi, 1977, PP. 323-339.  32.  Vickers, G.W., Plumtree, A., Sowerby, R. and Duncan, J.L., "Simulation of the Heading Process", ASME Paper No. 74-WA/ Prod-19.  33.  C o l l i n s , I.F., "The Upper Bound Theorem f o r R i g i d / P l a s t i c Solids Generalized to Include Coloumb F r i c t i o n " , J . Mech. Phys. Solids, Vol. 17, 1969, pp. 323-338.  34.  Sauerwine, F. and Avitzer, B. , "Limit Analysis of Hollow Disk Forgingy, Part I: Upper Bound", ASME Paper No. 77WA/Prod-2, 1977.  35.  Sauerwine, F. and Avitzer, B., "Limit Analysis of Hollow Disk Forging, Part 2: Lower Bound", ASME Paper No. 77-WA/ Prod-3, 1977.  36.  Dwivedi, S.N., Sharan, R. and Mishra, C.B., " P l a s t i c Deformation of Polygonal Disks under High Velocity Impact", Presented i n 8 t h Conference of U.S. National Congress of Applied Mechanics, June, 1978.  37.  Juneja, B.L. and Prakash, R., "An Analysis f o r Drawing and Extrusion of Polygonal Sections", Int. J . Mach. Tool Des. Res., Vol. 15, 1975, PP- 1-30.  38.  Nagpal, N., "Analysis of Plane-Strain Extrusion through A r b i t r a r i l y Shaped Dies using Flow Function", J . Engg. Ind.,  Vol. 99, 1977, PP. 544-548.  39.  Nagpal, N., "General Kinematically Admissible V e l o c i t y F i e l d f o r some Axisymmetric Metal Forming Problems, " J . Engg. Ind., Vol. 96, 1974.  40.  Zienkiewicz, O.C., The F i n i t e Element Method, Jrd Edition, McGraw-Hill, 1977.  92  41.  Lee, C H . and Kobayashi, S., "New S o l u t i o n s to R i g i d - P l a s t i c Deformation Problems u s i n g a M a t r i x Method," Trans. ASME, J . of Engrg. f o r Ind. V o l . 95, 1973. pp. 865-873.  42.  Godbole, P.H. and Z i e n k i e w i c z , O.G., "A P e n a l t y F u n c t i o n Approach to Problems of P l a s t i c Flow of M e t a l s w i t h Large Surface Deformation", J . S t r a i n A n a l y s i s , V o l . 10, 1975. pp. I 8 O - I 8 3 .  43.  Huebner, K.H., The F i n i t e Element Method f o r E n g i n e e r s , John Wiley & Sons, 1975.  44.  Strang, G. and F i x , G.J., An A n a l y s i s of the F i n i t e Element Method, P r e n t i c e - H a l l , Englewood C l i f f s , N.J., 1973.  45.  Rowe, G.W., "Recent Developments i n the Theory and P r a c t i c e of Metal Forming", Proc. T h i r d North American Metal-Working Research Conference, Carnegie-Melton U n i v e r s i t y , 1975. pp. 2-25.  46.  Alexander, J.M. Metal Forming", Research, 1977.  47.  Lee, E.H., M a l l e t , R.L. and Yang, W.H., " S t r e s s and Deformation A n a l y s i s o f the Metal E x t r u s i o n P r o c e s s " , Computer Methods i n A p p l i e d Mechanics and Engg., V o l . 10,  1977.  PP.  and P r i c e , " F i n i t e Element A n a l y s i s of Hot Proc. 18th I n t . Conf. Mach. T o o l Design pp. 267-274.  339-353.  48.  Lee, E. H., M a l l e t , R.L., and McMeeking, R.M., " S t r e s s and Deformation A n a l y s i s o f Metal Forming Processes, "Numerical M o d e l l i n g of Manufacturing Processes", ASME S p e c i a l P u b l i c a t i o n PVP-PB-025, 1977, pp. 19-34.  49.  O d e l l , E . I . , "A Study of W a l l I r o n i n g by the F i n i t e Technique", ASME Paper No. 77-WA/Prod-8.  50.  W i f i , A.S., "An Incremental Complete S o l u t i o n of the S t r e t c h Forming and Deep-Drawing of C i r c u l a r Blank u s i n g a H e m i s p h e r i c a l Punch", I n t . J . Mech. S c i e n c e , V o l . 18, 1976,  Element  pp. 23-31.  51.  Wang, N.M. and Budiansky, B., " A n a l y s i s of Sheet Metal Stamping by F i n i t e Element Method", General Motors Research P u b l i c a t i o n GMR-2423, Sept. 1977.  52.  Shah, S.N. and Kobayashi, S., "A Theory of Metal Flow i n Axisymmetric P i e r c i n g and E x t r u s i o n " , J . Rud. Engrg., V o l .  1977.  53.  PP. 73-103.  Kobayashi, S., " R i g i d - P l a s t i c F i n i t e Element A n a l y s i s of Axisymmetric Metal Forming Processes, Numerical M o d e l l i n g of Manufacturing P r o c e s s e s " , ASME S p e c i a l P u b l i c a t i o n PVP-PB-025, 1977, PP. 4 9 - 6 5 .  1,  93  54.  Kobayashi, S. and Matsumoto, H., "A note on the M a t r i x Method f o r R i g i d - P l a s t i c A n a l y s i s of Ring Compression", Proc. 18th MT DR Conference, London, 1977, pp. 3 - 9 .  55.  Kobayashi, S. and Chen, C.H., "Deformation A n a l y s i s of M u l t i - P a s s Bar Drawing and E x t r u s i o n " , CIRP Annalus, 1978.  £6.  P r i c e , J.W.H. and Alexander, J.M., "A study of Isothermal Forming or Creep Forming of a T i t a n i u m A l l o y " , 4 t h North American Metal-Working Research Conference, B a t t e l l e , Columbus, Ohio, 1976, pp. 46-57.  57.  Lung, M. and Mahrenholtz, "A F i n i t e Element Procedure f o r A n a l y s i s of Metal Forming P r o c e s s e s " , Trans, of CSME,  V o l . 2, No. 1,  1973-74.  58.  W i l k i n s , M.L., " C a l c u l a t i o n of E l a s t i c - P l a s t i c Flow", Lawrence R a d i a t i o n L a b o r a t o r y Report, UCRL-7322, Rev. Livermore, U n i v e r s i t y of C a l i f o r n i a , 1969.  59.  Gordon, P. and Karpp, R., " A p p l i c a t i o n of a New F i n i t e D i f f e r e n c e Method of M e t a l Forming", Proc. F o u r t h NAMRC, B a t t e l l e , Columbus Lab., 1976, pp. 72-79.  60.  Shabaik, A.H., "Computer S i m u l a t i o n of Metal Flow During E x t r u s i o n " , Proc. 1976 I n t . Conference on Computer S i m u l a t i o n f o r M a t e r i a l s A p p l i c a t i o n s , Nuclear M e t a l l u r g y , V o l . 20, P a r t 2, 1976, pp. 752-765-  61.  Woo, D.M., "On the Complete S o l u t i o n of the Deep Drawing Problem", I n t . J . Mech. S c i . , 10, 1968, pp. 8 3 - 9 4 .  62.  Wang, N.M. and Shammamy, M.R., "On the P l a s t i c B u l g i n g of a C i r c u l a r Diaphragm by H y d r o s t a t i c P r e s s u r e " , J . Mech. Phys. S o l i d s 17, 1969, pp. 4 3 - 6 1 .  63.  Wang, N.M., "Large P l a s t i c Deformation of a C i r c u l a r Sheet Caused by Punch S t r e t c h i n g " , J . Appl. Mech. 1970, pp. 431-440.  64.  Yamada, Y. and Yokochi, G., " E l a s t i c - P l a s t i c A n a l y s i s of the H y d r a u l i c Bulge T e s t by the Membrane Theory", Manf. Res. 21, 1969.  65.  Thomsen, E.G. and L a p s l l e y , J.T., "Experimental S t r e s s Determination w i t h i n a metal d u r i n g P l a s t i c Flow", Proc. Soc. E x p t l . S t r e s s A n a l y s i s . , 11, No. 2, 5 9 - 6 8 , 1954.  66.  Thomsen, E.G., Yang, C.T., and Bierbower, J.B., "An Experimental I n v e s t i g a t i o n of the Mechanics of P l a s t i c Deformation of M e t a l s " . U n i v e r s i t y of C a l i f o r n i a P r e s s , B e r k e l e y and Los Angeles, 1954, pp. 89-144.  67.  Thomsen, E.G., " V i s i o p l a s t i c i t y " , CIRP Conference, September, 1963.  1,  94  68.  Shabaik, A. A l t a n , T. and Thomsen, E.G., " V i s i o p l a s t i c i t y " , F i n a l Report Prepared f o r the U.S. Navy, Bureau of Naval Weapons, February, 1965.  69.  Shabaik, A. and Kobayashi, S., " I n v e s t i g a t i o n of the A p p l i c a t i o n of V i s i o p l a s t i c i t y Methods o f A n a l y s i s t o Metal Deformation Processing",- F i n a l Report prepared f o r the U.S. Navy, Bureau of Naval Weapons, Feb. 1966.  70.  Shabaik, A. and Thomsen, E.G., " I n v e s t i g a t i o n o f the A p p l i c a t i o n o f V i s i o p l a s t i c i t y Methods o f A n a l y s i s to Metal Deformation P r o c e s s i n g " , F i n a l Report prepared f o r the U. S.Navy, Bureau o f Naval Weapons, Feb. 1967.  71.  Shabaik, A. and Kobayashi, S., "Computer A p p l i c a t i o n to the V i s i o p l a s t i c i t y Method", J o u r n a l of E n g i n e e r i n g f o r I n d u s t r y , Trans. ASME, S e r i e s B., V o l . 8 9 , No. 2, May, 1967,  PP. 339-346.  72.  Lee, C.H., and S h i r o Kubayashi, "Matrix Method of A n a l y s i s f o r P l a s t i c Deformation Mechanics and i t s A p p l i c a t i o n to V i s i o p l a s t i c i t y " , Ann. CIRP, Volume 21, No. 1, pp. 71-72, 1972.  73-  Medrano, R.E. and G i l l i s , P.P., " V i s i o p l a s t i c i t y Techniques i n Axisymmetric E x t r u s i o n " , J o u r n a l of S t r a i n A n a l y s i s ,  V o l . 7, No. 3, 1972,'pp. 170-177.  74.  Shabaik, A.H., "Computer A i d e d V i s i o p l a s t i c i t y S o l u t i o n to Axisymmetric E x t r u s i o n Through Curved Boundaries", J o u r n a l of E n g i n e e r i n g f o r I n d u s t r y , Nov., 1972, pp. 1225-1231.  75-  Medrano, R.E., G i l l i e s , P.P., Conrad, H. and H i n e s l e y , C P . , "Use of M i c r o s t r u c t u r e f o r V i s i o p l a s t i c i t y A n a l y s i s " , J o u r n a l S t r a i n A n a l y s i s , Volume 9, Number 3, pp. 146-151, J u l y , 1974.  76.  Brown, J.H. and P.F. Thomson, "Mechanics of Deformation During C o l d R o l l i n g " , A u s t r a l i a ' s Conference on the Mechanics of S t r u c t u r e s and M a t e r i a l s , pp. 415-422, 1977-  77«  Robinson, J.N. and Shabaik, A.H., "Determination o f the R e l a t i o n s h i p between S t r a i n and Microhardness by means of V i s i o p l a s t i c i t y " , Machine Design, Volume 4, Number 9, pp. 2091-2095, Sept. 1973.  78.  Mohamed, S.A. and Tetelman, A.S., " A p p l i c a t i o n of the V i s i o p l a s t i c i t y Technique f o r D e r i v a t i o n o f the C r i t e r i o n f o r D u c t i l e Rapture I n i t i a t i o n i n F u l l y P l a s t i c Notched Bars", E n g i n e e r i n g F r a c t u r e Mechanics, V o l . 7, No. 4 ,  pp. 63l-'640, 1975.  79.  Green,A.P.,"The Use o f P l a s t i c i n e Models to Simlate the P l a s t i c Flow of M e t a l s , " P h i l . Mag. Ser. 7.,Vol. 42,Page 365-373,] 95? .  95  A P PEN D I X  COMPUTER PROGRAM FOR VISIOPLASTICITY 1 2 3 4  IMPLICIT SEAL*8 (A-H,0-.Z) COMMON XX (15) ,XI,YI,CC,CN,CA,Y(40G,2) ,EBR (400,20) COMMON/C1/EHO,YOID{15,2) ,DT DIMENSION XYT (2 , 20,20,2) ,0(20,20) ,V (20,20) ,XY(2,80,30 )  5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 46 47 48 49 50 51 52 53 54 55 56  DIMENSION X (15, 40 0) , ERBOR ( 2) , I SUB (8 0 ,8 0) REAL*4 Z1 (90,90) ,Z2 (90,90) ,TAUXY (80 ,80) LOGICAL REFINE DIMENSION EBBDT(80,80) ,EXDT (80,80) , EYDT (80,80) 1,GAMDT (80,80) , LANDA (80,80) BEAL*8 DEBB(80,80) BEAL*8 2LS(10,800) REAL*8 LANDA HEAL*8 X1 (20,20 ,2) ,Y1 (20,20, 2) DIMENSION AINT2(80) C C C C C C C Y C C OOBD C C 10 20  IX=NG PTS IN X IY=NO PTS IN Y IT=NO OF TIME STEPS DT=TIME INTERVAL BETWEEN TIME STEPS CC,CN ARE CONSTANTS WHERE SIGB=CC*EB**CN CA IS LOWER INTERVAL OF INTEGRATION FOR SIG SIGYOA IS CONSTANT ADDED TO SIGY XYT(L,I,J,K) CONTAINS: L= 1 X-COGED; L=2 Y-C FOB 1=1,11 J=1,IX L=1,2 READ(5,10) IX,IY,IT FORMAT (312) READ (5, 20) DT FORMAT(3F10.0) READ(5,20) CC,CN , CA, SIGYOA  C C NIX CONTAINS # GRID PTS IN X DIRECTION, NIY IN Y DIRECTION FOR PLOTS C READ (5,10) NIX,NIY BEAD (5, 20) BHO C C X S Y COORD READ FOB TIME=0 C BEAD(4,30) ( ( (XYT (K, I, J , 1^ ,K= 1 ^2) ,1 = 1 ,1Y) , J= 1 ,IX) DO 25 11=1,IX DO 2 5 IJ=1,IY IF ( I I . EQ. 1). XYT ( 1 , I J , I I , 1)=0.D0 XYT{2,IJ,II,1)=XYT(2,IJ,II,1)-1.0D0 25 CONTINUE 30 FORM AT(5X,2F6.3,1X, 2F 6. 3, 1 X, 2F 6. 3 ,1X,2F6.3,1X,2F6.3) LTM=2 NTM= 1 IDIM=20 IDIM.P=80 IXY=.IX*IY CALL AX.IS(0.,0., X»,-3,lO-,0.,0.,.2) CALL AXIS (0. ,0. , *Y', 1 ,10. ,90. ,0. ,.2) CALL PLOT (XYT (1 ,IY„ 1,1) *5» , XYT (2,IY, 1, 1) *5, ,3) CALLPLOT (XYT (1,IY,IX, 1) *5.,XYTf2,IY,IX,1) *5. ,2) J  COMPUTER PROGRAM FOR 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 81.6 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 1 10 111 112 113 114 115 116  35 C C  40  VISIOPLASTICITY  CALL PLOT (XYT (1 , 1,IX, 1) *5, ,XYT (2, 1, IX, 1) *5- ,2) DO 35 1=2,IX,2 DO 35 J= .2,11,3 X1(I,J,1)= XYT(1,J,I, 1) Y1 (I,J,1)=XYTf2*J rI#D i  INITIALIZE EBR TO ZERO DO 40 .1=1,20 DO 40 3=1,400 EBR (J,I) =0. DO CONTINUE FACT= 1- DO  C  c  FOR EACH TIME STEP EBB IS ACCUMULATED IT1=IT-1 DO 80 K= 1,IT1  c c c  45  X&Y COORD ARE READ FOR  NTM= 3-NTM 1TH=3-1TM READ {4, 30) {{ (XYT (KK,I, J,NTM) ,KK=1,2) ,1=1 ,IY) , J=1 DO 45 1=1,IX DO 45 J=1,IY IF (I.EQ. 1) XYT (1 ,J,.I,NTM) =0. DO XYT(2,J,I,NTH)=XYTi2,J,I,NTH)-1-0D0 CONTINUE .  c c c c  48  50 C C C C C  NEXT TIME STEP  U,V CALCULATED FOR THIS TIKE STEP U MOST BE >0, V MUST BE <0 DO 50 J=1,IX DO 50 1=1,IT U (I,J)=- (XYT(1,I,J,N.TM) -XYT { 1 ,1, J,1TM) ) /DT y {I , J)=- {XYT(2,I,J,NTM) -XYT (2,I,J,1TM^ )/DT IF (U (I, J) . GT. 0. DO) UtI,J)=0.D0 IF {V (I, J) . LT.O-DO) V(I,J)=0.D0 CONTINUE CURVE FITTING FOR U AN D V USING DLSQHS SET UP INDEPENDENT VARIABLES IN X DEPENDENT VARIABLES IN Y DO 60 3=1,IX 11=1Y* (J-1) DO 60 1=1,IY 1=1+11  CALL AUX (XYT (1,1,3, NTM) , XYT (2,1, J,N'TM) ,X(1 ,L) ) XLS{1 ,1) =-X{2,L) XLS (2,L)=-2.D0*X (10,L) XLS (3,1) =-.3-P0*X(11,L) XLS (4,1) =-4.D0*X{12,l) XLS (5,1) =-. 5D0*X (3,1) XLS (6,1) =-X (13,1) XLS (7,1) =-1. 5D0*X (14,1) X1S(8,1) =-1.D0/3.D0**X (4,1) X1S(9,1) =-2.D0/3.D0*X (15,1) X1S(10,1)=-. 25D0*X (5,1)  COMPUTER PROGRAM FOR VISIOPLASTICITY 1 17 118 119 120 121  DO 55 IK=1, 10 XLS {IK,1XY+ L) =X (I K+5 , L) Y(L,1)=U(I,J) Y (I XY+L , 1) = V (I, J) CALL DLSQHS (Y,XLS, 2*IXY, 10 , 1 , 803 , 10 ,ERROR  55 60  98  FALSE. ,IER  ,&200) 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149. 150 151 152 153 154 155 156 157 158 159  62 63  64  C C C  DO 62 IK=1,5 Y (IK,2} =0.00 DO 63 IK = 6,15 Y (IK,2) = Y(IK-5,1) Y (1, 1)=0. DO Y (2 ,1)=-Y (6 ,2) Y (3, 1) =-. 5DG*Y{ 10,2) Y(4, 1)=-1.DO/3.D0*Y (13,2) Y (5, 1)=-. 25D0*Y {15,2) DO 64 IK=6,9 Y (IK, 1) =0.D0 Y(10,1)=-2.D0*Y(7,2) Y (11,1) =-3. DG*Y{8,2) Y(12, 1)=-4. D0*Y{9,2) Y (13, 1) =-Y (11 ,2) Y(14,1|=-1.5D0*Y(12,2) Y(15, 1) =-2-D0/3.D0*Y (14,2) U&V COEFF SAVED FOR DO/DT,DV/DT IF (K. NE. (IT1-1) ) GO TO 410 DO 400 1=1, 15 DO 400 J=1,2 YOLD (I, J) = Y ( I , J) CONTINUE  400 410 C C THE VALUES OF EDTX,EDTY,GAMXY,EBRDT,AND EBR ARE CALCULATED C AT EACH TIME STEP C DO 52 1=2,IX,2 DO 52 J=2,IY,3 X1 (I,J,NTM) = XYT (1 , J,I,LTM) -DT * AU X2 {XYT *( 1 , J , I , NTM) , 1 XYT (2, J , I , NTM) , Y .{1, 1) ,- TRUE. ) Y1 ( I , J , NTM) =X YT (2 , J, I , LTM) - DT*AUX2(XYT{1,J,1,NTM), 1 XYT (2,J, I,NTM) ,Y(1,2) , - FALSE. ) CALL PL0T(X1 (I,J,LTM) *5„ , Y 1 ( I , J, LTM) *5. ,3) CALL PLOT (X1 (I, J , NTM) *5. , Y 1 ( I , J, NTM) *5. ,2) CALL SYMBOL (X1 (I, J , NTM) *5. ,Y1 (I, J,NTM) *5. ,. 14,5,0. ,- 1 )  160 161 162 163 164 165 166 167 168 169 170  52  CONTINUE I F (K.NE.IT1) GO TO 53 CALL PLOT (0. ,XYT (2,IY, 1 , NTM) *5. , 3) CALL PLOT (XYT (1, IY,IX, NTM) *5.,XYT (2 ,1Y,IX, NTM) *5. ,2) DO 49 1=2,IY CALL PLOT (XYT (1 , IY-I +1,IX, NTM) *5-, XYT{ 2, IY-I+ 1 , IX ,NTH )*5.,2) 49 CONTINUE CALL PLOT (12. ,0. ,-3) 53 CONTINUE DO 70 J=1,IX IL=IY*(J-1)  COMPUTER PROGRAM FOR VISIOPLASTICITY 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202  .  DO 70 1=1,IY IJ=I+IL CALL AUX {XYT {1, I, J , NTM) , XYT (2,1, J , NTM) , XX) CALL DEBIV(Y(1,1) ,XYT (1,1,3,NTM) ,XYT (2 ,1, 3,NTM) ,DUDX, DUDY,3) CALL DERIV (Y (1,2) ,XYT (1,1, J,NTM) ,XYTi(2 ,I,3,NTM) ,DVDX, DVDY, 1) GAMXY=DUDY+ DVDX DPACT= (DUDX**2)/3-D0+(GAMXY) **2/12. DO 65 EBR DT (I, J) =2- D0*DSQRT (DFACT) IF (K.EQ.IT1) FACT=.5D0 EBR (13, 1) =EBR ( I J , 1) +FACT*DT*EBRDT (I , J) 70 CONTINUE 80 CONTINUE C C 4TH DEGREE POLY FIT TO EBR C CALL D1SQHS-{EBB, X,IX*IY, 15,2,400, 15, ERROR, .FALSE. ,IER ,S200) C C MASTER GRID IS SET FOR FINAL PLOTS C C XY (2,1,3) CONTAINS THE MASTER GRID ST XY(1, 1,3) IS C THE X-COORD, XY (2,1,3) IS THE Y-COOSD FOR 1=1,IX 3=1,IY XMAX = 0-D0 DO 90 1=1,IY IF (XYT ( 1,1, IX, NTM) - LT. XMAX) GO TO 90 XMAX=XYT(1,1,IX,NTM) IMY=I 90 CONTINUE YMAX=XYT(2,IY,IX,NTM) 100 CONTINUE DY=YMAX/(NIY-1) DO 110 1=1,NIX XY(2,I, 1)=0.D0 DO  XY (2,1,3) =XY (2,1,3-1) +DY CONTINUE DX=XMAX/(NIX-1) DO 120 1=1,NIY XY(1,1,I)=0.D0 DO 120 3=2,NIX XY(1,3,I)=XY (1,3-1,1) +DX CONTINUE  216 217 218 219 220 221 222  110  120 C C C  110  3=2,HIY  203  204 205 206 207 208 209 210 211 212 213 214 215  TO FIND ZERO FILLS ON PLOTS CALL FILL<XYT,IX,IY,ISUB,XY,NIX,NIY,IMY,NTa,IDIM,IDIM  C C C  TO PLOT U, V,EBSDT, AND TAUXY DO 125 3=1,NIY DO 125 1=1,NIX EBR DT {I, 3) =0. DO DEB1 (I, 3) =0 . DO  COMPUTES PROGRAM FOR VISIOPLASTICITY 223 224 225 226 227 228 229 230 231 232 233 2 34 235  125  126  TAUXY (I, J) = 0.1)0 100 NIY10=NIY+1Q NIX10=NIX+10 DO 126 J=1,NIY10 DO 126 1=1,NIX 10 Z1 ( I , J) =0.D0 Z2{I, J)=0.D0 DO 130 J=1,NIY DO 135 1=1,NIX IF (I SD3 {.I, J) . EQ. 0) GOTO 130 Z 1 (1*5,3 + 5) =-AUX2 (XY (1,I,J) , XY{2,I, J) ,Y {1 ,1) ,.TRUE. ) Z2(I + 5,J+5) =AUX2 (XY (1 ,I,J) ,XY{2,I,3) ,Y (1,2) ,.FALSE.) CALL DERIV (Y, XY (1,1, J) ,XY (2,1, J) , EX DT (I, J) ,G AMDT (I,J)  ,3) 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273  CALL DERIV {Y (1,2) , XY (1,1, J) , XY (2,1, J) , EBRDT {I, J) , EYDT (I,J) ,3) GAMDT (I, J) = G AMDT (I, J) +EBSDT (I, J) DFACT=(3.DO*EXDT(I,J)** 2+.75DQ*G AMDT(I,J) **2) 131 EBR DT (I,J)=2-DO/3.DO *DSQRT(DFACT) DEBE (I,J)=AUX2{XY{1,1,J) ,XY (2,1,J),EBR,•FALSE.) LAN DA (I, J) = 1.5D0*EBRDT(I,J)/(CC*DEBR (I, J) **CN) TAUXY (I J) = GAMDT (I, J) /{2- DO*.L ANDA {I , J) ) 135 CONTINUE CONTINUE 130 C C PLOT U,V,EBRDT,TAUXY C DXY = XY(2,NIX,NIY)/XY (1,NIX,NIY) CALL PERS{Z1,IOIflP+10,NIX+10,NIY+10,DXY,.333,45. ,45. , 10., 10.) CALL PLOT (12. ,0. ,-3) CALL PERS(Z2,IDIMP+10,NIX+10,NIY+10,DXY,.333,45. , 45. , 10., 10.) CALL PLOT {12. ,0. ,-3) DO 137 J=1,NIY10 DO 137 1=1,NIX10 Z1 (I,J) =0. 137 Z2 (I,J) =0. DO 138 J=1,N.IY DO 138 1=1,NIX 138 Z1 (1 + 5, J + 5) =DEBR (I, J) CALL PEP,S{Z1,IDIMP+10,NIX+10,NIY + 10 ,DXY,. 333,45. ,45., 10., 10.) CALL PLOT(12. ,0. ,-3) DO 139 J=1,NIY10 DO 139 I=1,NIX10 139 Z 1 {I, J) =0. . DO 140 3=1,NIY DO 140 1=1,NIX 21 {1 + 5, J + 5) =EB1DT{I,J) 140 Z2 (1 + 5, J+5) =-TAUXY{I, J) CALL PERS{Z1,IDIMP+10,NIX+10,NIY + 10,DXY,. 333,45.,45. , 10., 10. ) CALL PLOT (12. ,0.,-3) CALL PERS {Z2,ID.IM? + 10,NIX+1G,NIY + 10 ,DXY,. 333,45. ,45. , 10., 10.) CALL PLOT {12. ,0. ,-3) C ;r  COMPUTER PROGRAM FOR VISIOPLASTICITY 274 275 276 277 278 279 280 281 282 2 83 284 285 286 2 87 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 3 13 314 315 316 317 318 319 3 20 321 322 323 324 325 326 327  C C  CALC SIGY  101  EXTERNAL DF1,DF2 DO 141 J=1,N.IY10 DO 141 1=1,NIX 10 Z1 (I,J)=0. 141 Z2 {I, J) =0. DGRID=(XY(2, 1,2) -XY (2,1, 1) ) /2 JGRID=1 1000 IF (CA,. GT. XY (2 , 1 , JGRID) + DGRID) GO TO 1010 GO TO 1020 1010 JGRID=JGRID*1 IF (JGRID. EQ. NIY) GO TO 10 20 GO TO 1000 1020 AINT2(1) =0. DO DO 1030 1=2, NIX AINT2 (I) =AINT2 (1-1) IF (ISOB (I, JGRID) . EQ. 0) GO TO 10 30 AINT2 (I) =AI13T2 (I) +DQUANK (DF2, XY( 1 ,1-1,JGRID) 1 XY{ 1 ,1, JGRID) 00 1 DO ,TOL,FIFTH) 1030 CONTINUE DO 150 J=1,NIY DO 150 1=1,NIX IF (ISUB(I,J)-EQ.0) GO TO 150 XI= XY (1 , I , J) YI=XY (2,1, J) AINT1=0.D0 IF {JGRID .EQ. J) GO TO 1060 JADD=0 IF (ISUB ( I , JGRID) . EQ. 1) GO TO 1050 JINC=1 I F (JGRID. GT.J) JINC=-1 1040 JADD=JADD+J.INC IF (JGRID+JADD. EQ. J) GO TO 1060 IF (ISUB ( I , JGRID+JADD). EQ. 1) GO TO 1050 GO TO 1040 1050 AINT1=DQUANK (DF1 ,XY (2 ,1, JGRID+JADD) ,XY (2,1, J) 001DO, TOL, FIFTH) 1060 Z1 (1*5, JS-5) =SIGYOA-AINT1-AINT2 (I) Z2(I + 5, J + 5) =Z1 (1+5, J+5) + (EXDT (I, J ) - EYDT (I,J) ) /LANDA (I ,J) 150 CONTINUE CALL PERS(Z1,IDIMP+10,NIX*10,NIY+10,DXY,.333,45 .,45. , 10., 10.) CALL PLOT (12. ,0. ,-3) CALL PERS (Z2,IDIMP+10,NIX+10,NIY +10,DXY,- 333,45 -,45., 10., 10-) CALL PLOTND STOP 200 STOP 1 END SOBROUTINE AUX(X,Y,XX) IMPLICIT REAL*8 (A-H,0-Z) DIMENSION XX (15) XX{ 1) =1- DO XX (2)=X XX(3) =X*X XX{4)=X*XX(3)  COMPUTES PROGRAM FOR VISIOPLASTICITY 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 3 66 367 368 369 370 371 372 373 374 375 376 377 378 379 3 80 380-5 382 3 83 384  XX(5)=X*XX(4) XX(6)=Y XX (7) =Y*Y XX(8)=Y*XX(7) XX{9) =Y*XX (8) XX( 10) = X*Y XX{11)=XX (10) *Y XX{ 12}=XX (1 1) *Y XX (13) = XX (10) *X XX (11) =XX(13) *Y XX {15) =XX (1 3) *X RETURN END FUNCTION AUX2 (XX,YY,P,LL) IMPLICIT REAL*8(A-H,Q-Z) C C C C  EVALUATE THE FITTED FUNCTION AT XX,YY P CONTAINS THE FITTED PARAMETERS LL IS TRUE IF THE INDEPENDENT VARIABLE MUST BE EVALUATED BY AUX  C  10  C C  LOGICAL LL COMMON X DIMENSION X (15) , P (1) IF (LL) CALL AUX (XX,YY,X (1) ) AUX2=0.D0 DO 10 1=1,15 AUX2=AUX2+P (I) *X<I) CONTINUE RETURN END FUNCTION DF 1 (YC) IMPLICIT REAL*8(A—H,0-Z) COMMON XX(15) ,XI,YI,C1,C2,CA,Y (4 00, 2} , PEBR (20, 20) COMflON/C1/RHO,YOLD(15,2),DT THE INTEGRAND DTAU/DX IS EVALUATED  c  EB1=AUX2 (XT ,YC,PEBR,.TRUE. ) CALL DERIV (Y{1, 1) , XI, YC,EXDT, DUDY,3) CALL DERIV(Y(1,2) , XI , YC, DV DX, DVDY,3 ) GAM DT=DU DY + DVDX EBRDT=2-D0/3-D0*DSQRT (3. DO *EXDT**2+. 75D0*GAMDT**2) CALL DERIV (PEB1;XI,YC,DEBRDX,DUM ,1) CALL DERIV2 (Y (1, 1) ,XI, YC,DUDXY,3) CALL DER.IV2 (Y (1,2) , XI,YC,DVDXX,1) DGAM DX= DUDX Y+DVDXX CALL DERIV2 (Y (1, 1) , XI, YC , D EXDX , 1 ) DEBRDT=(4.DO*EXDT*DEXDX+GA MDT*DGAMDX)/(3.D0*EBRDT) BF1=C1*EBR**C2/(3-D0*EBRDT)*(C2*DEBRDX*GA MDT/EBR 1 +DGAMDX - DEBSDT*GAMDT/EBRDT) U=AUX2(XI,YC,Y(1, 1) ,. FALSE.) V=aUX2(XI,YC,Y(1,2) FALSE.) VOLD=AUX2 (XI,YC,YOLD (1,2) ,-FALSE.) FACT=RHO* {(V-VOLD)/DT) DF1=DF1+FACT RETURN END  COMPUTER  PROGRAM FOR VISIOPLASTICITY 103  3 85 386 387 388 3 89 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413  FUNCTION DF2{X) IMPLICIT REAL*8(A-H,0-Z) COMMON XX{15) ,XI,YI,C1,C2,CA,Y (400, 2) , PEBR (20, 20) COM MON/C1/R HO,YO.L.D(15,2) ,DT REAL*8 LAMBDA C c c  THE INTEGBAND DTAU/DY IS EVALUATED  431  EBR=AUX2{X,CA,PEBB,.TRUE.) CALL DERIV (Y (1, 1) , X, CA ,EXDT, DO DY , 3) CALL DERIV {Y{1, 2) , X,CA,DVDX,EYDT,3) GAMDT=DUDY+DVDX EBR DT=2.DO/3.D0*DSQBT<3.D0*EXDT**2+.75D0*GAMDT**2) CALL DEHIV(PEBR,X,CA,DUM,DEBRDY,2) CALL DEB.IV2 (Y (1, 1) ,X,CA,DUDYY,2) CALL DEBIV2 (Y (1,2) , X ,CA, DVDXY ,3) DGAMDY=DUDY Y+DVDXY CALL DEBIV2 (Y (1, 1) ,X,CA,DEXDY,3) YEBRDT={4-DO*EXDT*DEXDY+GAMDT*DGAMDY)/{3.D0*EBBDT) CALL DER IV 2 {Y (1, 1) ,X,CA,DEXDX,1) CALL DEF.IV2 {Y {1, 2) ,X,CA,DEYDX,3) LAMBDA= (2. D0*C1*EBB**C2) / ( 3 . DO*EBRDT) CALL DERIV(? EBB,X,CA,DEB RDX> DUM,1) CALL DEBIV2{Y (1,2) ,X,CA,DVDXX, 1) DGAMDX=DEXDY+DVDXX DEBR DT=(4. DO*EXDT*DEXDX+GAMDT*DGAMDX) /(3.D0*EBBDT) DF2=LAM BDA*({(EYDT-EXDT)*{-C2*DEBBDX/EBB 1 *DEB8DT/EBBDT) +DEXDX-DEYD X) 1 +.5D0*(C2*DEBBDY*GAMDT/EBB + DGAMDY - YEBBPT*GAMDT/E BBDT) ) U=AUX2{X,CA,Y(1,1) ,. FALSE. ) V= A U X2 f X , C A , Y { 1 , 2) ,,. FALSE. ) UOLD=AUX2 (X,€A,YOLD (1 ,1) ,. FALSE. ) FACT=BHO* { (U-UOLD) /DT) DF2=DF2 + FACT BETURN END SUBROUTINE DERIV{A,X,Y,DUDX,DUBY,N) C C EVALUATE DERIV WBT X AND Y c IF N=1 DUDX, I F N=2 DUDY, OTHERWISE DUDX AN D DUDY C IMPLICIT REAL*8 (A-H,0-Z) COMMON XX {15) DIMENSION A{1) IF (N . EQ. 2) GO TO 10 DUDX=A (2)+2.D0*A {3) *X + 3. DQ* A (4) *XX {3) + 4.D0*A(5)*X X(4) 1 *A (10) *Y+A (11) *XX{7) + A(12)*XX(8) * 2. B0*A { 1 3) *XX { 1  432 433 434  1 +2.D0*A {14)*XX (11) + 3.D0*A (15)*XX (13) IF (N. EQ. 1) BETURN DUDY=A(6) + 2.D0*A (7) *Y +3.D0*A{8) *XX (7) + 4.D0*A{9)*  4 14 4 15 416 416.5 418 4 19 420 421 422 423 424 425 4 26 4 27 428 429 430  435 436 437  10 XX <8)  1 +A{10)*X + 2.P0*A (1 1) *XX{ 10) +3.D0*A (12) *XX {1 1) 1 + A(13)*XX(3) +2.D0*A (14) *XX (13) + A(15)*XX(4) BETURN  COMPUTER PROGRAM FOR VISIOPLASTICITY 438 439 440 441 4 42 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 4 62 463 464 465 466 467 468 469 470 471 4 72 473 4 74 4 75 476 477 478 479 480 481 482 483 484 485 4 86 487  END SUBROUTINE DERIV2(A,X,Y,DUDD,N) IMPLICIT REAI*8 (A-H,0-Z) C C  EVALUATE 2ND ORDER DERIV WRT X AND Y N=1 DXDX; N=2 DYDY; N=3 DXDY  c c  DIMENSION A{1) GO TO (10,20,30),N 10 DUDD=2~ D0*A (3) + 6.D0*A(4)*X + 1 2-D0*A (5) *X*X 1 +2. D0*A (13) *Y +2.D0*A(143 *Y*Y * 6. D0*A (1 5) *X*Y RETURN 20 DUDD=2. DQ*A (7) + 6.D0*A(8)*Y + 12. D0*A (9) *Y*Y + 1 2. D0*A(11)*X * 6..D0*A (12) *X*Y + 2. D0*A ( 14) *X*X RETURN 30 DUDD=A(10) + 2. D0*A(11)*Y *• 3.D0*A{12)*Y*Y 1 + 2.D0*A(13)*X + 4. DG*A (14) *X*Y + 3. D0*A (15) *X*X RETURN END SUBROUTINE FILL(XYT,IX,IY,ISUB,XY,NIX,NIY,IMY,NTM,IDI M,IDIMP) IMPLICIT REAL*8 (A-H,0-Z) DIMENSION XYT(2,IDIM,IDIM, 2) , ISUB (I DIMP, 1) ,XY (2,IDIMP 1\  » ';  C C C C  CONTAINS 1 WHERE A FUNCTION VALUE IS PLOTTED, 0 IF OUTSIDE BOUNDARY  ISOB  10  15  20  30  /  488 489 490 491 492  104  DO 10 1=1, NIX DO 10 J=1 ,NIY ISUB(I,J) = 1 XMIN=XYT(1,IMY,IX,NTM) DO 15 1=1, IY IF (XYT (1,1,IX, NTH) .GT. XMIN) GO TO 15 XMIN=XYT ( 1 , 1 , IX, NTM) IMYMIN=I CONTINUE DO 110 J=1,NIY DO 100 1=1,NIX IF {XY(1,I,J).LT.X¥T(1,IMYMIN,IX,NTM)) GO TO 100 IF (XY(1,I,J).GT.XYT(1,IMY,1X,NTM)) GO TO 80 LB=1 NB=2 IF{XY(2,I,J).LT-XYT(2,NB,IX,NTM) ) GO TO 30 LB=LB+1 NB=NB*1 IF (NB. LT. IY) GO TO 20 IF (XY (1 ,I,J) .LT. XYT (1 ,LB,IX, NTM) . AND. 1 XY(1,I,J) . LT.XYT(1,NB,IX,NTM) ) GO TO 100 IF (XY(1 ,1,J).GT.X¥T (1, LB,IX , NTM) .AND. . 1 XY(1,I,J) .GT.XYT(1,NB,IX,NTM)) GO TO 80 YN=XYT(2 ,LB,IX, NTM) + (XYT (2 , NB,IX , NTM) -XYT (2 ,LB, IX, NTM  )  1 /{XYT(1,NB,IX,NTM)-XYT(1,LB,IX,NTM) ) * 1 (XY { 1,1, J)-XYT (1,LB,IX, NTM) ) IF (YN.GT.XY ( 2 , 1 , J) . AND. XYT (1 ,LB,IX, NTM) . GE. 1 XYT(1,NB,IX,NTM}) GO TO 100 IF (YN.LT. XY (2,I,J) . AND. XYT (1 ,LB,IX,NTM) . LE-  COMPUTES PBOGBAM FOB VISIOPLASTICITY 493 494 495 496 497 4 98 4 99 500  105 1 80 90 100 110  XYT (1,NB,IX,NTM)) GO TO 100 DO 9 0 II=I,NIX ISUB(II,J)=0 GO TO 110 CONTINUE CONTINUE BETURN END  PBOGBAM FOB TWO  DIMENSIONAL PLOT 106 IMPLICIT BEAL*8 (A-H,0-Z) COMMON XX (1 5) ,XI,YI,CC,CN,CA,Y (400, 2) , EBR (400, 20) COMMON/C1/RHO,YOLD{15,2) ,DT DIMENSION XYT {2,20,20,2),0 (20,20),V(20,20),XY{2,80,80  1 2 2.5 3 ) 4 5 6 7 8 8.05 8.07 8.2 8.6 9 10 11 12 13 14 15 16 17 19 20 21 22 23 24 25 25.1 25.2 25.3 25.4  DIMENSION X (15,400) , ERROE (2) , ISO B (80, 80) REAL*4 Z 1 (90,90) ,Z2{ 90,90) ,TAUXY (80 ,80) LOGICAL BEF.INE , DIMENSION EBRDT (80, 80) , EXDT (80,80) , EYDT (8 0,80) 1,GAMDT(80,80),LANDA(80,80) REAL*8 DEBR (80,80) BEA1*8 XLS( 10,800) BEA.L*8 LAND A DIMENSION AINT2(80) C C C C C C C Y C C OORD C C 10 20  SIGYOA IS CONSTANT ADDED TO SIGY XYTfL,I,J,K) CONTAINS: L= 1 X-COOBD: L=2 Y-C FOR 1=1,IY J=1,IX L=1,2 BEAD(5, 10) IX,IY,IT FORMAT (312) BEAD (5, 20) DT FORMAT {8F10.0) READ(5,20) CC,CN ,CA , SIG YOA ,CB  C C NIX CONTAINS # GRID PTS IN X DIRECTION, NIY IN Y DIRECTION FOB PLOTS C BEAD (5,10) NIX,NIY HEAD (5, 20)  25-7  26 27 28 29 30 31 31.2 31.4 32 33 34 35 35.2 36 37 38 38.2 38.6 39 40  IX=NO PTS IN X IY=NO PTS IN Y IT=NO OF TIME STEPS DT=TIME INTERVAL BETWEEN TIME STEPS CC,CN ARE CONSTANTS SHEBE SIGB=CC*EB**CN CA IS LOSES INTEBV AL OF I NTEGBAT.ION FOR SIG  C C C  25 30  C C  BHO  X S Y COORD BEAD FOR T.IME = 0 SEAD(4,30) ( ( (XYT {K, I , J , 1) ,K= 1 ,2) ,1 = 1 ,1Y) , J= 1 , IX) DO 25 11=1,IX DO 25 IJ = 1,IY I F (II.EQ. 1) XYT (1 , I J , I I , 1) =0. DO IF(IJ.EQ.I) XYT ( 2 , I J , I I , 1) =0.D0 DO 25 IK=1,2 IF (XYT {IK , IJ,11, 1) - LT. O.D 0) XYT (IK , I J , I I , 1) = 0. DO CONTINUE FORMAT(5X,2F6.3,1X,2F6.3,1X,2F6.3,1X,2F6.3,1X,2F6.3) CALL PLOTIT {XYT(1,1,1,1),IX,IY) LTM = 2 NTM = 1 IDIM=20 IDIMP=80 IXY=IX*IY INITIALIZE EBB TO ZEBO  PBOGBAM FOR 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 57. 2 57. 4 58 59 60 60. 2 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 84. 6 1 84. 62 84. 6 3 84. 64 84. 65 84. 6 6 84. 67 84. 68 84. 69 84. 7 84. 71 84. 72 84. 73 84. 74  TWO DIMENSION AL PLOT  40 C C  C C c  45 C C c c  50 C C C C c  55 60  DO 40 1= 1, 20 DO 4 0 J=1,400 EBR(J,I) =0. DO CONTINUE FACT= 1. DO FOB EACH TIME STEP EBB I S ACCUMULATED IT1=IT-1 DO 80 K=1,IT1 XSY COORD ARE READ FOR  NEXT TIME STEP  NTM= 3-NTM LTM=3~LTM BEAD(4,30) ( { (XYT ( K K , I , J ,NTM) ,KK=1,2) ,1=1 ,IY) ,J=1,IX) DO 45 .1=1 ,IX DO 45 J= 1,IY I F (I.EQ.1) XYT (1 , J , I,NTM) =0. DO I F (J.EQ.1) XYT(2,J,I,NTM)=0.DO DO 45 IK=1,2 I F (XYT(IK,J,I,NTM) .LT.O.) XYT (IK , J , I ,NTM) =0. DO CONTINUE CALL PLOTIT (XYT (1,1,1 , NTM) ,1X,1 Y) U,V CALCULATED FOR THIS TIME U MUST BE >0, V MUST BE <0  STEP  DO 50 J=1,IX DO 50 1=1,IY U (I , J) =— (XYT (1,1, J , NTM) -XYT(1 ,1, J , LTM) ) /DT V ( I , J ) =- (XYT ( 2, I , J , NTM) - XYT (2 ,1, J,LTM) )/DT IF (U ( I , J) .GT.O. DO) U(I,J)=0.D0 I F (V ( I , J) .LT.O.DO) V(I,J)=0.D0 CONTINUE CURVE F I T T I N G FOB U AND V USING DLSQHS SET UP INDEPENDENT VABIABLES IN X ' DEPENDENT VABIABLES IN Y DO 60 J=1,IX IL=IY*(J-1) DO 60 I=1,IY L=I+IL CALL AUX (XYT (1,1, J , NTM) ,XYT (2,1, J,NTM) ,X (1,L) ) XLS(1,L)=-X (2,L) XLS (2 ,L) =-2.D0*X(10,L) XLS (3,L) =-3.D0*X (11,1) XLS (4,L)=-4.D0*X(12,L) XLS (5,L)=-. 5D0*X (3,L) XLS(6,L) =-X{ 13,L) XLS (7,L)=-1.5D0*X(14,L) XLS (8, L) =- 1. DO/3 . D 0*X (4 ,L) XLS (9,L) =-2. DO/3. D0*X (15,L) XLS ( 10,L) =-.25D0*X (5,L) DO 55 IK=1, 10 XLS (IK,IXY + L) =X (IK+5,L) Y(L,1) = U ( I , J ) Y (IXY+L, 1) =V (I,3)  PROG.RAM FOR TWO  DIMENSIONAL  PLOT  . 108 CALL D1SQHS (Y,XLS,2*IXY,10,1,800,10,ERROR,.FALSE.,IER  84.75 ,&200) 84.76 84.77 84.78 84.79 84.8 84.8 1 84.82 84.83 84.84 84.85 84.86 84.87 84.88 84.89 84.9 84.91 84.92 85 85.1 85.2 85.3 85.4 85.5 85.6 85.7 85.8 86 87 88 89 90 91 92 92.2 93 94 95 96 96. 16 96.2 97 98 101 102 103 104 105 106 107 108 109  DO 62 IK=1,5 Y(IK,2)=0-D0 DO 63 IK=6,15 Y (IK, 2) =Y ( I K - 5 , 1) Y(1,1)=0.D0 Y (2, 1) =-Y (6, 2) Y{3 ,1)=-.5D0*Y(1Q,2) Y (4, 1) =- 1.D0/3.D0*Y (1 3,2) Y (5 ,1) =-. 25D0*Y (15,2) DO 64 IK = 6, 9 Y(IK,1)=0.D0 Y(10,1) =-2. D0*Y(7,2) Y(11,1)=-3.D0*Y(8,2) Y (1 2, 1) =-4. D0*Y (9,2) Y (13, 1) =-Y{ 1 1,2) Y(14,1)=-1.5D0*Y (12,2) Y{15, 1) =-2. D0/3.D0*Y (14,2)  62 63  64  C C C  U&V  400 4 10 C C ARE C C  COEFF SAVED FOR  I F (K.NE. (IT1-1) ) GO DO 400 1=1,15 DO 400 J = 1, 2 YOLD ( I , J)=Y ( I , J) CONTINUE THE VALUES  DU/DT,DV/DT  TO 410  OF EDTX,EDTY,GAMXY,EBEDT,AND  EBR  CALCULATED AT EACH TIME  STEP  DO 70 J=1,IX IL=IY*{J-1) DO 70 I=1,IY IJ=I+IL CALL AUX (XYT (1, I , J,NTM) , XYT (2 ,1, J,NTM) ,XX) CALL DERIVCY(1,1),XYT(1,I,J,NTM),XYT(2,1,J,NTM),DUDX, DUDY,3) CALL DERIV(Y(1,2) , XYT (1 ,1, J , NTM) ,XYT (2 ,1, J,NTM) ,DVDX, DVDY,1) GAMXY=DUDY+DVDX DFACT= (DUDX**2)/3- D0+ (GAMXY) **2/12. DO E B R D T d , J)=2.D0*DSQ8T (DFACT) 65 I F (K.EQ.IT1) FACT=. 5D0 EBR ( I J , 1 ) = E B R ( I J , 1| +FACT*DT*EBRDT ( I , J ) 70 CONTINUE 80 CONTINUE C 4TH DEGREE POLY F I T TO EBR C c CALL DLSQHS(EBR,X,IX*IY,15,2,400,15,ERROR,.FALSE.,IER ,S200) c MASTER GRID I S SET FOR FINAL PLOTS c c c X Y ( 2 , I , J ) CONTAINS THE MASTER GRID ST XY{1, I,J) IS v  PROGBAM  FOB TWO  110  DIMENSIONAL  C 1 = 1,IX  111 1 12 113 114 1 15 1 16 1 17 1 18 1 19 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134  90 100  110  120 c c c  109  PLOT  THE X-CQQBD, X Y ( 2 , I , J ) I S THE Y-COOBD FOB J=1 /IY XMAX=0.D0 DO 90 .1=1,1 Y I F (XYT{ 1,1,IX,NTM) .LT.XMAX) GO TO 90 XMAX=XYT ( 1 , 1 , I X , NTM) IMY=.I CONTINUE YMAX=XYT(2,IY,IX,NTM) CONTINUE DY=YMAX/(NIY-1) DO 110 1=1,NIX XY (2,1, 1) =0. DO DO 110 J=2,NIY XY{ 2,1,3) =XY ( 2 , 1 , J-1) +DY CONTINUE DX=XMAX/ (NIX- 1) DO 120 1=1,NIY XY ( 1, 1, I) =0.D0 DO 120 J=2,NIX X Y ( 1 , J , I ) = X Y ( 1 , J - 1 , I ) +DX CONTINUE TO FIND ZEBO F I L L S ON PLOTS CALL F 1 1 L ( X Y T , I X , I Y , I S U B,X Y,NIX,NIY,IM Y,N TH,IDIM,IDIM  P)  135 136 137 138 138. 138. 139 139. 139. 140 140. 140. 141. 141. 144 145 146 147 148 149  c c c 4 8 2 6 4 8 2 6  150 151 152 152. 8 153 154 155 155. 2 156  125  126  TO PLOT U,y,EBBDT, AND  TAUXY  DO 125 J = 1,NIY DO 125 1=1,NIX EBRDT {I, J) =0.D0 DEBR ( 1 , 3) =0. DO TAUXY ( I , J) =O.D0 NIY10=NIY+1 0 NIX10=NIX+10 DO 1.26 3=1,NIY10 DO 126 1=1,NIX 10 Z1 ( I , J) =0. DO Z 2 ( I , J) =0.D0 DO 130 3=1,NIY DO 135 1=1,NIX I F (ISUB ( I , J) - EQ. 0) G O T O 130 Z1 (1+5,3 + 5) =-AUX2(XY (1,1,3) , X Y ( 2 , I , J) ,Y{1 , 1) ,.TBUE.) Z 2 ( I + 5, J + 5) =AUX2(XY ( 1 , I , J ) ,XY ( 2 , I , J ) ,Y (1 ,2) FALSE. ) CALL DE B l V (Y,IY (1,1, J) ,XY{2,I,J) ,EXDT(I,J) ,GAMDT(I,J)  CALL DEBIV (Y (1 ,2) , XY (1,1, J) ,XY (2,1, J) , EBRDT ( I , J) , EYDT d , J ) ,3) GAMDT (I , J) = G AMDT ( I , J) +EBRDT ( I , J) DFACT=(3. D0*EXDT ( I , J) **2+. 75D0*GAMDT (I ,J) **2) 131 EBRDT ( I , J) =2. DO/3. DO *DSQRT (DFACT) DEBR ( I , J) =AUX2 (XY ( 1 , I , J ) , X Y ( 2 , 1 , J) , EBR, . F ALSE. ) LAN DA ( I , J) = 1. 5D0*EBRDT(I,J) /(CC*DEBR (I ,J) **CN) TAUXY(I,J) = GAMDT ( I , J) / (2. D0*L AND A ( I , J) ) 135 CONTINUE 130 CONTINUE  PROGBAM FOE T»0 DIMENSIONAL PLOT 157 158 159 159.4 159.6 160 161 161.2 162 162-5 162-55 162.6 162.65 162.7 162-73 162.76 162.79 162.82 162-85 162.88 162.91 162.94 163 164 165 165-2 166 166.5 167 167.5 168 169 170 171 171.2 171.4 171.6  171.8 172 172-1 172.2 172-3 172-35 172-4 172.45 172.5 172.55 172-6 172-65 172.7 172.75 172.8 172.85  C C C  PLOT U,V,EBHDT,TAUXY  CALL PLOT 2 (TAUXY,XY,IX,IY, IS OB, NIX, NIY) CALL PLOT( 1 2.,0. ,-3) DXY=XY{2,NIX,NIY)/XY (1,NIX,NIY) C CALL PERS(Zl,IDIMP+l0,NIX+10,NIY+10,DXY,- 333,45.,45 10.,10.) C CALL PLOT {12. ,0. ,^-3) C CALL PERS(Z2„IDIMP + 10,NIX+10,NIY + 10,DXY,. 333.45.,45 10., 10.) C CALL PLOT {12. ,0. ,-3) DO 137 J=1,NIY10 DO 137 I=1,NIX10 Z1(I,J)=0. 137 Z2(I,J)=0. DO 138 J=1,HIY DO 138 1=1,NIX 138 Z 1(1*5, J + 5) =DEBR (I,J) C CALL PERS(Z1,IDIMP*10,NIX+10,NIY +10,DXY, . 333,45-,45 10.,10-) C CALL PLOT (12. ,0. , — 3) DO 139 J=1 NIYlQ DO 139 I=1,NIX10 139 Z1(I,J)=G. DO 140 J=1,NIY DO 140 1=1,NIX Z1(1 + 5, J+5) = EBRDT(I,J) 140 Z2(I*5, J + 5) =-TA0XY (I, J) C CALL ?ERS(Zl,IDIMP+10,NIX+10,NIY+10,DXY,.333,45-,45 10., 10.) C CALL PLOT (12. ,0. ,-3) C CALL PEBS(Z2,IDIMP + 10,NIX+10,NIY+10 ,DXY, . 333 ,45. ,45 10.,10.) C CALL PLOT (12-, 0. ,-3) C C CALC SIGY C • EXTERNAL DF1,DF2 DO 14 1 J=1,NIY10 DO 141 I=1,NIX10 21 ( I , J) = 0. 141 Z2(I,J)=0. DGRID= (XY (2,1,2)-XY (2, 1, 1) ) /2 JGRID=1 1000 IF (CA-GT. XY (2 ,1 , JGRID) +DGRID) GO TO 1010 GO TO 1018 1010 JGRID=JGBID+1 I F (J GRID. EQ-NIY) GO TO 1018 GO TO 1000 1018 IGBID=1 1 015 I F (ISOB (IGRID+1 , JGBID) - EQ.O) GOTO 1020 IGRID=IGBID+1 I F (IGBID. EQ. NIX) GO TO 1020 GO TO 1015 1020 DO 1030 1=1,NIX XI=XY( 1,1,1) I F (XI.LE.CB) GOTO 1025 f  .PROGRAM FOR TWO 172.9 172.95 173 173.05 173. 1 173.15 173.2 173.4 173.5 173.6 173.7 1 73- 8 173-9 174 174.1 174.2 174.3 174.4 174.5 174.6 174.7 174-8 174.9 175 176.4 177  DIMENSIONAL PLOT  111  IF (XI -GT. CB -AND. ISUB (I, JGRID) -NE. .0) GO TO 1025 XI=XY (1 ,IGRID, JGRID) 1025 IF (XI - EQ. CB) GO TO 1027 A.INT2 (I) =DQUANK (DF2,CB,XI, . GO 1 DO ,TOL, FIFTH) GO TO 10 30 1027 AINT2 (I) =0. DO 1030 CONTINUE DO 150 J=1,NIY DO 150 1=1,NIX IF (ISUB (I, J).EQ.O) GOTO 150 XI=XY(1,I,J) YI=XY(2,I,J) AINT1=0.D0 IF (JGRID .EQ. J) GO TO 1060 JADD=0 IF (ISUB ( I , JGSID) . EQ. 1) GO TO 1050 JINC=1 IF (JGRI D. GT. J) JINC=-1 1040 J ADD=JADD+J INC IF (JGRID + JADD. EQ. J) GO TO 1060 IF (ISUB (I, JGRID+JADD) .EQ. 1) GOTO 1050 GO TO 1040 ,XY (2,1, J) ,. 001 DO, 1050 AINT 1=DQ UANK (DF 1, XY (2,1, JGRID + JADD) TOL,FIFTH) 1060 Z1 (1+5, J+5) =SIGYOA—AIN.T1-AINT2 (I) T AUX Y (I , J) =Z1 (1 + 5, J+5) Z2 (1 + 5,J+5)=Z1 (1+5,J + 5)*(EXDT(I,J) E -YDT(I,J) ) /LANDA (I # «i  178 178.2 178.4 178.6 179 179.5 180 180.2 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200  CONTINUE CALL PLOT2(TAUXY,XY,IX,IY,ISUB,NIX,NIY) CALL PLOTND STOP CALL PERS (Zi ,IDIHP+10,NIX+10,NIY + 10,DXY,.333,4 5.,45., 10., 10.) CALL PLOT (12. ,0. ,-3) CALL PERS{Z2,.IDIMP + 10,NIX+10,NIY+10 ,DXY,. 333 ,45.,45. , 10., 10. ) CALL PLOTND STOP 200 STOP 1 END SUBROUTINE AUX(X,Y,XX) IMPLICIT REAL*8 (A—H,0—Z) DIMENSION XX { 15) XX(1)=1.D0 XX(2)=X XX(3)=X*X XX(4) =X*XX(3) XX (5) =X*XX (4) XX(6)=Y XX (7) =Y*Y XX (8) =Y*XX(7) XX (9)=Y*XX (8) XX ( 10) =X*Y XX(1 1)=XX{10) *Y XX ( 12)=XX( 1 1) *Y XX (13)=XX (1 0) *X XX { 14)=XX(13) *Y  150  PROGRAM  F O R T S O DIMENSIONAL  PLOT 112  201 202 2 03 204 205 206 207 208 2 09  XX ( 1 5 ) = X X £13) * X RETURN END FUNCTION AUX2(XX,YY,P,LL) I M P L I C I T REA1*8 (A-H,0-Z) C C  c c BE  210 211 211- 2 2 12 213 214 215 216 217 218 219 2 20 221 2 2 2 222. 5 223 224 2 2 5 226 227 228 229 230 231 232 233 234 235 236 2 3 7 238 238. 1 238. 2 238. 3 238. 4 238. 5 239 240 241 242 243 2 4 3 .1 243. 2 244 245 246 247 248  EVALUATED  EVALUATE T H E F I T T E D FUNCTION A T XX,YY P CONTAINS T H E F I T T E D PARAMETERS LL I S TRUE I F T H E INDEPENDENT V A R I A B L E B Y A U X  MUST  C  10  c c c  L O G I C A L L L COMMON X DIMENSION X { 15) , P ( 1 ) I F (LL) C A L L A U X(XX,YY,X (1)) AUX 2=0. D O DO 1 0 1=1 ,15 A U X 2 = A U X 2 + P ( I )* X ( I ) CONTINUE RETURN END FUNCTION DF1(YC) I M P L I C I T REAL*8 (A—H,Q—Z) COMMON X X ( 1 5 ) , X I , Y I , C 1 , C 2 , C A , Y ( 4 0 0 , 2 ) , P E B R ( 2 0 , 2 0 ) COMMON/C1/RHO,YOLD{15,2) ,D T THE  INTEGRAND  DTAU/DX  I S  EVALUATED  EBR=AUX2(XI,YC,PEBR,-TRUE. ) C A L L D E R I V ( Y ( 1 , 1 ) ,XI,YC,EXDT,DUDY,3) CALL D E R I V ( Y { 1 , 2 ) ,XI,YC,DVDX,DVDY,3) GAMDT=DUDY* DVDX EBR DT=2. D0/3.D0*DSQRT(3.D0*EXDT**2+-75D0*GAMDT**2) C A L L DERIV(PEBR,XI,YC,DEBRDX,DUM,1) C A L L D E R I V 2 ( Y(1 , 1 ),X I , Y C , D U D X Y , 3 ) C A L L D E R I V 2 ( Y ( 1 , 2 ) ,X I , Y C , D VD X X , 1 ) DGAMDX=DUDX Y+ DVDXX C A L L D E R I V 2 ( Y { 1 , 1), X I , Y C ,D E X D X , 1 ) DEBRDT=(4.DO*EXDT*DEXDX+GAMDT*DGAMDX)/(3.D0*EBRDT) D F 1 = C 1 * E B R * * C 2 / ( 3 .D 0 * E B R D T ) * ( C 2 * D E B R D X * G A MDT/EBR 1 +DGAMDX - DEBRDT*GAMDT/EBRDT) U=A U X 2 ( X I , Y C , Y ( 1 , 1 ) , . F A L S E . ) V = A U X 2 ( X I , Y C , Y ( 1 ,2) , . F A L S E . ) VOLD=AUX2 (XI,YC,YOLD ( 1 , 2 ) , . F A L S E . ) F A C T = R H O * , ' ( V - V O L D ) / D T ^ K T " ~ ,:' 1 '\ DF1=DF1 +FACT ~ " " ' • . RETUSN END FUNCTION DF2(X) I M P L I C I T REAL*8(A-H,0-Z) COMMON X X { 1 5 ) , X I , Y I , C 1 , C 2 , C A , Y ( 4 0 0 , 2 ) , P E B R ( 2 0 , 2 0 ) COMMON/C1/RHO,YOLD(15,2) , D T REAL*8 LAMBDA  c  c c  T H E  INTEGRAND  DTAU/DY  I S  EBR=AUX2 (X„CA,PEBR,.TRUE. ) CALL D E R I V ( Y ( 1 , 1 ), X , C A , E X D T , D U D Y , 3 )  EVALUATED  PROGRAM FOR TWO DIMENSIONAL PLOT 249 250 251 252 253 254 255 256 257 258 259 260 261 263 264 265 266 267 268 268. 1 268.2 268.3 268-4 268.5 269 270 271 272 273 274 275 276 277 2 78 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297  113  CALL DERIV (Y (1,2) , X,CA, DVDX.EYDT , 3) GAM DT=DU DY*DVDX EBBDT=2.D0/3.DO*DSQRT(3. DO *EXDT**2 + . 75D0*GAMDT** 2) CALL DERIV {PEBR , X,CA ,DUM ,DEBRD Y, 2) CALL DERIV2 (Y (1, 1) ,X ,CA, DUDYY, 2) CALL DERIV2 (Y {1,2) , X ,CA, DV DXY , 3) DGAM DY=DUDYY +DVDXY CALL DERIV2 (Y (1 , 1) ,X,CA, DEXDY,3) YE8RDT= ( 4. DO*EXDT*DEXDY+GAMDT*DGAMD Y) / (3. D0*EBRDT) CALL DERIV2 {Y (1 , 1) ,X,CA,DEXDX , 1) CALL DERIV2 (Y (1,2) , X ,CA , DE YDX , 3) LAMBDA= { 2. D0*C 1 *EBR**C2) / (3-D0*EBR0T) CALL DERIV (PEBR,X,CA, DEBRDX, DOM, 1) CALL DEBIV2{Y(1,2),X,CA,DVDXX,1) DGA MDX=DEXDY + DVDXX DEBR DT={4.DO*EXDT*DEXDX+GAMDT*DGAMD X)/(3.D0*EBRDT) DF2=LAMBDA*({(EYDT-EXDT)*{—C 2 * DE BSD X/E BR 1 +DEBRDT/EBRDT)+DEXDX-DEYDX) 1 +.5D0*(C2*DEBRDY*GAMDT/EBR + DGAMD Y - YEBSDT*GAMDT/E SRDT) ) O=A0X2(X,CA,Y{1, 1) ,.FALSE.) V=A0X2 (X,CA,Y (1 ,2} ,. FALSE. ) OOLD=AOX2(X,CA,YOLD(1,1) ,.FALSER . F ACT=RHO*\f<{0—UOLD) /DT i*>U* " * Y> VJ DF2=DF2 • FACT ' ' " ~ -------- - a - - " ' RETURN END SUBROUTINE DERIV (A, X, Y, DU DX , DU DY , N) C C EVALUATE DERIV WRTX AND Y C I F N=1 DODX, IF N=2 DODY, OTHERWISE DUDX &N D DUD Y C IMPLICIT REAL*8(A-H,0~Z) COMMON XX(15) DIMENSION A {1) IF (N .EQ. 2) GO TO 10 DUDX=A(2) +2.D0*A (3) *X + 3-D0*A (4) *XX (3) + 4.D0*A(5)*X X{4) 1 +A (10) *Y+A (1 1) *XX(7) + A(12)*XX(8) + 2.D0*A (13) *XX (1 0) 1 +2. D0*A ( 14) *XX(11) + 3-D0*A (15) *XX (13) I F (N.EQ. 1) BETUBN 10 DUDY=A(6) + 2.D0*A (7) *Y + 3.D0*A (8) *XX (7) + 4.D0*A(9)* XX (8) 1 • A ( 1 0 ) - * X + 2. D0*A (1 1)*XX (10) + 3. DO *A { 12) *XX (11) 1 + A ( 1 3 ) * X X ( 3 ) + 2.D0*A (14) *XX {13) + A(15)*XX(4) BETUBN END SUBROUTINE DERIV2(A,X,Y,DODD,N) IMPLICIT REAL*8(A-H,0-Z) C C EVALUATE 2ND OBDEB -DERIV WRT X AND Y C N=1 DIDX; N=2 DY DY; N=3 DXDY C DIMENSION A (1) GO TO (10,20,30) ,N 10 DUDD=2. D0*A (3) + 6.D0*A(4)*X • 1 2. D0*A (5) *X*X t  v  PBOGRAM FOR TWO 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337  DIMENSION A l PLOT  1 +2. D0*A(13)*Y +2.D0*A{14) *Y*Y + 6. D0*A(15)*X*Y RETURN 20 DUDD=2.D0*A(7) • 6. D0*A (8) *Y + 12. D0*A (9) *Y*Y + 1 2.DQ*A(11)*X • 6. D0*A(12) *X*Y + 2. DO*A (1 4) *X*X RETURN 30 DUDD=A(1Q) + 2. D0*A(11)*Y • 3. DO* A (12) *Y* Y 1 + 2.D0*A(13)*X *• 4. D0*A (1 4) *X*Y + 3. DO *A (15) *X*X RETURN END SUBROUTINE Fill(XYT,IX,IY,ISUB,XY,NIX,NIY,IHY,NTM,ID I M,IDIMP) IHPLICIT REAL*8 (A-H,0-Z) DIMENSION X YT (2,IDXM , IDIM, 2) , I SUB (I DI MP , 1) ,XY(2,IDIMP ,D C C ISUB CONTAINS 1 IBE.RE A FUNCTION VALUE IS C PLOTTED, 0 IF OUTSIDE BOUNDARY C DO 10 1=1,NIX DO 10 J=1,NIY 10 ISUB(I,J)=1 XMIN=XYT( 1,IMY,IX,NTM) DO 15 1=1,IY IF (XYT ( 1,1, IX, NTM) .GT. XMIN) GO TO 15 XMIN=XYT (1,1,IX, NTM) IHYMIN=I 15 CONTINUE DO 1 10 J=1,NIY DO 100 1=1,NIX I F (XY(1,I,J) -IT. XYT (1,IMYMIN,IX ,NTM) ) GO TO 100 IF (XY(1,1,J).GT.XYT { 1,IMY,IX,NTM)) GO TO 80 IB=1 NB=2 20 IF(XY(2,I,J) -IT. XYT (2,NB ,IX,NTM) ) GO TO 30 1B=1B+1 N B= N B + 1 IF (NB.IT.IY) GO TO 20 30 IF(XY(1,I,J) .IT. XYT (1 ,1B ,1 X, NTM) .AND. 1 XY(1,I,J) .LT. XYT (1 ,NB,IX,NTM) ) GO TO 100 IF (XY(1,I,J) .GT-XYT(1,LB,IX,NTM).A ND. 1 XY(1,I,J).GT.XYT (1 , NB, IX, NT M) ) GO TO 80 YN=XYT(2,LB,IX,NTH) + (XYT(2,NB,IX,NTM)-XYT(2,LB,IX,NTM ) )  338 339 340 341 342 343 351 352 353 354 355 356 357 358 359  114  80 90 100 110  1 /(XYT(1,NB,IX,NTM)-XYT(1,LB,IX,NTM) ) * 1 (XY (1,1, J) -XYT (1, LB, IX, NTM) ) IF (YN.GT.XY(2,I,J).AND.XYT{1,LB,IX,NTM).GE. 1 XYT(1»NB,IX,N3M)) GOTO 100 IF (YN.IT. XY (2,1, J) .AND. XYT (1 ,IB,I X,NTM)-IE. 1 XYT (1,NB,IX,NTM) ) GOTO 100 DO 90 11=1, NIX ISUB (II , J) =0 GO TO 110 CONTINUE CONTINUE RETURN END SUBROUTINE P10.TIT (XYT ,IX,IY) IMP1ICIT REAL*8 (A-H,0-Z)  PROGRAM FOR TWO DIMENSIONAL PLOT  115 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 3 80 381 382 383 384 385 386 387 388 389 390 3 90-2 391 392 393 394 395 396 397 398 399 400 401 403 4 04 4 05 407  DIMENSION XYT (2,20,20) CALL AXIS (0.,0.,'X ,- 1,10- ,0.,0.,.2) CALL AXIS (0, ,0- , ' Y • , 1 , 10. ,90- ,0- ,-2) DO 10 J=1,IY CALL PLOT (XYT { 1, J , 1) *5- , XYT (2,J,1)*5.,3) CALL SYMBOL (XYT(1,J,1)*5.,XYT ( 2 , J , 1 ) * 5 . 1 4 , 4 , 0-,-1) DO 10 1=2, IX CALL PLOT (XYT (1 , J,I) *5- , XYT (2, J, I) * 5. , 2) 0. ,-1) CALL SYMBOL (XYT (1, J , I) *5. , XYT (2 , J,I) *5.,-14,4, CONTINUE DO 20 J=1,IX CALL PLOT(XYT(1, 1,J) *5. , XYT (2 , 1 , J) *5 . , 3) DO 20 1=1,IY CALL PLOT(XYT(1,I,J)*5./XYT (2 , 1 , J) *5.,2) CONTINUE CALL PLOT (12.,0. ,-3) RETURN END SUBROUTINE PLOT 2 (TX Y , XY, IX , IY , ISUB, NIX , NIY) REAL*8 XY(2,80,80) DIMENSION TXY (80,80) ,ISUB(80,80) /TAUXY (80,80), X (80) DO 5 1=1,80 DO 5 J=1 ,80 TAUXY (I, J) =0.D0 DO 10 1=1,NIX DO 10 J=1,NIY TAUXY (I, J) =TXY ( I , J) CALL SCALE(TAUXY,6400, 1Q.,YMIN,DY,1 ) CALL AXIS (0.,10., »X»,1,10.,0.,0.,.2) CALL AXIS (0. ,0. * ,5, 10. ,90. ,YMIN,DY) NY=NIY/IY*3 KK=0 DO 50 J=1,NIY,NY L=NIX IF(ISUB (L, J) . EQ. 1) GO TO 20 L=L-1 GO TO 15 DO 30 1= 1 ,L X(I) =XY ( 1/1, J) *5. CALL LINE (X ,TAUXY ( 1 , J) , L , 1) NX=NIX/IX KK=KK+1 DO 40 K=1,L,NX CALL SYMBOL (X (K) ,TAUXY(K, J) ,. 14, KK, 0.,-1) CONTINUE CONTINUE END 1  10  20  5 10  15 20 30  40 50  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0100192/manifest

Comment

Related Items