nun ‘ .v mu: v ~71 ind! "(4‘60 "ucoI-\- w. 5.; {THESIS h LIB R in KY I Michigan: State 1 University This is to certify that the thesis entitled COUPLING OF THE FINITE ELEMENT AND BOUNDARY ELEMENT METHODS FOR THE SOLUTION OF PLANE PROBLEMS OF LINEAR ELASTICITY presented by Jose Jamie Garcia has been accepted towards fulfillment of the requirements for M.S. Mechanics degree in Major prof ssor Date 12/7/84 0-7639 MS U is an Affirmative Action/Equal Opportunity Institution MSU LIBRARIES 1—P- RETURNING MATERIALS: Place in book drop to remove this checkout from your record. FINES wiIl be charged if book is returned after the date stamped below. -M" COUPLING OF THE FINITE ELEMENT AND BOUNDARY ELEMENT METHODS FOR THE SOLUTION OF PLANE PROBLEMS OF LINEAR ELASTICITY By Jose Jaime Garcia A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Metallurgy, Mechanics and Materials Science 1984 ABSTRACT COUPLING OF THE FINITE ELEMENT AND BOUNDARY ELEMENT METHODS FOR THE SOLUTION OF PLANE PROBLEMS OF LINEAR ELASTICITY By Jose Jaime Garcia A. A hybrid finite element-boundary element formulation for the solution of elastostatics plane problems is presented. Only isotropic materials are considered. A finite element program, FEM, which uses triangular and quadrilateral isoparametric elements, is initially 'developed. Next, a direct boundary element program, BEM, is created using a piecewise straight boundary discretization. Constant tractions and linear displacements are assumed on each boundary element. This model allows traction discontinuities at the nodes and is compatible with the elements used in FEM. Finally, a mixed formulation is obtained in the program FEMBOU, which uses standard finite elements and multi-node elements or "superelements", whose symmetric stiffness matrices, created using the boundary integral equations, are incorporated into a standard finite element program, similar to FEM. Results of the three programs are compared for two examples. ACKNOWLEDGEMENTS The author wishes to express his gratitude and deep appreciation to his advisor, Professor Nicholas J. Altiero, for his able guidance and sincere encouragement during the course of this thesis research, and for his spending many hours with the author in the preparation of this thesis. Special thanks are due to General Dynamics Land Systems Division, Warren,- Michigan (Contract No. DEY-600484, Print. Investigator: N. Altiero) for its support of this research. The author also wishes to thank The Universidad del Valle, Cali, Colombia and The Fulbright Commission for their support of his studies in the United States during the last year. Finally, the author wishes to acknowledge ‘his parents for’ their continuous encouragement and ‘moral support in the past several years. ii rd (W LIST OF TABLES TABLE OF CONTENTS LIST OF FIGURES O O I O O 0 O O O O O O O O O O O O 0 CHAPTER I CHAPTER II CHAPTER III CHAPTER IV CHAPTER V APPENDIX BIBLIOGRAPHY. INTRODUCT I ON 0 O O I O O O O O O O O O FINITE ELEMENT MODELLING . . . . . . . II.1 DESCRIPTION OF THE FINITE ELEMENT METHOD II.2 II.3 BOUNDARY ELEMENT MODELLING . . . . . ISOPARAMETRIC QUADRILATERAL ELEMENTS . PLANE STRESS FINITE ELEMENT PROGRAM III.1 DESCRIPTION OF THE BOUNDARY ELEMENT METHOD III.2 PIECEWISE LINEAR ELEMENTS. . . III.3 PLANE STRESS BOUNDARY ELEMENT PROGRAM. . COUPLING THE FINITE ELEMENT AND BOUNDARY ELEMENT METHODS O O O O O O O O O O O O C O O 0 IV 0 l PROCEDURE 0 O O O O I I O O IV.2 PLANE STRESS HYBRID PROGRAM. EXAMPLES AND DISCUSSION. . . . . . . . V.I EXAMPLES . . . . . . . . . v.2 DISCUSSION OF RESULTS. . . . v.3 RECOMMENDATION FOR FUTURE STUDY. COMPUTER PROGRAM LISTS . . . . . . . iii .11 .14 .44 .44 .48 .58 .72 .72 .75 .83 .83 100 102 104 151 I..- ‘ a... V? 7'? ct. TABLE II.l III.1 V.l v.2 V.3 LIST OF TABLES Array NTBC(I,K) before and after PROFIL . . . . . . . Singular Integrals. . . . . . . . . . . . . . . . Vertical deflection of A and CPU time for cantilever beam prOblem. C O O O O O O O O O O O O O O O O O O 0 CPU time for problems of Figure V.5 . . . . . . . . . Number of equations and components of coefficient matrices for problems of Figure V.5 . . . . . . . . . iv .88 .99 .99 Figure 1.1 11.1 11.2 11.3 11.4 11.5 11.6 11.7 11.8 11.9 11.10 11.11 11.12 11.13 11.14 111.1 111.2 111.3 111.4 111.5 LIST OF FIGURES Mixed boundary value problem. . . . . . . Domain discretization for Finite Element Analysis . . Local and global coordinates for a four-node isoparametric element . . . . . . . . . . . Flow chart of Program FEM . . . . . . . . Example for PROFIL. . . . . . . . . . . . . Matrices [K]. [K] for the example ofl Figure 11.4.. . . . . . Arrays B(J) and IDIAG(J) for the example of Flow chart of PFORM . . . . . . . . . . . . Flow chart of ELEMENT . . . . . . . . . . . Location of elements of [K] in A(I) . . . . Location of elements of [K]12 in B(I) . . . Illustration of the first part of SOLVE . . Flow chart for first part of SOLVE. . . . Flow chart for second part of SOLVE . . . . Example problem to demonstrate input data fo Plane boundary value problem. . . . . . . . . Angles<1andcu of eq. (III.2). . . . . . . . Boundary discretization . . . . . . . . . . “(m) Definition of a1, bi’ x . . . . . . Flow chart of program BEM . . . . . . . . . and arrays A(J) and JDIAG (J) Figure II. r FEM . .12 .15 .20 .22 .25 .26 .28 .30 .32 .34 .36 .39 .42 .45 .45 .49 .49 .59 Figure 111.6 111.7 111.8 111.9 111.10 111.11 1V.1 1V.2 IV.3 1V.4 V.l v.2 v.3 v.4 v.5 V.6 V.8 V.9 v.10 v.11 v.12 v.13 Flow chart of part 4 of the program . . . . . . Singular contributions in [UC] and [UR] . . Non-singular contributions in [UC] and [UR] Flow chart for the modification of [UR] . . Flow chart for rearranging equations into final form. Example problem to demonstrate input data for BEM . Combined discretization . . . . . . Interaction between BOUN and FEMBOU . . . Flow chart for BOUN . . . . . . . . . . . . Example problem to demonstrate input data for FEMBOU. Cantilever beam . . . . . . . . .1. . . . . l3 node boundary discretization of cantilever beam. . . 19 node boundary discretization of cantilever beam. Stress distribution on section a—a, located free end 0 O O O C O O D O O O O O O O O O 0 Example No. 2 . . . . . . . . . . . . . . . BEM model for R/w = 0.5. . . . . . . . . FEM model for R/w = 0.5. . . . . . . . . . FEMBOU model for R/w a 0.5 . . . . . . . . BEM model for R/w = 0.2. . . . . . . . . . FEM model for R/w a 0.2. . . . . . . . . . FEMBOU model for R/w = 0.2 . . . . . Stresses 022 on section a—a for example 2, R/w Stresses 022 on section a-a for example 2, R/w vi 2.833 from .64 .67 .68 .70 .73 .76 .77 .80 .84 .85 .87 .89 .90 .91 .92 .93 .94 .95 .96 .97 CHAPTER 1 INTRODUCTION The classical mixed boundary value problem of linear elastostatics, shown in Figure 1.1, consists of an elastic body, R, loaded by specified tractions, t on part Bt of the boundary and specified displacements, i, u , on the remainder of the boundary Bu' Body forces are neglected here. i The displacement and stress fields everywhere in R, subject to the boundary conditions, are sought. Exact solutions have been obtained for many problems in which simple geometry is involved, such as rectangular, circular, elliptical bodies with particular boundary conditions. Nevertheless, it is quite difficult to obtain the exact solution for a problem in which either the body shape is more complicated or the load is more general. Under these circumstances, numerical techniques have been developed to approximate 'the solution of any general elastostatics problem. The most notable of these procedures are the Finite Element Method and the Boundary Element Method, which have become practical only after the development of the high-speed digital computer. The Finite Element Method [1,2,3],* first used for solving structural problems in the 1960's, has become a powerful engineering tool which can be applied to a wide range of problems such as fluid flow, electromagnetic fields, etc. The method consists of dividing the continuum into a set of imaginary discrete elements or "finite elements", which are defined to interact with the neighboring elements only at 1 * Brackets refer to references listed in BIBLIOGRAPHY Figure 1.1. / Mixed boundary value problem 3 specified points on their boundaries, called "nodes". Governing field equations are approximated within the elements by prescribed interpolating functions .or "shape functions". Initial engineering intuition was that such a model should produce results similar to those of the real problem, a fact which has since been amply proven by mathematicians. In the case of a three dimensional problem, it is necessary to handle a large amount of data for specification of elemental and nodal information of the discretized model. This is a major drawback of the method, especially when a large number of small elements is required to accurately predict high stress gradients. The Boundary Element Method [4,5], introduced for elasticity problems in the late 1960's and develOped through the 1970's, solves _governing integral equations by dividing the boundary of the body into a set of small segments or "boundary elements", on which certain properties (tractions, displacements, etc) are approximated by interpolating func- tions or "shape functions". In the case of the "direct" method for elastostatics, Somigliana's integral equation is discretized. At least one exact singular solution must be known in order to implement the approximate equations. The dimensional discretization is one order less than that for the Finite Element Method, this being an important advantage in particular when solving three dimensional problems. Moreover, field equations are satisfied exactly within the body, allowing accurate calculation of field properties, even if there are high gradients involv- ed. Furthermore, the Boundary Element ‘Method can. be implemented to handle problems over infinite or semi-infinite domains. Due to the fact that boundary data is discretized, problems arise for the method when solving cases in which the ratio of boundary area to domain volume becomes large, i.e. when most of the domain lies in the proximity of the boundary. In the hope that the difficulties of the two methods can be avoided, a mixed Finite Element - Boundary Element formulation has often been suggested and actually implemented using various approaches, [3,4,6,7,8,9,10]. In this dissertation, programs based on the Finite Element and Boundary Element methods are implemented, with the goal of using them to produce a third "combined" formulation, which is subsequently presented. Additionally, results of the three programs are compared for some example problems. The Finite Element Program, FEM, described in CHAPTER II, follows standard procedures for saving the global stiffness matrix, solving the equations, etc., as recommended by various Finite Element authors, [1,2,3,10]. Numerical integration is used to create elemental stiffness matrices. The extension to the three dimensional problem can be easily accomplished. The Boundary Element program, BEM, described in CHAPTER III, is based on the direct approach. The boundary discretization consists of straight segments in which displacements vary linearly and tractions are constant. This model allows traction discontinuities and is compatible with the finite elements used in FEM. The hybrid program, FEMBOU, described in CHAPTER IV, consists of the program, FEM, plus a subroutine, BOUN, which creates the stiffness matrix of a "superelement", in which the boundary integral equations are valid. In the final chapter of the work, results of the three programs are presented and compared with regard to accuracy and computational speed. 5 Some conclusions can be drawn from these results but no definitive statement can be made as to which is the "best" method. This work is just a beginning toward the goal of implementation of the methods to the general three dimensional anisotropic case. CHAPTER II FINITE ELEMENT MODELLING 11.1 DESCRIPTION OF THE FINITE ELEMENT METHOD The displacement approach, or direct formulation, of the Finite Element Method, consists of dividing the body into imaginary "finite elements" as shown for a plane problem in Figure 11.1. Each element is connected to neighboring elements at specified points called nodes. The components u 1, u2 of the displacement at a point within an element e is given approximately by n (i) n (1) u (x)- Z N (x)u , u (x)= r N (x)u (11.1) 1 i 1 2 H i 2 i=1 i=1 where ulu), uzu) are the components of the displacement of node i, N1(x) is a specified function corresponding to that node, or. "shape function" of the node i, and n is the number of element nodes. There are as many shape functions as nodes in the element. These shape functions must satisfy certain requirements, i.e., they are equal to l at the corresponding node and are equal to 0 at all other nodes of the element. Also, they must be able to reproduce exactly a state of constant stress [2]. Figure 11.1. Domain discretization for Finite Element Analysis 8 If the nodal displacements are known, the stress and strain can be calculated in the element. For a plane problem, the components of strain are \ E ' (5'3 gullaxl le/‘ax1 0 ...... aNn/3x1 0 61(1)\ u2(1) {a}: 522 =8u2/3x2 I 0 3N1/3x2......0 BNn/ax2 . < fl . (n) u1(n) 2B,: 311 /3x +311 /3x 3N /8x 3N /3x ....3N /3x 3N /3x u \JK122QL1211 n2n1]\2J or {e} =[B B B ]{u}(e) (II 2) 1 2.... n . where aNi/ax1 0 [B - 0 3N1/3x2 BNi/sz BNi/Bx1 9 )(e) and { u is the vector of nodal coordinates of the element. With the elastic constants of the material, Young's Modulus E and Poisson's Ratio u, the components of stress are (011W - 1 u 0 1 (£111 {01- 'a == E u 1 o e L = IDIIel=IDllBl {u}(e) i 22 L .1752. y 22 012 0 0 (1:3) 2e: L J L 2 J L 12; (11.3) in the case of plane stress. If the potential energy of the element, n(e)-%IVeioT{cidv-{F}(°)Tiu1(e) . or n‘e)-a{u}‘e)T IVEIBJTID][Bldv{u}(e)-IF}(e)TIU}(e) (II 4) (e) is minimized, a relationship between the nodal forces ,{F} }(e) nodal displacements {u is obtained. This relationship is [K] (8){ u} (e)..{P} (e) , and the ' 10 ](e) where the matrix [K is the stiffness matrix of the element e,i.e. [xi(e)=.ge[B]T[o][B]dv (11.5) The global stiffness matrix which relates the global displacement vector {u} to the global force vector {F}, is also obtained by minimizing the total potential energy of the body. The final system of equations can be written as, [K]{u}={F} (11.6) where {u} is a vector containing all nodal displacements, {F} is a vector containing externally applied nodal forces, and [K] is syumretric and banded. The entries of [K] are composed of contributions of the elemental stiffness matrices. The final system of equations can be solved for the 'unknown. displacements, and subsequently the state of stress and strain can be calculated at all points of the body using eqs. (11.3) and (11.2), respectively. The general algorithm of a finite element program is as follows: - Read the data of the discretized model. - Form the global stiffness matrix from the contributions of all elemental stiffness matrices. - Modify the system of equations for specified displacements. - Calculate displacements and stresses at prescribed points in the body. 11 \ 11.2 ISOPARAMETRIC QUADRILATERAL ELEMENTS In the case of an isoparametric element, the relation between the local and the global coordinate system is given by, n n , (i) . (i) where Ni(€,n) is the shape function of node i. In the case of a four-node quadrilateral, these can be represented by the equation, Ni(€ .n)=lr(1+€€i)(1+nni) (II.8) where i.endrli are the local coordinates of node i, as shown in Figure 11.2. 'For example, the local coordinates of node 2 are 72-1 and n2” -1 ‘9 so that Nzta.n)=k(1+€)(1-d) Values of the displacements in the element are related to those at the nodes by the same shape functions (thus the description "isoparametric") so that n u1(x)= z Ni(i,n)u1(i) , u2(x)= N.(S,Tl)u “’ (11.9) 181 1 l 2 II (‘13 1 where n is the number of nodes in the element. 12 X2 l ‘ n acx‘f’. x? I g 1'l . l 4. 'J > (11.13) [“21 J “('22.- 3“)” \{F}nsj .l h where {u}us is a vector of unknown displacements, {u}s is a vector of specified displacements (zero and non-zero), {F}s is a vector of known external forces, {F} ns is a vector of unknown external forces or reactions, NEQ is the number of equations to be solved and NKD is the number of specified displacements. The number of equationsis calculated as NEQ - (No. of nodes)x(No. of degrees of freedom)- NKD and the equations to be solved can be written as [KIllIuIns‘IFIS-[K]12{UIS=IF};’ (11.14) If the specified displacements are all zero, {F}S* = {F}S and it is not necessary to save [K] It is noted that the submatrix [K]11 is 12' symmetric and banded. 19 The objective of the subroutine PROFIL is to calculate the dimensions of two vectors, A(J) and B(J), in. which parts of [K]11 and [K]12 respectively are to be saved by columns, and to set two arrays, JDIAG(I) and IDIAG(1), that serve to identify positions of the elements of [K]11 and [K]12 in the vectors A(J) and B(J). Space is reserved in A(J) for all column entries in [K]11 from the first non-zero entry to the diagonal. With this organization, considered by Taylor in [3], less space is required for saving the matrix. On the other hand, some additional time must be spent in order to set the arrays JDIAG(I) and IDIAG(I). 11.3.2.1 FIRST PART OF PROFIL For each specified force condition, NTBC(1,K)-0, an. equation is counted and the number of that equation is assigned to NTBC(I,K). For each specified non-zero displacement, NTBC(I,K)-1, the specified non-zero displacement is counted and the negative of that number is assigned to NTBC(I,K). For each specified zero displacement, NTBC(I,K)=-1, NTBC(I,K) is set equal to zero. For the problem of Figure 11.4, the values of NTBC(I,K) before and after PROFIL are shown in Table 11.1. In this example there are 2 specified non-zero displacements and 3 specified zero displacements. 11.3.2.2 SECOND PART OF PROFIL For the same example of Figure 11.4, the form of the complete stiffness matrix is shown in Figure II.5(a), where the positions above the upper profile line are zero. The form of the submatrix [K]11 is shown in Figure II.5(b) and the positions marked as x are the ones to be saved in A(J) as is shown in the same figure. 20 F Specmed nonzero displacement I---o 1 1 2 ( ) (2) 4 3: =4 "‘ (a) Figure 11.4. Example for PROFIL 21 Table 11.1 Array NTBC(I,K) before and after PROFIL Node No. Degree of NTBC(1,K) NTBC(I,K) 1 Freedom, K Before PROFIL After PROFIL 1 0 1 1 2 -1 0 1 0 2 2 2 0 3 1 0 4 3 2 1 -1 1 1 -2 4 2 0 5 1 -1 0 5 22 NtIc(I.K)—’.(1) (0) (2) (3) (4) ('1) ('2) (5) (0) (0) uracum) NOD! women—«- 1 2 3 4 5 Moosnumsa DIRECTION —.- 1 2 1 2 1 2 1 2 1 2 DIRECTION J x x x * X 1 (1) 1 2 (0) t x x x * * x 1 (2) 2 x x * * x 2 (3) x * * x 1 '3 (4) 2 ('1) 1 (-2) 4 at: * x 2 (5) 1 (o) 5 2 (O) (a) A 1 Kn 2 K: 3» K”! i 1 x >< x x 5 “he __ 1 2 >< x x x 6 'Q” l 3 7 K24 I 6 3 x x x 8 K,‘ 9 4 x x 9‘ x. .L_F: 14 5 10 K. L x 11 K.. 13 K4: 1 14 1(5 (b) Figure 11.5. Matrices [K], [K]ll and arrays A(J) and JDIAG(J) for the example of figure 11.4 23 The array JDIAG(I) relates the entries of A(J) to those of [K]11. There are as many components in JDIAG(1) as there are equations. Each entry JDIAG(I) indicates the location in A(J) of the diagonal element of column I. The number of non-zero positions in a column I of [K]11 is equal to the maximum difference between the equation number I and the equation numbers corresponding to the degrees of freedom of the nodes adjacent to 1, plus 1, i.e. the number of elements in column 4 of [K]11 is equal to max {I4-3I , I4-2|}+ 1-2+1-3. The algorithm to set the array JDIAG(1) is as follows. - Do a loop on the elements. Do a loop on the nodes of the element. - Do a loop on the degrees of freedom. - Find the maximum difference between the equation I of that degree of freedom and the equation numbers. of the other degrees of freedom of the same element. Assign that value to JDIAG(I). - When the same equation number is found in another element, assign to JDIAG(I) the maximum 'value between the current maximum difference and the value already saved in JDIAG(I). - After finishing these 100ps, the entry JDIAG(I) is equal to the height of column 1 minus 1. - To obtain the final JDIAG(I) array, another loop is made to add the heights of all the columns. 111 the array B(J) are saved the columns of [K] corresponding to the Specified non-zero displacements. For the example of Figure 11.4, these are marked as asterisks in Figure 11.5. The number of entries in each Column is equal to the difference between the highest and the lowest e e quaticul number of the degrees of freedom with entries in that column, plus 1.. i.e. the number of entries in column No. 6 is (5-2)+1=4- There .. 24 are two entries in IDIAG(I) for each Specified non-zero displacement I. The first, IDIAG(ZI-l), indicates the equation number that corresponds to the first element of column I, and the second indicates the position in the array B(J) of the last element of column I. In Figure 11.6, are shown the arrays B(J) and IDIAG(I) for the example of Figure 11.4. The algorithm to calculate IDIAG(I) is as follows. - Do a loop on the elements. - Do a loop on the nodes of the element. - Do a loop on the degrees of freedom. - For each negative NTBC(I,K), find the highest and the lowest values of the equation numbers of the other positive NTBC(1,J) corresponding to the other degrees of freedom of the same element. Save these values in IDIAG(ZI-l) and IDIAG(ZI). - Add the heights of the columns to obtain the final array. 11.3.3 SUBROUTINE PFORM This subroutine fills the arrays A(J) and B(J) in which the submatrices [K]11 and [K]12 are saved. The flow chart for PFORM is shown in Figure 11.7. In the first part, PFORM calculates the entries of the matrix [D], which relates Stress to strain. These entries are saved in the arrays 12(1), PR(I) and G(I). TPhe second part of PFORM consists of a loop on the elements. For each element 1, the material properties are localized in D(l), D(2), D(3). the coordinates of the nodes in XE(1,L) and the boundary condition types Eat the nodes in NBC (1,L). Also the local stiffness matrix S(I,L) is initialized. With the arrays D(I) and XE(I,L), the stiffness of the element is calculated by the subroutine ELEMENT. With the arrays 25 B 1 “2(4) 2 “3(4) IDIAG 3 x...) 2 4 K 4..) 4 5 K1(2) 1 6 K212) 9 7 Kat-2) a K112) 9 K5“) Figure 11.6. Arrays 8(3) and IDIAG(J) for example of Figure 11.4 26 (: START > Loop on Materials j>——-——- I = l,NMAT E(I) H m PR(I) II rn II m C(I) Loop on elements :>_____T N = l,NEL i Localize: D(L), XE(I,L), NBC(I,L). Initialize: S(I,L) Calculate S(I,L) CALL ELEMENT I Add S(I,L) to A or] 8. CALL ADDSTIF (:CUNTINUE::>A C Li. D Figure 11.7. Flowchart of PFORM 27 NBC(1,L), JDIAG(I), IDIAG(1), the entries of the local stiffness matrix are located in A(J) or B(J) by the subroutine ADDSTIF. I I . 3 . 3 . 1 SUBROUTINE ELEMENT This subroutine calculates the stiffness matrix of three-node triangular and four-node quadrilateral elements by numerical integration. The flow chart of ELEMENT is shown in Figure 11.8. In the first part, ELEMENT calculates the entries W11, W12, W21, W22 of eq. (11.11). Each of these entries is approximately equal to n 3'" " W12 [Ve(3Nj/3xl)(3Ni/3x2)dv L 1|(BNj/axl)(aNi/3x2)detlJldtdn N2 . . kil WkkdetlJl(BNj/oxl)(8N1/8x2)(ck,nk) (11.15) where N is the number of points of integration in each direction, Wkk are the Gauss weights and Ek’ nk are the local coordinates of the Gauss POints. The subroutine PGAUSS calculates Wkk’ 5k, nk, depending on the Value of N. The subroutine SHAPE calculates the derivatives of the shape functions at the points of integration. The summation from 1 to N is Performed in ELEMENT. In the second part of ELEMENT, the entries Wll, W12, W21, W22 are rearranged and multiplied by the material constants to obtain the components of eq. (11.12). Subroutine SHAPE returns the values of the global derivatives 3Ni/8xl, aNi/ax2 in the array SHP(I,I(), where 1 indicates the shape func- tion number and K the derivative,i.e., K-1->3Ni/8x1, K=2n8Ni/3x2, K-3->Ni. The Sequence of operations in SHAPE is essentially that of equation (II'IO). First the local derivatives are formed, then the Jacobian matrix [J] and its inverse [J]-'1 are calculated, and finally the -1 multiplication of [J] by local derivatives is performed. When it is 28 C... ) I Calculate Gauss points and weights. ' CALL PGAUSS <} L =1,LL2 ll ‘ Calculate 3N,/ax,, aN./ax, J at each point of integration CALL SHAPE l [ Form Wll, w12, w21, w22 J l \ (:EONTINUE#J/ l Rearrange and multiply by material constants i (:JRETURN::> Figure 11.8. Flowchart of ELEMENT 29 required to find the stiffness matrix of a triangular element, the shape functions N3 and N4 are added, and the same procedure is followed in all other steps. 11.3.3.2 SUBROUTINE ADDSTIF This subroutine stores the values of elemental stiffness matrices in either A(J) or B(J), depending on the boundary condition NBC(I,K) at each degree of freedom. The algorithm for ADDSTIF is, - Do a loop on the nodes of the element, 181, NEN(N). - Do a loop on the degrees of freedom, K=1, NDOF. - For each degree of freedom: a) If NBC(I,K) is zero, continue to the other degrees of freedom. b) If NBC(I,K) is negative, store the column(NDOF*(2*I-l)+K) of S(L,M) is B(J), except for the entries with an NBC(L,M) value also negative or zero. c) If NBC(I,K) is positive, store the first part of the column in A(J), except for those entries with an NBC(L,M) value negative or zero. Only the first part of the column (up to the diagonal term) is stored since the matrix is symmetric. The way each entry of [K] is located in A(I) is illustrated in Figure 11.9. For a given address (M,J) in [K], the number (JDIAG(J)-J) is calculated. This number gives the positions in A(I) in which the first entry of the column would be located if it were different from zero. The proper position is obtained by adding the row number M. 30 JDMGLH _ --- ‘ ./ /’ / /’ i\&l §Z>CEQEJ \ /’ J/’ & Figure 11.9. Location of elements of [K] in A(I) 31 The location procedure in B(I) is illustrated in Figure 11.10. For each NBC(L,N) negative, the position of the first element of the column J is located in B(I), this position being given by the number IDIAG(2*(J-1)+1). Next, a number of positions equal to the difference between the equation number M, and the equation number corresponding to the first non-zero element in the column IDIAG(2*J-1) must be added to set the correct location for the entry. 11.3.4 SUBROUTINE SOLVE This subroutine uses a variation of the Gauss elimination method called Crout algorithm. This algorithm is thoroughly discussed by R.L. Taylor in [3]. In this section, the important facts of the method and the subroutine are summarized. The first part of SOLVE, which works when the logical arguments AFAC and BACK are true and false respectively, factors the stiffness matrix [K]11 stored in A(I) into the form . 1 0.......0 ull ulz.....u1n 121 100000000 0 Uzz00000u2n {Inna- . . -- [Luv] 0 l 0 0 :n11n2”'°'1~ LO 0unfl The procedure which calculates the entries of [L] and [U] consists of performing the matrix: multiplication and equating the respective terms. Following this procedure, it is found that 32 J IDIAG(2(J-l))+l IDIAG(ZJ-I) Ll ‘ :1) M-IDIAG (2.1-!) J l-———g——-1 \\\\\ L\\ RQQN NEQ Figure II.IO. Location of elements of [KJl2 in B(I) 33 ulngkln lni’kni/“n 1-1 u. -k -21 u.9i#1 1:] ij m-l in mg j-I kij-gillimumj 11j’ ujj (1*j.j#1) 1 / , if [K]11 is symmetric. 11’ “11 “11 In the subroutine, the sequence in which the entries are calculated and saved is shown in the example of Figure 11.11. It can be shown that both [L] and [U] have the same profile as [K]ll. If k n-O, u 0, u =-1 -0 also; and, in general, 1 1 2n' 2n 21“1n n-O since all the elements umn of the same column, above u =0; if k n if kin-O, u are 0. i in Since [K]11 is stored in A(I) by columns, addresses must be found using the array JDIAG(1). Also, since only the profile part of [K]11 is saved, some attention must be paid to the performance of the operations, i.e. the summation -1 Z l m=1 LJ. imumj is on the elements of column 1 and j so that, to do this operation, it is necessary to know which is the shortest column. This is done with the function MINO, which finds the minor value between two integers. _ 1 kn k1: k): u -k 1‘22 km 11- n L k3: "11 '2) km :11 -:- '21 12 "n : “22 R23 “szn'hiuu k3: _ J l31=u13+ull hfuzfi‘uzz “33:19:51" Issue“: “:7 Figure 11.11. Illustration of 34 u" k!) k1: “12"“: R22 1‘23 ———- _ [- un‘kia [.1H '21 uszza‘hnuuL * “22 1 “11 '2) In, t’22 '32 “a: d tne first 4. part “n ":2 km 1‘22 1‘21 km d "13 “23 k3: of SOLVE 35 A summary of the algorithm for the first part of SOLVE is as follows, - Do a loop on the columns. - If the height of the column equals 1, continue with other columns. - If the height of the column equals 2, calculate the diagonal element. - If the height is greater than two, calculate i-l “ij'kij’millimumj from J82 through the diagonal. - On the diagonal, calculate j-l 11 fun “‘11 ’ “Jj‘kn‘m_111m“mj A flow chart of this algorithm is shown in Fig. 11.12. The second part of SOLVE, which works when the arguments AFAC and BACK are false and true respectively, performs two operations. It solves the systems [L]{y}-{F}* and [U][x}={y} by forward and backsubstitution respectively. In these systems, {F}* is the force vector formed in MODIFY, {x} is the vector of unknown displacements, and {y} is an intermediate vector. 36 Loop on equations J = l,NEQ Calculate JH 1 Loop on hei ht of elements umn , of column Compute diagonal term L.__CZONTINUE J..._ ~ BACK _ true i-I false y, =E’-Zl.. Y... en] < RETURN > Figure II.l2. Flowchart for first part of SOLVE 37 For a 3x3 matrix the first operation is as follows, 21 31 In general 32 } =1 F* y1 1 y2'F2 '121y1 * ys'F3 '(131y11132Y2)' The elements of {y} are saved in the same vector {F} . For the same 3x3 matrix, the back substitution is u11 u12 “13- rx1] [V17 0 u22 '123 i x2 L3 y2 L’ .0 ° “33. Ba» Us) In general 1+1 y - E u x 1 manro 1” m 1+1 x a = y./u - Z i u 1 i m=NEQ ii x3’y3/“33 X27(y2'“23x3)/“22 x1'(y1'(“12x2+“13x3))/“11 (u /u im )x ii m 38 or 1+1 2 1 x X's/u . m-NEQ Note that the values of x1 are in terms of uii and 1mi(m#i), so that it is not necessary to calculate umi for mfli. A flow chart of this procedure is shown in Fig. 11.13. 11.3.5 SUBROUTINE MODIFY This subroutine has two parts. The first part reads and prints specified non-zero loads or displacements, and the second part modifies the load vector if there are specified non-zero displacements. In the first part the subroutine reads (in free format) the node number, the direction and the value of each non-zero force or displacement. The loads are stored in the first NEQ positions of F(I), and the displacements in the remaining NSD positions of F(I). No ordering is necessary to specify this data. The subroutine stops reading when a zero value is given for the node number. The second part of MODIFY performs the operation‘ IF}S*={F}S-[K112{u}s (11.16) 39 Loop on equations I=l, NEQ Loop on equations I = l, NEQ Y 301 Xi Y; ‘Zlmi X... M =NEO T RETURN:) Figure 11.13. Flowchart for second part of SOLVE 40 where the dimensions of the matrix [K]12, stored in the array B(J), are NEQxNSD. The algorithm to do this modification is as follows. - Do a loop on the equations. - Do a loop on the specified non-zero displacements. - Locate each element in B(J), and perform the product. The array IDIAG(I) is used to locate positions in B(J). 11.3.6 SUBROUTINE STRESS This subroutine calculates and prints the stresses in each element and the reactions at the points of specified displacement. The algorithm of the subroutine is as follows. - Do a loop on the elements. - For each element, localize D(K),XE(I,K). - Call SHAPE to calculate the global. derivatives of the shape functions. - Compute stresses following eq. (11.3). In order to organize the reactions, changes are made in some positions of the array NTBC(I,K). A negative value is assigned at those positions where NTBC(I,K) is equal to zero (corresponding to a specified zero displacement) starting from -NSD-1. Thus, for each negative NTBC(I,K), a reaction is computed. This change in NTBC(I,K) is actually made in PFORM, and it permits organization of the reactions in the last part of the array F(I). 41 Each reaction is obtained by calculating and adding elemental nodal forces at each location in which NTBC(I,K) is negative. This procedure requires knowledge of the stiffness matrices of those elements adjacent to the reactions. These matrices can be saved in disk storage during the execution. of 'PFORM; or can ‘be calculated again, as is done in the program. Although this procedure requires a larger number of operations, considerable storage is saved in central memory since only a reduced part of the global stiffness matrix needs to be saved. This is an important feature when treating large problems. 11.3.7 EXAMPLE OF INPUT DATA For the example of Figure 11.14, the input data is as follows, (TITLE, 1-80) EXAMPLE OF INPUT DATA, FEM (NMAT,NEL,NNO,NDOF,NLO,Free Format) 1 4 9 2 l (B(l),PR(1),E20.12,F20.12) 2.5E00 0.25 (NMAT(1),NEN(I),(NN(1,K),K81,NEN(I)),Free Format) HHD—‘H b-L‘bb Ult-‘Nr—o mum» «3me O‘UILAN (XE(I,1),NTBC(I.1),XE(I,2),NTBC(1,2),F15.5,15,F15.S,15) 0.0 -1 0.0 -1 0.0 -1 1.0 0.0 -1 2.0 2.0 0.0 2.0 1.0 2.0 2.0 x2 42 ll 3 9 .J ‘ (2) (a) 2 a :2. 1- ‘” (» AW 7 ,1. :X‘ 2. 2 . Figure 11.14. Example problem to demonstrate input data for FEM 43 40 40 4 COO NHO COO (LL, Free Format) 3 (N,I,F(2*(N-1)+I, Free Format) 7 1 1.0 8 1 2.0 9 1 1.0 0 0 0.0 (R, Free Format) 0.5 (NE1,NE2, Free Format) 1 3 0 204‘.— ote that stresses are to be calculated in elements 1,3 and 4 at four points defined by the local coordinate R-0.5. For the element 1 those points are located in (0.5,0.25), (1.5, 0.25). (1.5, 0.75) and (0.5, 0.75). CHAPTER III BOUNDARY ELEMENT MODELLING 111.1 DESCRIPTION OF THE BOUNDARY ELEMENT METHOD For the plane boundarydvalue problem of linear elasticity illustrated in Figure 111.1, the displacement at a point x on B is related to the displacements and tractions at all other points on B by Somigliana's identity, i.e. (x)u (x)+fB(UC) (x,§)uj(§)d§=fB(UR)i.j(x,§)tj(§)d§ (111.1) C111 3 i-j where the integral on the left hand side is interpreted in the Cauchy principal-value sense. The function (UC)i j(xy;) is the displacement, ui(x), due to a unit displacement discontinuity, C (i), applied in the j infinite elastic plane and (UR) (x,§) is the displacement, ui(x) due to 1.3 a unit force, Rj(§), applied in the infinite elastic plane. The coefficients a1 (x) are defined as follows for plane stress: 3 w/2§+(1+u)(sin2(a+w)-sin2a)/8w,i=j aij(x)= (111.2) (1+U)(5102(fl+w)-sin2a)/4n ,ifj 44 45 X' (X,.X;) _Xn Figure 111.1. Plane boundary value problem Figure 111.2. Angles a and u>of eq. (111.2) 46 where czandcu are shown in Figure 111.2. For a smooth boundary at x, (1)317 and aij-Li G‘ij. At a point x in R, the displacements and stresses can be calculated from the equations ui(X)-fB(UR) (x.§)cj(§)d§-IB(UC) i.j i.j(x,x)uj(x)ds . oik(x)-IB(UR)ik.j(x,§)tj(§)d§-IB(CC)ik.j(x.§)uj(§)d§ (111.3) where the influence functions (0R) (x,§) and (dC)1k (x,§) give the ik.j .j stress component Oik(X) due to a unit force, RJ(§), and a unit displacement discontinuity, Cj(§), respectively, applied in the infinite plane. At each point x on B and in each direction, either uj(x) or t.(x) is known. Therefore, eq. (111.1) can be used to solve for the unknown values of uj (x) and t1 (x), thus giving complete boundary information. The displacements and stresses at any internal point can then be determined by integration using eq. (111.3). It can be shown that, for plane stress the influence funtions of eq. (111.1) and (111.3) are given by (UR)1.k'[-(3-U)5iklog;3+ (1+u)qiqk]/(8NG) (UC)1.1‘[2(1+U)(51q13-52q23)+(1‘U)filq1+(3+U)52q2]/(4fip) (UC)1.2=[2(1+U)(-52q13-51q23)+(1+3u)52q1+(3+v)51q21/(4no) 47 (UC)2.1'[2(1+U)(~52q13—Elq23)+(3+u)32q1+(1+3u)31q2]/(4np) (UC)2.2'[2(1+U){-alq13*fizqz3)+(l-U)52q2+(3+u)51q11/(4Wo) (0R)11.1-[-2(1+U)q13-(1-U)q1]/(4flo) (OR)12.1'[2(1+U)q23-(3+u)q21/(4no) (0R)22.1=[2(1+U)q13-(1+3U)q1]/(4wo) 03R)11.2'[2(1+U)q23-(1+3u)q2]/(4np) awn)12.2-12<1+u>q13-<3+u>q11/<4no) 63R)22.2'[-2(l+u)q23-(1-u)q21/(4wo) are) -c(1«1)1<1+4 2-3 “)5 +2 <1-4 2>E 1/(2n 2) 11.1 q1 q1 1 q1‘12 q1 2 o «rc1 =c<1+ >1<1-8 2 2); +2 (1-4 2); 1/(2 2) 12.1 U q1 q2 2 q1‘12 q1 1 “o 2 - 2 - 2 (0C)22.1=G(1+U)[(1-8q1 q22)n1+2q1q2(1-4q2 )nzl/(Zwo ) (0°)11.2'(°C)12.1 (Ge)12 =(0C) .2 22.1 (0C) .G(1+u)[(1+4q22-8q24)32+2q1q2(1-4q22)fil]/(2n02) (111.4) 22.2 % where o=[(xl-§1)2+(x2_;2)2] q1=(x1-§1)/o qz-(xl-Ez)/o (111.5) and 31,52 are the components of the normal unit vector at a point on the boundary. 48 111.2 PIECEWISE LINEAR ELEMENTS Eq. (111.1) can be solved numerically if the boundary B is approximated by N straight segments, as shown in Figure 111.3. For this model, eq. (111.1) can be written as N (x(“))+ 2 fm(UC) m-l (x(“),§)u (;)d§ (n) (x )u 1.j j “11 1 N = Z [m(UR) m-l (x(n) ,1): (§)d§ (111.6) 1-3 J The displacements and tractions in each segment m can be approximated using shape functions so that (m-l) (m) N1(€)+uj N2(€) - . (m) (111.7) tj(x) tj where N1<€>=a<1-a) . N2(a)=35(1+::) (m-l) (m) i=N1(5)x +N2(5)x ds=%(sm-sm_1)dg=%gsmdg (111.8) Figure 111.3. Boundary discretization x1!!!) a; a b b2 "_"' 1 X (n) x‘m-” Figure 111.4. Definition of a,, bi’ §(m) 50 where 5 is a local coordinate for the segment m with value -1 at node m-l, value 0 at the center of the segment, and value 1 at node m. Note that the order of differenciation of t (x) in the interval is 3 less that that of u (x). This model allows discontinuities of t (x) on J i the boundary and it is consistent with the Finite Element triangular element, which has linear displacements and constant tractions on its sides. If eqs. (111.7) and (111.8) are substituted into eq. (111.6) the following is obtained N 2a1j(x(n))uj(n)+ Z Asm[fm(UC)i.j(x(n) 9€)N1(E)dguj (“1-1) m-l N +1me)i j(x‘n’,a)N2(e>dguj(m’1= 2 rm(un)i.j(x‘n),s>dgasmtj‘m) m=1 or N N (n) (n) (m.n) (m-l) (m.n) (m) (m.n)‘ (m) 2aij uj +mE1[Ai'j uj +Bi.j uj ]41/G;:1Ci j Fj n=1, ...... ,N (111.9) where aij(n)=aij(x(n)) u (n)=u.(x(n)) J (m.n)_ (n) , . E ; iii -A3 fm(UC)i.j(X 9g)1\1(g)da 51 (m.n) (n) r 1' 3L1 aosm fm(UC)1 j(x ,§)N2(>)ds (m.n). (n) , Cia G fm(UR)i.j(x ,9)dE ‘ (m). (m) Fj Asmtj (111.10) Note that the functions (UC)i j' (UR)1 j are singular when the point of coordinate x approaches the node n. This is the case when the element m ends or begins with the node n, i.e. when m-n, m=n+1 (n-m,m-1). The resulting singular integrals for these positions must be calculated analytically. A summary of those calculations is shown in Table 111.1, where b is defined in Figure (111.4). Note that the terms A1 2(n+1.n) i 1 2(n.n), which contribute to the coefficient of “2(n)’ produce a net result and B (n+l.n)+B (n.n)=(1-u)log(Asn/Asn+1)/(2n) A1.2 1.2 The same applies for A2 1(n+l.n) + B2 1(n.n). The nonsingular integrals are evaluated by numerical integration using the Gauss quadrature formula. If it is defined that §(m)=%[x(m)+x(m-1)] . g. m-.. m . =X.-;.gm> i i i i 2 R=aiai-2aibi;+bibi; (111.11) 52 Table 111.1 Singular Integrals n:node number m:segment number Condition Coefficients Value (m.n) n-m-1,m Ci.j [(3-u)(1-1og(2b))&i +(1mbibj/b21/(zm) j 3—. A (m.n):B (m.n) “Hm 1.1 1.1 o (m.n) (m.n) A2.2 ”32.2 (m.n) (m.n) (1- )/(2") “'m A1.2 "A2.1 U (m.n) (m.n) -(1-u)/(2n) n'm'l B1.2 “Bz.1 (m°n)=-A (m.n) (l-u)[2+21im[log(e/Asm)]]/(4n) n=m-l A 1.2 2.1 8+0 n-m B1 2(m.n)=_Bq 1(m.n) (l-U)[-2-211m[log(€/Asm)]]/(Aw) . ‘- 0 €+O Net contri- A1 2(n+1-n)+31 2(n.n) buti°n ' ' (l-U)log(Asn/Asn+1)/(2n) (n+1.n) (n.n) °f 5 A2.1 +32.1 and 6 53 as shown in Figure 111.4, the nonsingular integrals can be calculated for m#n,n+1 W - - 1.1(m'n) '[2b2(1+u)L‘(1+5)[(al-blg)3/R21dg+2b1(1+u)£:(1+€)[(32_b25)3/R2]dg 1.1(m.n:‘+b2(l-U)£:(1:5)[(al-b1€)/R]dS-b1(3+u)I:(1:5)[(a2-b25)/R]dg]/(4w) A A1.2(m'“7‘-[2b1(1+u>£1(115)[(al-blg)3/R21dg-2b2(1+u)£}(11g)[(a2-b25)3/Rz]dg ? (m.n) 1.2 -b1(1+3U)£:(1:5)[(al-b1€)/R]dg+b2(3+u)£=(1:5)[(az-b2€)/R]dg]/(4fl) L A2;1(m°ni)'[2b1(1+u)£f(l:€)[(al-bla)3/R21dS-2b2(1+u)£:(1:5)[(az-b25)3/R2]d€ > 2.1(m'ni) ‘b1<3+v>£1(135)I(al-bla)/R1da+b2(1+3u>1:(1:5)[(aZ—bzg)/R1dg]/(4,) \ _ _ (m-n) .[-2b2(1+U)L}(1+5)[(al-b1€)3/R2]dg-2bl(1+u)£:(1+€)[(82-b2€)3/R2]d€ 1 82.20”“: +b2(3+u)f,} (1:5,)[(al-b1€)/R]d€-b1(1-u)f.}(1:5)[(aZ-bZSVRldEJ/(M) 2.2 1.1(m'n)‘[‘5(3‘”>411°8d€+<1+u)fl [(31-b15)2/R]d€]/(8n) 1.2(mon)=C2. 1(m.n)=[ (141))1')‘ [(31-1315) (32-1325.) /R]dc:] /(8T") c2 2(m°“)=1-s(3-u)1}103(R)da+(1+u)1:[(az-b,:>2/R]da]/(81) (111.12) result groble nodal 54 Eq. (111.9) can be also written in matrix form as A [UC]{Gu}=[UR]{F} (111.13) This is a system of equations relating nodal displacements to resultant segment forces. In order to solve a well-posed elasticity problem, it is necessary to re-pose this system of equations in terms of nodal forces. Therefore, a transformation {F}=[r]{F+ (111.14) relating the vector of nodal forces {F} to the vector of segment forces A {F}, is introduced into the system (111.13). The simplest physical interpretation of the transformation is to replace the segment forces by nodal forces equal to the average of the segment forces adjacent to each node, or (n) " (n) " (n+1) Fi =8lFi +11 1 (111.15) The form of [F] for this transformation is '1 1 o o..o‘ o 1 1 o..o [1]=la o 0 1 1..o (111.16) 0 o o 1..0 _1 o o o..1d where I is a 2x2 identity matrix. 55 For an odd number of nodes, the inverse of [F] is - 1 1 -1 1 -1...1 _1 1 1 -1 1..-1 [r] . -1 1 1 -I...1 I -1 1 1..-1 b-I 1 -1 1 1 and eq. (111.13) becomes [UCJIGu}-[UR][1]71{F}. (111.17) (111.18) It should be noted that, for an even number of nodes, [F] has no inverse. A small perturbation is introduced into the system of equations (111.18) in order to enforce the equilibrium conditions N N N . z F1(i)=0 z F2(i)=0 2 (x1(1)1:.2(i)_x2(i)5.1(1)):0 i=1 i=1 181 or in matrix form [Q]{F}=0 where _ 1 0 1 0 1 0 ' [01- o 1 o 1 .... o 1 .‘x2(1) x1(1) _x (2) x1(2) _X2(N> x1(N)‘ (111.18) and(111.19): UC UR* QT 11} {Gu}= 0 Q o {A} 1UC11Gu1=IUR*1 {F1+101T{1} o=101(r}. Coupling eqs. where [UR*] = [UR][F]-1. In partioned form (111.19) (111.20) 56 The vector {)}(3 )1 1) contains three operators which introduce a perturbation into the system. Note that the higher the number of nodes in the model, the smaller the components of{x} . After solving eq. (111.20),. stresses and displacements can be calculated anywhere in the body using a discretized form of eq. (111.3) which can be written as ul(x> ( > ( > u * u ( ) . [UR 1(2,2n){F}-[UC 1(2,2n){u} 112 X 011(X) ( > * ( ) s , s “12(X) where [UR(“)]*=IUR(“)1111‘1.1UR(S)1*=IUR(S)1111‘1. All the entries of the matrices as [UR(u)], [UC(U)], [UR(S)], [UC(8)] are calculated by numerical integration. Eq. (111.12) are used (11) to compute [UR(U)] and [UC ]. For the entries of [UR(S)] and [UC(S)] the following equations apply (UR(S’) 1)-1-2(1+u>11 ((al-bla)3/121dg-(1-6>13 [(al-blg)/R]dg]/(8w) 13(2n- (UR(9’>1,(2n)=12(1+u)13 1(a2-62513/R21dg-(1+36)13 [(az-b25)/R]dg]/(81) (UR(S)) (UR(8)) (UR(5)) (UR(3)) (UC(8)) (UC(8)) (UC(8)) (UC(5)) (UC(S)) (UC(S)) (UC(S)) (UC‘S)) 2,(2n-1) 2.(2n)’[ 3,(2n-1) 3.(2n>'[ 1,(2n-3)) 1,(2n-1)/ \ 1,(2n-2) 57 =I+2(1+u)f3 1(al-bla)3/R21ds-(1+3u>fj [(al-b1€)/R]d€]/(81) -2(1+6)x; 1(a2-b2g13/R21dg-(1-6)13 [(aZ-b25)/R]dg]/(8fi) -12(1+u>11 [(az-bza)3/R2]da-(3+u)fj [(az-b2€)/R]d€]/(81) 2(1+u)1: 1(a1-bla)3/Rzlds-(3+u)rl [Cal-bla)/R]d£]/(81) -c(1+u>12b2(1:((1151/R1da+411(111)[(a1 -bls>2/R21da > -8£§<1¥a>1(a1 -bla>4/R3Idal -4bl[rj(1ia>[(a1-bla)(az—bza)/R2][1-a(a1-bla)2/R 1d611/(an) -c(1+u>1-2blit((1la>/R111-8(a1-bla>2(a2-bza>2/R21da > 1.(2n) J 2,(2n-3; +4b2 {Kl:i)[(a1-bla)(az-b2€)/R2][1-4(a1-b1£)2/R]d€1/(8n) ---c(1+u>1+2b2 1:1(IIa>/R111-8(a1-blg>2(a2-bzg)Z/Rzlda ? 1 -_ P _ P 2 _ - P 2 2,(2n-1)j "b1£u(1*9)[(al-blg>1-2b111:1(1$a>/11da+4 11(Ila)1(a2-b25>2/R21ds -81:(115)1(a2-b2;>4/R31d11 +4b.1 1:(115)1(a1-bla)(a.-b21)/R21[1-4(a,-b2:>2/R1d:31/(8.> 58 (3)) (3)) 3,(2n-3)"(UC 1,(2n-2) 3,(2n-1) 1,(2n) (S) (3)) (111.22) (3)) (uc(5)) =(Uc (UC -(Uc 3,(2n-2) )2,(2n-3) 3,(2n) 2,(2n-1) where ai’bi and R are defined by eq. (111.11). 111.3 PLANE STRESS BOUNDARY ELEMENT PROGRAM The boundary element program described in this section can solve plane stress elastostatics problems using linear interpolation for the displacements and constant segment tractions. The flow chart for the program is shown in Figure (111.5). In the first part, the data of the problem is read. Then the weights and points for the numerical integration of the non-singular integrals are assigned to the arrays W(L) and R(L). Next, the augmented matrices [UR] and [UC] are initialized. In part 4 of the program, the entries of the matrices [UC] and [UR] are calculated. In part 5, the operation [UR] [F I"1 is performed, and in part 6, the system of equations is rearranged in order to have all the unknowns on the same side of the system of equations. Subroutine GAUSS solves the equations, and the unknown are printed in part 8 of the program. Finally, stresses and displacements are calculated for the field points. 59 C START ) / I Input Data] Assign weights (and points for numerical integration I Initialize matrices UR, UC I Calculate entries of UR and UC 1 $0 Modify UR l Rearrange equations 1 Solve equations. CALL GAUSS I Output of displacements andfgrces 1111 Calculations for field points Figure 111.5. Flowchart of program BEM 111 60 111.3.1 INPUTTING THE PROBLEM DATA. The order and format in which the variables and arrays must be submitted to the program are as follows: TITLE X(I). Y(I) J,K,F(2*J+K—2) -The title of the problem. This must be written on one line in columns 1 through 80. -The number of nodes on the boundary. Must be odd number. Enter in free format. -Material properties where E is Young's modulus and PR is Poisson's ratio. Enter in free format. -Coordinates of nodes 1 through N, entered lcounter-clockwise in free format. -Nodes, J, and corresponding directions, K, at which displacements are specified (zero and non-zero). Enter in free format and end by inputting 0,0. -Nodes, J, and corresponding directions, K, at which non-zero forces are specified, and specified values F. Enter in free format and end by inputting 0,0,0.. 61 LL -No. of points used for numerical integration. LL maybe 3,4,5 or 6. I,XF,YF -Coordinates of each internal point at which displacements and stresses are to be computed. Enter in free format, one point per line. Always begin line with 181. Enter Isinteger other than one to stop. 111.3.2 CALCULATION OF COEFFICIENT MATRICES. A more detailed flow chart for part 4 of the program is shown in Figure 111.6. The order of the loops was chosen in order to avoid repeating operations. For each value of J, associated with the element J, parts of the (2x2) singular submatrices, shown in Figure 111.7, are calculated. Submatrices 1 and 3 correspond to the cases 3 and 4 respectively of Table 111.1 Submatrix 2 contains contributions from c: of eq. (111.2) and 13 from the net contributions of cases 5 and 6 of Tables 111.1, i.e. case 7 of that table. Submatrices 4 and 5 correspond to case 1 of Table 111.1. In the other loop, I, all of the nonsingular integrals are calculated by numerical integration. The positions in the matrices are illustrated in Figure 111.8. The functions for the numerical integration are calculated using eq. (111.12). 62 Loop on Jcolumns Calculate I \:/ singular entries Loop on rows\ 1 = 1, NM j/l Loop on points of Integration L = 1, LL 1 Calculate nonsingular rals Figure 111.6. Flowchart of part 4 of the program 63 J \ \ \ J \ \ . \\ [UC] \ \ ' 1 2 3 \ \ L . a J —-TR '1 \ \ \ , \ \ [on] \ \\ T]? "T \ \ L E Figure 111.7. Singular contributions in [UC] and [UR] 64 -r_J< \ ZS \ 2S J \ a V5 1 §§\\ .. 2.: \s \\ § [UR] \\ \ H J I \ s \ A \__.. \ _ S \_. I o . ‘0 - R W Figure 111 8. Non-singular contributions in [UC] and [U J 65 III 3.3 TRANSFORMATION OF COEFFICIENT MATRIX UR This part of the program performs the operation * -1 [UR] =[URJIF] in which [1"]-1 is defined in eq.(III.17) The most efficient way to obtain the matrix [UR]* is not to create [1']-1 and postmultiply by [UR], but to modify [UR] itself. The algorithm for this modification is better understood by examining the example of a 6 x 6 matrix: * * * [URll] [UR12] [UR13] UR11 UR12 UR13 1 -I 1 UR * UR * UR * UR UR UR 1 1 [ 21] [ 22] [ 23] ' 21 22 23 '1 UR * UR * UR * UR UR UR 1 I 1 .[ 31] [ 32] [' 3311 L. 31 32 33- L’ . ’ in which the entries [URZi] Care 2‘: * [UR21] ‘ [URzlI +'[UR22] ' [UR23] a— a- — '7 [UR221 [UR21] + [URZZI + [UR23] [[URZIl + [UR22] [UR23]]+~[UR221 --[UR21] +2[UR22] * it [UR23] a-[UR22] +2[UR23] In general, it can be shown that * [LRkl] =[URk1]+[LRk2]-[LRR31+[LRk4]-' . . .+[LRkn] (111.23) and [URki]*=-[URk(i_l)]*+2[URki] , i=2,...,N. (111.24) 66 The flow chart for this modification is illustrated in Figure 111.9. In the first loop, the first two columns of [UR]* are obtained using eq. (111.23). In the second loop, the other entries of [UR]* are generated using eq. (111.24). 111.3.4 CONSOLIDATION OF KNOWN AND UNKNOWNS BOUNDARY VALUES This part of the program consists of a reorganization of the system of eqs. (111.20) so that all unknowns go to the left hand side and all knowns go to the right hand side of the equations. This procedure is illustrated with the flow chart of Figure (111.10). The array NTBC (K,I) -which specifies the boundary condition at node K, directions I, is used to decide when to interchange columns, i.e. if NTBC(K,I,)=0, there is a specified force, and no interchange is necessary. The final vector of knowns is multiplied by the matrix of coefficients to produce the right-hand side of the final system, of equations. The final system of equations is then solved by the subroutine GAUSS, which uses the well known GAUSS elimination procedure. 111.3.5 CALCULATION OF STRESSES AND DISPLACEMENTS AT INTERNAL POINTS In this part of the program, eq. (111.21) are evaluated. Initially the matrices [UR(u)], [UC(u)], [UR(S)], [UC(S)], defined by eqs. (111.12) and (111.22), are calculated by numerical integration using the Gauss quadrature formula. These matrices are saved in the first five rows of the arrays [UC] and [UR]. Next [UR(U)] and [UR(S)] are modified to (u) produce [UR ]* and [UR(S)]*, respectively. Finally, the matrices are 67 3’ Loop on 1\\\ I = l, N 1 Calculate [UR]: eq. (111.22) Loop on K‘\\ 2f Calculate [UR]: eq. (111.23) "——“‘ 1 CONTINUE 7\\— C J J, Figure 111.9. Flowchart for the modification of [UR] 68 Loop on degrees of freedom I = l, 2 , l [KK=NTBC(K,I) I II C) KK KK KK=l Interchange columns , [2*(K-l)+1] of [UC] and [UR] with signs changed I @765} Multiply matrix at RHS by.vector of nowns, save in SAV(1) ,L Figure 111.10. Flowchart for rearranging equations into final form 69 multiplied by the vectors {F} and {u} to produce displacements and stresses at each chosen field point. 111.3.6 EXAMPLE OF INPUT DATA For the example of Figure 111.11, the input data is as follows (TITLE, 1-80) EXAMPLE OF INPUT DATA - BEM (N, Free format) 7 (E, PR, Free format) 2.5E00 .25 (X(1), Y(1), Free format) 0.0 ' 0.0 2.0 0.0 4.0 0.0 4.0 2.0 2.0 2.0 0.0 2.0 0.0 1.0 (J,K, Free format) 1 OO‘VH (3 id 1— to rd (J,K,F(2*J+K-2) Free format) 3 1 1.0 4 1 1.0 0 0 0.0 70 3‘2 c :1 6 5 4 1. 7 Co B. 1 A. I 2 3 I. x; J, :—-—-—p 2 2- Figure 111.11. Example problem to demonstrate input data for BEM 71 (LL, Free format) 5 (I,XF,YF, Free format) 1 2.0 0.5 l 2.0 0.7 1 1.0 1.0 O 0.0 0.0 Displacements and stresses are calculated at points A,B and C of Figure 111.11. CHAPTER 1V COUPLING THE FINITE ELEMENT AND BOUNDARY ELEMENT METHODS 1V.1 PROCEDURE The procedure outlined here for combining the two methods has been suggested. by several authors, .[3,4,6,7,8,9]. The method consists of dividing the body as shown in Figure IV.1. Elements defined by more than four nodes, such as the shaded element, are treated as "superelements", whose symmetric stiffness matrices are obtained by using the boundary integral equations. These matrices are incorporated into the global stiffness matrix for the entire body and the equations are solved as was done in the program FEM described in CHAPTER II. The field equations are satisfied exactLy in each superelement. Compatibility with the neighboring finite elements is guaranteed since the shape functions for the superelement boundaries have the same form as those for the isoparametric quadrilateral and triangular elements. Equilibrium is satisfied in an average sense. The procedure used to obtain a symmetric stiffness matrix for a superelement begins with solution of the system (111.20) for the set of displacement vectors in which only one of the components is l and all others are 0. The solution vectors can be arranged in a matrix [E] so that 72 73 X: X. Figure 1V. 1. Combined discretization 74 1E11u1(38)=111(83) . (1v.1) (39) and {F}(Se) are the vectors of nodal displacements and where {u} nodal forces, respectively, for the superelement. The same result may be reached if the inverse of the augmented matrix [UR]* is multiplied by [UC], but this procedure was shown to be less efficient for the problems examined. but this condition Note that the entry E should be equal to E 13 31' is, in general, only attained exactly as the number of nodes goes to infinity. However, since the boundary element model has exactly the same interpolating functions on the elements as the finite element triangular element, the matrix [E] of eq. (1V.1) is exactly equal to the stiffness matrix for that triangular element. It is important to note that the rows and columns of [E] each sum to zero, since equilibrium equations were enforced in the superelement. The necessary symmetric matrix of the superelement is obtained by minimizing the functional n(se)=s fV{€}T{o}dv-{u}(se)T{F}(Se) (1V.2) which can be written for the discretized model as H 75 where [K(Se)]=k{[E]+[E]T}. (1v.5) As the number of nodes in the superelement increases, [E] approaches [K]. 1V.2 PLANE STRESS HYBRID PROGRAM 1V.2.l GENERAL DESCRIPTION This combined program is the same finite element program FEM described in CHAPTER II, with the addition of the subroutine BOUN, which creates the stiffness matrix of the superelement and calculates stresses and displacements at the field points within the superelement. A few variations are also introduced into the subroutines PFORM and STRESS, allowing the program to call BOUN when required, see Figure 1V.2. The format and order for the input data is exactly the same as was described in CHAPTER IV, with the addition of the coordinates of the field points in the superelement as described in section 111.3.1. An example of the input data for a problem is given in section 1V.2.4. 1V.2.2 SUBROUTINE BOUN The flow chart for BOUN is shown in Figure IV.3. The first part of the subroutine, which works when the argument 1K is l, calculates the stiffness matrix of the superelement. The second part of BOUN computes the stresses and displacements at the field points when the argument IR is 2. Steps 1,2 and 3 of BOUN are the same as those already explained for the program BEM in CHAPTER III. 76 fNIN>4 '[ N1N£4 CALL SUBROUTINE CALL ELEMENT BOUN <.....1) - mm. - l fNIN>4 ’TJ CALL SUBROUTINE f BOUN (....,2)> , STRESS ,L NIN: No of element nodes Figure 1V.2. Interaction between BOUN and FEMBOU Calculate entries of UR and UC { Modify UR 1 9 9°- Solve system of equations for unit displacement vectors. A§ymmetrize 77 L _____w {I Calculate stresses and ‘displacements at field points (:RETURN > Figure IV.3. Flowchart for BOUN 78 In part 4 of the subroutine, the system of equations * 1 uc UR Q {F} 4' (1v.6 ) 0 Q o A} (2N+3,2N) (2N13’§N*3) is solved for the set of unit prescribed displacements of the form (0 m -< 0H0- (0 This is the same as solving the system UR* T Q {x} = {b} Q 0 for vectors {b} equal to each of the columsof[DC] . The method for 0 doing this is to perform the factorization * UR QT = [L][U] . Q 0 where [L] and [U] are lower and upper triangular matrices respectively, and then to solve by back and forward substitution for each column of [U:]. An additional array IPIV(I) contains row permutations if any are rgquired. Each solution vector is saved in the same matrix [UC] so that, when the procedure is completed, the first 2N rows of this matrix are equal to the matrix [E] of eq. (IV.1). A more detailed explanation of this method can be found in [12]. 79 1V.2.3 INPUT PROCEDURE The. input data procedure is almost the same as that for the FEM program, although few variations are introduced as follows. The number of points of integration is fixed into the program as 3 for the finite elements and 6 for the boundary elements, therefore these numbers do not need to be given to the program. The points in which stresses are to be calculated within the finite elements are defined by the local coordinate R as in the FEM program. The intervals of finite elements in which the stresses are to be calculated are given in ascending order, ending with zeros. If no stress calculation is required into the finite elements, two zeros must be input. The field points within each superelement must be given after the interval of elements corresponding to that particular superelement, ending with zeros as in the BEM program. If no stress calculation is to be made into the finite elements, the field points within superelements must be given in order from the first to the last superelement in the problem, each list ending with zeros as in the BEM program. An aclaratory example is presented in the next section. 1V.2.4 EXAMPLE OF INPUT DATA The input data for the problem of Figure 1V.4 is as follows. (TITLE, 1-80) EXAMPLE OF INPUT DATA, FEMBOU (NMAT,NEL,NNO,NDOF,NLO, Free format) 1 5 12 2 1 80 Figure 1v.4. input data for FEMBOU X2 1] “I. J 2 s a To 12 ' (4) (3) A 8 ° 4 7 - (5) (2) (I) 1 3 6 9 11 hx) I / 1. 1 1. T. 1 Example problem to demonstrate 81 (E(1),PR(1),E20.12,F20.12) 2.5E00 0.25 (MP(I),NEN(I).(NN(1,J),J=1,NEN(I)), Free format) 1 4 3 6 7 4 1 5 1 3 4 5 2 1 4 4 7 8 5 l 5 6 9 10 8 7 l 4 9 11 12 10 XE(I,1),NTBC(I,1),XE(I,2),NTBC(I.2).F15.5,IS.F15.5,15) -1 0.0 -1 NONONHONHON 00000000000 ( 0 0 l 1 1 2 2 2 3 3 4 4 (N,1,F(2*(N-1)+1), Free format) 12 1 1.0 8 2 1.0 0 0 0.0 Three different alternatives for calculation of stresses are as follows. (R, Free format) 0.5 (NEl, NE2, Free format) 1 4 0 0 (I, XF, YF, Free format) 1 0.5 1.0 0 0.0 0.0 1 2.5 1.0 0 0.0 0.0 Stresses are to be calculated into the finite elements 1 and 3, at four points defined by the local coordinate R=O.5; and at points A and B of superelements 2 and 4 respectively. (R, Free format) 0.5 82 (NEl, NE2, Free format) 0 0 (I, XF, YF, Free format) 0 0.0 0.0 l 2.5 1.0 0 0.0 0.0 Stresses are to be calculated at point B of superelement 4. (R, Free format) 0.5 (NEl, NE2, Free format) (I, XF, YF, Free format) 1 2 1 0.5 1.0 0 0.0 0.0 4 5 1 2.5 1.0 0 0.0 0.0 0 0 Stresses are to be calculated into the finite elements 1 and 5, at four points defined by the local coordinate R=0.5; and at points A and B of superelements 2 and 5 respectively. CHAPTER V EXAMPLES AND DISCUSSION V.1 EXAMPLES In first example, the cantilever beam shown in Figure V.1 was modelled using 13 and 19 nodes on the boundary, as illustrated in Figures v.2 and v.3. The CPU time on a PRIME 750 minicomputer and the computed vertical deflection of the point A are shown in Table V.1. The vertical deflection is compared to the exact solution for the same problem with the true parabolic load distribution on its end, i.e. 14.8 taken from [13]. Normal stresses, 011, on the section a-a are shown in Figure v.3 for the 19-node discretization. Note that the results of <31 given by BEM and FEMBOU are not reliable near the boundary. In the second example, half of the plate of Figure v.5 was modelled using 49 and 55 nodes on the boundary for the ratios R/w=0.5 and R/w-O.2 respectively, see Figures v.6, v.7, V.8, v.9, V.10 and v.11. Results for the normal stress, 022, on section a-a are plotted in Figures V.12 and V.13, as well as the exact solution for the problem in which the notch is a perfect semicircle, as taken from [14,15]. Note that the ratio of perimeter to area is 1.46 1/UL for the plate with R/w=0.5, and 0.53 1/UL for the plate with R/w=0.2, where UL is a length unit. 83 84 X2 1 A 1J67 1 j I / 2 l ” - - - l]. A ° ] / l 2 1 / l C 1 Figure v.1. Cantilever beam X1 85 if ._ A 013625 1 1 0.6875 ‘ 1015625 (b) ‘1? {If Figure v.2. 13 of (c) node boundary discretization cantilever beam 86 [0.074074 ’ 10.425926 ‘x’ftffi ) ) 10.62 592‘ ; e A e : 10 .074074 (a) (b) (c) Figure v.3. 1? node boundary discretization of cantilever beam I)... 1.61 1.4). u. 19 Figure V.4. Stress d'S 87 DEMO 'EMD 'EMIOU A J bution on section a-a, from free and k)! (Dr-I” located 88 Table V.1. Vertical deflection of A and CPU time for cantilever beam problem. Model Program Deflection of A Zdifference(l4.8) CPU Time(sec) 13 BEM -13.829 6.57 2.57 FEM -13.358 9.74 3.52 FEMBOU -14.194 4.09 3.41 19 BEM -14.230 3.85 4.52 FEM -l4.149 4.39 6.13 FEMBOU -14.292 3.43 5.81 89 °TW Figure v.5. Example No. 2 90 Figure v.6. BEM model for R/w = O. 91 Figure v.7. FEM model for R/w = 0.5 92 Figure V.8. FEMBOU model for R/w = 0.5 Figure V.9. 93 BEM model for R/w O .2 94 Figure V.10. FEM model for R/w = O 2 95 H]! “H Figure v.11. FEMBOU model for R/w = 0.2 BEN O FEM O FEMBOU A i.” ) Figure v.12. 96 U).- i n a-a for C) on secti 0.5 Stresses On example 2, R/W = 97 3.0 4 DEMO FEMO FEMBOU A 2.04) Suponlomom lFinin EloMoMt A N Figure v.13. Stresses oh on sec example 2, R/w = O 98 CPU time in seconds, on a PRIME 750 minicomputer is listed in Table v.2. The number of equations solved in each procedure and the number of components of the coefficient matrices are in Table V.3. 99 Table V.2 CPU time for problems of Figure v.5 CPU CPU CPU TOTAL Relation R/W Program READ FORM SOLVE,STRESS CPU 0.5 BEM 1.03 8.17 17.13 26.33 (49 nodes FEM 3.01 8.46 4.71 16.18 boundary) FEMBOU 2.60 8.25 6.00 16.85 0.2 BEM 1.27 12.19 27.55* 41.01 (55 nodes FEM 4.53 14.74 9.00 28.27 boundary) FEMBOU 3.19 17.23 8.87 29.29 * 27.55-18.52(solving)+0.62 (writting)+8.41 (stresses) Table V.3. Number of equations and components of coefficient matrices for problems of Figure V.5. Model Program No. of equations No. of comp. in coefficient matrices 49 nodes BEM 101 10201 FEM 225 3127 FEMBOU 197 3147 55 nodes BEM 113 12769 FEM 317 5808 FEMBOU 241 5681 100 V.2 DISCUSSION OF RESULTS Some conclusions can be drawn from the results of the two examples. As expected, it is much easier to prepare and input the data for BEM than for FEM. This factor can be quantified by comparing the time that each program spends reading the data. This advantage of. BEM ‘will undoubtedly be even higher when solving three dimensional problems. The time spent calculating the coefficient matrices [UC] and [UR] by BEM, is a little bit lower than the time spent creating the global stiffness matrices by FEM and FEMBOU. The comparison is made using 6 points of integration on each boundary element and 9 points of integration in each quadrilateral finite element. The difference is not substantial however. On the other hand, there appears to be a substantial difference in the time required for creation of the stiffness matrix of the superelement when the number of nodes goes beyond a certain value between 19 and 29. Recall that the computation of this stiffness matrix requires a matrix inversion process, the expense of which increases dramatically as the number of nodes increases. For the same number of boundary nodes, fewer equations are solved in the BEM program, but since the BEM matrix of coefficients is neither symmetric nor banded, a great effort is necessary for solving the equations. Conversely, the high number of equations in the FEM method are solved quite efficiently, due to the symmetric and handed form of the global stiffness matrix. This difference can be appreciated by looking at Table v.3 and at Column 5 of Table v.2. On the other hand, the time for the combined program is almost the same as for FEM. 101 Due to the boundary discretization of the boundary integral equations, field calculations are not reliable near the boundary, particularly when non-zero tractions are applied on the boundary as is the case in Figure V.13 near the line joining the finite elements to the superelement. Good accuracy is expected at 1 1/2 mesh length from the border. In the case of relatively thin elements, such as the plate for which R/w-O.5, the accuracy of the BEM program is not good compared to the accuracy of the other two programs working with analogous modells. On the other hand, the combined program produces the most accurate results at interior points far enough from the borders. This is an encouraging result in anticipation of future three dimensional programs. In general, the FEM program proved to be the most accurate and most efficient for the problems solved. However a considerable amount of time was spent preparing the data for each problem. Very good results were obtained with the FEMBOU program at a speed similar to that for FEM and with less time spent creating the input data. Some attention must be paid to improving the boundary element stress calculations near or at the borders and to improving the effectiveness of some of the algorithms in the boundary element programs, i.e. solving equations in BEM and creating the stiffness matrix of the superelement in FEMBOU. Some suggestions are offered below. 102 v.3 RECOMMENDATIONS FOR FUTURE STUDY The following suggestions or considerations are made regarding future work in this area. It appears to be possible to calculate stresses near or on the free-traction portion of the boundary in the programs BEM nad FEMBOU, this calculation would require exact evaluation of the singularities of [UC(U)] of eq. (111.21), which involves a procedure similar to that for the calculation of the singularities of [UC] as explained in Chapter 111. All the terms of [UR(S) ] of eq. (111.21) can be computed by numerical integration. Anoter technique, such as extrapolation or a Finite-Boundary eiement combination, must be devised for the stress calculation on the loaded boundaries. . The efficiency of the FEMBOU program can be improved a great deal in the cases where the superelement is not completely imbedded within the finite elements. For this case, only the "stiffness" matrix, corresponding to the nodes in the common finite-boundary border must be calculated, which implies that eqs. (111.20) must be solved only for the unit displacement vectors involving the common boundary-finite nodes. It would be wise to examine other methods for solving the system of equations in BEM. Should a good algorithm be found for this operation, the general efficiency of the program will be improved a great deal, [16]. Since the analytical calculation of singular entries is sometimes quite complicated, it might be advisable to use rigid body motion arguments to compute those entries, i.e. rigid body displacements in x 1 or x2 directions must produce a zero stress field so that all the even 103 entries of the same row of [UC] must add to zero (the same applies for the odd entries). This argument, useful here for checking intermediate results, should be used to prepare the three dimensional program. It is advisible to prepare a subroutine to generate and check graphically the subdivisions for a given Finite Element problem. This applies not only for the plane problem but particularly for the three dimensional case. Many articles about this topic can be found in the literature, [17,18,19] Few changes are necessary to extend FEM to the three dimensional case. Another part must be added to the subroutine ELEMENT for the calculation of the stiffness matrix of the the three dimensional element. Small changes are required in some dimensions and formats of the program. APPENDIX COMPUTER PROGRAM LISTS 104 00001:C (.0002: Convoifiooo'too o§9§ooooooonoovo o§§§§§§999§§0999«womboctocéootittntnnncofin 00003: 00004: PROGRAM BEM 00005:C 00006:C THIS PROGRAM EMPLOVS THE DIRECT BOUNDARY ELEMENT METHOD 00007:C (USING STRAIGHT BOUNDARY ELEMENTS CHARACTERIZED BY LINEAR 00008:C DISPLACEMENTS AND CONSTANT TRACTIONS) T0 SOLVE PLANE STRESS 00009:C PROBLEMS OF LINEAR ELASTICITY. UNIT THICKNESS IS ASSUMED. 00010:C 00011:C WE DEFINE THE FOLLOWING INPUT PARAMETERS, VARIABLES 00012:C AND ARRAYS: 00013ac 00014:C TITLE I THE TITLE OF THE PROBLEM. THIS MUST BE WRITTEN 00015:C ON ONE LINE IN COLUMNS I THROUGH 60. 00016:C N I THE NUMBER OF NODES ON THE BOUNDARY. MUST BE 00017:C ODD NUMBER. ENTER IN FREE FORMAT. 00018:C E,PR I MATERIAL PROPERTIES (E IS YOUNG'S MODULUS 00019:C AND PR IS POISSON'S RATIO). ENTER IN FREE FORMAT. 00020:C 00021:C X(I),Y(I) 00022:C I COORDINATES OF NODES 1,N ENTERED COUNTER-CLOCKWISE 000231C IN FREE FORMAT. 00024zC 00025:C J,K I NODES, J, AND CORRESPONDING DIRECTIONS,K, 00026:C AT WHICH DISPLACEMENTS ARE SPECIFIED(ZERO AND 00027:C NON-ZERO). ENTER IN FREE FORMAT AND END BY INPUTTING 00028:C 0, 0. 00029zC 00030:C J,K,F(2*J+K-2) ‘ 000313C I NODES,J, AND CORRESPONDING DIRECTIONS, K, AT WHICH 000323C NON-ZERO FORCES ARE SPECIFIED, AND SPECIFIED VALUES 000333C F. ENTER IN FREE FORMAT AND END BY INPUTTING 0, O, 0. 00034:C 000353C LL I NO. OF POINTS USED FOR NUMERICAL INTEGRATION. 000361C LL MAY BE 3,4,5, OR 6. 00037:C 00036:C I,XF,YF 000391C I COORDINATES OF EACH INTERNAL POINT AT WHICH DIS- 00040:C PLACEMENTS AND STRESSES ARE TO BE COMPUTED. 000418C ENTER IN FREE FORMAT, ONE POINT PER LINE. ALWAYS 000428C BEGIN LINE WITH III. ENTER IIINTEGER OTHER THAN ONE 000433C TO STOP. 00044:C 00045:C 00046:C§§§§§**§§*§**§§§§****§****§§§§****§§§*§§§Q§**§§*§*§*§*§§*§*§§*§§§G**§Q§ 00047:C 0004B: COMMON UR(203,200),NTBC(200.2) 00049: COMMON R(6),W(b),Wl(6),W2(6) 00050: COMMON ICGAUSS/ UC(203,203),SAV(203) 00051: DIMENSION X(100),Y(100),F(200) 00052: DIMENSION TITLE(20) 000531C 00054: CALL CLOCK(START) 00055: OPENN FORMAT(/,5X,'NUMBER OF NODES - WRITE(10,9)E,PR ',3x,15) FORMAT(/,5X,'YOUNGS MODULUS I ‘,EIZ.6.// +,5X,'POISSONS RATIO I ',F10.4) WRITE(10,IO) FORMAT(//,3X,'NODE NUMBER',6X,‘NODAL COORDINATES AND BOUNDARY COND +ITIONS') DO 11 1-1,~ WRITE(10,12)I,X(I),NTBC(I,1),Y(I),NTBC(I,2) FORMAT(/,3X,I5,5X,3(F15.5,15)) POINTS AND WEIGHTS FOR NUMERICAL INTEGRATION WRITE(10,13)LL FORMAT(//,‘POINTS OF INTEGRATION',3K,I5> CALL CLOCK THEN ‘ R(1)=-0.86113631159405 R(2)I-0.33998104358485 1,101 3:}: 00131: 00132: 00133: 00134: 00135: 00136:C 00137: 00138: 00139: 00140: 00141: 00142: 00143: 00144: 00145: 00146: 00147: 00148:C 00149: 00150: 00151: 00152: 00153: 00154: 00155: 00156: 00157: 00158: 00159: 00160: '00161: 00162: 00163: 00164: 00165:14 00166:C 00167:C 00168:C 00169: 00170: 00171: 00172: 00173: 00174:15 00175:C 00176: 00177: 00178: 00179: 00180:16 00181:C 00182: 00183: 00184: 00185: 00186:17 00187:C 00188:C 00189:C 00190: 00191: 00192: 00193: 00194:C 00195:C 106 F\(:-)'-'I\(;') 5(4)=’F(1) k(1)-0.3478548451?745 W(2)IQ.652I4515485254 ”(3)3W(2) “(4)3WII) ELSE 1F(LL.EQ.5) THEN R(1)--0.90617984593866 R(2)--0.53846931010568 R(3)-0.0 Rt4>c-R(2) W(1)-0.23692688505618 W(2)IO.47862867049936 W(3)-0.56888888888888 W(4)IW(2) W(5)IW(1) R(5)s-R(1) ELSE IF(LL.EO.6) THEN R(1)I-0.93246951420315 R(2)I-0.66120938646626 R(3)I-0.23861918608319 R(4)I-R(3) R(5)I-R(2) R(6)I-R(1) W(1)=0.17132449237917 W(2)I0.36076157304813 W‘3)I0.46791393457269 w<4)-w(3> W(5)IW(2) N(6)IN(I) END IF DO 14 LI1,LL H1(L)II.-R(L) W2(L)II.*R(L) INITIALIZE MATRICES DO 15 II1,NNP3 DO :5 JI1,NNP3 UC(J,I)I0.0 IF(I.GT.NN) GO TO :5 URtJ,I)I0.O CONTINUE DO 16 1-1,N UC(2*I-1,NNP1)I-1.O uc<2:x ,N~92>--1.o UR1NNP1,2*I-1)I1.0 UR(NNP2,2II >-1.o DO 17 1-2,~ uc<2*1-1.~~P3)-v<1>-v<1> uc<2*1 ,NNP3)IX(I)-X(1) UR(NNPZ,2*I-1)IY(1)-Y(I) UR(NNP3,2*I )-x<1>-x<1) SOME CONSTANTS PR2-2.#PR PR$=3.:PR F7I-(3.-PR)/2. PI-ACOS(—1.) LOOP ON COLUMNS 00196:C 00197: 00198: 00199: 00:00: 00201: 00202: 00203: 00204: 00205: 00206: 00207: 00208: 00209: 00210: 00211: 00212: 00213:C 00214: 00215: 00216: 00217: 00218:C 00219: 00220: 00221: 00222: 00223: 00224: 00225: 00226: 00227: 00228: 00229: 00230: 00231: 00232: 00233: 00234: 00235:C 00236: 00237: 00238: 00239:C 00240: 00241: 00242: 00243: 00244:C 00245: 00246: 00247: 00248: 0024928 00250: 00251: 00252: 00253: 00254:C 00255: 00256: 00257: 0025 : 00259:C 00260: 00261: 107 D0 19 J-1,N JJ:2*J IF(J.EO.1) THEN JH1-N JJN2=NN ELSE JN1=J-1 JJH2-JJ-2 END IF IF(J.E0.N) THEN JP1'1 JJPl-l ELSE JPl-J+1 JJPI8JJ+1 END IF UC(JJ-1,JJH2)-UC(JJ-1,JJM2)+2.-PR2 UC(JJ ,JJN2-1)-UC(JJ ,JJH2-1)-2.+PR2 UC(JJ-1,JJP1+1)-UC(JJ-1,JJP1+1)-2.+PR2 UC(JJ ,JJP1)8UC(JJ .JJP1)+2.-PR2 xn-X(J)-X(JH1) YM=Y(J)-Y(JH1) SF-SORT(XM*XH+YM*YM) CTF-YH/SF STFs-XM/SF - XH-X(JP1)-X(J) YM=Y(JP1)-Y(J) SSBSQRT(XH§XH+YH*YM) CTSIYH/SS 5T5--XM/SS ARG=CT84CTF+STS*STF IF(ARG.GT.1.0) THEN OH-PI44. ELSE 0M=(FI-ACOS(ARB))44. END IF AA:(2.+PR2)*(STS*CTS-STF*CTF) BB‘(2.+PR2)*(CTS*CTS-CTF§CTF) CCC(2.-PR2)*LOG(SF/SS) UC(JJ-1.JJ-1)-OM+AA UC(JJ,JJ-l)--BB-CC UC(JJ-1,JJ)--BB¢CC UC(JJ,JJ)-OH-AA AAII.-LDG(SF) BB-STF98TF CCI~CTF§STF DD'CTFiCTF UR(JJ-1,JJ-1)-AA*(3.-PR)+BB*(1.+PR) UR(JJ,JJ-1)=CC*(1.+PR) UR(JJ-1,JJ)-UR(JJ,JJ-1) URtJJ,JJ)-AA*(3.-PR)+DD§(1.+PR) AA-1.-LOG(SS) EB=STS*STS CCS-CTS*STS DD=CTS*CTS URtJJ-l,JJP1)=AA*(3.-PR)*BB*(1.+PR) UR(JJ.JJP1)=CC5(1.+FR) 00262: 00263: 00264:C 00265: 00266: 00267: 00268: 00269: 00270: 00271: 00272: 00273: 00274: 00275: 00276: 00277: 00278:C 00279:C 00280:C 00281: 00282: 00283: 00284: 00285: 00286: 00287: 00288:C 00289:C 00290:C 00291: 00292: 00293: 00294: 00295: 00296: 00297: 00298: 00299: 00300: 00301: 00302: 00303:C 00304: 00305: 00306: 00307:C 00308: 00309: 00310: 00311:C 00312: 00313: 00314: 00315:C 00316: 00317: 00318: 00319:C 00320: 00321: 00322: -- ’ C’\_.‘.-u~)._-: 00:24:18 00:25:19 00326: 00327:C 108 URtJJ~1.JJF1+1)=UH(JJ,JJF1) UR1JJ.JJP1+1)CHH*(S.-PR)+DU‘(1.+PR) xn-(X(J)+X(JH1))/:. YM=(Y(J)+Y(JM1))/2. BI=X(J)-xn E2=Y(J)-YM BB-E::BI+B2¢B2 P8-829(2.+PR2) P9-51.(2.+PR2) P10-824(1.-PR) P11=81*(3.+FR) P12=B1*(1.+PR3) F13=82v(3.+PR) P14-52.(1.+PR3) P15-B1*(1.-PR) LOOP ON RONS DO 19 1'1,N IF(I.EQ.J.OR.I.EO.JH1) GO TO 19 II-2*I AI'X(I)-XH A2=Y(I)-YH AASAlfiA1+A29A2 A882.*(A1*BI+A2*BZ) LOOP ON POINTS OF INTEGRATION DO 18 L-1,LL RR-AA-AB*R(L)+BB*R(L)0R(L) Axax-Az-azofitL) UF3-A181/RR UF68A131*UF3 UFl-UF6*UF3 AlBl-A2-82*R(L) UF73A181§UF3 UF4-A181/RR UF8-A1819UF4 UF2-UF49UFB UF5=LOG(RR) EE.(PB*UF1+P9§UF2+P109UF3-P11*UF4)*W(L) UC(II-!,JJH2-1)-UC(II-1,JJH2-1)+H1(L)*EE UC(II-1,JJ-1)-UC(II-1.JJ-1)+W2(L)OEE EE-(P9*UF1-P8*UF2-P12*UF3+P13*UF4)§N(L) UC(II-1.JJH21-UC(II-1,JJH2)+N1(L)*EE UC¢EE EE‘(P9*UF1-P9*UF2-P11‘UF3+PI4*UF4)*W(L) UC(II,JJH2-1).UC(II,JJH2-1)+N1(L)*EE UC(II,JJ-1)-UC(II,JJ‘1)+N2(L)*EE EE-(-P8*UF1-P9#UF2+813*UF3-P15GUF4)*W(L) UC(II,JJH2)-UC(II,JJH2>+EE9N1(L) UC(II,JJ)=UC(II,JJ)+w2(L)*EE URtII-1,JJ-1)-UR(11—1.JJ-1)+(P7*UF5+(1.+PR)*UF6>¢N(L)/2.0 UR(11-1.JJ)-UR(II-1,JJ)+(1.+PR)9UF7*W(L)/2. UR(11.JJ—1)-UR(II-1,JJ) UR(II.JJ)-UR(II,JJ)+(P7«UF5+(1.+PR)*UFB)4H(L)x2.0 CONTINUE CONTINUE CALL CLUCH(F1N2) 109 003:8:L MUDIFY Uh 00:29zc 00330: DU 22 1:1,N 00331: II=2¢I 00332: 11fl1811-1 00333: A=1.0 00334: DO 20 L82,N 00335: UR(IIH1,1)-UR(IIH1.1)+UR(IIN1.2*F-I)§A 00336: UR1 CALL CLOC}(F1N4) CALCULATIONS FOR FIELD POINT‘. NEITE(10.39) FORHAT(/,5x,‘RESULTS FOR FIELD POINTS'.//,6x,'COORDINATES'. +15X,'DISFLACEHENTS'.13X,'STRESSES',//) READ<8.*)IA,XF,YF IF(IA.NE.1) GO TO 38 DO 31 1-1,NN DO 31 k-1,5 UC(K.I)-0.0 UR(k,I>-0.0 DO 32 J-1,N JJ-2*J 1F(J.EO.1) THEN JJH2-NN JJHS-NN-I JHI-N ELSE JJH2-JJ-2 JJM3-JJ—3 Jnl-J-I END IF xn-(X(J)+X(JMI))/2. YH-(Y(J)+Y(JH1))/2. 81-X(J)-XH 82-Y(J)-YM 88-81*B1+82*B2 P8-82*(2.+PR2) P9-81*<2.+PR2) P10-82*(1.-PR) P11881*(3.+PR) P12-81*(1.+PR3) P13-B2*(3.+PR) P14-82¢(1.+PR3) PIS-819(1.-PR) Al-XF-XM A2-YF-YH AA-A1*A1+A2*A2 AB'2.*(A1*BI+A2*BZ) LOOP ON POINTS OF INTEGRATION DO 32 L-I.LL RR-AA-AB*R(L)*BBOR(L)*R(L) A191-A1-81*R(L) UF3 OAIBIIRR UFb DAlalfiUFS UFI 'UFbiUF3 AIBl-A2-32§R(L) UF7 IAIBI'UFS UF4 CAIBI/RR UFO IAIBIQUF4 UFZ 'UF4fiUFB UFS ILOG(RR) UF9SUF39UF3 UFIOIUF4§UF4 UF1181.0/RF\' UF12‘UF6'UFB*UF11 UF13=UF69UF11 UF14-UF89UF11 00460:L 00461: 00462: 00463: 00464:C 00465: 00466: 00467: 00468:c 00469: 00470: 00471: 00472:C 00473: 00474: 00475: 00476:C 00477: 00478: 00479: 00480: 00481:C 00482: 00483: 00484: 00485: 00486:C 00487: 00488: 00489: 00490: 00491:C 00492:C 00493: 00494: 00495: 00496: 00497:C 00498: 00499: 00500: 00501: 00502: 00503:C 00504: 00505: 00506: 00507: 00508: 00509:C 00510:C 00511: 00512: 00513: 00514: 00515:C 00516832 00517:C 00518:C 00519:C 00520: 00521: 00522: 00523: 00524 : 111 EE=(F6‘UF1+F9iUF:*FlD*UF3'?!ICUF4)ON(L) UC(1,JJH2-1)=UL(1,JJM2-l)¢wl(L)OEE UC‘1.JJ’1)‘UC(1.JJ-l)*N2(L)OEE E£=(F9iUF1-FSGUF2-P12iUF3+P13*UF4)¢w(L) UC(1.JJH2)-UC(1.JJM2)+w1(L)§EE UC(1.JJ)-UC(1,JJ)+N2(L)*EE EE‘(P9¢UF1-PB*UF2-F11*UF3+P14¢UF4)GW(L) UC(2,JJH2-1)-UC(2,JJM2-1)+N1(L)*EE UC(2,JJ-1)-UC(2,JJ-1)+WZ(L)*EE EEB(-P8§UF1-P9:UF2+P13*UF3-P15*UF4)9w(L) UC(2,JJH2)IUC(2.JJM2)+EE*N1(L) UC(2,JJ)-UC(2,JJ)+N2(L)§EE UR<1.JJ-11-UR(1,JJ-1)+(P7GUF5+(1.+PR)*UF6)cw(L)/2.0 UR(1.JJ)-UR(1,JJ)+(1.+PR)OUF7QN(L)l2.0 UR(2,JJ-1)-UR(1.JJ) UR(2,JJ)-UR(2,JJ)+(P79UF5+(1.*PR)*UF8)§W(L)/2.0 EE=(2.0+PR2)*UF1 UR(3,JJ-1)IUR(3,JJ-1)-(EE+(1.O-PR)*UF3)QW(L) UR(4,JJ-1).UR(4,JJ-1)+(EE-(l.0+3.0§PR)*UF3)*N(L) UR(5,JJ )IUR(5,JJ )+(EE-(3.0+PR)§UF3)*W(L) EE=(2.0+PR2)*UF2 UR(3,JJ )8UR(3,JJ )+(EE-(l.0+3.0*PR)*UF4)*H(L) UR(4.JJ )IUR(4,JJ )—(EE+(1.0-PR)§UF4)§N(L) UR(5,JJ-1)IUR(5,JJ-1)+(EE-(3.O+PR)§UF4)§N(L) EE-(2.*82*(UF11+4.*UF13-8.*UF1*UF3)“4.981*UF7*(UF11-4.*UF9)) +*N(L} UC(3,JJM3)IUC(3,JJH3)+EE*N1(L) UC(3,JJ-1)IUC(3,JJ-1)+EE*N2(L) EEI(4.*B2*UF7§(UF11-4.0*UF9)-2.*81*(UF1l-B.O*UF12))*N(L) UC(5,JJH3)-UC(5,JJH3)+EE*W1(L) UC(5,JJ'1)-UC(5,JJ-1)+EE§N2(L) UC(3,JJH2)-UC(3.0JH2)+EE*N1(L) UC(3,JJ )IUC(3,JJ )+EE*N2(L) EE-((UF11-8.0*UF12)*2.*82-UF7*(UF11-4.0*UF10)*4.§81)*W(L) UC(4,JJH3)-UC(4,JJM3)+EE*W1(L) UC(4.JJ-1)-UC(4,JJ-1)+EE*N2(L) UC(5,JJH2)-UC(5,JJH2)+EEfiw1(L) UC(5,JJ )-uc<5,aa )+EE*H2(L) EE-(4.9829UF7*(UF11-4.*UF10)-2.*Bl4(UF11+4.*UF14-B.‘UF2*UF4)) +*W(L) UC(4,JJM2)-UC(4,JJH2)+EE*N1(L) UC(4,JJ )IUC(4,JJ )+EE*N2(L) CONTINUE MODIFY UR AND OUTPUT FIELD RESULTS. A81.0 DO 34 K=2,N l.¥=2*li LLHI‘kL-l DU 33 L31.5 UR (L. l)=UF'< (L. 1 )+UH(L.K!M1 )OA 00526:33 005:7:34 00518: 00529: 00530: 00531: 00532: 00533:35 00534:C 00535: 00536: 00537: 00538: 00539: 00540:C 00541: 00542: 00543: 00544: 00545: 00546:36 00547: 00548: 00549: 00550: 00551: 00552: 00553:37 00554:C 00555: 00556:38 00557: 00558: 00559: 00560: 00561: 00562: 00563: 00564: 00565: 00566:C 112 UH:L.2)=UH+2.0§UR(L.K:H1) UR(L,KK )--UR(L,KL-2)+2.09UR(L,LK ) DISIIO.O DISZBQ.O SIGIIO.O SIGZ‘0.0 516380.0 DD 36 I-1,NN DISlsDIS1-UC(1,I)*SAV(I)+UR(1,I)*F(I) D152=DIS2-UC(2.I)48AV(1)+UR(2,I)*F(I) 5161-8161-(1.0+PR)*UC(3.I)*SAV(I)+UR(3,I)*F(I) 5182-5182-(1.0+PR)*UC(4,1)OSAV(I)+UR(4,I)*F(1) SIG3=SIG3-(1.0+PR)OUC(5,I)§SAV(I)+UR(5,I)*F(I) DISl-DIS1/(8.0*PI) 0152-0182/(8.0*PI) 818188181/(8.0*PI) SIGZsSIGZ/(8.0*FI) 5183-5163/(8.0*PI) wRITE(10,37)XF,YF,DIS1,DISZ.SIG1,SIG2,5163 FORMAT(7<3X,F10.5),/) IF(IA.EQ.1) GO TO 30 CALL CLOCK,NEN(I> IS THE NUMBER OF NODES OF THE ELEMENT, AND (NN(1,J),J-1.NEN(I)) ARE THE NUMBERS OF THE NODES NHICH FORM THE ELEMENT(LISTED COUNTER- CLOCKNISE).THESE PROPERTIES ARE READ IN FREE FORMAT. X(I,1),NT8C(I.1),X(I,2).NTBC(I,2) LL NODE PROPERTIES. X(I,1) IS THE COORDINATE x1 OF NODE I, NTBC(I.1) IS THE TYPE OF BOUNDARY CONDITION IN THE x1 DIRECTION AT NODE I, X(I.2) IS THE COORDINATE x2 OF NODE I. NTBC(I,2) IS THE TYPE OF BOUNDARY CONDITION IN THE X: DIRECTION AT NODE I. ENTER FOR I=1,NNO, ONE SET OF NODE FROFERTIES FER LINE. ENTER x<1,1) IN CO- LUMNS 1 THROUGH 15 IN F FORMAT, NTBC(I,1) IN COLUMNS 16 THROUGH 20 IN I FORMAT. X(I.2) IN COLUMNS 21 THROUGH 35 IN F FORMAT AND NTBC(I,2) IN COLUMNS 36 THROUGH 40 IN I FORnAT. THE BOUNDARY CONDITION CODE FOR NTBC(I.J) Is AS FOLLONS NTBC(I,J) --1 IF SPECIFIED ZERO DISPLACEMENT AT NODE I IN DIRECTION J. NTBC(I.J) - 0 RECTION J. NTBC(I.J) c 1 IF SPECIFIED NON-ZERO DISPLACEMENT AT NODE I IN DIRECTION K. THE DEFAULT VALUE FOR NTBC(I,J) IF SPECIFIED FORCE AT NODE I IN DI-. IS ZERO. NUMBER OF POINTS IN EACH DIRECTION USED FOR NUMERICAL INTEGRATION. LL MAY BE 2 OR 3. 00065:C fluubbgc 000:7:C 00068:C 00069:C 00070:C 00071:C 00072:C 00073:C 00074:C 00075:C 00076IC 00077:C 00078:C 00079:C 00080:C 00081:C 00082:C 00083:C 00084IC ALSO, I,II.R = R I: NEI,NE2= 115 I’Uh EHLH LOITL' LASE,I THROUGH NLO: THL NODE NUMBER.DIRECTION AND VALUE OF EACH SPECIFIED NON-ZERO FORCE OR DISPLACEMENT. ENTER IN FREE FORMAT AND END BY INFUTTING 0.0.0.0. LOCAL COORDINATE FOR STRESS SOLUTION IN ELEMENTS TO BE DESIGNATED BELOW. R MUST LIE BETWEEN 0. AND 1. ENTER IN FREE FORMAT. NOTE THAT R80. IMPLIES THAT STRESSES ARE TO BE COMPUTED AT CENTER OF ELEMENTS ONLY WHEREAS R‘I. IMPLIES THAT STRESSES WILL BE COMPUTED AT THE FOUR CORNERS. INTERVALS OF ELEMENTS IN WHICH STRESSES ARE TO BE COMPUTED IN ASCENDING ORDER. STRESSES WILL BE COMPUTED IN ELEMENTS NE! THROUGH NEZ. N51 THROUGH NE2,ETC. ENTER IN FREE FORMAT AND END BY INPUTTING 0,0. 000853Cfififi§fifififiififiGO§§§§§§§§§§*§§§§§fifiii§§§&*§§G§§§*§*O*§§§O§O§§*O§§§§§§§§OQQ 00086:C 00087: 00088: 00089: 00090: 00091: 00092:C 00093: 00094: 00095:C 00096: 00097: 00098: 00099:C 00100: 00101: 0010220 00103: 00104: 00105xc 00106: 00107: 00108:C 00109: 00110: 00111: 00112: 00113:10 00114: 00115:C 00116: 00117: 00118: 00119: 00120: 00121:C 00122: 00123: 00124:C 00125: 00126: 001273C COMMON MED,MSD COMMON / CFROFIL/ NMAT.NEL.NSD,NR,NNO,NDOF,NEO,NEN(200) +.NN(200,4),NTBC(200,2) COMMON /DIAG/JDIAG(500>,IDIAG(100) COMMON / CRED/E(3),PR(3),G(3),MP(200),X(200,2) OPEN(UNIT=B,FILE='DATA1') OPEN(UNIT-10,FILE¢'RESULTI') CALL CALL CALL CALL CALL CALL CALL CALL CALL CLOCK(START) RED1(NLO.LL) CLOCK(FINI) PROFIL! CLOCK(FIN2) PFORM CLOCK(FIN4) DO 10 I-I.NLO CALL MODIFY(I) CALL SOLVE(.FALSE.,.TRUE.,I) CALL STRESS CONTINUE CALL CLOCK(FINS) NRITE(10,§)'CPU TIME AFTER READ 8 '.FIN1-START NRITE(10,§)’CPU TIME AFTER PROFIL 3 '.FIN2-START WRITE(IO,*)'CPU TIME AFTER PFORM = '.FIN3-START NRITE(10,*)'CPU TIME AFTER SOLVEI I '.FIN4-START WRITE(10,*)'TOTAL CPU TIME - ',FINS-START CLOSE(UNITIIO) CLOSE(UNIT=8) CALL EXIT END 00128:CROOOQDDQODOSAOoa#9GGQ...§Q9§§§O§§§§.§§.Oo9o§§.o§§o.ooooG§O§¢O#o.§§§ 00129:C 001301C SUBROUTINE RED) DUIBIIC 116 Uu1323(.099...a...O90.QODOOOOOOOOOOOOOQOOOOOaa9.99.«O’DooOOOQQOOOOOOOQQOQO 00133:C 00134: 00135:C 00136: 00137: 00138: 00139: 00140: 00141:C 00142: 00143:1 00144:C 00145: 00146:C 00147: 00148:2 00149:C 00150: 00151: 00IS2I3 00153:C 00154: 00155:4 00156:5 00157:C 00158: 00159:C 00160: 00161:6 00162:C 00163: 00164:7 00165: 00166: 00167IC 00168: 00169:8 00170: 00171:C 00172: 00173I9 00174:C 00175: 00176:10 00177: 00178:C 00179:C 00180: 00181: 00182:11 00183: 00184: 00185:12 00186:13 00187:C 00188: 00189:14 00190: 00191:C 00192: 00193:IS 00194:IO 0019S:C 0019c: SUBROUTINE RED1(NLD.LL) COMMON /CFROFIL/ NMAT.NEL.NSD,NR,NNO,NDOF,NEO,NEN<2OO> +,NN(200,4).NTBC(200.2) COMMON /DIAG/ JDIAG(500),IDIAG(100) COMMON / CRED/ E(3),PR(3),G(3),MP(200),X(200,2) DIMENSION TITLE(20) READ(8,1) TITLE FORMAT<20A4) READ(8,*) NMAT,NEL,NNO,NDOF,NLO READ(8,2)(E(I),PR(I).I-1,NMAT) FORMAT(E20.12,F20.12) DO 3 1-1.NEL READ(8,G)MP(I),NEN(I),(NN(I,J),J-1,NEN(I)) CONTINUE DO 4 I-1.NNO READ:S,S) FORMAT FORMAT,I-I.NMAT) FORMAT(//,15,13X,E14.7,F20.6) NRITE(10,10> FORMAT(//,3X,‘ELEMENT NUMBER',4X,'MATERIAL NUMBER',4X, +'NUMBER OF NODES',4X,'NODE NUMBERS') DO 13 I-1,NEL wRITE(10,11>I,MF(I),NENII),(NNII.J),J-1,NEN(I>> FORMAT(/,4X,IS.2(I5X,15),7X,4(15)) IF(NEN(I).LE.4) GO TO 13 NRITEIIO,12>(NNII,J).J-5,NEN(I)) FORMAT CONTINUE NRITE(IO,14> FORMAT(/,3X,'NODE NUMBER‘,6X, +‘NODAL COORDINATES AND BOUNDARY CONDITIONS') DO 15 I-I.NNO wRITEIIO,IG)I,; WFITE(IU.I7)LL 00197:17 LII :1 98: 0019¢:C 00200: 00201: 00202:: C 117 FORMAT(/I,'NO. OF POINTS OF INTEGRATION IN EACH DIRECTION +,15,//) RETURN END 00:03;C§CQ§§§O§§Al96§§§§§§0§0¢9000Ofififit0.00GCQQOGOOGOOGGOGO§§§0OOfifitibOOOCOQO 00204:C 00205:C 00206:C SUBROUTINE PROFILI 00207;CQOOOOQGOROOOOOOOOOOGGOOOO9999*.OOOOOOOOOOOOOOOOGGOOOOOOOOODQOOOOOOQOOO 00208:C 00209: 00210:C 00211: 00212: 00213: 00214:C 00215:C 00216:C 00217: 00218: 00219: 00220: 0 221: 00222:C 00223: 00224: 00225:C 00226: 00227: 00228: 00229:C 00230: 00231: 00232: 00233: 00234: 00235: 00236: 00237:1 00238:C 00239:C 00240:C 00241: 00242: 00243: 00244: 00245: 00246:C 00247: 00248: 00249:C 00250:2 00251: 00252: 00253: 00254: 00255: 00256: 00257: 00258: 00259:3 00260:4 00261: 002OZ:C SUBROUTINE PROFILI COMMON /CFROFIL/ NMAT,NEL,NSD,NR,NNO,NDDF,NEQ,NEN(200) +,NN(200,4),NTBC(200,2) COMMON /DIAG/JDIAG(500),IDIAG(100) COUNT EQUATIONS, INITIALIZE JDIAG, IDIAG NEO'O NSD-0 DO I N'I,NNO DO 1 I-I.NDOF J8NTBC(N,I) IF(J.LT.0) THEN NTBC(N,I)-0 ELSE IF(J.E0.0) THEN NEQ-NEO+1 NTBC(N.I)-NEO ELSE NSD‘NSD-I NTBC(N.I)-NSD ND-Z'IABS(NSD) IDIAG(ND-1)-20000 IDIAG(ND)-O END IF CONTINUE COMPUTE COLUMN LENGTHS OF STIFFNESS VECTOR DO 10 N-I,NEL DO 9 1-1,NEN(N1 II-NN(N,I) IF(II.EQ.O) GO TO 9 DO 8 K-1.NDOF rK-NTGC(II,K> IF(KL)2.8,5 kV-IABS(KK) DO 4 J-I,NEN(1> JJ-NN(N,J) IF(JJ.EQ.0) GO TO 4 DO 3 L-I,NDOF LLsNT8C(JJ,L) IF(LL.LE.0) GO TO 3 IDIAG(2RHK—I)-HINO(IDIAG(Ztrr-I),LL) IDIAG(ZGKF)IMAYO(IDIAG(2RVL),LL) CONT INUE CONTINUE GO TO 8 0020?:5 002b4: 00265: 00260:C 00257: 00268: 00269:C 00270: 00271: 002 2: 00273: 00274:C 00275: 00276: 00277: 00278: 00279: 00280: 00281:6 00282:7 00283:8 00284:9 00285:10 00286:C 00287:C 00288:C 00289: 00290: 00291: 00292: 00293:11 00294: 00295: 00296:14 00297: 00298:12 00299:C 00300: 00301: 00302: 00303: 00304:13 00305: 00306:C 00307: 00308 : 00309:C 00310:C 00311:C 00312: 00313: 00314: 00315: 00316: 00317: 00318: 00319: 00320:C 00321: 00322: 00323: 00324: 00325: 00326: 00327: 00328: 118 D0 7 J=1.NFN(N) JJ=NN(N.J> IFIJJ.EU.0) 60 T0 7 DO 6 L81,NDOF LL-NTBC(JJ,L) IF(LL.LT.0) THEN LL-IASSILL> IDIAG<2«LL-II-MINOIIDIAGI2.LL-1),RL) IDIAG(2*LL1-MAX0(IDIAG(2FLL),Kh) ELSE IF(LL.GT.0) THEN M=MAXO+1 NAD-JDIAG(NEQ) NRITE(10,14)NEQ,NAD FORMAT(/,3X,'NO. OF EQUATIONS',IS,3X,‘COMPONENTS OF STIFFNESS'. +15) IFINSD.EQ.O)RETURN NED-1 NSD-IABS(NSD) IF(NSD.EQ.1) RETURN DO 13 N‘2,NSD IDIAG(2’N)IIDIAG(2*N-2)+1+IDIAG(2*N)-IDIAG(2*N-1) NED'IDIAG(2*NSD) RETURN END FUNCTIONS MINO AND MAXO FUNCTION MIN011,K) IF(I.LE.K) THEN MINOsI ELSE MINOBK END IF RETURN END FUNCTION MAXO(1,K) IF(I.GE.F) THEN MAXO=I ELSE MA¥0=1 END IF RETURN END 00329:C 00:30;COCO§OOCOOG§§IQCOIGOQGfibi§§*9§‘§filfii§§ifii§§fiifilift!GOiO‘.§QQIOiii0.00990 00331:C ‘:’(’3. s" o 00333:C 00:34;COO...5*.OOOOOofiaOOOOOOGOOOGOOO969.6.99.50GOO...OOOOOO§96§¢OOOOOOOOCQ§OO 00335:C 00336: 00337:C 00338: 00339: 00340: 00341: 00342: 00343: 00344:C 00345:C 00346:C 00347: 00348:2 00349: 00350: 00351:3 00352:C 00353:C 003548C 00355:100 00356: 00357: 00358: 00359:4 00360:C 00361:C 00362:C 00363: 00364:C 00365:C 00366:C 00367: 00368: 00369: 00370: 00371:C 00372: 00373: 00374: 00375: 00376: 00377:5 00378: 00379:6 00380: 00381: 00382:7 00383:8 00384:C 00385:C 00386:C 00387: 00388: 00389:9 00390:C 00391:C 00392:C 00393: 00394:C 119 SUBROUTINE PFOHN SUBROUTINE PFORM(LL) COMMON ICPROFIL/ NMAT.NEL.NSD.NR,NNO,NDOF,NEO,NEN(200) +,NN(200,4),NTBC(200.2) COMMON IDIAG/ JDIAG(500),IDIAG(100) CONMON lABST/ A(7000).B(500>,F(400) COMMON /CRED/ E(3),PR(3),G(3),MP(200),X(200,2) DIMENSION D(3),XE(4.2),NBCC4,2).S(8,8) INITIALIZE A AND 8 DO 2 1‘1,JDIAG(NEO) A(I)'0.0 ‘ IF(NSD.EQ.0) GO TO 100 D0 3 1'1,IDIAG(2*NSD) B(I)30.O CALCULATE MATERIAL PROPERTIES DO 4 N=I,NMAT NI-RR-E(N)/(1.0—PR(N)*§2) PR(N)-E(N)*N1 B(N)-E(N)§(1.0-N1)/2.0 LOOP ON THE ELEMENTS DO 10 N81,NEL LDCALIZE PROPERTIES IN ARRAYS xE HMIHP(N) D(1)'E(MH) D(2)-PR(HM) D(3)IG(HM) DO 8 1-1,NEN(N) .II-NN(N,I) IF<11.NE.0) GO TO 6 DO 5 L-:,NDOF xE(I,L>-o.o N8C(I,L)-0 GO TO 9 DO 7 L-1,NDOF XE(I,L)-X(II.L) N8C 00395:C 00396:C 00397: 00398:10 003998C 00400:C 00401:C 00402: 00403: 00404: 00405: 00406: 00407: 00408: 00409 : I I 00410:C 00411: 00412: 00413:C 120 ADD THE VALUES OF LOCAL STIFFNESS INTO GLOBAL CALL ADDETIF(S,NLC,NLN(N>,NDOF> CONTINUE REARRANGE THE ARRAY NTBC FOR REACTIONS NR=0 DO 11 1-1,NNO DO 11 K=1,NDOF KK-NTBC S(J1+1,K1)-S(J1+1,Kl)+W12#SHP(K,1) S(J1+1,K1+1)-S(J1+1,K1+1)+N12*SHP(K.2) K18K1+NDOF J18J1+NDOF NSL-NINiNDOF REARRANGE PRODUCTS INTO STIFFNESS MATRIX DO 4 J=I.N9L.NDOF Do 4 r=J.NSL.NDOF 00461: 00462: LND46T4 00464: 00465:C 00466: C104 ’57 3 00468: 00469: 00470:C 004718 00472; 00473: 00474:4 00475:C 00476: 00477: 00478:C 121 ”1185:.1,1 ) H12=S(J,k+1) NL1-S\J+I.D) N22=S(J41,L+1) $1.] ,h)-=I’>(1)iN11+D(3)§NZ‘2 S(J,k91)=D(2)¢N12+D(3)9N21 S(J+1,k)sD(2)Aw21+D:3)§wI2 S(J+1,L+1)-D(1)ow22+D(3)Ow11 S(K,J)-S(J,P) S(h,J+1)=S(J+1,K) S(K+1.J)-S(J,k+1) S(K+1,J+1)-S(J+1.K+1) RETURN END 00479:Ctifiic.0999.OQOOOOOOOOOOO§94§O0.90449.49OAOAOQCOOOAOOOAuoOAO§§§§OOAOAO§O 00480:C 00481:C 00482:C SUBROUTINE SHAPE 00483:C§O§**§§§§O*§*§O§§Q§{#0.0..§§§§§§§§§§¥*§§§§§GO§§§§O§§§§§O§§0§O6fififi§0909* 00484:C 00485: 00486:C 00487: 00488: 00489:C 00490: 00491:C 00492: 00493:C 00494: 00495: 00496: 00497:1 00498: 00499:C 00500:C 00501:C 00502: 00503:2 00504:C 00505:C 00506:C 00507:3 00508: 00509: 00510:C 00511: 00512:4 00513:C 00514: 00515:C 00516: 00517: 00518: 00519: 00520:C 00521:C 00522:C 00523: 00524: . . r -.=' ‘.M..—"-’: 00526:5 SUBROUTINE SHAPE(SS,TT,XE,SHP,XSJ,NDOF.NIN) DIMENSION SHP(4,3),XE(4,2),S(4),T(4),XS(2,2) +,SX(2,2) DATA S/-065'005'005g-00S/QT/—OOSQ-OOS,OOSQOOS/ FORM 4 NODE OUADRILATERAL SHAPE FUNCTIONS DO 1 1.1.4 ‘ SHF’(I ,3)8(O.5+S(I)*SS)§(O.5+T(I)*TT) SHP(I,1)'S(I)*(O.5+T(I)*TT) SHP(I,2)-T(I)*(0.5+S(I)*SS) IF(NIN.GE.4) GO TO 3 FORM TRIANGLE SHARE FUNCTIONS BY ADDING THIRD AND FOURTH TOGETHER DO 2 1:1,: SHP(3,I)-SHP(3,I)+SHP(4,I) CONSTRUCT JACOBIAN AND ITS INVERSE DO 4 I-I,NDOF DO 4 J-I,2 XS(I,J)-0.0 DO 4 K-I,NIN XS(I,J)-XS(I,J)+XE(K,I)*SHP(K,J) xsa-x5(1,1>.x5(:,2>-xs<1.2).x5(2,1) SX(1.1)-XS(2,2)/XSJ SX(2.2)-XS(1,1)/XSJ SX(1.2)--XS(1,2)IXSJ SX(2,1)--XS(2,1)IXSJ FORM GLOBAL DERIVATIVES DO 5 181,N1N TR:5HR(I.1).SI(I,I>+SHP*SY(2.1) 5H;Ix‘g):SHE(1,I)asx:1,2)+SHF(I,2)9SX(2»3) SHF(I.1)'TF 00527:C 00528: 00519: 005301C 005:1gCiOOOOQOAOOOOOOOOOOOOO§QC§QOA6999.6...o.QOAOOOOOOOOOOOO9656996§§9¢O§9§a9 00532:C 00533:C 00534:C 00535:C9§¢O§OO¢OiOOOOQOOQQOOOOOQQOOOOO96.9..9...09906696966949;9.6969999.It... 00536:C 00537: 00538:C 00539: 00540: 00541: 00542: 00543:C 00544:C 00545:C 00546: 00547:C O0548:C 00549:C 00550:2 00551: 00552: 00553: 00554:3 00555: 00556:C 00557:C 00558:C 00559:4 00560: 00561: 00562: 00563: 00564:5 00565:C 00566: 00567: 00568:C 00569;C6**6.Q«#4....§§§Ofiuii9000*DOOGOOODOODA65.666906...90*O§OOOOOOOOGOAO§§ 00570:C 00571:C 00572:C 00573;COGQGO§§§§§§*§**O§§§O§§§O*§*§§*‘§*§§§O§*§G§§§§*§§§O§§§§§§§§G§§§§§OQ§§QOG 00574:C 00575: 00576:C 00577: 00578: 00579: 00580:C 00581: 00582: 00583: 00584: 00585: 00586: 00587:C 00588:C 00589:C (30591:: g (.1059 1 3 00592: 122 RE I URN EIJL' SUBROUTINE GAUSS SUBROUTINE PGAUSS(L,R,Z,N) DIMENSION LR(9),LZ(9),LN(9) +,R(9),z(9>,w(9) DATA LR/-1.1,1,-1,0.1.0,-1.0/,LZ/-1,-1,1,1,-1.0,1,0.0/ DATA LW/4i25.4*40,64/ GAUSS POINTS AND HEIGHTS FOR TNO DIMENSIONS GO TO (2,4).L 2’2 INTEGRATION G-1./SORT(3.) DO 3 LI-l,4 R(LI)-G*LR(LI) Z(LI)=G*LZ(LI) N(LI)'I. RETURN 3&3 INTEGRATION G'SORT(O.6) H-1./81. DO 5 LI*1.9 R(LI)-G*LR(LI) Z(LI)'G*LZ(LI) N(LI)IH*LN(LI) RETURN END SUBROUTINE ADDSTIF SUBROUTINE ADDSTIF(S,NBC,NIN,NDOF) COMMON IABST/ A(7000),B(500),F(400) COMMON /DIAG/ JDIAG(500),IDIAG(100) DIMENSION S(S,S).NBC(4,2) JCIC) DO 3 J=1,NIN DO 3 JJ-1,NDOF JCIJC+I k-N8C(J,JJ) IF(K.LT.0) THEN LOCAL STIFFNESS IN 8 FOR SFECIFIED NON-ZERO DISPLACEMENTS k-IAFS(1) IF(F.EO.I) THEN L=1 0:359: : C'Utv‘r 4‘ 1,1115ng 00596: 00597: 00598: 00590; (:QbCu‘n 00601 : 00602: 00603: 00604: 00605: 00606:1 00607:C 00608:C 00609:C 00610: 00611: 00612: 00613: 00614: 00615: 00616: 00617: 00618: 00619: 00620:2 00621: 00622: 00623:C 00624: 00625: 00626:C 123 ELSL L‘IDIAG(2O(L—1))*1 END IF LL=IDIAE(2.F-1) L‘L-LL IC=0 DO 1 III,NIN DO I II'I,NDOF IC=IC*1 M8NBC(I.II) IF(M.LE.O) GO TO I M=L+M B(M)=B(M)+S(IC,JC) CONTINUE LOCAL STIFFNESS IN A FOR SPECIFIED FORCE ELSE IF(K.GT.0) THEN L-JDIAG(K)-K 18-0 DO 2 I-I,NIN DO 2 II-I,NDOF ICaIC+I M=N8C(I,II) IF(M.BT.K.OR.M.LE.0) GO TO 2 N=M+L , A(M)-A(M)+S(IC,JC) CONTINUE END IF CONTINUE RETURN END 00527:C§§§§§§§§*9.9§§§94§¢§§§§6*§6§9§§O6§§§§§§§§¢¢*§§*§§O§¢§§§§§6§§§OOOOOOOOOO 00628:C 00629:C 00630:C SUBROUTINE SOLVE 00531gC9696460699.666046490096060§066600§06666466O666660666606669§66666966666. 00632:C 00633: 00634:C 00635: 00636: 00637: 00638: 00639: 00640: 00641:C 00642: 00643:C 00644: 00645: 00646:C 00647:C 00648:C 00649: 00650: 00651: 00652: 00653: 00654:1 00655: 00656: 00657:C 00658:C SUBROUTINE SOLVE(AFAC,BACK.LL) LOGICAL AFAC,BACK COMMON IABST/ A(7000).C(500),B(400) COMMON IDIAG/ JDIAG(500),IDIAG(100) COMMON ICPROFIL/ NNAT,NEL,NSD,NR.NNO,NDOF,NEO,NEN(200) +,NN(200,4),NTBC(200,2) DIMENSION R(2) PUT STIFFNESS A(I) INTO TRIANGULAR FORM AENGY'O. C’ JRsO LOOP ON COLUMNS DO 6 J-1,NEO JD-JDIAG(J) JH-JD-JR ISzJ-JH+2 1F(JH-2)6.3,1 IF(.NOT.AFAC) GO TO 5 IEcJ-I r=JR+2 LOOP ON ROWE Unbquc 00551:): 00661: 00662: 00663: 00664: 00665:: 00666:3 00667: 00668: 00669: 00670: 00671: 00672: 00673: 00674: 00675: 00676:4 00677:5 00678:6 00679: 00680:C 00681:C 00682:C 00683: 00684: 00685:' 00686:7 00687: 00688: 00689:8 00690: 00691: 00692: 00693: 00694: 00695: 00696: 00697:9 00698:10 00699: 00700:11 00701: 00702:12 00703: 00704: 00705: 00706:C 00707:C 00708:C 00709: 00710: 00711: 00712: 00713: 00714: 00715: 00716: 00717: 00718: 00719:13 00720:!4 00721:15 00722:C 00723: 00724: 124 D0 2 1315.16 Ifi=JDIAG(1-1) ID=JDIA6<11 IH=HINO(ID-IR-1.I-IS+1) 1F(1H.GT.0)A(I)=H(r>-DOT(A(K—IH).A(ID-1H>,IH) K2L41 IF(.NOT.AFAC) GO TO 5 IR-JR+1 IE=JD-1 KSJ‘JD DU 4 I-IR,IE ID-JDIAG(K+I) IF(A(ID).E0.0.0)GO TO 4 D-Atl) A(I)OA(I)/A(ID) A(JD)-A(JD)-D§A(I) CONTINUE IF(BACK)B(J)-B(J)-DOT(A(JR+1),8(IS-1),JH-1) JR-JD IF(.NOT.BACK> RETURN SOLVE FOR LOAD VECTOR IN 8 DO 7 1-1,NEO ID-JDIAGII) IF(A(ID).NE.0.0) B(I)-S(I)/A(ID) AENGY-AENGY+8(I)§B(I)§A(ID) J-NEO JD-JDIAG(J) D-B(J) J-J-1 1F(J.LE.0) GO TO 11 JR-JDIAG(J) IF(JD-JR.LE.1) GO TO 10 IS-J-JD+JR+2 K-JR-IS+1 DO 9 I-IS.J B(I)-B(I)-A(1+K)*D JD-JR GO TO 8 CONTINUE NRITE(10,12)LL FORMAT(//.3X,'FINAL DISPLACEMENT FOR LOAD NO.'; +2X,IS,//,3X,'NODE NUMBER',10X,‘DISPLACEMENT X',8X, +‘DISPLACEMENT Y‘) CALL CLOCK(FINISH) PRINT DISPLACEMENTS DO :4 J-I.NND DO 13 JJ-I,NDOF JR-NTBC(J,JJ) IF(JR.LT.-NSD) THEN R-o.o ELSE IF(JR.GT.0) THEN R(JJ)IB(JR) ELSE R(JJ)-E(NEO+IABS(JR)) END IF CONTINUE NEITEtIo.Is>J.R,R<2) rORNAT(/.5x,15.I:x,FI:.7,9x,F13.7; RETURN EIJD 00725:C 00726:C 007271C 00728: (11.:729: 00730: 00731: 00732:1 00753: 00734: 00735:C 00735;Ca§o¢99a9.9ab§d§9§9§§09§§§§6¢OOOO§566§OQO§¢§¢§99099oo.§99.§99§§§99§§§09. 00737:C 00738:C 00739:C 00740;C§§090§§§§§§§§§§§O§§§§§§O§§O§§fiffififi0§§§§§*§§§G§O§§§§§§§§§&*§§*§§§§.§§§§§ 00741:C 00742: 00743:C 00744: 00745: 00746: 00747: 0074838 00749:C 00750:C 00751: 00752:1 00753: 0075 :2 00755: 00756: 00757:3 00758: 00759: 00760:C 00761: 00762: '00763: 00764:4 00765:C 00766: 00767: 00768:C 00769: 00770: 00771:' 00772: 00773:5 00774: 00775: 00776:C 00777:C 00778:C 00779:C 00780:6 00781: 00782: 00783: 00784: 00785: 00786: 00787: 00786: 00789: Cnfi7903 125 FUNCIION DOT FUNCTION DOT(A,b,N) DIMENSION A(N),B(N) DOT=0.0 DO 1 1=1,N DDT:DOT+A(1>OB(I) RETURN END SUBROUTINE MODIFY SUBROUTINE HODIFY(LL) COMMON IABST/ A(7000),B(5001,F(400) COHHON IDIAG/ JDIAG(500),1DIAG(100) COMMON ICPROFIL/ NMAT,NEL,NSD,NR,NNO,NDOF,NEO,NEN(200) +,NN(200,4),NT8C(200,2) READ AND PRINT NON-ZERO SPECIFIED FORCES AND DISPLACEMENTS DO 1 K-1.(NEO+2*NSD+NR) F(k)-0.0 NRITE(10,2)LL FORMATtIH,//,5X.'LOAD CASE NO.’,2X.15,//, +5X,‘NON-ZERO SPECIFIED DISPLACEMENTS OR FORCES‘, +//,5X,‘NODE NUMBER'.2X,'COORDINATE NO.’,5X,'COND1TION AND VALUE') READ(8.*) 1.11.8 IF(I.EQ.0) GO TO 6 1K8NTBC(I,II> IF(IK.GT.0) THEN F(IK)-R HRITE(10,4)I,II.R FORMAT(//,8X,15,8X,I5.10X,'FORCE'.12X.F10.5 ELSE IFtIK.LT.-NSD) THEN WRITE(10,§)'ERROR NODE ',I,' ',II,‘ 15 ZERO DISPLACEMENT‘ ELSE 1KIIABS(1K) F(IK+NEQ)-R NRITE(10,5)I,II,R FORMAT(//,8X,15.8x,15,10X,‘DISPLACEHENT'.5X.F10.5) END IF GO TO 3 MODIFY RIGHT HAND SIDE OF EQUATION FOR SPECIFIED NON-ZERO DISPLACEMENTS IF(NSD.EQ.O)RETURN DO 8 III.NEO JR-O DO 7 J81,NSD JJ=IDIAG(2‘J-1) LiIIDIAG(2*J) JH=11~JR J1=JJ-I IF(JL.GI.O) GO TO 7 IF(JK.LT.0.AND.I.BT.(JJ+JH‘1)) GO TO 7 IF‘J1.ED.0) THEN (“20791 : [NJ-79:: 00793: 00794: 00795:? 00796:8 (113797: C 00798: 00799: 00800:C 126 F(1)=F(I)-F(NLU+J)¢8(JH*1) ELSE F<1)=F(1)-F(NEO+J)¢E(JR+1+l-JJ) END IF JRBLL CONTINUE RETURN END 00801gC§§§§§§*§§§§§O§§§I§G§§§§§i§§§§§§§§i9}!Qi§§§§§§§§§b§§§§§§i§§§§*§*§Q§§§§§f 00802 : C 00803:C 00804:C SUBROUTINE STRESS 00805; C§§§§§§§§i§§§§*§O§§*§§§§§§§§§§§§§§§§§fifl§§§*§§§§¢§§O§§§§§O§§O§§§§§§§§§§§ 00806:C 00807: 00808:C 00809: 00810: 00811: 00812: 00813: 00814: 00815:C 00816:C 00817:C 00818:C 00819:C 00820: 00821: 00822: 00823: 00824: 00825: 00826: 00827: 00828: 00829: 00830: 00831: 00832:1 00833:C 00834: 00835: 00836: 00837: 00838: 00839: 00840: 00841: 00842: 00843: 00844: 00845: 00846: 00847:C 00848: 00849: 00850: 00851: 00852: 00853: 00854: 00855: 00856: SUBROUTINE STRESS COMMON / CPROFIL/ NHAT,NEL,NSD,NR,NNO,NDOF,NEO.NEN(200) +,NN(200,4),NTBC(200,2) COMMON ICRED/ E(3),PR(3).G(3),HF(200),X(200,2) CONHON IABST/ A(7000),8(500),F(400) DIMENSION D(3),XE(4,2).UE(8),STR(3),SIG(3) +,SG(4),TG(4),SHP(4,3),S(8,8) LOCALIZE MATERIAL PROPERTIES.NODAL COORDINATES AND NODAL DISPLACEMENTS FOR EACH ELEMENT 1K-0 READ(8,*)R READ(8,*)NEl,NE2 SBtI)=-R SG(2>-R SG(3)-R 58(4)--R TG(1)--R TG(2)--R TG(3)-R T8(4)-R NR1TE(10,:) FORMAT(//.'ELEMENT NO.',3X,'COORDINATES',20X,'STRESSES') DO 10 I-1,NEL IF(IK.EO.1) THEN READ(8,§)NE1,NE2 IK-O ELSE CONTINUE END IF IF(NE2.EO.I)IK-1 n-HE<1) D(1)-E(M) D(2)-PR(H) D(3)-G(H) N1N=NEN(I) DO 2 J81.NIN JJ-NN(1,J) J2=2¢(J-1) DO 2 K-1,NDOF J2=J2+1 1F(JJ.E0.0) THEN XE(J,1.)=0.0 UE-x(JJ,rJ IL-N18C(0J.t) IF(11.LT.-NSD) THEN UE(J:)-0.0 ELSE IF(1:.GT.0) THEN UE(J2)-F(Kr) ELSE UE(J2:-F(NEO+1ABS(kk)) END IF END IF CONTINUE CALCULATE REACTIONS "2(2). DO 3 LL=1,8 DO 3 LLL‘I,B S(LL,LLL)-0.0 DO 5 J81,NEN(I) JJ’NN (I ,J) 0232*(J-1) DO 5 K-1,NDOF J2=J2+1 KkiNTBC(JJ,K) IF(LK.ST.0) 60 TO 5 IF(H.ED.O)CALL ELEMENT(0,XE,S,NEN(I),NDOF,3) M=1 HCCNEQ+NSD+IAES(KK) DO 4 L-I,8 F(MC)-F(HC)+S(J2,L)*UE(L) CONTINUE CALCULATE STRAINS AND STRESSES IF(NE1.E0.0) GO TO 10 IFtI.LT.NEI.OR.I.GT.NE2) GO TO 10 LL-4 IF(NIN.EO.3.0R.R.E0.0.0) LL-I DO 8 L'I,LL CALCULATE SHAPE FUNCTIONS AND THEIR DERIVATIVES CALL SHAFE(SG(L),TB(L),XE,SHP,XSJ,NDOF,NIN) DO 6 K¢1,3 ' STR(K)=0.0 XC80.0 YC-0.0 DO 7 J81,NEN(I) J2-2*J J1-J2-1 XC-XC+SHP(J,3)GXE(J,I) YC'YC+SHF(J.3)*XE(J.2) STR(1)'STR11)+SHF'(J,1)::UE(J1) STR(2)-STR(2)+SHP(J,2)*UE(J2) STR(3)=STR(3)+SHR(J,1)*UE(J2)+SHP(J,2)*UE(J1) SIG(1)-D(I)*STR(1)*D(2)*STR(2) SIG(2)8D(2)§STR(I)+D(I)*STR(2) SIG(3)ID(3)*STR(3) WRITE(10,9)I.XC.YC.SIB(1),SIG(2).SIG(Z) CONTINUE OUTPUT OF REACTIONS 128 00923: WRITE(10.11) 00924311 FORMAT(//,3X,'NODE NO.'.3x.'DIRECTION‘,31.'REACTION',/) 00925: DO 13 l-1,~ND 0092b: DO 13 :81.NDOF 00927: kk-NTBC(I.:; 00928: 1F(K|.GT.0) GO TO 13 00929: HCsNEO+NSD+IABStKL) 00930: R-F(HC) 00931: NRITE(10.12)I,h,R 00932x12 FORMAT(3X,15,6X,15,3X,F12.6,/) 00933:13 CONTINUE 00934:C 00935: RETURN 00936: 'END BOTTOM S ("91;") 1 3 C CIUCICIZ ; C 129 00003;C9§O§§§§§O§§§§§G§i§§§§i§§§§iffil¥fii§§§9§§§§§Oii§§9§§90§§§ilhiiiiil§§§ 000043C 00005: 00006:C 00007:C 00008:C 00009; C 00010xc 00011:C 00012:c 00013:C 00014zc 00015:C 00016:C 00017:C 00018:C 00019:C 00020:C 00021:C 00022:C 00023:C 00024:C 00025:C 00026:: 00027:C 00028:C 00029:C 00030:C 00031:C 00032:C 00033:C 00034:C 00035:C 00036:C 00037:C 00038:C 00039:C 00040:C 00041:C 00042:C 00043:C 00044:C 0004S:C 0004b:C 00047:C 00048:C 00049:C 00050:C 00051:C 00052:C 00053:C 00054:C 00055: 0005¢:C 00057:C 00058:C 00059:C (3006C: : C 000:1:0 00062:C 00063:C PROGRAM FEMBOU THIS PROGRAM EMPLOYS A COMBINATION OF THE FINITE ELEMENT AND DIRECT BOUNDARY ELEMENT METHODS FOR THE SOLUTION OF PLANE STRESS PROBLEMS OF LINEAR ELASTICITY. IS ASSUMED. UNIT THICYNESS FINITE ELEMENTS MAY BE DEFINED BY THREE OR FOUR NODES AND BEM-GENERATED ELEMENTS MAY BE DEFINED BY ANY ODD NUMBER OF NODES GREATER THAN FOUR. NE DEFINE THE FOLLOWING INPUT FARAMETERS,VARIABLES AND ARRAYS: TITLE NMAT NNO NDOF NEL NLO THE TITLE OF THE PROBLEM. THIS MUST BE WRITTEN ON ONE LINE IN COLUMNS I THROUGH 80. THE NUMBER OF MATERIALS OF WHICH THE BODY IS COMPOSED. ENTER IN FREE FORMAT. THE TOTAL NUMBER OF NODES IN THE MODEL. FREE FORMAT. THE NUMBER OF DEGREES OF FREEDOM OF THE PROBLEM. THIS PARAMETER IS EQUAL TO 2 FOR PLANE PROBLEMS.ENTER IN FREE FORMAT. THE NUMBER OF ELEMENTS USED. ENTER IN FREE FORMAT. THE NUMBER OF LOAD CASES TO BE SOLVED FOR THE SAME GEOMETRY. ENTER IN FREE FORMAT. ENTER IN E(I),FR,NE~(I),(NN(I,J),J-I,NEN'cru TIME WRITEI10,*)’CPU TIME WRITE(10..)'CPU TIME AFTER READ c'.FIN1-START AFTER PROFIL t',FINZ-STAFT AFTER PFORM =’.FIN?-START AFTER SOLVE(I) c'.FIN4-START 00130: 00131:C 00132: 00133: 00134: 00135: 00136:C 00137:C999{1090*}*9§§*O§§§O§OO§§§9.0§§§O§§O§§O§§OOO§O§OOO§§OO§O§OO§OOQQ§O§ 00138:C 00139:C 00140xc C101 41: C§§§§O§§§§§§*§§Q§§D***§§O§§§O§*fii...’96*‘G.’*§§O§O§§.§QQ‘I‘Iii‘li-CNIO*5. 00142:C 00143: 00144:C 00145: 00146: 00147: 00148: 00149: 00150:C 00151: 00152:1 00153:C 00154: 00155:C 00156: 00157:2 00158:C 00159: 00160: 00161:3 00162:C 00163: 00164:4 00165:5 00166:C 00167:C 00168:C 00169: 00170:6 00171:C 00172: 00173:8 00174: 00175: 00176:9 00177:C 00178: 00179:10 00180: 00181:C 00182: 00183: 00184: 00185: 00186: 00187:11 00188: 00189: 00190: 00191:12 00192:13 00193:C 00194; 00195:14 131 NAITETID,.I TOTAL CFO TIME .-,r1N5_51gR1 CLOSE CLOSETUNITxa) CALL EXIT END SUBROUTINE REDl SUBROUTINE REDltNLO) COMMON ICFROFIL/ NMAT,NEL,NSD.NR,NNO.NDOF,NEO.NEN(100) +,NN(100,30),NTBC(200,2) COMMON IDIAG/ JDIAG(300),IDIAG(100) COMMON / CRED/ E(3),PR(3),G(3),MP(100),X(200,2 DIMENSION TITLE(20) READ(8,1) TITLE FORMAT(20A4) READ(9,§) NMAT,NEL,NNO,NDOF,NLO READ(8,2)(E(1),PR(1),I-1,NMAT) FORMAT(E20.12,F20.12) DO 3 I-1,NEL READ(S,*)MP(I),NEN(I),(NN(I,J),J=1,NEN(I)) CONTINUE DO 4 1-1,NNO READIS,S)) FORMAT(/,56X,415) CONTINUE WPITE(IO,I4) FORMAT(/.3X.'NODE NUMBER'.6X. LILI19¢ : (ILI197; 0019S:Ib 00199:16 001T"):C 00201: 00202 3 00203:C 132 *‘NODHL COORDINATLS AND BOUNDARY CONDITIONS ) DU 15 I=I,NNU WRITE(10,161I,(XtI,J),NTBC(I,J),J=1,NDUF) FORMAT(/,3x,15,5x,3tF15.5,15)) RETURN END 002043COO...OGOOGOOOGDOOOOOfiiOOOOOQOOO999.90.00.99...GOO99900099...O‘DOQOO‘OO 00205:C 00206:C 00207:C SUBROUTINE PROFILI 00208:C§O§Q§§§iflffiifiifilfifiifiiiiiii.‘§§§§i§§§§*0ifi’OGGOOOIOQ§GI§§fi*§§ii§§li§§§& 00209:C 00210: 00211:C 002 2: 00213: 00214: 00215:C 00216:C 00217:C 00218: 00219: 00220: 00221: 00222: 00223:C 00224: 00225: 00226:C 00227: 00228: 00229: 00230:C 00231: 00232: 00233: 00234: 00235: 00236: 00237: 00238:1 00239:C 00240:C 00241:C 00242: 00243: 00244: 00245: 00246: 00247:C 00248: 00249: 00250:C 00251:2 00252: 00253: 00254: 00255: 00256: 00257: 00258: 00259: ("3260: 3 00261:4 SUBROUTINE PROFILI COMMON /CPROFIL/ NMAT,NEL,NSD,NR,NNO,NDOF,NEO,NEN(100) +,NN(100,30),NTBC(200,2) COMMON /DIAG/JDIAG(300),IDIA8(100) COUNT EQUATIONS, INITIALIZE JDIAG, IDIAG NEO=O NSDSO DO 1 N-1,NNO DO 1 I-1,NDOF J=NTBC(N,I) IF(J.LT.0) THEN NTBC(N,I)-0 ELSE IF(J.E0.0) THEN NEO‘NEQ+1 NTBC(N,I)INEQ ELSE NSDBNSD-I NTBC(N.I)-NSD ND-2*IABS(NSD) IDIAG(ND‘11820000 IDIAG(ND)-O END IF CONTINUE COMPUTE COLUMN LENGTHS OF STIFFNESS VECTOR DO 10 N-I,NEL DO 9 1-1,NEN(N) II-NN(N,I) IF(II.E0.0) GO TO 9 DO 9 K-1,NDOF Kk-NTBC(11,K) kk-IABS(KK) DO 4 J=I,NEN(1) JJ-NN(N.J) IF(JJ.EQ.0) GO To 4 DO 3 L-1,NDOF LL-NTBC(JJ,L) IF(LL.LE.0) GO TO 3 IDIAG(2SFk-1)=M1NO(IDIAG(ZoKk-I).LL) IDIAG(29LLI3HAIO(IDIAG(2*FV),LL) CONT INUE CONTINUE 00262: 00263:C 00264:5 00265: 00255: 00267:C 00268 : 00269: 00270: C 00271: 00272: 00273: 00274: 0027':C 00276: 00277: 00278: 00279: 00280: 00281: 00282:6 00283:7 00284:8 00285:9 00286:10 00287:C 00288:C 00289:C 00290: 00291: 00292: 00293: 00294:11 00295: 00296: 00297:14 00298: 00299:12 00300:C 00301: 00302: 00303: 00304 : 00305:13 00306: 00307:C 00308: 00309: 00310:C 00311:C 00312:C 00313: 00314: 00315: 00316: 00317: 00318: 00319: 00320: 00321:C 00322: 00323: 00324: 00325: 00326: 00327: 133 EL: TO 8 D0 7 J=1,NEN(N) JJ=NN(N,J) Ifth.ED.0) DO TO 7 DO a L:I,NDOF LL-NTSCIJJ,L) IF(LL.LT.0) THEN LL=IABS(LL) IDIAG(29LL*1)8MINO(IDIAG(2*LL-1),KF) IDIAG(2*LL)¢MAXO(IDIAG(2*LL),KK) ELSE IF(LL.GT.0) THEN M=MAXO=MAXO(JDIAG(MI,IABS(KK-LL>) ELSE CONTINUE END IF CONTINUE CONTINUE CONTINUE CONTINUE CONTINUE COMPUTE TOTAL LENGTH OF STIFFNESS VECTOR NAD=1 JDIAG(1)-1 IF(NEQ.EO.1) GO TO 12 DO 11 N’2,NEQ JDIAG(N)-JDIAG(N)+JDIAG(N-I)+1 NADBJDIAG(NEQ) WRITE(10,141 NEQ,JDIAG(NED) FORMAT(/,'NO. OF EQUATIONS',IS,3X,'COMPONENTS OF STIFFNESS', +15) IF(NSD.E0.0)RETURN NED-1 NSD-IABS(NSD) IFO=I END IF 00328: (M 177.19 : OUT-3H! C 00:318C.999...OOQOOQOOObfiffiififffiQOifiifi§iihfiifiO999.9.fiOQOOOOOQOG.OQOOOOOOOOQOOO 00332:C 00333:C 00334:C 00:35:C§§§§I§*§§C§.O§§§O9.0.0.999G9.9....G990.§§§O§§§OOOOOO§§O§OOQ駧§O§§§OQGQ 00336:C 00337: 00338:C 00339: 0034CH 00341: 00342: 00343: 00344: 00345: 00346:C 00347: 00348: 00349:1 00350 : 00351:C 00352:C 00353:C 00354: 00355:2 00356: 00357: 00358:3 00359:C 00360:C 00361:C ‘00362x4 00363: 00364: 00365: 00366:5 00367:C 00368:C 00369:C 00370: 00371:C 00372: 00373:C 00374:C 00375:C 00376: 00377: 00378: 00379: 00380:C 00381: 00382: 00383: 00384: 00385: 00386: 00387: 00388: 0038°:6 CHJ39CH 00391: 0039:: 00393:8 134 RETURN END SUBROUTINE PFORM SUBROUTINE PFORM COMMON /CFROFIL/ NMAT,NEL.NSO,NP,NNO,NOOF,NEO.NEN(100) +,NN(100,30),NTBC(200,2) COMMON IDIAG/ JDIAG(300>,IDIAG(100) COHMON lABST/ A(7000),B(500),F(500) COMMON ICRED/ 5(3),PRtS),6(3),MP(100),X(200,2) COMMON ICBOUN/ 88(3364) DIMENSION D(3),XE(58),NBC(5B).S(8,8> LL'3 WRITE(10,1)LL FORMAT(//,'NO. OF POINTS OF INTEGRATION ON EACH DIRECTION', *15,//) INITIALIZE A AND 8 DO 2 1-1.JDIAB(NEO) A(I)-0.0 1F(NSO.E0.0) GO TO 4 DO 3 1-1.IDIAG(2*NSD) B(I)-0.0 CALCULATE HATERIAL PROPERTIES DO 5 N-1,NNAT WISPR(N) E‘N)'E(N)/(1.0-PR(N)§*2) PR(N)-E(N)*W1 B(N1-E(N)*(1.0-w1)l2. LOOP ON THE ELENENTs OPEN(UNIT-9,FILE-'LOCAL') DD 10 N-1,NEL LOCALIZE PROPERTIES IN LOCAL ARRAYS XE(I),NBC(I),D(L) HH-HP(N) D(l)-E(HH) D(2)-PR(HM) D(3)‘G(HN) DO 8 I-1,NEN(N) II-NN(N,I) IF(II.NE.0) 80 TO 6 NBC(I)-O NBC(1+NEN(N))-0 XE(I)-0.0 XE(I+NEN(N))-0.0 GO TO 8 XE(I)=X(II,1) XE(I+NEN(N))8X(II,2) NBC(1)-NTBC(II,1) NBC(I+NEN(N:18NTEC(II.2) CONTINUE 00394:C 00395:C 00396:C 00297: 00398: 00399: 00400x9 00401:C 00402:C 00403:C 00404: 00405: 00406:C 00407: 00408: 00409: 00410:C 00411:C 00412:C 00413: 00414: 00415: 00416: 00417:10 00418:C 00419:C 00420:C 00421:C 00422: 00423: 00424: 00425: 00426: 00427: 00428: 00429:11 00430:C 00431: 00432: 00433:C 135 INITIALIZE LOCAL STIFFNEES MATRIA 1F(NCN(N).LE.4) THEN OO 9 1:1,8 DO 9 LPI,8 S(k.1)=0.0 CALCULATE STIFFNESS OF FINITE ELEMENTS CALL ELENENT(D.XE.S,NEN(N),NDOF,LL) CALL ADDSTIF(S,NBC,NEN(N),NDOF,B) ELSE NN1-2*NEN(N) NNIZ'NNI’NNI CALCULATE STIFFNESS OF SUPER-ELEMENTS CALL BOUN(NEN(N>.NN1.D.XE.SB.E.F,I) CALL ADDST1F(SB,NBC,NEN(N),NDOF,NN1) NRITE(9,§)(SB(I),I-1,NN12) END IF CONTINUE REARRANBE THE ARRAY NTBC FOR REACTIONS NR-O DO 11 1-1,NNO DO 11 Ks1,NDOF KK-NTBC(I,K) 1F(KK.NE.0) GO TO 11 NR-NR+1 NTBC(1,K)--NSD-NR CONTINUE RETURN END 00434:C**§§§§§*§§§.§§§§§§§***§§§§OIO§§i9§§§fi§**§§*§§§§§§*§§i§§§*i§§§§§‘*§§*§§§ 00435:C 00436:C 00437:C SUBROUTINE ELEMENT 0043a:C9509OCCCGCCCOOCCG§CC§§§§§iii§§§9§§QCOOCQOOOOOGOC*0.§§§§§9§§**§§99§§¢§§§ 00439:C 00440: 00441:C 00442: 00443: 00444:C 00445:C 00446:C 00447: 00449: 00449: 00450:C 00451:C 00452:C 00453:C 00454: 00455:C 00456:C 00457:C 00458: 0045°:C SUBROUTINE ELEMENT.SHP(4.3) +,s LINTILL'LL DO 2 L-I,LINT CALCULATE SHAPE FUNCTIONS AND THEIR DERIVATIVES AT THE POINTS OF INTEGRATION CALL SHAPE(SG(L),TB(L),XE,SHF,XSJ,NDOF,NIN) CORRECT ELEMENT OF AREA XSJBXSJ*NG(L) 004:0:L 004OI:L ("34 b: : 00463: 00404: 00465: 00466: 00467:C 00468: 00469: 00470: 00471: 00472: 00473:1 00474:2 00475: 00476:C 00477:C 00478:C 00479: 00480: 00481: 00482: 00483: 00484: 00485:C 00486: 00487: 00488: 00489: 00490:C 00491: 00492: 00493: 00494:4 00495:C 00496: 00497: 00498:C 136 LFEATE PRODUCT OF DERIVATIVES AND STIFFNESS JI‘I DU 2 J‘IgNIN NII'SHF(J,I)OX§J NIZ‘SHP(J,2)§XSJ II'JI DO 1 ksJ,NIN S(JI,k1)-S(J1.k1)+N11*SHP(K.1) S(JI,h1*1)-S(J1,K1+1)+N11*SHP(K,2‘ S(J1+1.F1)=S(J1+1,K1)+N12*SHP=C(:1*N21+D(3)*w12 S(J+1,K+I)=D(1>*N22+D(3)&N11 S(K.J)-S(J.K) S(K,J+1)-S(J+1,K) S(K+1,J)-S(J,K+1) S(K+1,J+1)-S(J+1,K+1) RETURN END 00499; Club.*§§*OO‘Ififi'flfiNfifi‘lO‘IO'IQQiii‘NI-‘D‘Ii.**§§§*§§*§§§§§§O§§‘§*§§*§*i§*§fiifi#:ININI". 00500:C 00501:C 00502:C SUBROUTINE SHAPE 00503:C§§§§§§00§§9§§§§§§§§0§9§§§§§§5§§§§40§o*§*§*0§4§§*§*§§§§§§§§90549§§§O§§§§ 00504:C 00505: 00506:C 00507: 00508: 00509: 00510:C 00511:C 00512:C 00513: 00514: 00515g 00516:1 00517: 00518:C 00519:C 00520:C 00521: 00522:2 00523:C 00524:C 0052::C SUBROUTINE SHAPE(SS,TT,XE,SHP,XSJ,NDDF,NIN) DIMENSION SHP(4,3),XE(NIN,2),S(4),T(4),X8(2,2) *.SX(2,2) DATA S/-005g005’005,-005/,T/-005'-005‘0359005/ FORM 4 NODE OUADRILATERAL SHAPE FUNCTIONS DO 1 1-1.4 SHP(I,3)-(0.5#S=o.o DO 4 K-I,NIN ‘ x5(I.J)-XS(1,J)+XE(r.1)cSHP(L.J> xsa-xs<1,1>«xs<2.2)-xs<1.2)«x5(2.1) SX(1.1)-XS(2,2)/XSJ SX(2.2)=XS(1.1)/X5J 5x:1,2)--xs\1,2>/x5J SX(2,1)--XS(2,1)/XSJ FORM GLOBAL DERIVATIVES DO 5 I-1,NIN TP-SHP(I.1)§SX(1,1>+SHF(I,2)§SX(2.1) SHP(I.2)-SHP(I,1)*SX(1,2)+SHP(I,2)*SX(2,2) SHP(I.1)-TP RETURN END 00550:C§§§¢§§§Qfi§§*.§9§§§*¢§§§*§§&§§*4§§4§§§§§§*§§¢§666464§§fifi¢¢§§66644;566446 00551:C 00552: 00553:C SUBROUTINE PGAUSS 00554;C§§§*§6*.£5§§§6§§§§94466§§§6§§4*966**§§#66§66§6§§§§§§¢§§4§¢§§6*§.§§66§§o 00555:C 00556: 00557:C 00558:. 00559: 00560: 00561: 00562:C 00563:C 00564:C 00565: 00566: 00567: 00568: 00569: 00570:5 00571:C 00572: 00573: 00574:C SUBROUTINE PGAUSS(L.R,Z,N) DIMENSION LR(9),LZ(9),LN(9) +,R(9),Z(9),U(9> DATA LR/‘igl‘lg-IQC’gI'Og-l .0/,LZ/’1,-1 ’1'11-1.0’l '0‘(:)/ DATA LN/4*25,4§40,b4/ 3.3 INTEGRATION G-SORT(O.6) HiI./81. DO 5 LI-1,9 R(LI).BOLR(LI) ZCLIIIG§L2(LI) N(LI)=H*LN(LI) RETURN END 00575:COffifiifififififififi§§O§§§*§§§i§§§f§§§§fi§0§§§§§§O*§§*OiifiGGQOGOfifii‘OGGGO.900! 00576:C 00577:C 00578:C SUBROUTINE ADDSTIF 005793CO§§§**§§§Q9*.§§§§§§§§§il§§§§O§§OOQ§§O*§§*OO§¢O{#00*§§6§§§O§§§O§§§§§§§§§ 00580:C 00581: 005 2:8 00583: ' 00584: 00585: 00586:C 00587: 00588: 00589: 00590: 00591: SUBROUTINE ADDSTIF(S,NEC,NIN,NDOF,NN1> COMMON IABST/ A(7000).8(500),F(500) CONMON /DIAB/ JDIAG(300).IDIAG(100) DIMENSION S(NNI,NN1).NBC(NIN,2) 3:80 DO 3 J=I,NIN L3 3 JJ=I.NDOF JC=JC+I L=NBC(J.JJ) 00592: 00593: 00594:L 00595:C 00596: 00597: 00598: 00599: 00600 : 00601 : 00602: 00603 : 00604: 00605: 00606: 00607: 00608: 00609: 00610: 00611: 00612:1 00613:C 00614:C 00615:C 00616: 00617: 00618: 00619: 00620: 00621: 00622: 00623: 00624: 00625: 00626:2 00627: 00628:3 00629:C 00630: 00631: 00632:C 138 IF IF(M.LE.0) GO TO 1 HIL+H B(H)-8(H)+S(IC,JC) CONTINUE LOCAL STIFFNESS IN A FOR SPECIFIED FORCE ELSE IF(K.GT.0) THEN L'JDIAG(K)-H 1831:) DO 2 I81,NIN DO 2 II-1,NDOF IC‘IC+I HINBC(I,II) IF(H.GT.K.OR.H.LE.O) GO TO 2 H-H+L A(H).R(H)+S(IC,JC) CONTINUE, END IF CONTINUE RETURN END 00633;Cf’fiifif’*9.§§*§*§§9§§§§§§§§9§§§Q§§§§§§§§fififilfifl0§§§§0§0¢§§§f§§§*§*§§§§*§ 00634:C 00635:C 00636:C SUBROUTINE SOLVE 006373C§§§§§§§§§QO§ON{OOQ§§0§§§****9*§§§§*§9§§***§*§§*§§*§'6444649646644446644 00638:C 00639: 00640:C 00641: 00642: 00643: 00644: 00645: 00646: 00647:C 00648:C 00649:C 00650: 00651: 00652:C 00¢53:C 00654:C 00655: 00656: 00657: SUBROUTINE SOLVE(AFAC,EACK,LL) LOGICAL AFAC.8ACK COMMON IABST/ A(7000),C(500),8(500) COMMON IDIAG/ JDIAG<300),IDIAG(100) COMMON ICPROFIL/ NHAT,NEL,NSD,NR,NNO.NDOF,NEQ.NEN(100) +,NN(100,30).NTBC(200,2) DIHENSION R(2) PUT STIFFNESS A IN TRIANGULAR FORM AENGY=0.0 JR30 LOOP ON COLUMNS DU 6 J=1,NEO JD=JDIAG(J) JH=JD-JR (3.3558 g 0065~: 00660:1 00661: 00662: 00663:C 00664:C 00665:C 00666: 00667: 00668: 00669: 00670: 00671x2 00672:3 00673: 00674: 00675: 00676: 00677: 00678: 00679: 00680: 00681: 00682:4 00683:5 00684:6 00685: 00686:C 00687:C 00688:C 00689: 00690: 00691: 00692:7 00693: 00694: 00695:8 00696: 00697: 00698: 00699: 00700: 00701: 00702: 00703:9 00704810 00705: 00706:11 00707: 00708:12 00709: 00710: 00711: 00712:C 00713:C 00714:C 00715: 00716: 00717: 00718: 00719: 00720: (H.177 I.“ I : 00722: (”'17: T. - - - - . 139 IS=J-JH+2 IF(JH-:)b,3.1 IF(.NUT.AFAC) IE=J-1 K'JR92 60 TO 5 LOOP ON RUNS D0 2 I-I8,IE IR=JDIAG(I-1) ID=JDIAG(I) IH=HINO,IDIAG(100) CONMON /CPROF1L/ NMAT,NEL,NSD,NR,NNO.NDDF,NEO.NEN(100) +.NN(100,30),NTBC(200,2) READ AND PRINT NON-ZERO SPECIFIED FORCES AND DISPLACEMENTS DO I K-1,(NEQ+2*NSD+NR) F(K)I0.0 NRITE:10,2)LL I FORMAT(1H,//,5X,'LOAD CASE NO.',2x,15,//, +5X,’NON-ZERO SPECIFIED DISPLACEHENTS OR FORCES'. +//,5X,'NODE NUMBER',2X,'COORDINATE NO.’,5X,'CONDITION AND VALUE') READ(B,*) I,II,R IF(I.E0.0) GO TO 6 IKONTBC(I.II) IFIIK.OT.0) THEN F(IK)-R wRITE(:o,4)I.II,R FORMAT(/l,8x,15,8x,15.10X,'FDRCE',12X,F10.5) ELSE 1F(1K.LT.-NSD) THEN WRITEt10,*)'ERROR NODE '.1,' '.II,' IS ZERO DISPLACEHENT‘ ELSE Ik-IA88(IV) F(IK+NEQ)-R NRITE(10,5)1,11,R FORMAT(//,8X.IS.8X,15.10X,'DISPLACEMENT‘,5X,F10.5) END IF GO TO 3 MODIFY RIGHT HAND SIDE OF EQUATION IF SPECIFIED NON-ZERO DISFLACEMENTS IF(NSD.E0.0)EETURN DO 8 1'1,NEO JR=0 DU 7 J=1.NSD 141 00790: JJ‘IDIAO(ZCJ°1) 00791: ILFIDIAC(2¢J> 00792: JH=IL-Jh 00793: JTFJJ-I 00794: 1F(JL.GT.0) 60 T0 7 00795: IF(JL.LT.0.AND.1.GT.(JJ+JH-1)) GO TO 7 00796: IF(JN.E0.0) THEN 00797: F(I)-F(I)-F(NEO+J)IB(JR+1) 00798: ELSE 00799: F(I)3F(1)-F(NED+J)OB(JR+1+1-JJ> 0.3800: END IF 00801:? JR=LK 00802:8 CONTINUE 00803:C 00804: RETURN 00805: END 00806:C 00807;C§§§§§§§§¢§§6*.96.09.66.699.90.96§§§§¢§§o§6§§§0§050666¢§§§§S§6§§§0§§§9§Q 00808:C 00809:C 00810:C 00811:C§§§0§§§§ibiQGG§Q§§§§§§§§§§§§*§§§O§§§§II§§Q.0*.i§§§¥§§6§§OO¥§§§O§§O§§§§§ 00812:C SUBROUTINE STRESS 00813: SUBROUTINE STRESS 00814:C 00815: COMMON / CPROFIL/ NMAT,NEL,NSD.NR,NNO,NDOF,NEO,NEN<100) 00816: +,NN(100,30),NT8C(200,2) 00817: COMMON ICRED/ E(3),PR(3),E(3),MP(100),X(200,2) 00818: COMMON IABST/ A(7000),8(500),F(500) 00819: COMMON ICBOUN/ 88(3364) 00820: DIMENSION D(3),XE(58),UE(58),FE(58).STR(3).SIG(3) 00821: +,SB(4),TB(4),SHP(4,3),S(8,8) 00822:C 00823:C 00824:C LOCALIZE MATERIAL PROPERTIES.NODAL COORDINATES 00825:C AND NODAL DISPLACEMENTS FOR EACH ELEMENT 00826:C 00827:1 FORMAT(//,'ELEHENT ND.',3X,'CODRDINATES',20X,‘STRESSES') 00828:C 00829: RENIND(9) 00830:C 00831: READ(8,§)R 00832: READ(8,*)NE1,NEZ 00833: SG(1)I-R 00834: 88(21-8 00835: 88(31-R 00836: SG(4)--R 00837: TG(1)--R 00838: TG(2)--R 00839: TG(3)-R 00840: TG(4)-R 00841:C 00842: DO 12 I-1.NEL 00843:C 00844: NN12-NEN(I)§NEN(I)¢4 00845: IF(NEN(I).GT.4)READ(9,§)(SB(K),H-1.NN12) 00846: 1F(IK.EQ.1) THEN 00847: READ(8,*)NE1,NE2 00848: 1k=0 00849: ELSE 00850: CONTINUE 00851: END IF 00852: 00853: IF:NE2.EO.I) 11=1 00854: M=HF 00855: D(II=E(M) 00856: 00987: 00858: 00859:C 00860: 00861: 00862: 00863: 00864: 00865:c 00866: 00867: 00868: 00869: 00870: 00871: 00872: 00873: 00874: 00875: 00876:2 00877:C 00878:C 00879:C 00880: 00881: 00882: 00883:3 00884:C 00885: 00886: 00887: 00888: 00889: 00890: 00891: 00892: 00893: 00894: 00895: 00896: 00897:4 00898:C 00899: 00900: 00901:5 00902: 00903:6 00904:C 00905:C 00906:C 00907: 00908: 00909:C 00910:C 00911:C 00912: 00913: 0091 4 : 00915:7 00916:C 00917:C 00918:C 00919:C 00920: 00921: 142 DIC)‘PR(H) L(3)=8(H) NINrNLNII) DD 2 J=1.NIN JJaNNtI,J) J2=2¢(J-1) XE(J)-X(JJ.1) XE(J+NIN)-X(JJ,2) DO 2 KsI,NOOF J2=J2+1 KL-NTBC(JJ.K) IF(KK.LT.-NSD) THEN UE(J2)-0.0 ELSE IFIKK.GT.O) THEN UE(J2)-F(kK) ELSE UE(J2)-F(NED+IABS(KK)) END IF CONTINUE CALCULATE REACTIONS H80 DO 3 J-1,8 DO 3 K-1,8 S(k,J)-0.0 DO 6 J-1,NEN(I) JJ-NN(I.J) J2-2*(J-1) DO 6 K-1,NDOF J2-J2+1 KK-NT8C(JJ.K) IF(KK.BT.O) 80 TO 6 IF(M.E0.0.ANO.NIN.EO.4>CALL ELEMENT(D,XE,S,NEN(I>,NDDF,3) M-I HCINEQ+NSD+IABS(KK) IF(NIN.EQ.4) THEN DO 4 L-1.8 F(HC)-F(HC)+S(J2.L)§UE(L) ELSE DO 5 L-l,2*NIN F(HC)-F(NC)+SB(2*NIN*(L-1)+J2)*UE(L) END IF CONTINUE CALCULATE STRAINS AND STRESSES IF(NIN.GT.4) THEN 1132*NIN CALCULATE FE(I), THE VECTOR OF NODAL FORCES IN THE S. ELEMENTS DO 7 L31.II FE(L)¢0.0 DO 7 h-1,11 FE(L)=FE(L)+SB(II*(K-1)+L)*UE(K) BOUN(...,2) CALCULATES STRESSES AND DISPLACEMENTS IN THE S. ELEMENTS CALL EDUN(NIN.II.D,XE.SB,UE.FE.Z; 1F(I.NE.NEL) NRITE(10.1) 009221C UUQQZg 00924: 00925: 00926: 00927: 00928: 00929: 00930:C 009311C 00934. 3 C 00 33: 00934: 00935:8 00936: 00937: 00938: 00939: 00940: 00941: 00942: 00943: 00944: 00945:9 00946:C 00947: 00948: 00949: 00950:C 00951:10 00952:11 00953: 00954:12 00955: 00956:C 00957:C 00958:C 00959: 00960:13 00961: 00962: 00963: 00964: 00965: 00966: 00967: 00968:14 00969:15 00970:C 00971: 00972: 00973:C 143 ELSL IF(NEI.E0.0) GO TO I: IF11.LT.NEI.DR.I.GT.NE2) GO TO 12 IF (I.EC‘.I)NF;ITE(IO,I) LL84 IF(NIN.ED.3.0R.R.E0.0.0) LL‘I DO IO L31,LL CALCULATE SHAPE FUNCTIONS AND THEIR DERIVATIVES CALL SHAPEISG(L),TG(L).XE.SHF,XSJ.NDDF.NIN) DO 8 r=1,3 STR(K)=0.0 XC=0.0 YC=0.0 DD 9 J-1,NEN(I) J2=2§J J1-J2-1 XC=XC*8HP(J.3)*XE(J) YCsYC+SHP(J,3)§XE(J+NIN) STR(1)-STR(1)+SHP(J.1)GUE(J1) STRtZ)-STR(2)+SHP(J,2)*UE(J2) STR(3)-STR(3)+SHF(J,1)*UE(J2)+SHP(J.2)*UE(JII SIB(I)=D(I)*STR(I)+D(2)*STR(2) SIG(2)-D(2)*STR(I)+D(I)*STR(2) SIE(3)ID(3I*STR(3) WRITE(10,11)I,XC,YC,SIG(1).SIG(2),SIB(3) FORMATt/,15,5(2X,F12.6)) END IF CONTINUE CLOSE OUTPUT OF REACTIONS NEITE(10,13) FORMAT(//,3x,'NODE NO.',3x,'DIRECTION'.3x,'REACTION',/) DO 15 I-1.NNO DO 15 K-1,NDOF KKsNTBCtI,K> IF(KK.GT.0) 80 TD 15 MC~NEO+NSD+IAES RaFIMC) NRITE(10.14)I,K,R FORHAT(3X,IS,6X,IS,3X,F12.6,/) CONTINUE RETURN END 00974:00000509000490.009600005099096490060000000060000*.9.006000094405000966600 00975:C 00976:C 00977:C SUBROUTINE BOUN(N,NN,D,XE,UC.SAV,F,IE) 00978:C§§§Q§0§§§O§§§§§CO§§099.990..00.666.000.00000000006000560006600600660000 00979:C 00980: 00981:C 00982: 00983: 00984: 00985: 00986: 00987: SUBROUTINE BOUN(N,NN,D.XE,UC.SAV,F,IE) DIMENSION UR(61,61) DIMENSION RC6),N(6),W1(6),N216) DIMENSION UC(NN,NN) DIMENSION XE(N.2) *.IFIV(61),D(?) +.x;<1>.YF(I>.800(58).F(58> “0988: 00989:C 0099”: 00991: 00982: ODS”): : 00994: 00995; 00996: 00997: 00998: 00999:C 010001C 01001:C 01002: 01003: 01004:1 01005:C 01006: 01007: 01008: 01009: 01010: 01011: 01012: 01013: 01014: 01015: 01016: 01017: 01018: 01019: 0102012 01021: 01022:C 01023:C 010241C 01 25: 01026: 01027: 01028: 01029: 01030:3 01031:C 01032: 01033: 01034: 01 35: 01036:4 01037:C 01038: 01039: 01040: 01041: 01042:5 01043:C 010441C 01045:C 01046:C 01047: 01048: (1104 9 : 01050: 01051: 01052: 01053: 144 9,3181). Fk=1.0-2.0OD(3)/011) E-0111011.0-FR«FR) P1=AC081-1.0) PR2=:.§FR PRZ=:.9PR F7t-(3.-PR>/2. NNF' 1 II»'NN+1 NNP2=NN+2 NNF3=NN+3 POINTS AND WEIGHTS FOR NUMERICAL INTEGRATION LL86 URITE(10,1)LL FORMAT(//,'POINTS OF INTEGRATION IN SUPERELEMENT',3X,I R<1)--0.93246951420315 R(2)--0.66120938646626 R(3)--0.23861918608319 R<41n-R<3) R(5)I-R(2) R(6)--R(1) w(1)-O.17132449237917 w(2)-0.36076157304813 N(3)-0.46791393457269 w<4)-w(3) w(s>-w(2; N(6)IN(1) DD 2 L-1,LL NIIL)-1.—R(L) w2(L>-1.+R(L) IF(18.ED.2) GO TO 38 INITIALIZE MATRICES DO 3 1.1,NNF3 DO 3 J'I,NNP$ UR(J,I)-0.0 IF(I.GT.NN.OR.J.GT.NN) GO TO 3 UC(J,I).0.0 CONTINUE DO 4 I'I.N UR(2¢I ,NNP2)I1.O UR(NNPI,2*I-1).I.O UR(NNP2,2§I 1.1.0 DO 5 1-2,N UR(2*I-1,NNP3)-XE(1,2)-XE(I.2) UR(2¢1 ,NNF3)-XE(I.1)-XE(1,1) UR(NNF3,2*I-1)-XE(1,2)-XE(I,2) URINNP3,2cI )-XE(I,1)-XE(1,1) LOOP ON COLUMNS DO 7 J-1,N aJ-2«J IF(J.EO.1) THEN JHI‘N JJM2=NN ELSE amz=J-1 c .J ) 01054: 0105:: 01056: 01057: 01058: 01059: 01060: 01061: 01062: 01063:C 01064: 01065: 01066: 01067: 01068:C 01069: 01070: 01071: 01072: 01073: 01074: 01075: 01076: 01077: 01078: 01079: 01080: 01081: 01082: 01083: 01084: 01085:C 01086: 01087: 01088: 01089:C 01090: 01091: 01092: 01093: 01094:C 01095: 01096: 01097: 01098: 01099:C 01100: 01101: 01102: 01103: 01104:C 01105: 01106: 01107: 01108: 01109:C 01110: 01111: 01112: 01113: 01114:C 01115: 01116: 01117: 01118: 01119: 145 JJHE'8JJ-2 END IF Ir (J . ED. N) THEN JFI‘I JJPI”! ELSE JPI‘J+I JJPI‘JJ+I END IF UC(JJ-1,JJH2)-UC(JJ-1,JJH2142.-PRZ UC+BB*(1.+PR) UR UC<11-1,JJH2)-UC(II-1,JJH21+N1(L14EE UC(11-1.JJ1-UC(Il—1,JJ)+N2(L)*EE EE-(P9*UF1-P8GUF2-P11*UF3+P14*UF4)#N(L) UC(II,JJH2-1)-UC(II,JJH2-1)+NI(L)*EE UC(II,JJ-1)¢UC(II,JJ-1)+N2(L)*EE EE-(-P8*UF1-P94UF2+P13*UF3-P154UF4)*w(L) UC(II,JJM2)-UC-UR4UF61~w/2.o UR(II-1,JJ)IUR(II-1.JJ)+(1.+PR)*UF7*w(L)l2.0 URtII,JJ-1)-UR(II-1,JJ) ' UR-0.0 B(NNF3)30.O DD 19 I31,NNR2 IF(IPIV(I).EU.I) GO TO 17 TEN=B£IFIV(I)) EHIF‘IV(I))'E*(I) 01252: 01253:17 01254:18 01255:C 01256: 01257: 01258: 01259: 01260: 01261:19 01262: 01263:20 01264:C 01265: 01266:41 01267:21 01268:C 01269:C 01270:C 01271: . 01272: 01273: 01274:22 01275: 01276:C 01277:C 01278:C 01279:38 01280:39 01281: 01282:30 01283: 01284: 01285: 01286: 01287:31 01288: 01289: 01290: 01291: 01292: 01293: 01294: 01295: 01296: 01297: 01298: 01299:C 01300: 01301: 01302: 01303: 01304: 01305: 01306: 01307: 01308: 01309: 01310: 01311: 01312: 01313:C 01314: 01315: 01316: 01317: 148 8(1)=TEM DO 18 L-I+1,NNP3 81L)-B-UC(I,K) - IF118.ED.1) RETURN CALCULATIONS FDR FIELD POINTS NRITE(10,39) FORMAT(/,5X,'RE8ULTS FOR FIELD PDINTS'.//,6X,'COORDINATES‘, +15X,'DISPLACEMENTS',13X,'STRESSES'.//) READ(8,§)1A,XF(1),YF(1) IF(IA.NE.1) RETURN DD 31 I-I,NN DO 31 K-I,s UC(K,I)-0.0 UR(K,1)-0.0 DO 32 J-1,N JJ-ZGJ IF(J.EO.1Y THEN JJM2-NN JJM3-NN-1 JM1-N ELSE JJn2-JJ-2 JJM3-JJ-3 JH1-J-1 END IF xn-/2.o YM-(XE(J,2)+XE(JM1,2>)/2.0 Bl-XE(J,1)-XM 828XE(J,2)-YM 88-81i81+82*82 P8882*(2.+PR2) P9-81*(2.+PR2) PlO¢82§(1.-PR) P11-81913.+PR) 812-81*(1.+PR3) P13-82*(3.+PR) RI4-32»(I.+RR3> P15-81*(1.-PR1 AI'XF(1)-XM A2‘YF(I)-YM AA-A1*A1+A29A2 A832.*(A1981+A2*82) 01318:C 01319:: 01320:C 01321: 01322: 01323: 01324: 01325: 01326: 01327: 01328: 01329: 01330: 01331: 01332: 01333: 01334: 01335: 01336: 01337: 01338: 01339:C 01340: 01341: 01342: 01343:C 01344: 01345: 01346: 01347:C 01348: 01349: 01350: 01351:C 01352: 01353: 0 354: 01355:C 01356: 01357: 01358: 01359: 01360:C 01361: 01362: 01363: 01364: 01365:C 01366: 01367: 01368: 01369: 01370:C 01371:C 01372: 01373: 01374: 01375: 01376:C 01377: 01378: 01379: 01380: 01381: 01382:C 01383: 149 LOOP ON POINTS OF INTEGRATION DO 32 L91,LL RR-AA-AB¥R(L)+BB*R(L)*R(L) A1813AI-BI*R(L) UF3 'AIBI/RR UFb IAIBIGUF3 UFI 'UFb'UFS A1518A2-82*R(L) UF7 .AIBIOUFE UF4 BAIEI/RR UFB IAIBI'UF4 UF2 =UF4*UFB UF5 ILOG(RR) UF9-UF3*UF3 UFIOIUF4iUF4 UF1131.O/RR UFIZ-UF69UF99UFII UF13-UF6*UFII UFI4'UF8*UF11 EE'(P8*UF1+P9*UF2+PIO*UF3-R11*UF4)*N(L) UC¢1,JJM2-I)-UC(1,JJM2-1)+W1(L)*EE UC(1.JJ-1)-UC(1,JJ-1)+N2(L)*EE EE‘(P9*UF1-P8*UF2-P12*UF3+P13*UF4)*N(L) UC(I,JJM2)IUC(1,JJM2)+W1(L)*EE UC(1,JJ1-UC(1,JJ)+N2(L)*EE EE'(P9*UF1-P8*UF2-PI1*UF3+P14*UF4)§W(L) UC(2,JJM2-I)-UC(2,JJM2-11+N1(LJGEE UC(2.JJ-1)-UC(2,JJ-1)+H2(L)*EE EE-(-P8*UFI-P9iUF2+P13*UF3-P15*UF4)*N(L) UC(2,JJM2)'UC(2,JJM2)+EE*N1(L) UC(2,JJ)-UC(2,JJ)*W2(L)*EE UR(1,JJ-1)-UR(I,JJ-I1+(P7*UF5+(1.+PR)*UF6)*N(L)I2.0 UR(1,JJ)-UR(1,JJ)+(1.+PR)§UF7*N(L)/2.0 UR(2,JJ-1)-UR(1,JJ) UR(2,JJ)-UR(2,JJ)+(P7*UF5+(1.+PR)*UFB)§N(L)/2. EE-(2.0+PR21*UF1 UR(3,JJ-1)IUR(3,JJ-I)-(EE+(1.O-PR)*UF3)*N(L) UR(4,JJ-I)-UR(4,JJ-11+(EE-(1.0+3.0*PR)*UF3)*N(L) UR(5,JJ )IUR(5,JJ )+(EE-(3.0+PR)*UF3)§W(L) EE-(2.0+PR2)*UF2 UR(3,JJ )CUR(3,JJ 1+(EE-(1.0+3.0*PR)*UF4)*N(L) UR(4,JJ )‘UR(4,JJ )-(EE+(1.0-PR)*UF4)*N(L) UR(5,JJ-I)IUR(5,JJ-1)+(EE-(3.O+PR)lUF4)*W(L) EE-(2.*32*(UF11+4.*UF13-8.*UF1*UF3)-4.§BI*UF7*(UF11-4.*UF9)) +*N(L) UC(3,JJM3)-UC(3,JJM3)+EE*w1(L) UC(3.JJ-I)-UC(3,JJ-I)+EE*w2(L> EE'(4.*82*UF7*(UFI1-4.0*UF9)-2.*BI*(UF11-9.0*UF12))*N(L) UC(5.JJM3)-UC(5,JJM3)+EE*N1(L) UC(5,JJ-1)8UC(5,JJ-1)+EE*N2(L) UC(3,JJM2)=UC(3,JJM2)+EE*W1(L) UC(3,JJ )IUC(3,JJ )+EE*N2(L) EE=((UF11-8.0*UF12)*2.*82-UF74(UF11-4.0*UF10)*4.§81)¢N(L) 01384: 01385: 01386: 01387: 0138818 013891C 01390: 01391: 01392: 01393: 013941C 01395132 013961C 0139718 0139818 01399: 01400: 01401: 01402: 01403: 01404: 01405133 01406134 01407: 01408: 01409: 01410: 01411: 01412135 014131C 01414: 01415: 01416: 01417: 01418: 0141918 01420: 01421: 01422: 01423: 01424: 01425136 01426: 01427: 01428: 01429: 01430: 014311 01432137 0143318 01434: 014351 014361 014371C BOTTOM 8 150 UC(4,JJM3)=UC(4,JJM3)+EE¢W1(L) UC(4,JJ-1)-UC(4,JJ-1)+EE¢W2(L) UC(5,JJM2)-UC(5,JJM2)+EE*W1(L) U615.JJ >-UC(5,JJ )+EEow2(L) EE=(4.*BE*UF7*(UF11-4.*UFIO)-2.#819(UF11*4.IUF14-8.*UF2*UF4)1 +*N(L) UC(4,JJM21.UC(4.JJM2)+EE9W1(L) UC14,JJ )3UC(4,JJ )+EE*W2(L) CONTINUE NODIFY‘UR A-1.0 DO 34 K-2,N KK-ZfiK KKMIIKK-1 DO 33 L-1,5 UR(L.1)-UR(L.1)+UR(L,KKM1)*A URTL.2)-UR(L.2I+UR(L.KK 14A a--A DO 35 K-2,N KK-2iK KKM1-KK-1 no 35 L-1,5 UR(L,KKM1)I-UR(L,KK-3)+2.0*UR1L.KKM1) UR(L,KK )--UR+2.0*UR(L,KR ) DISl-0.0 DISZ'0.0 SIGI-0.0 SIG2-0.0 SIG3-0.0 no 36 I-I,NN 0181-0131-U811,I)*SAV(1)+UR(1,I)§F(I) DIsz-Dzsz-UCI2.I):EAV¢ETI> 8181-8161-(1.0+PR)§UC(3,I)*8AV(I)+UR(3,I)*F(I) 8182-8182-(1.0+PR)¢UC(4,I)!SAV(I)+UR(4,I)*F(I) 8183-8183-(1.0+PR)*UC(5,I)§SAV(I)+UR(5,I)&F(1) DIS1-DIS1/(8.0*PI1 DIS2¢DIS2I(8.0*PI1 8181-8181/18.0§PI) 8182-8182/(8.0*PI) SIG3-SIG3/18.0*P1) WRITE(10,37)XF(1),YF(1),DISI,DIS2,8181,8182,SIG3 FDRMAT(7(3X,F10.5),/) IF(IA.EQ.1) GO TO 30 RETURN END 2. 3. 10. 11. 12. 13. 14. 15. 151 BIBLIOGRAPHY Reddy, J.M., An Introductionggo the Finite Element Method, McGraw Hill Book Co., New York, (1984). Huebner, K.H., Thornton, E.A., The Finite Element Method for Engineers, John Wiley, Sons, New York, (1982). Zienkiewics, O.C., The Finite Element Method, 3rd edition, McGraw Hill Book Co., New York, (1977). Banerjee, P.K., Boundary Element Methods in Ehgineering Science, McGraw Hill Book Co., London, (1981). Brebbia, C.A. and Walker, 8., Boundary Element Techniques in Engineering, Newnes-Butterworths, London, (1980). Beer, 6., Finite Element, Boundary Element and Coupled Analysis of Unbounded Problems in Elastostatics, International Journal for Numerical Methods in Engineerin , 12, (1983). Mustoe, G.G.W., Volait, F. and Zienkiewics, O.C., A Symmetric Direct Boundary Integral Equation Method for Two Dimensional Elastostatics, Res Mechanica, g, (1982). Zienkiewics, O.C., Kelly, D.W. and Bettess, P., The Coupling of the Finite Element Method and Boundary Solution Procedures, International Journal for Numerical Methods in Engineering , _l_1_, (1977). Beer, 6., BEFE - a combined boundary element finite element computer program, Adv. Eng. Software, John Wiley, Sons, New York, (1982). Lapidus, L. and Pinder, 6., Numerical Solution of Partial Differential Equations in Science and Engineering, John Wiley, Sons, New York, (1982). Gupta, A.K., Efficient Numerical Integration of element stiffness matrices, International Journal for Numerical Methods in Engineering, 12, (1983). Johnson, L.W. and Riess, R.D., Numerical Analysis, Addison-Wesley Publishing Company, Reading, Massachusetts, (1982). Timoshenko, S.P. and Goodier, J.N., Theory of Elasticity, McGraw Hill Book Co., New York, (1953). Savin, G.N., Stress Concentration Around Holes, Pergamon Press, New York, (1961). Peterson, R.E., Stress Concentration Design Factors, John Wiley. Sons, New York, (1953). ‘ ..l.‘ 16. 17. 18. 19. 152 Crotty, J.M., A Block Equation Solver for large Unsymmetric Matrices arising in the Boundary Integral Equation Method, International Journal for Numerical Methods in Engineering, 1_8, (1982). Moscardini, A.O., Lewis, B.A. and Cross, M., AGTHOM—Automatic Generation of Triangular and Higher Order Meshes, International Journal for Numerical Methods in Engineering, 12, (1983). Nguyen-Van-Phai, Automatic. Mesh. Generator' with the Thetrahedrum Elements, International Journal for Numerical Methods in Engineering, 18, (1982). Perucchio, R., Ingraffea, A.R. and Abel, J.F., Interactive Computer Graphic Preprocessing for three-dimentional Finite Element Analysis, International Journal for Numerical Methods in Engineering, 18, (1982). "'TITI'ITmlfltLflflfiflflfilflfifllflffifljflflflfliflfliflfl'ES