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.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

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., University of Roork.ee;. 1971 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE i n THE FACULTY OF GRADUATE STUDIES (Department of Mechanical Engineering) We accept th i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA November, 1978 (c) Surendra Nath Dwivedi, 19 7 8 In presenting t h i s thesis i n p a r t i a l f u l f i l l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t freely available for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the Head of my Depart-ment or by his representatives. I t i s understood that copy-ing or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. (Surendra Nath Dwivedi) Department of Mechanical Engineering The University of B r i t i s h Columbia 20 75 Wesbrook Place Vancouver, B.C., Canada V6T 1W5 D a t e •  3 o [ f 11 7 g ABSTRACT The v i s i o p l a s t i c i t y approach i s developed to enable the complete stress history of any steady or non-steady, quasi-s t a t i c or impact, plane s t r a i n p l a s t i c deformation process to be determined from a record of the deformation pattern. The ve l o c i t y f i e l d i s determined experimentally and for dynamic conditions high speed photographs are taken of a g r i d pattern marked on the end surface of the specimen. D i g i t i z a t i o n of the instantaneous g r i d node positions allows the v e l o c i t y f i e l d s to be obtained at predetermined time i n t e r v a l s through-out the transient deformation period. Hence, the s t r a i n - r a t e , equivalent s t r a i n rate, equivalent 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 obtained. A three dimensional surface f i t t i n g procedure, using fourth order polynomials, i s used to smooth the scalar com-ponent of the experimentally determined v e l o c i t y f i e l d . The condition of continuity ( e-^=-e^ r^ • for plane strain), i s imposed on both surfaces thereby reducing the number of i n -dependent parameters from 30 to 10. Besides smoothing the experimental points t h i s procedure has the d i s t i n c t advantage that the polynomials can be read i l y d i f f e r e n t i a t e d for de-termining strain-rates and that deformation can be referred to a master reference g r i d that i s fixed with respect to time. Plane-strain upsetting tests, conducted at a speed of 0*02 ft/min give results that agree closely with the well docu-merited ' f r i c t i o n h i l l 1 type of normal stress d i s t r i b u t i o n for qua s i - s t a t i c rates of s t r a i n . However.,with the specimen deformed at a speed of 15;7 ft/sec the normal stress d i s t r i b u t i o n i s r a d i c a l l y d i f f e r e n t exhibiting a saddle type d i s t r i b u t i o n . The e f f e c t of s t r a i n rate on the interface and body stresses w i l l have s i g n i f i c a n t bearing on a number of metal forming operations. iv' ACKNOWLEDGEMENTS The author sincerely thanks the thesis supervisor, Professor G.W. Vickers for his unending words of encouragement and whose amiable attitude made thi s work a pleasant experience. Sincere thanks must also go to the other members of the thesis committee i . e . , Professor H. Ramsey, I.S. Gartshore, N.G. Ely, R.E. McKechnie and Norman Franz for t h e i r time to time f r u i t f u l discussions and expert opinions. The author i s most grateful to a l l the other members of the department for t h e i r i n s p i r i n g i n t e r e s t and assistance. Thanks are also due to my wife Shashi for some typing help and for her moral support. Above a l l , I thank God for giving me the opportunity to better myself academically and s p i r t u a l l y during the course of the research. 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 PLASTICITY IN METAL WORKING 3 1.1.1 Slab 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 F i n i t e Element Method 6 1.1.6 F i n i t e Difference Method . .. 7 1.2 VISIOPLASTICITY 8 CHAPTER TWO - NON-STEADY, PLANE-STRAIN, DYNAMIC AND QUASI-STATIC VISIOPLASTICITY 9 2.1 EQUATIONS FOR QUASI-STATIC VISIOPLASTICITY 10 2.1.1 Equations for Three Dimensional Non-steady State, Quasi-static, v i s i o -plas t i c i t y 10 2.1.2 Equations for Plane-Strain, Non-steady State, Quasi-static V i s i o p l a s t i c i t y 13 2.1.3 Determination of the Stress F i e l d from the Strain F i e l d 14 2.2 EQUATIONS FOR DYNAMIC VISIOPLASTICITY 16 2.2.1 Equations for Three Dimensional, Non-steady State, Dynamic, V i s i o p l a s t i c i t y . . . 16 v i Page 2.2.2 Equations for Plane Strain Non-Steady State Dynamic V i s i o p l a s t i c i t y 18 2.2.3 Determination of the Stress F i e l d from the Strain F i e l d 18 2.3 GENERAL PROCEDURE FOR THE STUDY OF DEFORMATION USING VISIOPLASTICITY 19 2.4 GENERALIZED EQUATIONS FOR PLANE-STRAIN VISIOPLASTICITY FOR COMPUTATIONAL PURPOSES 2 3 2.5 THE COMPUTER PROGRAM 2 9 CHAPTER THREE - VARIABLE SPEED, CONTROLLED VELOCITY PROFILE, SINGLE CYCLE IMPACTING PRESS.. 36 3.1 DESCRIPTION OF EQUIPMENT 37 3.1.1 Background Information 3 7 3.1.2 Description of the Press 39 3.1.3 Operation of the Press 39 3.2 KINEMATIC ANALYSIS OF THE MECHANISM 40 3.2.1 Derivation of Expressions for Velocity and Acceleration 4 0 3.2.2 Maximum Velocity and Acceleration of the Mechanism 4 3 3.3 CONTROL 47 3.3.1 Control for Single Cycle Operation 4 7 3.3.2 Control for Synchronizing the High Speed Camera 5 0 CHAPTER FOUR - EXPERIMENTAL PROCEDURE AND DISCUSSION.. . . 51 4.1 EXPERIMENTAL PROCEDURE 51 4.2 DISCUSSION 53 4.3 SOURCES OF ERRORS 84 v i i Page CONCLUSIONS 8 6 SUGGESTIONS FOR FURTHER, WORK 8 7 REFFERENCES-, 88 APPENDIX 95 v i i i LIST OF FIGURES Page 1. Elemental cube for derivation of the equation of motion 17 2'. Position of point P at " instant t n ' and't n+i . ... 2 0 3 .-- Grid-lines -on the specimen. .- : . V 31 4. Flow chart of computer program 32 5. Variable speed, controlled v e l o c i t y p r o f i l e , single cycle impacting press 38 6. Drive mechanism of impacting press 41 7. Motion of the c y c l o i d a l cam 44 8. Details of e l e c t r i c a l c i r c u i t for control 46 9. Plow diagram for relay sequence for single cycle control 4 8 10. C i r c u i t operation block diagram for high speed camera synchronization 49 11. Distortion of g r i d l i n e s during deformation for 0.02 ft./min. deformation speed (a) F i r s t step; and (b) second step 55 12. Distortion of g r i d l i n e s during deformation for 0.02 ft./min. deformation speed (c) Third step; and (d) fourth step.... 56 13. Distortion of g r i d lines during deformation for 15.7 ft . / s e c . deformation speed. ( F i r s t step) 57 14. Distortion of g r i d l i n e s during deformation for 15.7 f t . / s e c . deformation speed. (second step).... 58 15. Distortion of g r i d l i n e s during deformation for 15.7 f t . / s e c . deformation speed. (third step) 59 16. Distortion of g r i d l i n e s during deformation for 15.7 ft . / s e c . deformation speed. (fourth step).... 60 i x Page 17. Distortion of grid l i n e s during deformation for 15.7 f t . / s e c . deformation speed. ( f i f t h step).... 61 18. Distortion of g r i d l i n e s during deformation for 15.7 f t . / s e c . deformation speed. (sixth step).... 62 19. Distortion of g r i d l i n e s during deformation for 15.7 f t . / s e c . deformation speed. (seventh step).. 63 20. Grid node point movements during deformation for 0.02 ft./min. deformation speed 64 21. Grid node point movements during deformation for 15.7 f t . / s e c . deformation speed. 65 22. Three dimensional p l o t of horizontal v e l o c i t y (u) as a function of x and y for 0.02 ft./min. deformation speed 66 23. Three dimensional p l o t of v e r t i c a l v e l o c i t y ( v) as a function of x and y for 0.02 ft./min. deformation speed 67 24. Three dimensional plot of horizontal v e l o c i t y (u) as a function of x and y for 15.7 f t . / s e c . deformation speed 68 25. Three dimensional plot of v e r t i c a l v e l o c i t y (y0 as a function of x and y for 15.7 f t . / s e c . deformation speed 69 26. Three dimensional plot of e f f e c t i v e s t r a i n - r a t e as a function of x and y for 0.02 ft./min deformation speed 70 27. Three dimensional plot of e f f e c t i v e s t r a i n - r a t e as a function of x and y for 15.7 f t . / s e c deformation speed 71 28. Three dimensional plot of t o t a l e f f e c t i v e s t r a i n as a function of x and y for 0.02 ft./min deformation speed 72 29. Three dimensional plot of t o t a l e f f e c t i v e s t r a i n as a function of x and y for 15.7 f t . / s e c 73 X Page 30. Three dimensional p l o t of normal stress ( a v) as a function of x and y for 0.02 ft./min. deformation speed 7 4 31. Three dimensional plo t of normal stress ( a x) as a function of x and y for 0.02 ft./min. deformation speed 75 32. Three dimensional plot of shear stress ( T X y ) as a function of x and y for 0.02 ft./min. deformation speed 76 33. Three dimensional plo t of normal stress ( Oy) as a function of x and y for 15.7 f t . / s e c . deformation speed 77 34. Three dimensional p l o t of normal stress ( a x) as a function of x and y for 15.7 f t . / s e c . deformation speed. 7 8 35. Three dimensional p l o t of shear stress ( x x v ) as a function of x and y for, 15.7 f t . / s e c . deformation speed 7 ^ 36. Two dimensional plot of 1 o. y ' as a function of x for 0.02 ft./min. deformation speed... 8 0 37. Two dimensional p l o t of ' T X y ' as a function of x for 0.02 ft./min. deformation speed 81 38. Two dimensional p l o t of ' a- y' as a function of x for 15.7 ft./min. deformation speed 8 ^ 39. Two dimensional plo t of x X y as a function of x for 15.7 ft./min. deformation speed 83 NOMENCLATURE yt the component of v e l o c i t y i n x - d i r e c t i o n v the component of v e l o c i t y i n y - d i r e c t i o n W the component of v e l o c i t y i n z - d i r e c t i o n £ s t r a i n r a t e alone x - d i r e c t i o n x * £y . s t r a i n r a t e along y - d i r e c t i o n V ; shear s t r a i n - r a t e i x y £ e f f e c t i v e s t r a i n - r a t e £ t o t a l e f f e c t i v e s t r a i n 0 e f f e c t i v e s t r e s s 0" normal s t r e s s along 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 along z - d i r e c t i o n T shear s t r e s s xy p the den s i t y of the m a t e r i a l V the volume of the p l a s t i c a l l y deforming body £ vol u m e t r i c s t r a i n - r a t e v F the t r a c t i o n on the p a r t of surface S r 0 the v e l o c i t y presented on the remainder surface S. A Lagrange m u l t i p l i e r a the penalty f u n c t i o n 0 the angle made fcyrthe d r i v i n g arm w i t h v e r t i c a l C H A P T E R O N E INTRODUCTION & LITERATURE SURVEY 1 2 CHAPTER ONE  INTRODUCTION The mechanism of p l a s t i c deformation plays a v i t a l role i n many i n d u s t r i a l metal-working processes. However, i t has not proved possible to analyse completely many of these processes using the general basic equations "derived from the - theory - of p l a s t i c i t y . This i s primarily due to unclearly defined boun-dary conditions; for example, the actual f r i c t i o n a l conditions present at the metal-die interface are frequently unknown. Many s i m p l i f i e d alternative methods have been developed and used to study certain of the metal forming process. In these analyses certain assumptions and s i m p l i f i c a t i o n s are made regarding the processes and the behaviour of the materials during deformation. However, i n spite of these i d e a l i z a t i o n s the solutions often lack uniqueness and completeness. One approach 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 to determine the complete stress picture through-out the deformation zone i n certain steady-state extrusion and forging operations. I t requires that the v e l o c i t y f i e l d be determined experimentally and hence the strain-rates ] .;. A .v and f i n a l l y stress f i e l d s can a l l be obtained. This method has been shown to give r e a l i s t i c solutions and i t s application has been extended during 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 material deformation i n dynamic and non-steady condi-3 tions. The relevant equations and procedure have been embodied into a s p e c i a l l y developed computer program,. so that the com-plete stress history of any steady or non-steady, q u a s i - s t a t i c or impact plane s t r a i n , deformation can be determined from a record of the deformation pattern. Special attention has been given i n t h i s work to smoothing the experimentally determined v e l o c i t y f i e l d s , a point which has caused some d i f f i c u l t y i n the past. Results from t h i s work have been suitably compared with previous steady-state results for v e r i f i c a t i o n purposes. 1.1. PLASTICITY IN METAL-WORKING While the analyses of metal-working processes has been re-s t r i c t e d by the complexities involved, some approaches have been made. A number of these i n common use are the slab (or equilibrium) method, uniform deformation energy method, s l i p l i n e solutions, upper and lower bound solutions, f i n i t e difference and f i n i t e element methods. For completeness a b r i e f description of these common approaches i s given below. 1.1.1. Slab or Equilibrium Method. The method introduced by Sachs (1) i n 19 31, consists of i s o l a t i n g a small elemental volume of the material under goihrg deformation and observing the behaviour of thi s element as i t moves through the working zone. Since t h i s element i s an i n -tegral part of the material, i t should always be i n a state of equilibrium. The assumption i s made that stresses on a plane surface perpendicular to the dir e c t i o n of the flow are p r i n c i p a l stresses and that these do not vary on thi s plane. 4 Analysis of the equilibrium condition results i n one or more d i f f e r e n t i a l equations which together with the necessary boundary conditions, give the deformation stresses. Since the e f f e c t of redundancy, f r i c t i o n and pattern of flow are not considered, t h i s method gives an underestimate of the deformation stresses. However, the analysis i s s t r a i g h t -forward and i t has been widely used i n wire and tube drawing problems as well as hot and cold r o l l i n g of s t r i p and sheet (1). 1.1.2. Uniform Deformation Energy Method: Siebel (2) proposed t h i s approach i n 1932 i n which the amount of deformation i s determined by considering the shape of an element of material before and a f t e r deformation. I t hence gives only the average forming pressure as a function of s p e c i f i c i n t e r n a l energy and i s generally used for steady-state metal-working processes. 1.1.3. S l i p Line Method: Hencky(3) introduced the s l i p l i n e theory i n 1923. I t can be used for determining the l o c a l stress and v e l o c i t y d i s t r i b u t i o n during deformation, although i t i s r e s t r i c t e d to plane s t r a i n conditions and requires a predetermined pattern of flow. The s l i p l i n e solution consists of families of c u r v i l i n e a r or straight l i n e s , which are perpendicular to each other and correspond to the directions of maximum and minimum constant shear stress. These lines s a t i s f y the s t a t i c equilibrium con-d i t i o n , y i e l d condition and the pattern of flow everywhere i n the p l a s t i c zone of the material. These shear or s l i p l i n e s 5 are known as c h a r a c t e r i s t i c s of 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 tool and the material outside the s l i p l i n e are considered as r i g i d ( i . e . the metal ahead and behind the p l a s t i c zones and the tool material have an i n f i n i t e modulus of e l a s t i c i t y ) . The s l i p l i n e solution i s not optimal or unique and also gives values higher than the true solution. This method has been widely used for the study of many metal deformation processes (4-14), some of the l a t e s t work has involved s l i p l i n e solutions for anisotropic materials (15, 16) and has taken account of f r i c t i o n on the die-workpiece interface (17, 18). Also Ewing (19) and l a t e r C o l l i n s (20) have produced s l i p l i n e solution using numerical computation by power series and by matrix operational methods. 1.1.4. Limit Analysis: The mathematical model of l i m i t analysis places upper and lower estimates on the load required for deformation. This l i m i t analysis i s based on two extremum theorems put forward by Prager and Hodge (8), and Drucher and Prager (21). H i l l (7) gave the mathematical proof of these theorems, which are based on the assumption that the material i s r i g i d and p e r f e c t l y p l a s t i c . They can be stated as: a) Upper Bound Theorem. If a kinematically admissible v e l o c i t y f i e l d e x i s t s , the loads required to be applied to cause the v e l o c i t y f i e l d to operate constitute an upper bound solution. 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 xists such that the stresses are everywhere just below 6 those necessary to cause y i e l d i n g , then the loads associated with that f i e l d constitute a lower bound solution. These techniques have been used extensively (22-39) to study metal-working processes, such as forging, 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 of the most powerful tech-niques for solving two dimensional problems i n metal-working but at present has a li m i t e d p o t e n t i a l for complex problems due to economic constraints. This method was introduced by Argyris (40) i n 19 54. In t h i s approach, the deforming area or continuum i s subdivided into an equivalent system of elements, known as f i n i t e elements. The f i n i t e elements may be t r i a n g l e s , group of t riangles, q u a d r i l a t e r a l etc. for two dimensional studies and tetrahedra, rectangular prisms or hexahedra etc. for three dimensional studies. Once a displacement model i s selected, an element s t i f f n e s s matrix i s derived using v a r i a t i o n a l p r i n -c i p l e s . The algebraic equations for the whole continuum are then assembled and solutions for unknown displacements at the nodal points can be obtained. By use of the computed displace-ments and the stress and s t r a i n r e l a t i o n s , the stresses at the nodal points may be determined. The solution i s based on ex-tremum p r i n c i p l e according to which the actual solution mini-mizes the functional $ , where <f>= / a e d V r / I . y dS with the constraint that & = 0 v 7 where v = volume of the p l a s t i c a l l y deforming body cP = e f f e c t i v e stress e = e f f e c t i v e s t r a i n rate i = volumetric s t r a i n rate v F = the t r a c t i o n of the part of surface S p and ii = the v e l o c i t y prescribed on the remainder surface A modified functional has been given by Lee & Kobayashi (41) using the Lagrange m u l t i p l i e r so that cb = / CT e dV + 17 A e dV - F • udS v V O p where A = the Lagrange m u l t i p l i e r . While Godbole and Zienkie-wicz (42) have suggested that the functional be modified using a penalty function, a , as follows: cb = / a e d V + / -^ (e ) 2 dv - / F • u dS SF where a = very large number. By use of above functionals, many metal-working processes such as upsetting drawing, piercing etc. (40-57) have been studied. 1.1.6. F i n i t e Difference Method: This method i s one of the most recent techniques to be used for the study of metal-working processes. I t requires that the continuum be divided into a number of grids and that 'difference 1 ( i . e . f i n i t e ) quantities are substituted for d i f f e r e n t i a l quantities across the grids. Thus for a given d i f f e r e n t i a l equation with boundary conditions a set of simul-taneous equations can be substituted, which can be solved numerically using a computer. The size of the g r i d spacing determines the accuracy of the solution. The f i n e r the g r i d , 8 the better i s the accuracy obtained, but this at the expense of computer cost. Studies which describe the application of this technique to forging, extrusion and sheet-metal processes are given i n references (58-64). 1.2. VISIOPLASTICITY The v i s i o p l a s t i c i t y was introduced by Thomsen (56, 66, 6 7) and l a t e r developed and extended by Shabaik, Kobayashi et a l . (68, 69, 70, 71, 72). In t h i s method, the g r i d l i n e patterns are photographed for each incremental step of deformation and thus the movement of g r i d points can be determined. From en-larged photographs of consecutive g r i d patterns the instantan-eous v e l o c i t i e s of a l l g r i d node across the surfaces can be found. The s t r a i n s , s t r a i n rates, t o t a l e f f e c t i v e s t r a i n can thus be determined for a l l points i n the deformation region and f i n a l l y the stress 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 actual one and gives information of a l l strains and stresses over the entire deformation region. I t may be used for both work-hardening and non-workhardening materials. Details of the basic 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 applied to forging and axisymmetric and plane s t r a i n extrusion and r o l l i n g pro-cesses. (Reference 68 to 76). Recently i t has been used for investigating the relationship between s t r a i n and microhardness (77), crack propagation and for the derivation of c r i t e r i a for duc t i l e rupture of f u l l y p l a s t i c notched bars (78). 9 C H A P T E R T W O 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 QUASI-STATIC VISIOPLASTICITY 2.1.1 Equations for Three Dimensional Non Steady State Quasi-s t a t i c V i s i o p l a s t i c i t y  The following equations i n three dimensions are used to describe the mechanism of p l a s t i c deformation of an i s o t r o p i c s o l i d . • • • The s t r a i n rate £ , e . ,y can be given terms of x y xy v e l o c i t y components as follows: • _ 3u | 2 ; 3 V E = , 9 W x 3 x, y 3y, z 3 z ; Y '/ _ 3 u + 3 v Y _ 3w + 3 v and xy ~ 3y ax', yz ^y 3z 7Y. = _3_u + 3w zx 3z 3 x (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 equilibrium, neglecting a l l body forces are: 3 a 3 i : 3 T 3x 3y 9z 8 Txy + i i x + 5 T y z = o 3x g y 3 z 11 3T da 3 T x z + IJl + 2 = o 3 x 3 y 9 z The L e v y - V o n M i s e s s t r e s s a n d s t r a i n r a t e r e l a t i o n s h i p s ( o r f l o w r u l e ) i s g i v e n b y • • • • * • E x _ £ y _ £ z _ Y x y _ ^ y z _ ^zx _ \ a +p ""' a +p a +P ~ 2^ ~ 2 T ~ 21 x ^ y z x y y z z x ! ( 2 . 3 ) w h e r e p = - j ( CT v+ CT „ + °rr) / x 2 <-a e = e f f e c t i v e s t r a i n r a t e , a n d o = e f f e c t i v e s t r e s s . The Von M i s e s y i e l d c o n d i t i o n i s 2 2 2 2 2 2 -2 ( a - a ) z + ( a - C T ) ^ + ( a - a ) +6 ( r + T + T ) = 20 x y Y z z x xy z y z x — ( 2 . 4 ) The a b o v e e q u a t i o n ( 2 . 4 ) c a n be e x p r e s s e d i n t e r m s o f p r i n c i p a l s t r e s s e s a ^, c 2 a n d as ( V a 2 ) 2 + ( a 2 - a 3 ) 2 + ( a ^  o ^ ) 2 = 2 a 2 w h e r e a = e f f e c t i v e s t r e s s , w h i c h i s c o n s t a n t f o r p e r f e c t l y p l a s t i c m a t e r i a l s a = crCe) ^ f o r n o n - w o r k h a r d e n i n g p l a s t i c m a t e r i a l a = CT(£) f f o r w o r k h a r d e n i n g p l a s t i c m a t e r i a l a = ( e f e , T ) , f o r m a t e r i a l w h i c h i s a f f e c t e d b y s t r a i n r a t e , s t r a i n a n d t e m p e r a t u r e . From e q u a t i o n ( 2 . 3 ) , we may w r i t e 12 • 3— G 3e X 25 • e = * Y e z 2 a r 3~ Y _ t xy 2^ • • z" j •=-a Y = s i 1 . zx a X(a +p) x = ^ >a y +p) = . A ( a y +p) A ( a z +p) T ^ 2 * T xy , xy = 2 A x y z zx A T zx (2.5) Now subtracting the second equation (2.5) from f i r s t equation (2.5) gives e _ e = A r ( a - a ) x Y 2 a . x 7 and s i m i l a r l y for the other equations (2.5) e _ = 3~ ( a ~ a .) Y z? — Y z 2 a = 3 e ( a ~ a ) Z X Z X 2 a 3 Y ^ T X Y x ^ 2 a 3 y 3 x 2 f T Y Z = r z y z 2 a A Tzx = V- x 2 t z x ^ 7 ( 2 - 6 ) 2 a Squaring the l e f t hand side and right hand side of a l l the 7 equations (2.6) and adding, we get (e " e ) 2+( e - e ) 2+( e - e ) 2 + \( y 2 + Y 2 + Y 2 x y y z z x xy 7 Z , z x 13 /3 e' I(a - a ) 2 + ( a - a ) 2 + ( a -a ) 2+6T 2+6T 2+6T 2 ] L x y . y z z X x 7 *z z* J (2.7) Combining equation (2.4) and (2.7) gives . , 9 • • v 2 , • • >2 3 , '2 ' 2 - 2 , ( £ x - £ y ) 2 + (e - e z ) . 2 v + (e z-e x) + T ( Y + Y + Y ) J y xy yz zx f 0-2 9 - 2 |x 2a — 3 e So that the equivalent s t r a i n rate may be given as ! 2 '< ( e — e IF 1 x e ) 2 + ( e " e ) 2 + ( e - e ) ^ 3 ( y 2 + y ^ + - l \ [ y y. z z x "z\ x y y z . Y z x i f ; - e ) 2 + ( e - e ) 2 + ( e - e ) 2 > x y y z z x o . • 2 - 2 . 2 I ( Y + TJ- 'xy ( y + Y + Y ), yz zx ( e - e ) 2 + ( e - e ) 2 + ( I - 'e)2 y - z 3 ( - 2 + ' 2 + * 2 Y + Y + Y : '4* X Y Y z z x 1 2 (2.8) 2.1.2. Equations for Plane Strain Non-Steady State Quasi-Static V i s i o p l a s t i c i t y  For the plane s t r a i n condition, e , e , Y „ „ # y „„r Z 2 JL ZX Y , y , x and ' yz zx yz zx ~ ', also -CL = i ( f o r a r i g i d p e r f e c t l y p l a s t i c material). Hence the basic equations for plane s t r a i n are given as follows: The equation for s t a t i c equilibrium from equation (2.2) i s da 9T • x + . .xy ' 3x ",9y = 0 9T 3a + • 3x 9y Y _ (2.10) 14 The flow rule from equation (2.3) now becomes 8 x _ , E y _ Yxy _ • (2.11) a +p a +p 2 T x y c xy where A = 3 | e a Von Mises y i e l d c r i t e r i a from equation (2.4) becomes ( a -:,a ) 2+ \a - i ( a + a )l 2 + [ i ( a + a )~ a ] 2 x y 1 y 2 x y J L 2 u x u y xJ +6( x 2+ 0 + 0) = 2 a 2 xy or a'-= [ -nfa - a ) 2+3 T 2 ] ' 5 4* x y xy J (2.12) The e f f e c t i v e s t r a i n rate can be written from equation (2.8) as L 2 - " 2 3 * 2, ^  -E = T [ 3 e x + T Y ^ y ; J : (2.i3) 2.1.3. Determination of the stress f i e l d from the s t r a i n f i e l d Once a l l the normal and shear strains are known throughout the deforming region from equation (2.1) the stresses may be calculated. To determine the stress f i e l d from the s t r a i n f i e l d , the following steps are required. From equation (2.3) e e _ 3 P ( a +p^ a -p) 3e ( a ~ a ) x y - -5- x ^ y r = -- u x °y 2.a 2 a 15 2.0-a x ~ (e - e ) = x y a x £x-£y Di f f e r e n t i a t e equation (2.14) with respect to x gives (2.14) •3a, 3a x . 3x 3x e — e x -> y 3 X L A _J From the equilibrium equation (2.10) (2.15) 3a x 3T xy 3 X 3y Substituting i n equation (2.15) we get 3a Y 1 3T xy 3x 3y 3 T x E - e x . y (2.16) From equation (2.10) 3a, 9 T . xy 3y 3x (2.17) Using equation (2.16) and (2.17) with the known value of 0"(x,y) at x = 0 and y = bt i . e . a(o,a) we get, X - '3T V 3 T a (x,y) = a (0,a) - • •*( - xy dy - f y Y a ' av J , £ - £ xy+ 3 \ x y \ 3y 3x\ ^ / dx —1 y=a (2.18) From equation (2.14) x f y + ^ F ^ (2.19) and from equation (2.11) xy xy 2x (2.20) where A- = c T Z (2.21) 2.2 EQUATIONS FOR DYNAMIC VISIOPLASTICITY 2.2.1 Equations for Three Dimensional Non-Steady State Dynamic V i s i o p l a s t i c i t y  For the case of dynamic deformation the Levy-Von Mises stress and s t r a i n rate r e l a t i o n s h i p and Von Mises y i e l d c r i t e r i a are s i m i l a r to the three dimensional non-steady q u a s i - s t a t i c case. However, the s t a t i c e q u i l i b r i u m equation i s replaced by the equation of motion. The equation of motion i s obtained by considering a generic elemental cube subject to three normal and three independent shear stresses as shown i n F i g . (2.1). I f x,y are the c u r r e n t coordinates of a p a r t i c l e then x =• x ( x o ' V t } y = y ( x 0 ' Y o ' h ) where x 0/ y Q are the i n i t i a l coordinates at time t = 0 . The components of v e l o c i t y along the x, y, z axes are given by u , v , ' W r e s p e c t i v e l y . Writing.the equation of motion along x axis gives (c + x dx) dy dz - a dy dz + (x °- yx dy)dz dx x x ~ T i x y ay - T . , _ dz dx + ( T z x + -y^- ) dx dy - T z x d x dy yx (dx dy d 2 l - ^ { P u } ( 2 . 22 ) Dividing throughout by dx dy dz 3o 3 T 3T . • ,. / »> o o \ xx j. yx + . zx = £ _ ( p u ) 3~>: 3y 3 z * p<• S i m i l a r l y equation of motion along 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 FOR DERIVATION O F T H E EQUATION OF M O T I O N 18 Hy* + Jlxi + 2lx* = j | C p v ) (2.24) • 3x 3 y 9 z and - 3 T Z X 3 7 zy 3 a z z = J j (Pv.-) (2.25) 9x - 3y 9z ~ 9 i Considering p as constant, we get 9a 9x 3 T xx , yx . zx g:u • — + ->— = p — 3x 9y -,9 Z 9-t 9T 9a 3T yjs + . yy + .. yz = p • 9z 9t (2.26) 9- 9a zx , zy , .. zz 3.w 3 x ? y 9 z ' 2 . 2 . 2 . equations f o r Plane S t r a i n Non-Steady State Dynamic V i s i o p l a s t i c i t y  For the plane s t r a i n c o n d i t i o n T = T . = - — — 0. ^ zx zy 9t vso that the equation of motion ( 2 .26 ) can be w r i t t e n as 9c . 9T xx . . xy • 9 u + — K~TT~ = P 3 X 9 y y 91 . (2.27) 3T 3 a „ Vx , . yy 9 v J - + — = p — 9 x 3 y 91 . . 2.2.3. Determination of the Stress F i e l d from the S t r a i n F i e l d Proceeding i n the same way as in the case of 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, the equation (2.15) can be w r i t t e n as 3c 9 a ' I Y = X _ _2 c ,x- •- y • (2.15) 19 Now from the equation of motion (2.2 7) we get 3a x 3.T 3 X xy + dY p 3 t Substituting i n equation (2.15) gives 3a. .3x 3x 3y 3 u _ _ j _ [ £ x- £ y 3 t 3x (_ . J S i m i l a r l y from the equation of motion (2.2 7) ,8 a 3 T 3y , . 3 V 3x + pTt (2.28) (2.29) (2.30) Using equation (2.29) and (2.30) with the known values of o(x,y) at x=o and y=a, i . e . a(o,a) we get Gy(x,y) = ay(o,a) - / y ( 8 x x y 3 v> dy X " 3T x y . 3y 3 x I x- y + P 3u dx 3 t _ j y From equation (2.14) a = a + x y £ x-ey and from equation (2.li) 3s where A = — (2.31) (2.32) (2.33) (2.34) 2.3. 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 from ex-perimental data and thus the s t r a i n rates,*.'equivalent s t r a i n rates and f i n a l l y stresses can be determined. 20 FIG. 2.2 POSITIONS OF POINT P AT I N S T A N C E S T „ A N D T n + 1 21 Grid lines are marked on an end face of the specimen, which i s deformed at a predetermined speed. The deforming g r i d pattern i s photographed using a high speed camera. The g r i d points at a l l stages of deformation are d i g i t i z e d from enlarged photographs and the d i g i t a l p o s i t i o n a l data used as input to determine the instantaneous 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 Fig. (2.2), where a g r i d point formed by I row and J column has coordinates x , y ..at a p a r t i c u l a r instant of time t and grid-coordinates of x , , , y , , at the instant n ^ n+1 Jn+1 t , , . The instantaneous horizontal v e l o c i t y u and the v e r t i c a l n+1 J v e l o c i t y v i s then given by x , , -x u. . n+1 n : A X Since the above components of v e l o c i t y are obtained from the d i g i t i z e d coordinates of experimental g r i d points an e f f i c i -ent smoothing procedure i s required. The smoothing procedure mentioned by Shabaik (71) i s based on a simple averaging of the points. This procedure has caused d i f f i c u l t i e s i n the past in treating data which i s i l l - d e f i n e d and also tends . to be time consuming i n operation. Further for non-steady state conditions a reference g r i d i s needed, which can be fixed with respect to time (called a master g r i d ) . The simple averaging technique gives g r i d node positions that continually change with time. In order to surmount these d i f f i c u l t i e s a number of alternate methods were considered and f i n a l l y a three dimen-- t A t v. (2.35) 22 sional surface smoothing procedure was adopted, which treats x, y and u and also x, y and v as separate surfaces and f i t s a complete fourth order polynomial through the experimental points i . e . smoothing i s done in three dimensions to a surface formed from the scaler components of the vector f i e l d . The condition of continuity ( i . e . I * ^ *- or 3 u - 3 v for X ¥ 3 x 3 y plane strain) can also be imposed within the surface f i t t i n g procedure, thereby reducing the number of independent para-meters from 15 for each surface (a t o t a l of 30) to 10 for both surfaces. The smoothing procedure mentioned by Shabaik (71) does not take account of continuity and merely checks to see i f the error i s less than 15%. Besides f i t t i n g a smooth surface, the polynomials have the d i s t i n c t advantage that they can be readily d i f f e r e n t i a t e d for determining strains rates, and that the deformation can be automatically referred to a master reference g r i d . This means that strains and stresses can be determined for fi x e d points within the non-steady deformation zone. Also stresses can be evaluated d i r e c t l y at any position of x and y purely by sub-s t i t u t i n g the coordinates of a point (not necessarily a g r i d point) required, whereas the simple averaging technique requires an incremental evaluation of stresses over contiguous g r i d points u n t i l the required g r i d point i s reached. After c a l c u l a t i n g e x and e ^ using equation (2.1) , the e f f e c t i v e s t r a i n rates at a l l g r i d points for a l l instances of deformation can calculated from the equation (2.13). 23 In order to calculate T x^. (equation 2.33), (equation 2.31) and (equation 2.32) the value of A .is required. For non-workhardening materials the value of A A i s purely a function of e f f e c t i v e s t r a i n rate ( e ) as e f f e c t i v e stress ( a ) i s constant. For a workhardening material, the value of "5 must be obtained from experimental a vs e material data taken at the relevant s t r a i n rate conditions. I t i s normal to f i t a curve such as a = ^ e where c and n are material constants. Thus i f 7 i s known at a l l instances of time, the value of equivalent s t r a i n e at any deformation time t and hence a may be determined incrementally (assuming small intervals ©'f time) from the expression F . . = ' ft 7 dt (2.36) i D o e ± j where 7 ^ i s the e f f e c t i v e s t r a i n rate of a p a r t i c u l a r g r i d point as i t moves along i t s deformation path. 2.4. GENERALIZED EQUATIONS FOR PLANE-STRAIN VISIOPLASTICITY FOR COMPUTATIONAL PURPOSES The equations i n a form suitable for the development of the computer program are given below. The c a l c u l a t i o n of u and v i s done using equation (2.35) u. . _ Xn+1 - X n _ AX 1 3 -W 1 ^ A t y -Y v. . = ^n+1 n = Ay_ J t , . - t At n+1 n 24 Curve f i t t i n g of the v e l o c i t i e s i s done by a l i b r a r y sub-routine c a l l e d DLSQHS. This f i t s a complete fourth order polynomial i n x and y and determines the 15 parameters (con-stants) for u and v for the equations given below." The equa-tions used for computational purposes can be given as 2 3 4 u(x,y) = a^+a2x+a3x +a^x +a^x + a 6 Y + a 7 Y + a 8 Y + a 9 Y + a 1 0 X y 2 2 +a 1 1xy +a 1 2xy +a 1 3x y+a 1 4x y +a 1 5x y (2.37) S i m i l a r l y : 2 3 4 v(x,y) = b 1+b 2x+b 3x b 4 x +b5x 2 3 4 +b 6y+b yy +b gy +b gy +b 1 0xy+b 1 1xy 2+b 1 2xy 3 2 2 2 3 +b 1 3x y+b 1 4x y +t>15x y (2.38) Thus the v e l o c i t y component u anv v can be expressed • • • e e Y separately and the cal c u l a t i o n for s t r a i n rates x, y, ' xy and can be done as follows: "K = If = a 2 + 2 a 3 x + 3 a 4 x 2 + 4 a 5 x 3 + a 1 0 y +a 1 1y 2+a 1 2y 3+2a 1 3xy+2a 1 4xy 2 +3a 1 5x 2y (2.39) 25 '£y = If = V 2 b 7 y + 3 b 8 Y 2 + 4 b 9 y 3 + b 1 0 X 2 2 +2b 1 1xy+3b 1 2xy +k>13x +2b 1 4x 2y+b 1 5x 3 (2.40) and ^ = + |Z = a c+2a_y+3a Qy 2+4a Qy 3 'xy . 3y 3x . 6 7J S"7 9^  +a 1 Qx+2a 1 1xy+3a 1 2xy 2 2 2 3 + a 1 3 x + 2 a i 4 x y + a l 5 x +b0+2b0x+3b/.x2+4bnx3  + b 1 0 y + b l l y 2 + b 1 2 y 3 +2b 1 3xy+2b 1 4xy 2+3b 1 5x 2y (2.41) * * 3 u 3 v The condition of continuity ( e = - e or T — =- -— ) i s J X y 3x 3y imposed on the curve f i t t i n g by requiring the coefficients ;.to be related as follows. a 2 = " b6 2a 3= - b 1 Q 3a4= - b 1 3 • 4a 5= - b 1 5 a±1= - 3b 8 (2.42) 3 a I 5 =-2b 1 4 26 a 1 2 = " 4 b 9 a 1 0 = - 2 b 7 2 a 1 3 = " 2 b l l 2a 1 4= - 3b 1 2 The p a r t i a l derivative of -y with respect to x and y i s given by xy 3 * 2 Y x y = a 1 0+2a 1 1y+3a 1 2y +2a 1 3x+4a 1 4xy 3x 2 2 + 3a, cx +2b-+6b/1x+12b1-x 15 3 4 5 +2b 1 3y+2b 1 4y 2+6b 1 5xy (2.43) and hence 8Y xy T'iT a * = ( a 1 0 + 2 b 3 ) y + ( a l l + b 1 3 ) Y + ( 3 a 1 2 + 2 b 1 4 ) y + (2a 1 3+6b 4)xy. + ( 4a 1 4+6b 1 5) xy 2 + (3a 1 5+12b 5)x 2y (2.44) S i m i l a r l y 3 Y 2 ' xy = 2a 7+6a gy+12a gy +2a 1 1x+6a 1 2xy +2a 1 4x 2+b 1 0+2b i ; Ly+3b 1 2y 2 2 +2b 1 3x+4b 1 4xy+3b 1 5x (2.45) 27 and / _ i ] | y d x = ^ a ^ b ^ ) x + ( 6 a 8 + 2 b 1 1 ) x y 2 ( 2 a 1 4 + 3 b 1 5 } x 3  + ( a l l + b 1 3 } x + 3 + (12a g+3b 1 2) y x + (3a 1 2+2b 1 4)x-y (2.46) Further equation (2.36) can be expressed as e i+ .. »• ( £ i + £ i+1) A t -4- £. l - 2 K.. (2.47) The normal stress -<?y: can be calculated using equation (2.31) , i . e . Y - y(o,a)- y/ •. 3 Txy + p J _ v X/ a v a x a t ; y o xy. < 3 r x-"• y -, ..' 3 u 3y 3 x ^ ^ 3 t dx y = a (2.31) The terms for t h i s equation may be calculated separately. The f i r s t termo" y (o,a) i s determined experimentally at each inter-val of time. The second term can be calculated from the equa-tio n . Y T ,; = xy xy 2A (2.33) so that 3 T xy = = 1 3 3 3 x ( Y x y ^ 7 ^ ( _ 2 T ) - Y —n xy c e 3x Y xy 3 £ . a 28 = £ —a_ 3 9X Y _n xy fe c 3 c 3 Y gX^ ' 3Y r xy 3>^ 1 e _ Y xy 3S, I ax £ Y - n-1 „ •-n xy e 3 . 2 e n 8 Y „ „ MHY 3 X 7 + T .3x 3x S i m i l a r l y (2.48) 3 T XV _ c 'V " 3 n-Y — n-1 — _ n 3 Y /—<n Y /— \ xy F. 3 e + _s_ ^xy _(;£) xy 3 ^ e ^ sy ( i ) 2 3 Y . X y (2.49) Using equation (2.1), the t h i r d term of equation (2.31) can be computed as 3 ( £ x ^ <y) "^X A £ X £ y 3 A ,3 3 X + -r A 3 x • V b x t-yJ-( £ x £ y) 3 • (/J 7 ) . 1 E " e )TI .2 3 x - n + T- L - - * A 2Ce A 3 x x (2.50) But 3X 2c E = _ ! r t~ n e ( e ) 2c L t 2c gx -n-1 [ 7 ( 7 ) n] } 3e + ( e ) n_3_ ( e ) ] 3 x 3 x 3 J? 1 " 2c" n (f^T 3E + 3 (7 ) 3 x 2c ( i ) n ^  3x n A 3s + _A 3s_ - 3 X - ax •e E; A[_ IL_ 37 + 1 3s - 3 X - 3x e s :J . (2.51) d i f f e r e n t i a t i n g equation (2.13) with respect to x, we get 2 3s 3 x 3 _ r 2 ( 3 s l 2 3X x 4 x y = 1 [ (3 e 2 i 3 ; 2 ) 2 (6 e v___ ex 3 x 2y 5 yxy) 3 x 4 x y ' e 8x 4 x y ' 9x -h = (3 e 2 + 3 i 2 ) H e ^ x + 1 i _!_xy7T V x 4 x y / |_ x — 2 x y. 5 - T i r J (2.52) So that the t h i r d term of equation (2.31) , using equations (2.50-2.52) can be computed as 2 £ ~ £ ~K (e„ - e..){ — -g^ + - — > + — ^ -X y £ d d (2.53) Thus the normal stresses a and c?x can be computed using equations (2.31, 2.32, 2.48, 2.49, 2.50, 2.51, 2.52 and 2.53) with the known values of a(o,a). "2.5 COMPUTER PROGRAM The flow chart of the computer program developed for plane-s t r a i n dynamic and qu 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 Fig. (2.4) and a l i s t i n g of the program i s given i n Appendix I. The running instructions and the main steps i n the program are described below: (1) Input a l l data required for the calculations (a) Read IX = No. of experimental g r i d l i n e s p a r a l l e l to y axis IY = No. of experimental g r i d l i n e s p a r a l l e l to x axix The above are given in FORMAT (312) 30 IT = No. of the time steps DT = Time i n t e r v a l between two consecutive photographs (this need not be constant time interval) FORMAT (8F10.0) NIX = No. of gr i d lines p a r a l l e l to y axis i n master g r i d NIY = No. of g r i d l i n e s p a r a l l e l to x axis i n master gri d Both i n FORMAT (312) CA = Constant a used for a (o,a) CC = Constant used i n true stress and true s t r a i n r e l a t i o n CN = Constant (or index) for workhardening SIGOA = a (o, a), known value of y at x = o and y = a A l l above i n FORMAT (8F10.0) (b) Read, the instantaneous coordinates of the g r i d points of consecutive photograph. The format i s given as FORMAT (5X, 2F6.3, IX, 2F6.3, IX, 2F6.3, IX, 2F6.3, IX, 2F6.3) ) Calculate the values of horizontal 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 at each g r i d point using equation (2.35) ) F i t a 4th degree polynomial through the three dimensional curves for u and v as a function of x and y using the computer l i b r a r y routine c a l l e d 'DLSQHS'. This program provides a least square f i t to a l i n e a r func-ti o n of M parameters ( i . e . M independent variables) and N data points by the Householder transformation 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 IX 7 8 FIG. 2.3. GRID . LINES ON THE S P E C I M E N 32 FLOW CHART OF COMPUTER PROGRAM Read IX, IY, IT, DT, NIX, NIY, CA, CC, CN and (o,a) T Read (x,y), the g r i d points at each time i n t e r v a l . P l o t them Ca 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 U i j 3 n+1 n t - t n+1 n . V « = y -y (equ. 2 . 3 5 ) n+1 n t - t n+1 n E F i t the curves to the u and v values with complete fourth degree polynomial i n x and y (Subroutine DLSQHS) Calc u l a t e and store s t r a i n -r ates f o r a l l points, using subroutine DERIV x~ 3X y ay (equ. 2 . 1 , 2 . 1 3 , 2 . 3 8 , 2 . 3 9 , 2.40) Y J~ + 1 - v - t r 3 U xy ~sy o » V ax xy E C a l c u l a t e and store the t o t a l e f f e c t i v e s t r a i n at a l l s p e c i f i e d node points ( 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 c a l c u l a t i o n s ( i i ) F i l l zer.o f o r the points, outside the boundary of the deforming zone by . use of subroutine FILL X P l o t i (a) 00 u = f ( x , y ) , .v = f ( x , y ) , £, x = f(x,y), 6 y = f ( x , y ) , t = f(x,y), i n two dimensional form using subroutine PLOT, PLOTAND & SYMBOL i n three dimensional form using subroutine PERS Calculate 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 1 Plot:• ( i ) <7~ as a fun c t i o n of x and y y ( i i ) (J"^ as a fun c t i o n of x and (iii)T^* as a function of x xy and y (a) in two dimensional form using subroutine PLOT, PLOTAND & SYMBOL (b) in three dimensional form using subroutine PERS STOP FIG. 2.4 33 That i s i t minimises N M 2 (v. - a.x..) E 1 E 1 i = l j=l DLSQHS transforms the matrix X to an upper triangular form v i a Householder transformations, and then solves the system by backward substitution. I f the command REFINE i s defined by 'TRUE': a correction vector i s computed from the residual errors between the dependent variables and the f i t t e d values. Correction vectors are then applied to the solution and re-computed i t e r a t i v e l y u n t i l convergence i s obtained. DLSQHS i s most e f f e c t i v e for problems, where the correction of the matrices i s unknown and the scale of d i f f e r e n t variables varies widely. (4) impose the"condition of continuity i . e . e = - e or - x y au^ _§v • • -- - :- 1"- - ' ; (5) Calculate the strain-rates £ "/ "£ , shear stra i n - r a t e xy and ef f e c t i v e strain-rate using equations (2 . 3 9 ) , ( 2 . 4 0 ) , (2.41) and (2.73) ,are combined into a 1DERIV' ( d e t a i l s are given in 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 integrating along the path of p a r t i c u l a r p a r t i c l e or g r i d node (refer equa-tions 2.36 and 2 . 4 7 ) . (7) Calculate shear stress Txy using equations ( 2 . 33 , 2.34 and 2 . 4 1 ) . 34 (8) Calculate independently the terras from equations (2.31), (2.48), (2.49), (2.50), (2.51), (2.52), and (2.53) and thereby calculate a y using the computer l i b r a r y sub-routine 'DQUANK'. This subroutine integrates a function "3D f(x) when the l i m i t s a and b are given. i . e . I = ! j f(x) dx. I t i s b a s i c a l l y based on Simpon's three points i n -tegration & improved by using an adjustment term of f i f t h degree i n place of the three degree term. The absolute error can be limited to any a r b i t r a r i l y s p e c i f i e d value. Calculate stresses ° x using equations (2.31) and (2.31). (9) Plot the d i f f e r e n t quantities, u, v, ex, e y , Yxy, F , 7f C T y,° x, 1 xy, etc. for d i f f e r e n t x, y values. The p l o t t i n g vectors are taken at defined master g r i d node points. Any point within the master g r i d but outside the deformation zone are given zero values by the subroutine 'FILL'. (a) These values can be plotted i n two dimensions either as a function of x or a function of y. For t h i s purpose the subprogram 'PLOT' can be used. Thus subprogram i s the basic p l o t subroutine. I t generates the pen movement required to move the pen i n a straight l i n e from i t s pre-sent position to the position indicated i n the c a l l . I t i s also used to relocate the o r i g i n of the p l o t t e r coor-dinate, system i n the X di r e c t i o n . To ensure that p l o t t i n g i s complete, a second subprogram 'PLOTND' i s used. For 35 clear 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. This plots alphanumeric characters and spe c i a l symbols. (10) Three dimensional orthographic displays can be made using the subroutine 'PERS'. The above values are taken as 'Z' values and are plotted i n three dimensional form as a function of x, and y. 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 investigation of forging operations such as heading and upsetting, where st r a i n - r a t e , forming-speed and forming load are important, a device with special requirements i s needed. Obviously as wide a range of impact speed as possible i s required together With a v e l o c i t y p r o f i l e which i s sensibly i n -dependent of forming load and which may be adjusted to s u i t circumstances. The commercial alternatives that are available have been developed for s p e c i f i c applications and of necessity have a lim i t e d f l e x i b i l i t y . For example, the crank press, used i n many forging operations has a variable stroke and maximum vel o c i t y and has a well controlled v e l o c i t y p r o f i l e , which i s sensibly independent of forming load. However, the maximum vel o c i t y i s l i m i t e d to approximately 15 to 20 f t . / s e c . With the drop hammer and high v e l o c i t y forming machines (such as the petroforge) a higher v e l o c i t y can be reached (15-30 ft. / s e c . for the drophammer, 90-100 ft . / s e c . for the petroforge) but the stroke and v e l o c i t y p r o f i l e during deformation are determined purely by the resistance of the workpiece. That i s the high SOLENOID. FIG 3.1 VARIABLE S P E E D , C O N T R O L L E D VELOCITY PROFILE, S INGLE C Y C L E IMPACTING P R E S S 39 v e l o c i t y results i n high energy and with low strength targets very l i t t l e of the energy i s absorbed i n p l a s t i c work. Further the high v e l o c i t y devices cannot be used for low v e l o c i t y work. With these design c r i t e r i a in:mind a variable speed, con-t r o l l e d v e l o c i t y p r o f i l e , single cycle, impacting press was designed and b u i l t within the department. 3.1.2 Description of the Press: A diagram of the impacting press i s shown i n Fig . 3.1. The drive comprises a modified Whitworth quick-return mechanism consisting of a crank and a drive arm together with a variable speed D.C. motor, a flywheel, bearings etc. The end of the drive arm i s attached by a connecting rod to a c y c l o i d a l cam. In single cycle operation, the cam i s made to engage with an upper platen (or ram) which impacts the workpiece. The upper platen and cam are both mounted on multirod supports with l i n e a r b a l l bushings. A brake i s provided on the flywheel for emergency purposes. 3.1.3 Operation of the Press: The drive wheel i s rotated at a p a r t i c u l a r speed by ad-justment of the D.C. motor c o n t r o l l e r causing the drive arm to o s c i l l a t e about i t s lower pivot. The single cycle tripping.? mechanism connects the drive arm with the cam and the cam en-gages with the cam followers on the upper-platen. The upper-platen i s thus forced down towards the workpiece. The platen achieves a maximum v e l o c i t y when the drive arm i s i n a central position a f t e r which time the platen i s brought to rest. On the return stroke of the cam the platen i s returned to i t s 40 i n i t i a l p osition. The t r i p p i n g mechanism then disengages the connecting rod from the cam and the drive arm continues to o s c i l l a t e f r e e l y about i t s lower pivot. The stroke, v e l o c i t y and acceleration p r o f i l e of the upper-platen are determined solely by the cam contour and the speed setting of the D.C. motor. A c y c l o i d a l cam i s used for high v e l o c i t y work to minimise excessive wear, through shock and vibrat i o n . The lower platen height can be adjusted r e l a t i v e to the upper platen thereby determining the working portion of the stroke. 3.2 KINEMATIC ANALYSIS OF THE MECHANISM 3.2.1. Derivation of Expressions for Velocity and Acceleration: The drive mechanism i s shown schematically i n Fig. (3.2.). The centre l i n e of the drive arm i s indicated by l i n e AB, where point A denotes the fixed lower pivot. The path of B i s i n -dicated by the dotted l i n e with B 1 and B", showing the extreme points. The point C i s the centre of rotation of the drivewheel, a distance P above the fixed pivot of the drive arm. The eccen-t r i c i t y of point E, the cam follower, about C i s given by a distance Q. The l i n e CD i s a reference for angle e. The length of AB i s L. The path of point F denoting the cam po s i t i o n , i s shown by dotted l i n e . The length of connecting rod, BF = R and the distance of cam position F from l i n e AC (or AH) be X. The angle made by AB with l i n e AC i s 4> 42 Then from geometry = arc tan QCosQ P+QSinQ (3.1) D i f f e r e n t i a t e <\> with respect to t, we1 get ,= _ PQSinQ+Q2 d& d t '"~ " P2+2PQSin9+q2 ^ (2.3) D i f f e r e n t i a t i n g again, we get dt ( d t ) "> : (P 2+2PQSine+Q 2) 2 (3.3) Now from F i g . (3.2) HF = HG + GF = L Sin cb + GF = X but GF = [(BP' BG2)1 h = [R2 -{• L(l-Cos <f> ) } 2 ] h O O 1 Hence, X = L Sin <j) + [R - ' { L (1-Cos cb) } ] ^  D i f f e r e n t i a t i n g with respect to t, we get dx ± L d_j dt dt Cos* - L(l-Cos cb ) s i n * [ R 2 4 L ( 1-Cos c b ) } 2 ; ^ ^ = L dj> [Cos cb - W Sin <b ] dt V (3.4) where W = L (1-Cos cj> ) and V = [(R 2 - W 2)]^ r D i f f e r e n t i a t i n g again and rearranging the terms, we get d 2x = L d 2 cb (Cos <f> + W) 2 , 2 v dt dt' - L ,d<b. [ ( S i n <b + L Sin tb  cdt ; V (W + 1) + W Cos <b ) ] 2 V V (3.5) 43 3.2.2. Maximum Velocity and Acceleration of the Mechanism: The press has following dimensions:-P = 14 inch. Q = 6 inch. L = 24 inch. R = 26 inch. Ratio of driver to driven, pulley dia. = 0.573 For a flywheel speed of 250 RPM, ( i . e . a motor speed of 434.4 2 x 250 rpm) the angular v e l o c i t y of the flywheel = . 6 Q ,—r = 15 rad/sec. (3.6) The po 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 <j> = o or 0 = - 90' Hence from equation (3.2) we get d£ = _•• 'd9; PQ Sin e + Q2 dt dt p 2 + 2 PQsine+Q 2 = (- a) (14x6) (-1)+(6x6) (14x14) +(2x14x6 (-1) ) + (6x6) = 0. 75 a) (omega) From equation (3.4) dx _ d (j) (Cos (j) - W Sin cf> ) ; here cj> = o dt L d t V _ _ dA (1 - W x o ) -J Td(j) " L d t v " L d t 11 I t = 2 x .75 u = 1.5 u ( f t . / s e c . ) ' - ( 3 * 7 12 dt 44 45 Now the displacement of platen, y, i n r e l a t i o n to the c y c l o i d a l cam (refer. F i g . 3.3) i s given by h TT£_ _ 1 Sin 2ja_r (3.8) ^ — IT L ' 2 - L 1 if where I = distance t r a v e l l e d along x axis at a p a r t i c u l a r time L 1 = length of c y c l o i d a l cam i n x di r e c t i o n h = maximum distance t r a v e l l e d by upper-platen or follower i n the y d i r e c t i o n Velocity of follower i . e . platen i s given as (£4.) h v d t ; , ., COS2TT£ 1 (3.9) v = -jTT-f 1 - — - ) v = dt (3.10) max T , Si m i l a r l y the equation for acceleration of the upper-platen, a, i s given by (£-*) 2 3_ . . 2h TT dt ^ i . . . 2 v.I , h Cos2 -a d Z a- v = n S i n + — (1- ) —x- . . y 2 j i L 1 T I 2 (3.11) L ' L L dt For the present cam, h = 4.0 inch, and L ' = 8", so that vmax ,d£. , dx. 2x4.0x ldt ; _ 2x4.0x (dt' 8 8 2x4. 0x1.5 ui 8 = 1.5 a) Using the value of u =15 rad./sec. from equation (3.6) v max = 22.5 ft. / s e c . PUSH-BUTTON R E D -/Fuse Power METAL-CAS NO I MICRO-SWITCH • 9 3 500 P<wer i n RED 68 k 2 0u ORANGE JLACK o f f " i I L l l a t c h ORANGE ORANGE T 3 4 1 5 / 3 0 0 V JT n e t BLUE (LACK DC I IRL f S o l e n o i d i P o w e r k f f " E f t BRN l 9 9 l l 2 IRL, TTT >( BLACK 6 RL3 ' R e s e t ' 10 C BRN *ED 20 u 500V ^ C 4 £ l 5 0 6 7 t BRN BLUE 117V, 60A. BLUE + 170V-p Camer P o w e r 5K t From cam. C o n t r o l o I HJETAL-CHASIS ' 0 6 6 6 <£kND CAMERA L N 1 SYNCHRN. 110V AC(OUT) " j I - -L- + SOLENOID L_o O r - p — U r j — ° I 0 9 s J-100 u -M u 10V . L I, 1 T T To CAMERA SYNCHRONIZATION 0 0 1 I _o U N I T I— FIG. 3.4 ELECTRICAL CIRCUIT FOR SINGLE - CYCLE OPERATION CONTROL 47 3.3 CONTROL With single cycle operation i t i s es s e n t i a l that engagement of the drive arm with the cam occur at the correct point i n the cycle and have s u f f i c i e n t time to engage properly. This i s achieved by requiring that a sensing device placed at the extreme point of the drive arm motion be actuated i n conjunction with a push-button s t a r t switch before engagement occurs. A solenoid then retracts and allows the maximum time for the two parts to engage. Synchronization of the high speed camera with the specimen deformation i s also accomplished . using the signal from the remote sensing switch. A time delay device i s used to vary the s t a r t of filming so that d i f f e r e n t impact speeds can be accommodated. A further sensing device placed at the opposite extreme of the drive arm movement indicates when the event i s completed and triggers the magnetic camera brake. Details of the elec t r o n i c c i r c u i t r y for engagement and disengagement of the drive arm and for sychronization of the high speed camera are given i n the following sections. 3.3.1. Control for Single Cycle Operation: The e l e c t r i c a l c i r c u i t used for c o n t r o l l i n g the engagement of the driv i n g arm of the cam i s given i n Fig . (3.4) and the flow diagram of the relay sequence i s given i n Fig. (3.5). The sequence of events i s as follows: (i) Single cycle s t a r t switch (push button) i s pressed which sets the latch relay (RL2) through reset relay (RL3), faormally OPTICAL-SENSORS 2 1 AC JL HIGH SPEED CAMERA LOGIC HIGHSPEED CAMERA ^ SYNC. < CONTOL HIGH SPEED CAMERA RELAY ON RATCHET OFF TRIGGER ^ R L 1 PULSE POWER ON RESET 117 V AC POWER SUPPLY POWER ON PULSE FOF.MER 170V DC SOLENOID i AC POWER RELAY START PUSH BUTTON AC L4 LATCH RELAY RL 2 RESET RELAY LATCH CONTROL RL 3 PULS FORT E VtER i MICRO SWITCH ( A ) FIG. 3.5 FLOW D I A G R A M FOR R E L A Y S E Q U E N C E FOR SINGLE C Y C L E C O N T R O L OPTICAL SENSOR 1 CAMERA RELAY | SINGLE | CYCLE OPTICAL SAMPLE SENSOR 2 PULSE 1ms LAMP DRIVER CONTROL UNIT CAMERA RELAY LAMP INDICATOR FIG.3.6 CIRCUIT OPERATION BLOCK DIAGRAM FOR HIGH SPEED CAMERA SYNC. 50 closed contact p o s i t i o n ) . ( i i i ) The ratchet relay steps to 'ON' as the drive arm contacts the micro switch (A) at the bottom of the stroke. This puts a signal on the solenoid relay (RL4) through the la t c h relay (RL2), and energizes the solenoid and hence the t r i p p l i n g mechan-ism. (iv) When the drive arm contacts the micro switch (B) for the second time the ratchet relay steps to the 'OFF' position. The reset relay (RL) momentarily energizes and the la t c h relay (RL2) deenergizes v i a RL3. 3.3.2 Control for Synchronizing the High Speed Camera: The block diagram of the sychronization control of high speed camera i s given i n Fig. (3.6). I t consists of an o p t i c a l sensor microswitch (A) to i n i t i a t e a filming signal and an integrated c i r c u i t timers type 556 to give a variable time delay to the actuation of the camera relay and hence the camera. Reference pulses of 1 m.sec duration are compared with the pulses received from the o p t i c a l sensor no. 2. (in Fig. (3.6)) and when these correspond, a lamp indicator i s triggered showing correct synchronization. CHAPTER FOUR EXPERIMENTAL PROCEDURE AND DISCUSSIONS 4.1. E XPE RIMEN TAL P ROCE DURE A plane s t r a i n upset compression-test was made a t an im-pact speed 15.7 f t . / s e c . The p l a n e - s t r a i n specimen was made from p l a s t i c i n e u s i n g a metal mould. The specimen measuring l " x l - l / 2 " x 2 " was c u t with a f i n e wire and a l u b r i c a n t o f s i l i -cone grease was used to prevent the specimen from s t i c k i n g i n the mould. The end s u r f a c e o f the specimen was sprayed with b l a c k p a i n t and a square g r i d p a t t e r n was made by s p r a y i n g white p a i n t on the bl a c k s u r f a c e through a wire mesh g r i d (14 mesh, i n c h . ) . The use of bl a c k and white p a i n t gave good c o n t r a s t f o r h i g h speed photography. P l a n e - s t r a i n 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 the specimen on the lowest p l a t e n o f the hig h speed impacting p r e s s between two p a r a l l e l l u b r i c a t e d p e x i g l a s p l a t e s ( c a l l e d a Kudo appara t u s ) . The upper p l a t e n o f the ram was a d j u s t e d t o impact the specimen at the maximum v e l o c i t y i n the c y c l e and a l l con-t r o l s and f i l m i n g s y n c h r o n i z a t i o n (as d i s c u s s e d i n Chapter 3) a p p r o p r i a t e l y s e t . The h e i g h t of the camera was kept, such t h a t the o b j e c t i v e l e n s o f the camera was on the same h e i g h t as t h a t of the specimen and the plane of specimen was p a r a l l e l to the plane o f the l e n s . ( s e e F i g . A ) 52 FIG. (A) LOCATION OF THE SPECIMEN IN KUDO APPARATUS 53 4.2. DISCUSSION The g r i d points, for the dynamic plane s t r a i n upset com-pression test, were d i g i t i z e d at 0.00133 second time i n t e r v a l s (every 4th frame of the high speed f i l m ) . Also results from a q u a s i - s t a t i c plane-strain compression te s t done by Shabaik (71) (see Fig. 15) at a speed of 0.02 ft./min. were d i g i t i z e d . The d i g i t i z e d g r i d points are plotted for the four steps of qu a s i - s t a t i c deformation i n Fig. 4.1 and 4.2 and i n Figs. 4.3 to 4.9 for the seven steps of dynamic deformation. The movement of . certain grid-node points during deformation i s plotted i n Figs. 4.10 and 4.11. The smoothed horizontal 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) are given as a function of x and y in Figs. 4.12 to 4.15 for l a s t time i n t e r v a l in both s t a t i c and dynamic tests. Plots of e f f e c t i v e s t r a i n rate ( e~ ) are given in Figs. 4.16 and 4.17, while t o t a l e f f e c t i v e s t r a i n accumulated incre-mentally for a l l time in t e r v a l s for the 0.02 ft./min. and for 15.7 ft . / s e c . deformation speeds are shown i n Figs. 4.18 and 4.19. The normal and shear stresses ( 0 y , 0 x, T xy) are plotted as functions of x and y for the 0.02 ft./min. deformation speed in Figs. 4.20, 4.21 and 4.22,and for the 15.7 f t . / s e c . deforma-tion speed i n Figs. 4,23, 4.24 and 4.25. The values of a (o,a) for 0.02 ft./min. deformation speed was taken from reference (71) as 31000 p s i at the o r i g i n with the relationship between _ — o 25 e f f e c t i v e stress and s t r a i n as a = 31000 e " . The values of a (o,a„) for the p l a s t i c i n e specimen was taken as 20 p s i at the 54 o r i g i n from reference (79). The stresses ° x and-T'cxy are also plotted as a function of x for the CL 02 ft./min. deformation speed i n Figs. 4.26 and 4.27 and for 15.7 ft . / s e c . deformation speed i n Figs. 4.28 and 4.29. The results from the specimen deformed at 0.02 ft./min. agree clo s e l y with well documented results for 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 of normal ( a y ) i n t e r -face stress d i s t r i b u t i o n with maximum stress occurring at the central portion of the specimen. The normal stress d i s t r i b u t i o n for the specimen deformed at 15.7 ft/sec i s r a d i c a l l y d i f f e r e n t , showing a saddle type d i s t r i b u t i o n of normal stress, with the maximum stress occurring near the periphery of the contact zone. The interface shear stress d i s t r i b u t i o n s also change form with s t r a i n rate. The dramatic change of normal stress d i s t r i b u t i o n with s t r a i n rate i s t o t a l l y at variance with currently held views and furthermore i t occurs at quite moderate v e l o c i t i e s which are c e r t a i n l y well within the range of those encountered i n many metal forming operations. 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. > ) in ) > L n ) ) ) <=> N c \ ' <s V. 7 ' \ \ 7 ^ \ w 7 /• \ ' \' : — 7 . ' ^' f f t c f f \ 7 ' \ k 7 S \ R 7 /• \ k—7 < 7* f V * 7 > k 7 C > k ) f ^ V. 7 f \ ^ 7 ' \ \ -7 7 X < ) r \ r ' ^ 7 / V 7 ' \ \ 7 Y S s. 7 / } \ ) { \ \ 7 v 7 v 7 v- ' / \ \ 7 • N K 7 f \ K—~7 f >, f . \ \ —7 / \ s. ) / \ ^ 7 { \ \ 7 / \ < ) * V c > \ 7 ' -s \ -7 ' "s \ 7 f S k 7 k 7 i- i V / L > ^ / £ -A \ / \. / f \ \ 7 f 0 i. \ c > \ J / N / •> ; ^ 7 /_ > ^ 7 \ 7 /• > \ 7 K ' ( ^( > s > ; C— " v J f ^ 7 ^ c. > ) in a — j ~ ~ m a ( 7 / \ < i / 2 ' \ 7 \ / / > ^ V / > \ 3 «C 7\ -) 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 i n o 0.0 0.5 1.0 1.5 2.0 X 2.5 3.0 FIG. 4.6 DISTORTION OF GRID LINES DURING DEFOMATION FOR 1 5 . 7 FT./SEC. DEFORMATION SPEED. 61 62 63 FIG.4.9 DISTORTION OF GRID LINES DURING DEFOLIATION FOR 15.7 FT./SEC. DEFORMATION SPEED. 64 68 p r G 4 . 1 4 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 T H R E E 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 OFf7^ AS A FUNCTION OF X FOR 1 5 . 7 FT./SEC. DEFORMATION SPEED. 4.3. SOURCES OF ERROR 84 The quality of the g r i d node p o s i t i o n a l data used as input to the program i s the most important single item a f f e c t i n g the re s u l t s . Sparse or poorly d i g i t i z e d data i s l i k e l y to cause inconsistencies i n the output stress d i s t r i b u t i o n s . The surface f i t t i n g routines w i l l smooth certa i n i r r e g u l a r i t i e s but there i s a l i m i t to their c a p a b i l i t i e s . The v e l o c i t y of the upper plater varies during the cycle, and during specimen deformation, according to the equations 3.4 and 3.9. A further f l u c t u a t i o n of plateja v e l o c i t y may occur as the drivewheel speed changes during the cycle. I n i t i a l energy balance calculations show this change i s l i k e l y to be very small p a r t i c u l a r l y with low strength p r o j e c t i l e s . A prio r knowledge of the actual upper-platen v e l o c i t y p r o f i l e during deformation i s not required as this i s obtained auto-matically from the d i g i t i s e d displacement data and a knowledge of the time increment between frames of the high speed photo-graphs. Plane-strain deformation was achieved using a Kudo apparatus. While this assured plane-strain condition i t did introduce a f r i c t i o n a l drag on the end faces of the specimen. The e f f e c t was minimized using s i l i c o n grease as a lubricant and from examination of deformed specimen i t was concluded the e f f e c t was not important. In the analysis the material was assumed to be strain-rate 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 metal-forming m a t e r i a l s . I t i s p o s s i b l e to r e l a t e e f f e c t i v e s t r e s s , a, to both 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 to the 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 could 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 h a v e b e e n d e v e l o p e d f o r d y n a m i c a nd q u a s i - s t a t i c , s t e a d y o r n o n - s t e a d y d e f o r m a t i o n p r o c e s s e s . The e f f e c t o f i m p a c t v e l o c i t y o n t h e m e c h a n i s m o f d e f o r m a -t i o n d u r i n g d i f f e r e n t m e t a l - w o r k i n g p r o c e s s e s c a n b e s t u d i e d u s i n g t h i s w o r k . I t i s c l e a r f r o m t h e i n i t i a l s t u d y o f u p s e t t i n g , t h a t s t r a i n and s t r e s s d i s t r i b u t i o n v a r y s i g n i f i c a n t l y w i t h s t r a i n - r a t e . 6. SUGGESTIONS FOR FURTHER WORK 87 The method developed enables the stress d i s t r i b u t i o n s to be determined i n many dynamic metal-forming operations. A st a r t i n g point for thi s i s to determine the change of stress 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 density and surface geometry for plane-strain upsetting operations. Modifications can be made to the surface f i t t i n g routines to accommodate the constraints that the v e l o c i t y gradient i s zero along the y axis, and that the v e r t i c a l component of ve l o c i t y , v, i s equal to the platen v e l o c i t y for points on the platen-workpiece i n t e r f a c e . I t i s l i k e l y that a 5th order polynomial would then be needed for surface f i t t i n g . 88 R E F E R E N C E S References 89 1. Hoffmann, 0 . , and Sachs, G.,"Introduction to the Theory of P l a s t i c i t y f o r Engineers," McGraw H i l l , New York, 1953. 2. S i e b e l , E., "Die Formgebung i n bildsmen Zustande", Verlang 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 Gleichgewichts i n p l a s t i s c h e n Korpern", 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 Satz uber das Gleichgewicht", A. Angew. Math. Mech., _3_, 401-406, 1 9 2 3 . 5. Caratheodory, C , and Schmidt, E., "Uber die Hencky-Prandtlschen Kurven", Z. Angew. Math. Mech., 468-4 7 5 , 1 9 2 3 . 6. G e i r i n g e r , H., "Beitrang 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," Proc. T h i r d I n t e r n a t i o n a l Congress of A p p l i e d Mechanics, 1 . , 185-190, 1 9 3 0 . 7. H i l l , R., The Mathematical Theory of P l a s t i c i t y Clarendon Press, Oxford, 1950. 8. Prager, W. and Hodge, P.G. J r . , Theory of P e r f e c t l y  P l a s t i c S o l i d s . Wiley, New York, 1951 . 9. Lee, E.H., "The T h e o r e t i c a l A n a l y s i s of Metal Forming Problems i n Plane S t r a i n " , J . A p p l i e d Mechanics, 19 , 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 Mechanical  Engineers, D. Van Nustrand Co., L t d . , 1962. 11. Bishop, J.F., "On the E f f e c t of F r i c t i o n on Compression and In d e n t a t i o n Between F l a t Dies", 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 Plane 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 P l a n e - S t r a i n E x t r u s i o n 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 f o r C o n s t r u c t i n g A n i s o t r o p i c S l i p - L i n e F i e l d " , I n t . J . Mech. Sc., 1974, pp. 241-248 90 15. Rice, J.R.,-"Plane S t r a i n S l i p - L i n e Theory for Anisotropic R i g i d / P l a s t i c Materials", J . Mech. Phys. Solids, V21 N2, March, 1973, 63-74. 16. Booker, J.R. and Davis, E.I., "A General Treatment of P l a s t i c Anisotropy under conditions of Plane S t r a i n " , J. Mech. Phys. Solids, Vol. 20, 1972, pp. 239-250. 17. Shabaik, A., "Effect of 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 University, Canada, 1973, PP. 221-238. 18. Wilson, W.R.D., "Slip-Line Solutions 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 Conditions", Proc. 5 t h . NAMRC, SME, 1977, PP. 80-86. 19. Ewing, D.J.F., "A Series Method f o r Constructing P l a s t i c F i e l d s , " J . Mech. Phys. Solids, Vol. 15, 196?, pp. 105-114. 20. C o l l i n s , I.F., "The Algebraic-Geometry of Slip-Line F i e l d s with Applications to Boundary Value Problems", Proc. Roy. S o c , Ser. A, Vol. 303, 1968, pp. 317-338. 21. Drucker, D.C, Prager, W. and Greenberg, H.J., "Extended Limit Design Theorems for Continuous Media", Quart. Appl. Math., 2, 1952, pp. 381-389. 22. Avitzer, B., Metal Forming; Process and Analysis, McGraw H i l l , New York, 1968. 23. Johnson, W., "Estimate of Upper Bound Loads for Extrusion and Coining Operations", P r o c Inst. Mech. Engrs. (London), 173, pp. 61-72, 1957. 24. Kudo, H. , "An Upper Bound Approach to Plane S t r a i n Forging and Extrusion - I", Int. J . Mech. Science, 229-252, i 9 6 0 . 25. Kudo, H., "An Upper Bound Approach to Plane S t r a i n Forging and Extrusion - I I " , Int. J. Mech. S c i . , 1: 229-252, i 9 6 0 . 26. Kudo, H., "An Upper Bound Approach to Plane S t r a i n and Extrusion - I I I " , Int. J. Mech. S c i . , 1: 366-368, i 9 6 0 . 27. Kobayashi, S., "Upper Bound Solution of Axisymmetric Forming Problems - I", Presented at the Production Engineering Conference of ASMS, May, 1963. 28. Kobayashi, S., "Upper Bound Solutions of Axisymmetric Forming Problems I I " , Trans. ASME, Series 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 in Upset Forging of Rings", ASME Paper No. 77-WA/Prod.-35. 30. Lahoti, G.D. and Altan, T., "Prediction of Temperature Distributions in Tube Extrusion using a Velocity Field 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 in 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. Collins, I.F., "The Upper Bound Theorem for Rigid/Plastic Solids Generalized to Include Coloumb Friction", 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. 77-WA/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., "Plastic Deformation of Polygonal Disks under High Velocity Impact", Presented in 8 t h Conference of U.S. National Congress of Applied Mechanics, June, 1978. 37. Juneja, B.L. and Prakash, R., "An Analysis for 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 Arbitrarily Shaped Dies using Flow Function", J. Engg. Ind., Vol. 99, 1977, PP. 544-548. 39. Nagpal, N., "General Kinematically Admissible Velocity Field for some Axisymmetric Metal Forming Problems, "J. Engg. Ind., Vol. 96, 1974. 40. Zienkiewicz, O.C., The Finite Element Method, Jrd Edition, McGraw-Hill, 1977. 92 41. Lee, CH. and Kobayashi, S., "New Solutions to R i g i d - P l a s t i c Deformation Problems using a Matrix Method," Trans. ASME, J. of Engrg. for Ind. Vol. 95, 1973. pp. 865-873. 42. Godbole, P.H. and Zienkiewicz, O.G., "A Penalty Function Approach to Problems of P l a s t i c Flow of Metals with Large Surface Deformation", J. S t r a i n Analysis, Vol. 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 Engineers, John Wiley & Sons, 1975. 44. Strang, G. and Fix, G.J., An Analysis of the F i n i t e Element  Method, Prentice-Hall, Englewood C l i f f s , N.J., 1973. 45. Rowe, G.W., "Recent Developments i n the Theory and Practice of Metal Forming", Proc. Third North American Metal-Working Research Conference, Carnegie-Melton University, 1975. pp. 2-25. 46. Alexander, J.M. and Price, " F i n i t e Element Analysis of Hot Metal Forming", Proc. 18th Int. Conf. Mach. Tool Design Research, 1977. pp. 267-274. 47. Lee, E.H., Mallet, R.L. and Yang, W.H., "Stress and Deformation Analysis of the Metal Extrusion Process", Computer Methods i n Applied Mechanics and Engg., Vol. 10, 1977. PP. 339-353. 48. Lee, E. H., Mallet, R.L., and McMeeking, R.M., "Stress and Deformation Analysis of Metal Forming Processes, "Numerical Modelling of Manufacturing Processes", ASME Special Publication PVP-PB-025, 1977, pp. 19-34. 49. Odell, E.I., "A Study of Wall Ironing by the F i n i t e Element Technique", ASME Paper No. 77-WA/Prod-8. 50. W i f i , A.S., "An Incremental Complete Solution of the Stretch Forming and Deep-Drawing of C i r c u l a r Blank using a Hemispherical Punch", Int. J. Mech. Science, Vol. 18, 1976, pp. 23-31. 51. Wang, N.M. and Budiansky, B., "Analysis of Sheet Metal Stamping by F i n i t e Element Method", General Motors Research Publication GMR-2423, Sept. 1977. 52. Shah, S.N. and Kobayashi, S., "A Theory of Metal Flow i n Axisymmetric Piercing and Extrusion", J. Rud. Engrg., Vol. 1, 1977. PP. 73-103. 53. Kobayashi, S., "Rigid-Plastic F i n i t e Element Analysis of Axisymmetric Metal Forming Processes, Numerical Modelling of Manufacturing Processes", ASME Special Publication PVP-PB-025, 1977, PP. 49-65. 93 54. Kobayashi, S. and Matsumoto, H., "A note on the Matrix Method for R i g i d - P l a s t i c Analysis of Ring Compression", Proc. 18th MT DR Conference, London, 1977, pp. 3-9. 55. Kobayashi, S. and Chen, C.H., "Deformation Analysis of Multi-Pass Bar Drawing and Extrusion", CIRP Annalus, 1978. £6. Price, J.W.H. and Alexander, J.M., "A study of Isothermal Forming or Creep Forming of a Titanium 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 for Analysis of Metal Forming Processes", Trans, of CSME, Vol. 2, No. 1, 1973-74. 58. Wilkins, M.L., "Calculation of E l a s t i c - P l a s t i c Flow", Lawrence Radiation Laboratory Report, UCRL-7322, Rev. 1, Livermore, University of C a l i f o r n i a , 1969. 59. Gordon, P. and Karpp, R., "Application of a New F i n i t e Difference Method of Metal Forming", Proc. Fourth NAMRC, Ba t t e l l e , Columbus Lab., 1976, pp. 72-79. 60. Shabaik, A.H., "Computer Simulation of Metal Flow During Extrusion", Proc. 1976 Int. Conference on Computer Simulation f o r Materials Applications, Nuclear Metallurgy, Vol. 20, Part 2, 1976, pp. 752-765-61. Woo, D.M., "On the Complete Solution of the Deep Drawing Problem", Int. J. Mech. S c i . , 10, 1968, pp. 83-94. 62. Wang, N.M. and Shammamy, M.R., "On the P l a s t i c Bulging of a C i r c u l a r Diaphragm by Hydrostatic Pressure", J. Mech. Phys. Solids 17, 1969, pp. 43-61. 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 Stretching", 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 Analysis of the Hydraulic Bulge Test by the Membrane Theory", Manf. Res. 21, 1969. 65. Thomsen, E.G. and Lapslley, J.T., "Experimental Stress Determination within a metal during P l a s t i c Flow", Proc. Soc. Exptl. Stress Analysis., 11, No. 2, 59-68, 1954. 66. Thomsen, E.G., Yang, C.T., and Bierbower, J.B., "An Experimental Investigation of the Mechanics of P l a s t i c Deformation of Metals". University of C a l i f o r n i a Press, Berkeley 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. 94 68. Shabaik, A. Altan, 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., "Investigation of the Application of V i s i o p l a s t i c i t y Methods of Analysis to Metal Deformation Processing",- F i n a l Report prepared for the U.S. Navy, Bureau of Naval Weapons, Feb. 1966. 70. Shabaik, A. and Thomsen, E.G., "Investigation of the Application of V i s i o p l a s t i c i t y Methods of Analysis to Metal Deformation Processing", F i n a l Report prepared for the U. S.Navy, Bureau of Naval Weapons, Feb. 1967. 71. Shabaik, A. and Kobayashi, S., "Computer Application to the V i s i o p l a s t i c i t y Method", Journal of Engineering f o r Industry, Trans. ASME, Series B., Vol. 89, No. 2, May, 1967, PP. 339-346. 72. Lee, C.H., and Shiro Kubayashi, "Matrix Method of Analysis for P l a s t i c Deformation Mechanics and i t s Application 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 Extrusion", Journal of S t r a i n Analysis, Vol. 7, No. 3, 1972, 'pp. 170-177. 74. Shabaik, A.H., "Computer Aided V i s i o p l a s t i c i t y Solution to Axisymmetric Extrusion Through Curved Boundaries", Journal of Engineering for Industry, Nov., 1972, pp. 1225-1231. 75- Medrano, R.E., G i l l i e s , P.P., Conrad, H. and Hinesley, CP., "Use of Microstructure for V i s i o p l a s t i c i t y Analysis", Journal Str a i n Analysis, Volume 9, Number 3, pp. 146-151, July, 1974. 76. Brown, J.H. and P.F. Thomson, "Mechanics of Deformation During Cold Rolling", A u s t r a l i a ' s Conference on the Mechanics of Structures and Materials, pp. 415-422, 1977-77« Robinson, J.N. and Shabaik, A.H., "Determination of the Relationship between Stra 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., "Application of the V i s i o p l a s t i c i t y Technique f o r Derivation of the C r i t e r i o n for Ductile 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", Engineering Fracture Mechanics, Vol. 7, No. 4 , pp. 63l-'640, 1975. 79. Green,A.P.,"The Use of P l a s t i c i n e Models to Simlate the P l a s t i c Flow of Metals, " P h i l . Mag. Ser. 7.,Vol. 42,Page 365-373,] 95? . 95 A P P E N D I X COMPUTER PROGRAM FOR VISIOPLASTICITY 1 IMPLICIT SEAL*8 (A-H,0-.Z) 2 COMMON XX (15) ,XI,YI,CC,CN,CA,Y(40G,2) ,EBR (400,20) 3 COMMON/C1/EHO,YOID{15,2) ,DT 4 DIMENSION XYT (2 , 20,20,2) ,0(20,20) ,V (20,20) ,XY(2,80,30 ) 5 DIMENSION X (15, 40 0) , ERBOR ( 2) , I SUB (8 0 ,8 0) 6 REAL*4 Z1 (90,90) ,Z2 (90,90) ,TAUXY (80 ,80) 7 LOGICAL REFINE 8 DIMENSION EBBDT(80,80) ,EXDT (80,80) , EYDT (80,80) 9 1,GAMDT (80,80) , LANDA (80,80) 10 BEAL*8 DEBB(80,80) 11 BEAL*8 2LS(10,800) 12 REAL*8 LANDA 13 HEAL*8 X1 (20,20 ,2) ,Y1 (20,20, 2) 14 DIMENSION AINT2(80) 15 C 16 C IX=NG PTS IN X 17 C IY=NO PTS IN Y 18 C IT=NO OF TIME STEPS 19 C DT=TIME INTERVAL BETWEEN TIME STEPS 20 C CC,CN ARE CONSTANTS WHERE SIGB=CC*EB**CN 21 C CA IS LOWER INTERVAL OF INTEGRATION FOR SIG Y 22 C SIGYOA IS CONSTANT ADDED TO SIGY 23 C XYT(L,I,J,K) CONTAINS: L= 1 X-COGED; L=2 Y-C OOBD 24 C FOB 1=1,11 J=1,IX L=1,2 25 C 26 READ(5,10) IX,IY,IT 27 10 FORMAT (312) 28 READ (5, 20) DT 29 20 FORMAT(3F10.0) 30 READ(5,20) CC,CN , CA, SIGYOA 31 C 32 C NIX CONTAINS # GRID PTS IN X DIRECTION, NIY IN Y DIRECTION FOR PLOTS 33 C 34 READ (5,10) NIX,NIY 35 BEAD (5, 20) BHO 36 C 37 C X S Y COORD READ FOB TIME=0 38 C 39 BEAD(4,30) ( ( (XYT (K, I, J, 1^  ,K= 1 ^ 2) ,1 = 1 ,1Y) , J= 1 ,IX) 40 DO 25 11=1,IX 41 DO 2 5 IJ=1,IY 42 IF (II. EQ. 1). XYT (1,IJ,II, 1)=0.D0 43 XYT{2,IJ,II,1)=XYT(2,IJ,II,1)-1.0D0 46 25 CONTINUE 47 30 FORM AT(5X,2F6.3,1X, 2F 6. 3, 1 X, 2F 6. 3 ,1X,2F6.3,1X,2F6.3) 48 LTM=2 49 NTM= 1 50 IDIM=20 51 IDIM.P=80 52 IXY=.IX*IY 53 CALL AX.IS(0.,0.,JX»,-3,lO-,0.,0.,.2) 54 CALL AXIS (0. ,0. , *Y', 1 ,10. ,90. ,0. ,.2) 55 CALL PLOT (XYT (1 ,IY„ 1,1) *5» , XYT (2,IY, 1, 1) *5, ,3) 56 CALLPLOT (XYT (1,IY,IX, 1) *5.,XYTf2,IY,IX,1) *5. ,2) COMPUTER PROGRAM FOR VISIOPLASTICITY 57 CALL PLOT (XYT (1 , 1,IX, 1) *5, ,XYT (2, 1, IX, 1) *5- ,2) 58 DO 35 1=2,IX,2 59 DO 35 J= .2,11,3 60 X1(I,J,1)= XYT(1,J,I, 1) 61 35 Y1 (I,J,1)=XYTf2*J irI#D 62 C 63 C INITIALIZE EBR TO ZERO 64 DO 40 .1=1,20 65 DO 40 3=1,400 66 EBR (J,I) =0. DO 67 40 CONTINUE 68 FACT= 1- DO 69 C 70 c FOR EACH TIME STEP EBB IS ACCUMULATED 71 IT1=IT-1 72 DO 80 K= 1,IT1 73 c 74 c X&Y COORD ARE READ FOR NEXT TIME STEP 75 c 76 NTM= 3-NTM 77 1TH=3-1TM 78 READ {4, 30) {{ (XYT (KK,I, J,NTM) ,KK=1,2) ,1=1 ,IY) , J=1 79 DO 45 1=1,IX 80 DO 45 J=1,IY 81 IF (I.EQ. 1) XYT (1 ,J,.I,NTM) =0. DO 81.6 XYT(2,J,I,NTH)=XYTi2,J,I,NTH)-1-0D0 85 45 CONTINUE . 86 c 87 c U,V CALCULATED FOR THIS TIKE STEP 88 c U MOST BE >0, V MUST BE <0 89 c 90 48 DO 50 J=1,IX 91 DO 50 1=1,IT 92 U (I,J)=- (XYT(1,I,J,N.TM) -XYT { 1 ,1, J,1TM) ) /DT 93 y {I , J)=- {XYT(2,I,J,NTM) -XYT (2,I,J,1TM^ )/DT 94 IF (U (I, J) . GT. 0. DO) UtI,J)=0.D0 95 IF {V (I, J) . LT.O-DO) V(I,J)=0.D0 96 50 CONTINUE 97 C 98 C CURVE FITTING FOR U AN D V USING DLSQHS 99 C SET UP INDEPENDENT VARIABLES IN X 100 C DEPENDENT VARIABLES IN Y 101 C 102 DO 60 3=1,IX 103 11=1Y* (J-1) 104 DO 60 1=1,IY 105 1=1+11 106 CALL AUX (XYT (1,1,3, NTM) , XYT (2,1, J,N'TM) ,X(1 ,L) ) 107 XLS{1 ,1) =-X{2,L) 108 XLS (2,L)=-2.D0*X (10,L) 109 XLS (3,1) =-.3-P0*X(11,L) 1 10 XLS (4,1) =-4.D0*X{12,l) 111 XLS (5,1) =-. 5D0*X (3,1) 112 XLS (6,1) =-X (13,1) 113 XLS (7,1) =-1. 5D0*X (14,1) 114 X1S(8,1) =-1.D0/3.D0**X (4,1) 115 X1S(9,1) =-2.D0/3.D0*X (15,1) 116 X1S(10,1)=-. 25D0*X (5,1) COMPUTER PROGRAM FOR VISIOPLASTICITY 98 1 17 DO 55 IK=1, 10 118 55 XLS {IK,1XY+ L) =X (I K+5 , L) 119 Y(L,1)=U(I,J) 120 60 Y (I XY+L , 1) = V (I, J) 121 CALL DLSQHS (Y,XLS, 2*IXY, 10 , 1 , 803 , 10 ,ERROR FALSE. ,IER ,&200) 122 DO 62 IK=1,5 123 62 Y (IK,2} =0.00 124 DO 63 IK = 6,15 125 63 Y (IK,2) = Y(IK-5,1) 126 Y (1, 1)=0. DO 127 Y (2 ,1)=-Y (6 ,2) 128 Y (3, 1) =-. 5DG*Y{ 10,2) 129 Y(4, 1)=-1.DO/3.D0*Y (13,2) 130 Y (5, 1)=-. 25D0*Y {15,2) 131 DO 64 IK=6,9 132 64 Y (IK, 1) =0.D0 133 Y(10,1)=-2.D0*Y(7,2) 134 Y (11,1) =-3. DG*Y{8,2) 135 Y(12, 1)=-4. D0*Y{9,2) 136 Y (13, 1) =-Y (11 ,2) 137 Y(14,1|=-1.5D0*Y(12,2) 138 Y(15, 1) =-2-D0/3.D0*Y (14,2) 139 C 140 C U&V COEFF SAVED FOR DO/DT,DV/DT 141 C 142 IF (K. NE. (IT1-1) ) GO TO 410 143 DO 400 1=1, 15 144 DO 400 J=1,2 145 400 YOLD (I, J) = Y (I, J) 146 410 CONTINUE 147 C 148 C THE VALUES OF EDTX,EDTY,GAMXY,EBRDT,AND EBR ARE CALCULATED 149. C AT EACH TIME STEP 150 C 151 DO 52 1=2,IX,2 152 DO 52 J=2,IY,3 153 X1 (I,J,NTM) = XYT (1 , J,I,LTM) -DT * AU X2 {XYT *( 1 , J , I , NTM) , 154 1 XYT (2, J,I, NTM) , Y .{1, 1) ,- TRUE. ) 155 Y1 (I, J, NTM) =X YT (2 , J, I , LTM) - DT*AUX2(XYT{1,J,1,NTM), 156 1 XYT (2,J, I,NTM) ,Y(1,2) , - FALSE. ) 157 CALL PL0T(X1 (I,J,LTM) *5„ , Y 1 (I, J, LTM) *5. ,3) 158 CALL PLOT (X1 (I, J, NTM) *5. , Y 1 (I, J, NTM) *5. ,2) 159 CALL SYMBOL (X1 (I, J, NTM) *5. ,Y1 (I, J,NTM) *5. ,. 14,5,0. ,- 1 ) 160 52 CONTINUE 161 IF (K.NE.IT1) GO TO 53 162 CALL PLOT (0. ,XYT (2,IY, 1 , NTM) *5. , 3) 163 CALL PLOT (XYT (1, IY,IX, NTM) *5.,XYT (2 ,1Y,IX, NTM) *5. ,2) 164 DO 49 1=2,IY 165 CALL PLOT (XYT (1 , IY-I +1,IX, NTM) *5-, XYT{ 2, IY-I+ 1 , IX ,NTH )*5.,2) 166 49 CONTINUE 167 CALL PLOT (12. ,0. ,-3) 168 53 CONTINUE 169 DO 70 J=1,IX 170 IL=IY*(J-1) COMPUTER PROGRAM FOR VISIOPLASTICITY 171 DO 70 1=1,IY 172 IJ=I+IL 173 CALL AUX {XYT {1, I, J, NTM) , XYT (2,1, J, NTM) , XX) 174 CALL DEBIV(Y(1,1) ,XYT (1,1,3,NTM) ,XYT (2 ,1, 3,NTM) ,DUDX, DUDY,3) 175 CALL DERIV (Y (1,2) ,XYT (1,1, J,NTM) ,XYTi(2 ,I,3,NTM) ,DVDX, DVDY, 1) 176 GAMXY=DUDY+ DVDX 177 DPACT= (DUDX**2)/3-D0+(GAMXY) **2/12. DO 178 65 EBR DT (I, J) =2- D0*DSQRT (DFACT) 179 IF (K.EQ.IT1) FACT=.5D0 180 EBR (13, 1) =EBR (IJ, 1) +FACT*DT*EBRDT (I , J) 181 70 CONTINUE 182 80 CONTINUE 183 C 184 C 4TH DEGREE POLY FIT TO EBR 185 . C 186 CALL D1SQHS-{EBB, X,IX*IY, 15,2,400, 15, ERROR, .FALSE. ,IER ,S200) 187 C 188 C MASTER GRID IS SET FOR FINAL PLOTS 189 C 190 C XY (2,1,3) CONTAINS THE MASTER GRID ST XY(1, 1,3) IS 191 C THE X-COORD, XY (2,1,3) IS THE Y-COOSD FOR 1=1,IX 3=1,IY 192 XMAX = 0-D0 193 DO 90 1=1,IY 194 IF (XYT ( 1,1, IX, NTM) - LT. XMAX) GO TO 90 195 XMAX=XYT(1,1,IX,NTM) 196 IMY=I 197 90 CONTINUE 198 YMAX=XYT(2,IY,IX,NTM) 199 100 CONTINUE 200 DY=YMAX/(NIY-1) 201 DO 110 1=1,NIX 202 XY(2,I, 1)=0.D0 203 DO 110 3=2,HIY 204 XY (2,1,3) =XY (2,1,3-1) +DY 205 110 CONTINUE 206 DX=XMAX/(NIX-1) 207 DO 120 1=1,NIY 208 XY(1,1,I)=0.D0 209 DO 120 3=2,NIX 210 XY(1,3,I)=XY (1,3-1,1) +DX 211 120 CONTINUE 212 C 213 C TO FIND ZERO FILLS ON PLOTS 214 C 215 CALL FILL<XYT,IX,IY,ISUB,XY,NIX,NIY,IMY,NTa,IDIM,IDIM 216 C 217 C TO PLOT U, V,EBSDT, AND TAUXY 218 C 219 DO 125 3=1,NIY 220 DO 125 1=1,NIX 221 EBR DT {I, 3) =0. DO 222 DEB1 (I, 3) =0 . DO COMPUTES PROGRAM FOR VISIOPLASTICITY 223 224 125 TAUXY (I, J) = 0.1)0 NIY10=NIY+1Q 100 225 NIX10=NIX+10 226 DO 126 J=1,NIY10 227 DO 126 1=1,NIX 10 228 Z1 (I, J) =0.D0 229 126 Z2{I, J)=0.D0 230 DO 130 J=1,NIY 231 DO 135 1=1,NIX 232 IF (I SD3 {.I, J) . EQ. 0) GOTO 130 233 Z 1 (1*5,3 + 5) =-AUX2 (XY (1,I,J) , XY{2,I, J) ,Y {1 ,1) ,.TRUE. ) 2 34 Z2(I + 5,J+5) =AUX2 (XY (1 ,I,J) ,XY{2,I,3) ,Y (1,2) ,.FALSE.) 235 CALL DERIV (Y, XY (1,1, J) ,XY (2,1, J) , EX DT (I, J) ,G AMDT (I,J) ,3) 236 CALL DERIV {Y (1,2) , XY (1,1, J) , XY (2,1, J) , EBRDT {I, J) , EYDT (I,J) ,3) 237 GAMDT (I, J) = G AMDT (I, J) +EBSDT (I, J) 238 DFACT=(3.DO*EXDT(I,J)** 2+.75DQ*G AMDT(I,J) **2) 239 131 EBR DT (I,J)=2-DO/3.DO *DSQRT(DFACT) 240 DEBE (I,J)=AUX2{XY{1,1,J) ,XY (2,1,J),EBR,•FALSE.) 241 LAN DA (I, J) = 1.5D0*EBRDT(I,J)/(CC*DEBR (I, J) **CN) 242 TAUXY (I ; r J) = GAMDT (I, J) /{2- DO*.L ANDA {I , J) ) 243 135 CONTINUE 244 130 CONTINUE 245 C 246 C PLOT U,V,EBRDT,TAUXY 247 C 248 DXY = XY(2,NIX,NIY)/XY (1,NIX,NIY) 249 10., CALL PERS{Z1,IOIflP+10,NIX+10,NIY+10,DXY,.333,45. 10.) ,45. , 250 CALL PLOT (12. ,0. ,-3) 251 10., CALL PERS(Z2,IDIMP+10,NIX+10,NIY+10,DXY,.333,45. 10.) , 45. , 252 CALL PLOT {12. ,0. ,-3) 253 DO 137 J=1,NIY10 254 DO 137 1=1,NIX10 255 Z1 (I,J) =0. 256 137 Z2 (I,J) =0. 257 DO 138 J=1,N.IY 258 DO 138 1=1,NIX 259 138 Z1 (1 + 5, J + 5) =DEBR (I, J) 260 10., CALL PEP,S{Z1,IDIMP+10,NIX+10,NIY + 10 ,DXY,. 333,45. 10.) ,45., 261 CALL PLOT(12. ,0. ,-3) 262 DO 139 J=1,NIY10 263 DO 139 I=1,NIX10 264 139 Z 1 {I, J) =0. . 265 DO 140 3=1,NIY 266 DO 140 1=1,NIX 267 21 {1 + 5, J + 5) =EB1DT{I,J) 268 140 Z2 (1 + 5, J+5) =-TAUXY{I, J) 269 10., CALL PERS{Z1,IDIMP+10,NIX+10,NIY + 10,DXY,. 333,45. 10. ) ,45. , 270 CALL PLOT (12. ,0.,-3) 271 10., CALL PERS {Z2,ID.IM? + 10,NIX+1G,NIY + 10 ,DXY,. 333,45. 10.) ,45. , 272 CALL PLOT {12. ,0. ,-3) 273 C COMPUTER PROGRAM FOR VISIOPLASTICITY 274 C CALC SIGY 101 275 C 276 EXTERNAL DF1,DF2 277 DO 141 J=1,N.IY10 278 DO 141 1=1,NIX 10 279 Z1 (I,J)=0. 280 141 Z2 {I, J) =0. 281 DGRID=(XY(2, 1,2) -XY (2,1, 1) ) /2 282 JGRID=1 2 83 1000 IF (CA,. GT. XY (2 , 1 , JGRID) + DGRID) GO TO 1010 284 GO TO 1020 285 1010 JGRID=JGRID*1 286 IF (JGRID. EQ. NIY) GO TO 10 20 2 87 GO TO 1000 288 1020 AINT2(1) =0. DO 289 DO 1030 1=2, NIX 290 AINT2 (I) =AINT2 (1-1) 291 IF (ISOB (I, JGRID) . EQ. 0) GO TO 10 30 292 AINT2 (I) =AI13T2 (I) +DQUANK (DF2, XY( 1 ,1-1,JGRID) 293 1 XY{ 1 ,1, JGRID) 00 1 DO ,TOL,FIFTH) 294 1030 CONTINUE 295 DO 150 J=1,NIY 296 DO 150 1=1,NIX 297 IF (ISUB(I,J)-EQ.0) GO TO 150 298 XI= XY (1 , I , J) 299 YI=XY (2,1, J) 300 AINT1=0.D0 301 IF {JGRID .EQ. J) GO TO 1060 302 JADD=0 303 IF (ISUB (I, JGRID) . EQ. 1) GO TO 1050 304 JINC=1 305 IF (JGRID. GT.J) JINC=-1 306 1040 JADD=JADD+J.INC 307 IF (JGRID+JADD. EQ. J) GO TO 1060 308 IF (ISUB (I, JGRID+JADD). EQ. 1) GO TO 1050 309 GO TO 1040 310 1050 AINT1=DQUANK (DF1 ,XY (2 ,1, JGRID+JADD) ,XY (2,1, J) 001DO, TOL, FIFTH) 311 1060 Z1 (1*5, JS-5) =SIGYOA-AINT1-AINT2 (I) 312 Z2(I + 5, J + 5) =Z1 (1+5, J+5) + (EXDT (I, J ) - EYDT (I,J) ) /LANDA (I ,J) 3 13 150 CONTINUE 314 10., CALL PERS(Z1,IDIMP+10,NIX*10,NIY+10,DXY,.333,45 10.) .,45. , 315 CALL PLOT (12. ,0. ,-3) 316 10., CALL PERS (Z2,IDIMP+10,NIX+10,NIY +10,DXY,- 333,45 10-) -,45., 317 CALL PLOTND 318 STOP 319 200 STOP 1 3 20 END 321 SOBROUTINE AUX(X,Y,XX) 322 IMPLICIT REAL*8 (A-H,0-Z) 323 DIMENSION XX (15) 324 XX{ 1) =1- DO 325 XX (2)=X 326 XX(3) =X*X 327 XX{4)=X*XX(3) COMPUTES PROGRAM FOR VISIOPLASTICITY 328 XX(5)=X*XX(4) 329 XX(6)=Y 330 XX (7) =Y*Y 331 XX(8)=Y*XX(7) 332 XX{9) =Y*XX (8) 333 XX( 10) = X*Y 334 XX{11)=XX (10) *Y 335 XX{ 12}=XX (1 1) *Y 336 XX (13) = XX (10) *X 337 XX (11) =XX(13) *Y 338 XX {15) =XX (1 3) *X 339 RETURN 340 END 341 FUNCTION AUX2 (XX,YY,P,LL) 342 IMPLICIT REAL*8(A-H,Q-Z) 343 C 344 C EVALUATE THE FITTED FUNCTION AT XX,YY 345 C P CONTAINS THE FITTED PARAMETERS 346 C LL IS TRUE IF THE INDEPENDENT VARIABLE MUST BE EVALUATED BY AUX 347 C 348 LOGICAL LL 349 COMMON X 350 DIMENSION X (15) , P (1) 351 IF (LL) CALL AUX (XX,YY,X (1) ) 352 AUX2=0.D0 353 DO 10 1=1,15 354 AUX2=AUX2+P (I) *X<I) 355 10 CONTINUE 356 RETURN 357 END 358 FUNCTION DF 1 (YC) 359 IMPLICIT REAL*8(A—H,0-Z) 360 COMMON XX(15) ,XI,YI,C1,C2,CA,Y (4 00, 2} , PEBR (20, 20) 361 COMflON/C1/RHO,YOLD(15,2),DT 362 C 363 C THE INTEGRAND DTAU/DX IS EVALUATED 364 c 365 EB1=AUX2 (XT ,YC,PEBR,.TRUE. ) 3 66 CALL DERIV (Y{1, 1) , XI, YC,EXDT, DUDY,3) 367 CALL DERIV(Y(1,2) , XI , YC, DV DX, DVDY,3 ) 368 GAM DT=DU DY + DVDX 369 EBRDT=2-D0/3-D0*DSQRT (3. DO *EXDT**2+. 75D0*GAMDT**2) 370 CALL DERIV (PEB1;XI,YC,DEBRDX,DUM ,1) 371 CALL DERIV2 (Y (1, 1) ,XI, YC,DUDXY,3) 372 CALL DER.IV2 (Y (1,2) , XI,YC,DVDXX,1) 373 DGAM DX= DUDX Y+DVDXX 374 CALL DERIV2 (Y (1, 1) , XI, YC , D EXDX , 1 ) 375 DEBRDT=(4.DO*EXDT*DEXDX+GA MDT*DGAMDX)/(3.D0*EBRDT) 376 BF1=C1*EBR**C2/(3-D0*EBRDT)*(C2*DEBRDX*GA MDT/EBR 377 1 +DGAMDX - DEBSDT*GAMDT/EBRDT) 378 U=AUX2(XI,YC,Y(1, 1) ,. FALSE.) 379 V=aUX2(XI,YC,Y(1,2) FALSE.) 3 80 VOLD=AUX2 (XI,YC,YOLD (1,2) ,-FALSE.) 380-5 FACT=RHO* {(V-VOLD)/DT) 382 DF1=DF1+FACT 3 83 RETURN 384 END COMPUTER PROGRAM FOR VISIOPLASTICITY 103 3 85 FUNCTION DF2{X) 386 IMPLICIT REAL*8(A-H,0-Z) 387 COMMON XX{15) ,XI,YI,C1,C2,CA,Y (400, 2) , PEBR (20, 20) 388 COM MON/C1/R HO,YO.L.D(15,2) ,DT 3 89 REAL*8 LAMBDA 390 C 391 c THE INTEGBAND DTAU/DY IS EVALUATED 392 c 393 EBR=AUX2{X,CA,PEBB,.TRUE.) 394 CALL DERIV (Y (1, 1) , X, CA ,EXDT, DO DY , 3) 395 CALL DERIV {Y{1, 2) , X,CA,DVDX,EYDT,3) 396 GAMDT=DUDY+DVDX 397 EBR DT=2.DO/3.D0*DSQBT<3.D0*EXDT**2+.75D0*GAMDT**2) 398 CALL DEHIV(PEBR,X,CA,DUM,DEBRDY,2) 399 CALL DEB.IV2 (Y (1, 1) ,X,CA,DUDYY,2) 400 CALL DEBIV2 (Y (1,2) , X ,CA, DVDXY ,3) 401 DGAMDY=DUDY Y+DVDXY 402 CALL DEBIV2 (Y (1, 1) ,X,CA,DEXDY,3) 403 YEBRDT={4-DO*EXDT*DEXDY+GAMDT*DGAMDY)/{3.D0*EBBDT) 404 CALL DER IV 2 {Y (1, 1) ,X,CA,DEXDX,1) 405 CALL DEF.IV2 {Y {1, 2) ,X,CA,DEYDX,3) 406 LAMBDA= (2. D0*C1*EBB**C2) /(3. DO*EBRDT) 407 CALL DERIV(? EBB,X,CA,DEB RDX> DUM,1) 408 CALL DEBIV2{Y (1,2) ,X,CA,DVDXX, 1) 409 DGAMDX=DEXDY+DVDXX 410 DEBR DT=(4. DO*EXDT*DEXDX+GAMDT*DGAMDX) /(3.D0*EBBDT) 411 DF2=LAM BDA*({(EYDT-EXDT)*{-C2*DEBBDX/EBB 412 1 *DEB8DT/EBBDT) +DEXDX-DEYD X) 413 1 +.5D0*(C2*DEBBDY*GAMDT/EBB + DGAMDY - YEBBPT*GAMDT/E BBDT) ) 4 14 U=AUX2{X,CA,Y(1,1) ,. FALSE. ) 4 15 V= A U X2 f X , C A , Y { 1 , 2) ,,. FALSE. ) 416 UOLD=AUX2 (X,€A,YOLD (1 ,1) ,. FALSE. ) 416.5 FACT=BHO* { (U-UOLD) /DT) 418 DF2=DF2 + FACT 4 19 BETURN 420 END 421 SUBROUTINE DERIV{A,X,Y,DUDX,DUBY,N) 422 C 423 C EVALUATE DERIV WBT X AND Y 424 c IF N=1 DUDX, IF N=2 DUDY, OTHERWISE DUDX AN D DUDY 425 C 4 26 IMPLICIT REAL*8 (A-H,0-Z) 4 27 COMMON XX {15) 428 DIMENSION A{1) 429 IF (N . EQ. 2) GO TO 10 430 DUDX=A (2)+2.D0*A {3) *X + 3. DQ* A (4) *XX {3) + 4.D0*A(5)*X X(4) 431 1 *A (10) *Y+A (11) *XX{7) + A(12)*XX(8) * 2. B0*A { 1 3) *XX { 1 432 1 +2.D0*A {14)*XX (11) + 3.D0*A (15)*XX (13) 433 IF (N. EQ. 1) BETURN 434 10 DUDY=A(6) + 2.D0*A (7) *Y +3.D0*A{8) *XX (7) + 4.D0*A{9)* XX <8) 435 1 +A{10)*X + 2.P0*A (1 1) *XX{ 10) +3.D0*A (12) *XX {1 1) 436 1 + A(13)*XX(3) +2.D0*A (14) *XX (13) + A(15)*XX(4) 437 BETURN COMPUTER PROGRAM FOR VISIOPLASTICITY 104 438 END 439 SUBROUTINE DERIV2(A,X,Y,DUDD,N) 440 IMPLICIT REAI*8 (A-H,0-Z) 441 C 4 42 C EVALUATE 2ND ORDER DERIV WRT X AND Y 443 c N=1 DXDX; N=2 DYDY; N=3 DXDY 444 c 445 DIMENSION A{1) 446 GO TO (10,20,30),N 447 10 DUDD=2~ D0*A (3) + 6.D0*A(4)*X + 1 2-D0*A (5) *X*X 448 1 +2. D0*A (13) *Y +2.D0*A(143 *Y*Y * 6. D0*A (1 5) *X*Y 449 RETURN 450 20 DUDD=2. DQ*A (7) + 6.D0*A(8)*Y + 12. D0*A (9) *Y*Y + 451 1 2. D0*A(11)*X * 6..D0*A (12) *X*Y + 2. D0*A ( 14) *X*X 452 RETURN 453 30 DUDD=A(10) + 2. D0*A(11)*Y *• 3.D0*A{12)*Y*Y 454 1 + 2.D0*A(13)*X + 4. DG*A (14) *X*Y + 3. D0*A (15) *X*X 455 RETURN 456 END 457 SUBROUTINE FILL(XYT,IX,IY,ISUB,XY,NIX,NIY,IMY,NTM,IDI M,IDIMP) 458 IMPLICIT REAL*8 (A-H,0-Z) 459 1 \ DIMENSION XYT(2,IDIM,IDIM, 2) , ISUB (I DIMP, 1) ,XY (2,IDIMP 460 » '; C 461 C I S O B CONTAINS 1 WHERE A FUNCTION VALUE IS 4 62 C PLOTTED, 0 IF OUTSIDE BOUNDARY 463 C 464 DO 10 1=1, NIX 465 DO 10 J=1 ,NIY 466 10 ISUB(I,J) = 1 467 XMIN=XYT(1,IMY,IX,NTM) 468 DO 15 1=1, IY 469 IF (XYT (1 ,1,IX, NTH) .GT. XMIN) GO TO 15 470 XMIN=XYT (1 ,1 , IX, NTM) 471 IMYMIN=I 4 72 15 CONTINUE 473 DO 110 J=1,NIY 4 74 DO 100 1=1,NIX 4 75 IF {XY(1,I,J).LT.X¥T(1,IMYMIN,IX,NTM)) GO TO 100 476 IF (XY(1,I,J).GT.XYT(1,IMY,1X,NTM)) GO TO 80 477 LB=1 478 NB=2 479 20 IF{XY(2,I,J).LT-XYT(2,NB,IX,NTM) ) GO TO 30 480 LB=LB+1 481 NB=NB*1 482 IF (NB. LT. IY) GO TO 20 483 30 IF (XY (1 ,I,J) .LT. XYT (1 ,LB,IX, NTM) . AND. 484 1 XY(1,I,J) . LT.XYT(1,NB,IX,NTM) ) GO TO 100 485 IF (XY(1 ,1,J).GT.X¥T (1, LB,IX , NTM) .AND. . 4 86 1 XY(1,I,J) .GT.XYT(1,NB,IX,NTM)) GO TO 80 487 YN=XYT(2 ,LB,IX, NTM) + (XYT (2 , NB,IX , NTM) -XYT (2 ,LB, IX, NTM 488 / ) 1 /{XYT(1,NB,IX,NTM)-XYT(1,LB,IX,NTM) )* 489 1 (XY { 1,1, J)-XYT (1,LB,IX, NTM) ) 490 IF (YN.GT.XY (2 ,1 , J) . AND. XYT (1 ,LB,IX, NTM) . GE. 491 1 XYT(1,NB,IX,NTM}) GO TO 100 492 IF (YN.LT. XY (2,I,J) . AND. XYT (1 ,LB,IX,NTM) . LE-COMPUTES PBOGBAM FOB VISIOPLASTICITY 1 0 5 1 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 493 494 80 495 90 496 497 100 4 98 110 4 99 500 PBOGBAM FOB TWO DIMENSIONAL PLOT 106 1 IMPLICIT BEAL*8 (A-H,0-Z) 2 COMMON XX (1 5) ,XI,YI,CC,CN,CA,Y (400, 2) , EBR (400, 20) 2.5 COMMON/C1/RHO,YOLD{15,2) ,DT 3 DIMENSION XYT {2,20,20,2),0 (20,20),V(20,20),XY{2,80,80 ) 4 DIMENSION X (15,400) , ERROE (2) , ISO B (80, 80) 5 REAL*4 Z 1 (90,90) ,Z2{ 90,90) ,TAUXY (80 ,80) 6 LOGICAL BEF.INE , 7 DIMENSION EBRDT (80, 80) , EXDT (80,80) , EYDT (8 0,80) 8 1,GAMDT(80,80),LANDA(80,80) 8.05 REAL*8 DEBR (80,80) 8.07 BEA1*8 XLS( 10,800) 8.2 BEA.L*8 LAND A 8.6 DIMENSION AINT2(80) 9 C 10 C IX=NO PTS IN X 11 C IY=NO PTS IN Y 12 C IT=NO OF TIME STEPS 13 C DT=TIME INTERVAL BETWEEN TIME STEPS 14 C CC,CN ARE CONSTANTS SHEBE SIGB=CC*EB**CN 15 C CA IS LOSES INTEBV AL OF I NTEGBAT.ION FOR SIG Y 16 C SIGYOA IS CONSTANT ADDED TO SIGY 17 C XYTfL,I,J,K) CONTAINS: L= 1 X-COOBD: L=2 Y-C OORD 19 C FOR 1=1,IY J=1,IX L=1,2 20 C 21 BEAD(5, 10) IX,IY,IT 22 10 FORMAT (312) 23 BEAD (5, 20) DT 24 20 FORMAT {8F10.0) 25 READ(5,20) CC,CN ,CA , SIG YOA ,CB 25.1 C 25.2 C NIX CONTAINS # GRID PTS IN X DIRECTION, NIY IN Y DIRECTION FOB PLOTS 25.3 C 25.4 BEAD (5,10) NIX,NIY 25-7 HEAD (5, 20) BHO 26 C 27 C X S Y COORD BEAD FOR T.IME = 0 28 C 29 SEAD(4,30) ( ( (XYT {K, I , J, 1) ,K= 1 ,2) ,1 = 1 ,1Y) , J= 1 , IX) 30 DO 25 11=1,IX 31 DO 25 IJ = 1,IY 31.2 IF (II.EQ. 1) XYT (1 ,IJ,II, 1) =0. DO 31.4 IF(IJ.EQ.I) XYT (2,IJ,II, 1) =0.D0 32 DO 25 IK=1,2 33 IF (XYT {IK , IJ,11, 1) - LT. O.D 0) XYT (IK ,IJ,II , 1) = 0. DO 34 25 CONTINUE 35 30 FORMAT(5X,2F6.3,1X,2F6.3,1X,2F6.3,1X,2F6.3,1X,2F6.3) 35.2 CALL PLOTIT {XYT(1,1,1,1),IX,IY) 36 LTM = 2 37 NTM = 1 38 IDIM=20 38.2 IDIMP=80 38.6 IXY=IX*IY 39 C 40 C INITIALIZE EBB TO ZEBO PBOGBAM FOR TWO DIMENSION AL PLOT 41 DO 40 1= 1, 20 42 DO 4 0 J=1,400 43 EBR(J,I) =0. DO 44 40 CONTINUE 45 FACT= 1. DO 46 C 47 C FOB EACH TIME STEP EBB IS ACCUMULATED 48 IT1=IT-1 49 DO 80 K=1,IT1 50 C 51 C XSY COORD ARE READ FOR NEXT TIME STEP 52 c 53 NTM= 3-NTM 54 LTM=3~LTM 55 BEAD(4,30) ( { (XYT (KK,I,J ,NTM) ,KK=1,2) ,1=1 ,IY) ,J=1,IX) 56 DO 45 .1=1 ,IX 57 DO 45 J= 1,IY 57. 2 IF (I.EQ.1) XYT (1 , J , I,NTM) =0. DO 57. 4 IF (J.EQ.1) XYT(2,J,I,NTM)=0.DO 58 DO 45 IK=1,2 59 IF (XYT(IK,J,I,NTM) .LT.O.) XYT (IK , J , I ,NTM) =0. DO 60 45 CONTINUE 60. 2 CALL PLOTIT (XYT (1,1,1 , NTM) ,1X,1 Y) 61 C 62 C U,V CALCULATED FOR THIS TIME STEP 63 c U MUST BE >0, V MUST BE <0 64 c 65 DO 50 J=1,IX 66 DO 50 1=1,IY 67 U (I , J) =— (XYT (1,1, J , NTM) -XYT(1 ,1, J , LTM) ) /DT 68 V(I,J) =- (XYT ( 2, I, J , NTM) - XYT (2 ,1, J,LTM) )/DT 69 IF (U (I, J) .GT.O. DO) U(I,J)=0.D0 70 IF (V (I, J) .LT.O.DO) V(I,J)=0.D0 71 50 CONTINUE 72 C 73 C CURVE FITTING FOB U AND V USING DLSQHS 74 C SET UP INDEPENDENT VABIABLES IN X 75 C ' DEPENDENT VABIABLES IN Y 76 c 77 DO 60 J=1,IX 78 IL=IY*(J-1) 79 DO 60 I=1,IY 80 L=I+IL 81 CALL AUX (XYT (1,1, J , NTM) ,XYT (2,1, J,NTM) ,X (1,L) ) 84. 6 1 XLS(1,L)=-X (2,L) 84. 62 XLS (2 ,L) =-2.D0*X(10,L) 84. 6 3 XLS (3,L) =-3.D0*X (11,1) 84. 64 XLS (4,L)=-4.D0*X(12,L) 84. 65 XLS (5,L)=-. 5D0*X (3,L) 84. 6 6 XLS(6,L) =-X{ 13,L) 84. 67 XLS (7,L)=-1.5D0*X(14,L) 84. 68 XLS (8, L) =- 1. DO/3 . D 0*X (4 ,L) 84. 69 XLS (9,L) =-2. DO/3. D0*X (15,L) 84. 7 XLS ( 10,L) =-.25D0*X (5,L) 84. 71 DO 55 IK=1, 10 84. 72 55 XLS (IK,IXY + L) =X (IK+5,L) 84. 73 Y(L,1) = U(I,J) 84. 74 60 Y (IXY+L, 1) =V (I,3) PROG.RAM FOR TWO DIMENSIONAL PLOT . 108 84.75 CALL D1SQHS (Y,XLS,2*IXY,10,1,800,10,ERROR,.FALSE.,IER ,&200) 84.76 DO 62 IK=1,5 84.77 62 Y(IK,2)=0-D0 84.78 DO 63 IK=6,15 84.79 63 Y (IK, 2) =Y (IK-5, 1) 84.8 Y(1,1)=0.D0 84.8 1 Y (2, 1) =-Y (6, 2) 84.82 Y{3 ,1)=-.5D0*Y(1Q,2) 84.83 Y (4, 1) =- 1.D0/3.D0*Y (1 3,2) 84.84 Y (5 ,1) =-. 25D0*Y (15,2) 84.85 DO 64 IK = 6, 9 84.86 64 Y(IK,1)=0.D0 84.87 Y(10,1) =-2. D0*Y(7,2) 84.88 Y(11,1)=-3.D0*Y(8,2) 84.89 Y (1 2, 1) =-4. D0*Y (9,2) 84.9 Y (13, 1) =-Y{ 1 1,2) 84.91 Y(14,1)=-1.5D0*Y (12,2) 84.92 Y{15, 1) =-2. D0/3.D0*Y (14,2) 85 C 85.1 C U&V COEFF SAVED FOR DU/DT,DV/DT 85.2 C 85.3 IF (K.NE. (IT1-1) ) GO TO 410 85.4 DO 400 1=1,15 85.5 DO 400 J = 1, 2 85.6 400 YOLD (I, J)=Y ( I , J) 85.7 4 10 CONTINUE 85.8 C 86 C THE VALUES OF EDTX,EDTY,GAMXY,EBEDT,AND EBR ARE C A L C U L A T E D 87 C AT EACH TIME STEP 88 C 89 DO 70 J=1,IX 90 IL=IY*{J-1) 91 DO 70 I=1,IY 92 IJ=I+IL 92.2 CALL AUX (XYT (1, I, J,NTM) , XYT (2 ,1, J,NTM) ,XX) 93 CALL DERIVCY(1,1),XYT(1,I,J,NTM),XYT(2,1,J,NTM),DUDX, DUDY,3) 94 CALL DERIV(Y(1,2) , XYT (1 ,1, J , NTM) ,XYT (2 ,1, J,NTM) ,DVDX, DVDY,1) 95 GAMXY=DUDY+DVDX 96 DFACT= (DUDX**2)/3- D0+ (GAMXY) **2/12. DO 96. 16 65 v EBRDTd, J)=2.D0*DSQ8T (DFACT) 96.2 IF (K.EQ.IT1) FACT=. 5D0 97 EBR (IJ,1)=EBR(IJ, 1| +FACT*DT*EBRDT (I,J) 98 70 CONTINUE 101 80 CONTINUE 102 C 103 C 4TH DEGREE POLY FIT TO EBR 104 c 105 CALL DLSQHS(EBR,X,IX*IY,15,2,400,15,ERROR,.FALSE.,IER ,S200) 106 c 107 c MASTER GRID IS SET FOR FINAL PLOTS 108 c 109 c XY(2,I,J) CONTAINS THE MASTER GRID ST XY{1, I,J) IS PROGBAM FOB TWO DIMENSIONAL PLOT 109 110 C THE X-CQQBD, XY(2,I,J) I S THE Y-COOBD FOB 1 = 1, IX J=1 /IY 111 XMAX=0.D0 1 12 DO 90 .1=1,1 Y 113 IF (XYT{ 1,1,IX,NTM) .LT.XMAX) GO TO 90 114 XMAX=XYT (1,1,IX, NTM) 1 15 IMY=.I 1 16 90 CONTINUE 1 17 YMAX=XYT(2,IY,IX,NTM) 1 18 100 CONTINUE 1 19 DY=YMAX/(NIY-1) 120 DO 110 1=1,NIX 121 XY (2,1, 1) =0. DO 122 DO 110 J=2,NIY 123 XY{ 2,1,3) =XY (2,1, J-1) +DY 124 110 CONTINUE 125 DX=XMAX/ (NIX- 1) 126 DO 120 1=1,NIY 127 XY ( 1, 1, I) =0.D0 128 DO 120 J=2,NIX 129 XY(1,J,I)=XY(1,J-1,I) +DX 130 120 CONTINUE 131 c 132 c TO FIND ZEBO FILLS ON PLOTS 133 c 134 CALL F11L(X YT,IX,IY,IS U B,X Y,NIX,NIY,IM Y,N TH,IDIM,IDIM P) 135 c 136 c TO PLOT U,y,EBBDT, AND TAUXY 137 c 138 DO 125 J = 1,NIY 138. 4 DO 125 1=1,NIX 138. 8 EBRDT {I, J) =0.D0 139 DEBR (1, 3) =0. DO 139. 2 125 TAUXY (I, J) =O.D0 139. 6 NIY10=NIY+1 0 140 NIX10=NIX+10 140. 4 DO 1.26 3=1,NIY10 140. 8 DO 126 1=1,NIX 10 141. 2 Z1 ( I , J) =0. DO 141. 6 126 Z 2 ( I , J) =0.D0 144 DO 130 3=1,NIY 145 DO 135 1=1,NIX 146 IF (ISUB (I, J) - EQ. 0) GOTO 130 147 Z1 (1+5,3 + 5) =-AUX2(XY (1,1,3) ,XY(2,I, J) ,Y{1 , 1) ,.TBUE.) 148 Z2(I + 5, J + 5) =AUX2(XY (1,I,J) ,XY (2,I,J) ,Y (1 ,2) FALSE. ) 149 CALL DE B l V (Y,IY (1,1, J) ,XY{2,I,J) ,EXDT(I,J) ,GAMDT(I,J) 150 CALL DEBIV (Y (1 ,2) , XY (1,1, J) ,XY (2,1, J) , EBRDT (I, J) , EYDT d,J) ,3) 151 GAMDT (I , J) = G AMDT (I, J) +EBRDT ( I , J) 152 DFACT=(3. D0*EXDT (I, J) **2+. 75D0*GAMDT (I ,J) **2) 152. 8 131 EBRDT ( I , J) =2. DO/3. DO *DSQRT (DFACT) 153 DEBR (I, J) =AUX2 (XY (1,I,J) , X Y (2,1, J) , EBR, . F ALSE. ) 154 LAN DA (I, J) = 1. 5D0*EBRDT(I,J) /(CC*DEBR (I ,J) **CN) 155 TAUXY(I,J) = GAMDT (I, J) / (2. D0*L AND A ( I , J) ) 155. 2 135 CONTINUE 156 130 CONTINUE PROGBAM FOE T»0 DIMENSIONAL PLOT 157 C 158 C PLOT U,V,EBHDT,TAUXY 159 C 159.4 CALL PLOT 2 (TAUXY,XY,IX,IY, IS OB, NIX, NIY) 159.6 CALL PLOT( 1 2.,0. ,-3) 160 DXY=XY{2,NIX,NIY)/XY (1,NIX,NIY) 161 C CALL PERS(Zl,IDIMP+l0,NIX+10,NIY+10,DXY,- 333,45.,45 10.,10.) 161.2 C CALL PLOT {12. ,0. ,^ -3) 162 C CALL PERS(Z2„IDIMP + 10,NIX+10,NIY + 10,DXY,. 333.45.,45 10., 10.) 162-5 C CALL PLOT {12. ,0. ,-3) 162-55 DO 137 J=1,NIY10 162.6 DO 137 I=1,NIX10 162.65 Z1(I,J)=0. 162.7 137 Z2(I,J)=0. 162-73 DO 138 J=1,HIY 162.76 DO 138 1=1,NIX 162.79 138 Z 1 (1*5, J + 5) =DEBR (I,J) 162.82 C CALL PERS(Z1,IDIMP*10,NIX+10,NIY +10,DXY, . 333,45-,45 10.,10-) 162-85 C CALL PLOT (12. ,0. , — 3) 162.88 DO 139 J=1 fNIYlQ 162.91 DO 139 I=1,NIX10 162.94 139 Z1(I,J)=G. 163 DO 140 J=1,NIY 164 DO 140 1=1,NIX 165 Z1(1 + 5, J+5) = EBRDT(I,J) 165-2 140 Z2(I*5, J + 5) =-TA0XY (I, J) 166 C CALL ?ERS(Zl,IDIMP+10,NIX+10,NIY+10,DXY,.333,45-,45 10., 10.) 166.5 C CALL PLOT (12. ,0. ,-3) 167 C CALL PEBS(Z2,IDIMP + 10,NIX+10,NIY+10 ,DXY, . 333 ,45. ,45 10.,10.) 167.5 C CALL PLOT (12-, 0. ,-3) 168 C 169 C CALC SIGY 170 C 171 • EXTERNAL DF1,DF2 171.2 DO 14 1 J=1,NIY10 171.4 DO 141 I=1,NIX10 171.6 21 (I, J) = 0. 171.8 141 Z2(I,J)=0. 172 DGRID= (XY (2,1,2)-XY (2, 1, 1) ) /2 172-1 JGRID=1 172.2 1000 IF (CA-GT. XY (2 ,1 , JGRID) +DGRID) GO TO 1010 172-3 GO TO 1018 172-35 1010 JGRID=JGBID+1 172-4 IF (J GRID. EQ-NIY) GO TO 1018 172.45 GO TO 1000 172.5 1018 IGBID=1 172.55 1 015 IF (ISOB (IGRID+1 , JGBID) - EQ.O) GOTO 1020 172-6 IGRID=IGBID+1 172-65 IF (IGBID. EQ. NIX) GO TO 1020 172.7 GO TO 1015 172.75 1020 DO 1030 1=1,NIX 172.8 XI=XY( 1,1,1) 172.85 IF (XI.LE.CB) GOTO 1025 .PROGRAM FOR TWO DIMENSIONAL PLOT 111 172.9 IF (XI -GT. CB -AND. ISUB (I, JGRID) -NE. .0) GO TO 1025 172.95 XI=XY (1 ,IGRID, JGRID) 173 1025 IF (XI - EQ. CB) GO TO 1027 173.05 A.INT2 (I) =DQUANK (DF2,CB,XI, . GO 1 DO ,TOL, FIFTH) 173. 1 GO TO 10 30 173.15 1027 AINT2 (I) =0. DO 173.2 1030 CONTINUE 173.4 DO 150 J=1,NIY 173.5 DO 150 1=1,NIX 173.6 IF (ISUB (I, J).EQ.O) GOTO 150 173.7 XI=XY(1,I,J) 1 73- 8 YI=XY(2,I,J) 173-9 AINT1=0.D0 174 IF (JGRID .EQ. J) GO TO 1060 174.1 JADD=0 174.2 IF (ISUB (I, JGSID) . EQ. 1) GO TO 1050 174.3 JINC=1 174.4 IF (JGRI D. GT. J) JINC=-1 174.5 1040 J ADD=JADD+J INC 174.6 IF (JGRID + JADD. EQ. J) GO TO 1060 174.7 IF (ISUB (I, JGRID+JADD) .EQ. 1) GOTO 1050 174-8 GO TO 1040 174.9 1050 AINT 1=DQ UANK (DF 1, XY (2,1, JGRID + JADD) ,XY (2,1, J) ,. 001 DO, TOL,FIFTH) 175 1060 Z1 (1+5, J+5) =SIGYOA—AIN.T1-AINT2 (I) 176.4 T AUX Y (I , J) =Z1 (1 + 5, J+5) 177 Z2 (1 + 5,J+5)=Z1 (1+5,J + 5)*(EXDT(I,J) -EYDT(I,J) ) /LANDA (I 178 # « i 150 CONTINUE 178.2 CALL PLOT2(TAUXY,XY,IX,IY,ISUB,NIX,NIY) 178.4 CALL PLOTND 178.6 STOP 179 CALL PERS (Zi ,IDIHP+10,NIX+10,NIY + 10 ,DXY,.333,4 5. ,45., 10., 10.) 179.5 CALL PLOT (12. ,0. ,-3) 180 CALL PERS{Z2,.IDIMP + 10,NIX+10,NIY+10 ,DXY,. 333 ,45. ,45. , 10., 10. ) 180.2 CALL PLOTND 181 STOP 182 200 STOP 1 183 END 184 SUBROUTINE AUX(X,Y,XX) 185 IMPLICIT REAL*8 (A—H,0—Z) 186 DIMENSION XX { 15) 187 XX(1)=1.D0 188 XX(2)=X 189 XX(3)=X*X 190 XX(4) =X*XX(3) 191 XX (5) =X*XX (4) 192 XX(6)=Y 193 XX (7) =Y*Y 194 XX (8) =Y*XX(7) 195 XX (9)=Y*XX (8) 196 XX ( 10) =X*Y 197 XX(1 1)=XX{10) *Y 198 XX ( 12)=XX( 1 1) *Y 199 XX (13)=XX (1 0) *X 200 XX { 14)=XX(13) *Y P R O G R A M F O R T S O D I M E N S I O N A L P L O T 1 1 2 2 0 1 X X ( 1 5 ) = X X £ 1 3 ) * X 2 0 2 R E T U R N 2 0 3 E N D 2 0 4 F U N C T I O N A U X 2 ( X X , Y Y , P , L L ) 2 0 5 I M P L I C I T R E A 1 * 8 ( A - H , 0 - Z ) 2 0 6 C 2 0 7 C E V A L U A T E T H E F I T T E D F U N C T I O N A T X X , Y Y 2 0 8 c P C O N T A I N S T H E F I T T E D P A R A M E T E R S 2 0 9 c B E L L I S T R U E I F T H E I N D E P E N D E N T V A R I A B L E M U S T E V A L U A T E D B Y A U X 2 1 0 C 2 1 1 L O G I C A L L L 2 1 1 - 2 COMMON X 2 1 2 D I M E N S I O N X { 1 5 ) , P ( 1 ) 2 1 3 I F ( L L ) C A L L A U X ( X X , Y Y , X ( 1 ) ) 2 1 4 A U X 2 = 0 . D O 2 1 5 D O 1 0 1 = 1 , 1 5 2 1 6 A U X 2 = A U X 2 + P ( I ) * X ( I ) 2 1 7 1 0 C O N T I N U E 2 1 8 R E T U R N 2 1 9 E N D 2 2 0 F U N C T I O N D F 1 ( Y C ) 2 2 1 I M P L I C I T R E A L * 8 ( A — H , Q — Z ) 2 2 2 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 ) 2 2 2 . 5 C O M M O N / C 1 / R H O , Y O L D { 1 5 , 2 ) , D T 2 2 3 c 2 2 4 c T H E I N T E G R A N D D T A U / D X I S E V A L U A T E D 2 2 5 c 2 2 6 E B R = A U X 2 ( X I , Y C , P E B R , - T R U E . ) 2 2 7 C A L L D E R I V ( Y ( 1 , 1 ) , X I , Y C , E X D T , D U D Y , 3 ) 2 2 8 C A L L D E R I V ( Y { 1 , 2 ) , X I , Y C , D V D X , D V D Y , 3 ) 2 2 9 G A M D T = D U D Y * D V D X 2 3 0 E B R D T = 2 . D 0 / 3 . D 0 * D S Q R T ( 3 . D 0 * E X D T * * 2 + - 7 5 D 0 * G A M D T * * 2 ) 2 3 1 C A L L D E R I V ( P E B R , X I , Y C , D E B R D X , D U M , 1 ) 2 3 2 C A L L D E R I V 2 ( Y ( 1 , 1 ) , X I , Y C , D U D X Y , 3 ) 2 3 3 C A L L D E R I V 2 ( Y ( 1 , 2 ) , X I , Y C , D V D X X , 1 ) 2 3 4 D G A M D X = D U D X Y + D V D X X 2 3 5 C A L L D E R I V 2 ( Y { 1 , 1 ) , X I , Y C , D E X D X , 1 ) 2 3 6 D E B R D T = ( 4 . D O * E X D T * D E X D X + G A M D T * D G A M D X ) / ( 3 . D 0 * E B R D T ) 2 3 7 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 M D T / E B R 2 3 8 1 + D G A M D X - D E B R D T * G A M D T / E B R D T ) 2 3 8 . 1 U = A U X 2 ( X I , Y C , Y ( 1 , 1 ) , . F A L S E . ) 2 3 8 . 2 V = A U X 2 ( X I , Y C , Y ( 1 , 2 ) , . F A L S E . ) 2 3 8 . 3 V O L D = A U X 2 ( X I , Y C , Y O L D ( 1 , 2 ) , . F A L S E . ) 2 3 8 . 4 F A C T = R H O * , ' ( V - V O L D ) / D T ^ K T " ~ ,:' 1 '\ D F 1 = D F 1 + F A C T ~ " " ' • 2 3 8 . 5 2 3 9 . R E T U S N 2 4 0 E N D 2 4 1 F U N C T I O N D F 2 ( X ) 2 4 2 I M P L I C I T R E A L * 8 ( A - H , 0 - Z ) 2 4 3 C O M M O N 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 ) 2 4 3 . 1 C O M M O N / C 1 / R H O , Y O L D ( 1 5 , 2 ) , D T 2 4 3 . 2 R E A L * 8 L A M B D A 2 4 4 c 2 4 5 c T H E I N T E G R A N D D T A U / D Y I S E V A L U A T E D 2 4 6 c 2 4 7 E B R = A U X 2 ( X „ C A , P E B R , . T R U E . ) 2 4 8 C A L L D E R I V ( Y ( 1 , 1 ) , X , C A , E X D T , D U D Y , 3 ) PROGRAM FOR TWO DIMENSIONAL PLOT 113 249 CALL DERIV (Y (1,2) , X,CA, DVDX.EYDT , 3) 250 GAM DT=DU DY*DVDX 251 EBBDT=2.D0/3.DO*DSQRT(3. DO *EXDT**2 + . 75D0*GAMDT** 2) 252 CALL DERIV {PEBR , X,CA ,DUM ,DEBRD Y, 2) 253 CALL DERIV2 (Y (1, 1) ,X ,CA, DUDYY, 2) 254 CALL DERIV2 (Y {1,2) , X ,CA, DV DXY , 3) 255 DGAM DY=DUDYY +DVDXY 256 CALL DERIV2 (Y (1 , 1) ,X,CA, DEXDY,3) 257 YE8RDT= ( 4. DO*EXDT*DEXDY+GAMDT*DGAMD Y) / (3. D0*EBRDT) 258 CALL DERIV2 {Y (1 , 1) ,X,CA,DEXDX , 1) 259 CALL DERIV2 (Y (1,2) , X ,CA , DE YDX , 3) 260 LAMBDA= { 2. D0*C 1 *EBR**C2) / (3-D0*EBR0T) 261 CALL DERIV (PEBR,X,CA, DEBRDX, DOM, 1) 263 CALL DEBIV2{Y(1,2),X,CA,DVDXX,1) 264 DGA MDX=DEXDY + DVDXX 265 DEBR DT={4.DO*EXDT*DEXDX+GAMDT*DGAMD X)/(3.D0*EBRDT) 266 DF2=LAMBDA*({(EYDT-EXDT)*{—C 2 * DE BSD X/E BR 267 1 +DEBRDT/EBRDT)+DEXDX-DEYDX) 268 1 +.5D0*(C2*DEBRDY*GAMDT/EBR + DGAMD Y - YEBSDT*GAMDT/E SRDT) ) 268. 1 O=A0X2(X,CA,Y{1, 1) ,.FALSE.) 268.2 V=A0X2 (X,CA,Y (1 ,2} ,. FALSE. ) 268.3 OOLD=AOX2(X,CA,YOLD(1,1) ,.FALSER . 268-4 F ACT=RHO*\f<t{0—UOLD) /DT i*>U* " * Y> VJ 268.5 DF2=DF2 • FACT ' ' " v~ -------- - a - - " ' 269 RETURN 270 END 271 SUBROUTINE DERIV (A, X, Y, DU DX , DU DY , N) 272 C 273 C EVALUATE DERIV WRTX AND Y 274 C IF N=1 DODX, IF N=2 DODY, OTHERWISE DUDX &N D DUD Y 275 C 276 IMPLICIT REAL*8(A-H,0~Z) 277 COMMON XX(15) 2 78 DIMENSION A {1) 279 IF (N .EQ. 2) GO TO 10 280 DUDX=A(2) +2.D0*A (3) *X + 3-D0*A (4) *XX (3) + 4.D0*A(5)*X X{4) 281 1 +A (10) *Y+A (1 1) *XX(7) + A(12)*XX(8) + 2.D0*A (13) *XX (1 0) 282 1 +2. D0*A ( 14) *XX(11) + 3-D0*A (15) *XX (13) 283 IF (N.EQ. 1) BETUBN 284 10 DUDY=A(6) + 2.D0*A (7) *Y + 3.D0*A (8) *XX (7) + 4.D0*A(9)* XX (8) 285 1 •A(10)-*X + 2. D0*A (1 1)*XX (10) + 3. DO *A { 12) *XX (11) 286 1 + A(13)*XX(3) + 2.D0*A (14) *XX {13) + A(15)*XX(4) 287 BETUBN 288 END 289 SUBROUTINE DERIV2(A,X,Y,DODD,N) 290 IMPLICIT REAL*8(A-H,0-Z) 291 C 292 C EVALUATE 2ND OBDEB -DERIV WRT X AND Y 293 C N=1 DIDX; N=2 DY DY; N=3 DXDY 294 C 295 DIMENSION A (1) 296 GO TO (10,20,30) ,N 297 10 DUDD=2. D0*A (3) + 6.D0*A(4)*X • 1 2. D0*A (5) *X*X PBOGRAM FOR TWO DIMENSION Al PLOT 1 1 4 298 1 +2. D0*A(13)*Y +2.D0*A{14) *Y*Y + 6. D0*A(15)*X*Y 299 RETURN 300 20 DUDD=2.D0*A(7) • 6. D0*A (8) *Y + 12. D0*A (9) *Y*Y + 301 1 2.DQ*A(11)*X • 6. D0*A(12) *X*Y + 2. DO*A (1 4) *X*X 302 RETURN 303 30 DUDD=A(1Q) + 2. D0*A(11)*Y • 3. DO* A (12) *Y* Y 304 1 + 2.D0*A(13)*X *• 4. D0*A (1 4) *X*Y + 3. DO *A (15) *X*X 305 RETURN 306 END 307 SUBROUTINE Fill(XYT,IX,IY,ISUB,XY,NIX,NIY,IHY,NTM,ID I M,IDIMP) 308 IHPLICIT REAL*8 (A-H,0-Z) 309 DIMENSION X YT (2,IDXM , IDIM, 2) , I SUB (I DI MP , 1) ,XY(2,IDIMP ,D 310 C 311 C ISUB CONTAINS 1 IBE.RE A FUNCTION VALUE IS 312 C PLOTTED, 0 IF OUTSIDE BOUNDARY 313 C 314 DO 10 1=1,NIX 315 DO 10 J=1,NIY 316 10 ISUB(I,J)=1 317 XMIN=XYT( 1,IMY,IX,NTM) 318 DO 15 1=1,IY 319 IF (XYT ( 1,1, IX, NTM) .GT. XMIN) GO TO 15 320 XMIN=XYT (1,1,IX, NTM) 321 IHYMIN=I 322 15 CONTINUE 323 DO 1 10 J=1,NIY 324 DO 100 1=1,NIX 325 IF (XY(1,I,J) -IT. XYT (1,IMYMIN,IX ,NTM) ) GO TO 100 326 IF (XY(1,1,J).GT.XYT { 1,IMY,IX,NTM)) GO TO 80 327 IB=1 328 NB=2 329 20 IF(XY(2,I,J) -IT. XYT (2,NB ,IX,NTM) ) GO TO 30 330 1B=1B+1 331 N B= N B + 1 332 IF (NB.IT.IY) GO TO 20 333 30 IF(XY(1,I,J) .IT. XYT (1 ,1B ,1 X, NTM) .AND. 334 1 XY(1,I,J) .LT. XYT (1 ,NB,IX,NTM) ) GO TO 100 335 IF (XY(1,I,J) .GT-XYT(1,LB,IX,NTM).A ND. 336 1 XY(1,I,J).GT.XYT (1 , NB, IX, NT M) ) GO TO 80 337 YN=XYT(2,LB,IX,NTH) + (XYT(2,NB,IX,NTM)-XYT(2,LB,IX,NTM ) ) 338 1 /(XYT(1,NB,IX,NTM)-XYT(1,LB,IX,NTM) )* 339 1 (XY (1,1, J) -XYT (1, LB, IX, NTM) ) 340 IF (YN.GT.XY(2,I,J).AND.XYT{1,LB,IX,NTM).GE. 341 1 XYT(1»NB,IX,N3M)) GOTO 100 342 IF (YN.IT. XY (2,1, J) .AND. XYT (1 ,IB,I X,NTM)-IE. 343 1 XYT (1,NB,IX,NTM) ) GOTO 100 351 80 DO 90 11=1, NIX 352 90 ISUB (II , J) =0 353 GO TO 110 354 100 CONTINUE 355 110 CONTINUE 356 RETURN 357 END 358 SUBROUTINE P10.TIT (XYT ,IX,IY) 359 IMP1ICIT REAL*8 (A-H,0-Z) PROGRAM FOR TWO DIMENSIONAL PLOT 1 1 5 360 DIMENSION XYT (2,20,20) 361 CALL AXIS (0.,0.,'X1,- 1,10- ,0.,0.,.2) 362 CALL AXIS (0, ,0- , ' Y • , 1 , 10. ,90- ,0- ,-2) 363 DO 10 J=1,IY 364 CALL PLOT (XYT { 1, J, 1) *5- , XYT (2,J,1)*5.,3) 365 CALL SYMBOL (XYT(1,J,1)*5.,XYT (2,J,1)*5.14,4, 0-,-1) 366 DO 10 1=2, IX 367 CALL PLOT (XYT (1 , J,I) *5- , XYT (2, J, I) * 5. , 2) 368 CALL SYMBOL (XYT (1, J , I) *5. , XYT (2 , J,I) *5.,-14,4, 0. ,-1) 369 10 CONTINUE 370 DO 20 J=1,IX 371 CALL PLOT(XYT(1, 1,J) *5. , XYT (2 , 1 , J) *5 . , 3) 372 DO 20 1=1,IY 373 CALL PLOT(XYT(1,I,J)*5./XYT (2 , 1 , J) *5.,2) 374 20 CONTINUE 375 CALL PLOT (12.,0. ,-3) 376 RETURN 377 END 378 SUBROUTINE PLOT 2 (TX Y , XY, IX , IY , ISUB, NIX , NIY) 379 REAL*8 XY(2,80,80) 3 80 DIMENSION TXY (80,80) ,ISUB(80,80) /TAUXY (80,80), X (80) 381 DO 5 1=1,80 382 DO 5 J=1 ,80 383 5 TAUXY (I, J) =0.D0 384 DO 10 1=1,NIX 385 DO 10 J=1,NIY 386 10 TAUXY (I, J) =TXY (I, J) 387 CALL SCALE(TAUXY,6400, 1Q.,YMIN,DY,1 ) 388 CALL AXIS (0.,10., »X»,1,10.,0.,0.,.2) 389 CALL AXIS (0. ,0. * ,5, 10. ,90. ,YMIN,DY) 390 NY=NIY/IY*3 3 90-2 KK=0 391 DO 50 J=1,NIY,NY 392 L=NIX 393 15 IF(ISUB (L, J) . EQ. 1) GO TO 20 394 L=L-1 395 GO TO 15 396 20 DO 30 1= 1 ,L 397 30 X(I) =XY ( 1/1, J) *5. 398 CALL LINE (X ,TAUXY ( 1 , J) , L , 1) 399 NX=NIX/IX 400 KK=KK+1 401 DO 40 K=1,L,NX 403 CALL SYMBOL (X (K) ,TAUXY(K, J) ,. 14, KK, 0.,-1) 4 04 40 CONTINUE 4 05 50 CONTINUE 407 END 

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}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0100192/manifest

Comment

Related Items