Q) \ “5 LIBRARY Michigan State University This is to certify that the thesis entitled A PENALTY-BASED INTERFACE TECHNOLOGY USING COLLOCATION FOR CONNECTING INDEPENDENTLY MODELED SUBSTRUCTURES presented by Ramcharan Kumar Dhondi has been accepted towards fulfillment of the requirements for the MS degree in Mechanical Engineeriqu flew Major Professor’s Signature [fl/fl? /03 Date MSU is an Atfimrative Action/Equal Opportunity Institution PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 6/01 cJCIRC/DaIeDuepeS-p. 15 A PENALTY-BASED INTERFACE TECHNOLOGY USING COLLOCATION FOR CONNECTING INDEPENDENTLY MODELED SUB STRUCTURES By Ramcharan Kumar Dhondi A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Mechanical Engineering 2003 ABSTRACT .A PEN ALTY-BASED INTERFACE TECHNOLOGY USING COLLOCATION FOR CONNECTING INDEPENDENTLY MODELED SUBSTRUCTURES By Ramcharan Kumar Dhondi A computationally efficient and. effective method has been developed to connect independently modeled substructures whose finite element nodes do not coincide along the interface. This method exploits the computational ease of the collocation method to provide an efficient algorithm. The constraint formulation is based on the penalty parameter method. The meshes on either side can be non-uniform with respect to element size. The present algorithm can easily be programmed in any commercial finite element software that allows for user—defined elements. This method has been implemented as a user element in Abaqus Standard v5.8/6.2/6.3 and numerous simulations have been performed to establish its effectiveness. The most important consideration for obtaining good results using the collocation method is the location of collocation points. An effort was made to provide an easy method for selection of collocation points. After a thorough search of available literature and performing hand calculations, it was concluded that all the nodes along the interface should be used as collocation points to obtain a reliably accurate solution. The elements along the interface are constrained to the interface element, a cubic spline having equally spaced nodes, which may or may not coincide with nodes on either side. The results obtained are compared with those obtained without an interface element and with those using a previously formulated interface element. To my parents Mr. Subhash Dhondi and Mrs. Krishna Veni Subhash and my brothers. iii ACKNOWLEDGMENTS I would like to thank God for bringing me in this beautiful world and providing me with everything required to achieve my goals. I would always remain grateful to my ever-helpful parents, who have constantly helped, encouraged and blessed me with impeccable selflessness. I would like to thank my brothers for their support without which it would have been a difficult task to reach where I am now. I would like to thank my advisor Dr. Ronald C. Averill, a person with immense compassion and understanding. His unconditional support academically and otherwise really helped me get through difficult times. I would like to thank Dr. Antonio Pantano for all his help. Without his help the implementation of the collocation based interface element would have been a nigtmare for me. I would also like to thank my colleague Mr. Narasimhan for his patience. I would like to thank Dr. Dashin Liu and Dr. Farhang Pourboghrat for taking out some time from their busy schedule to be in my thesis committee. Dhondi Ramcharan Kumar May 2003 East Lansing, Michigan iv TABLE OF CONTENTS LIST OF TABLES .................................................................. vii LIST OF FIGURES ................................................................ viii CHAPTER 1 INTRODUCTION .............................................. 1 1.1 Preliminary Information .............................................. 1 1.2 Literature Review ...................................................... 3 1.2.1 Penalty Method ................................................ 3 1.2.2 Least Squares Method ......................................... 4 1.2.3 Collocation Method ............................................ 4 1.2.4 Interface Element Technology ............................... 5 1.2.5 Multipoint Constraint Method ................................ 6 1.2.6 Hybrid Interface Method ...................................... 6 1.2.7 Penalty hybrid Interface Method ............................. 7 1.3 Objective of Present Study ........................................... 7 1.4 Organization of the Thesis ............................................ 8 CHAPTER 2 FORMULATION OF INTERFACE ELEMENT ........ 10 2.1 Hybrid Interface Method ............................................ 11 2.2 Penalty Hybrid Interface Method .................................. 14 2.3 Penalty Hybrid Interface Element using Collocation ............ 15 2.4 Determination of Penalty Parameters .............................. 17 2.5 Automatic Round-Off Error Control .............................. 18 CHAPTER 3 PENALTY HYBRID INTERFACE ELEMENT USING 3.1 COLLOCATION ............................................. 21 Axial Load-Two Bar Elements ..................................... 21 3.1.1 Lagrange Multiplier Method ................................ 22 3.1.2 Penalty Method ................................................ 24 3.1.3 Collocation Based Interface Element ...................... 26 3.1.4 Comparison of the Three Methods ......................... 27 3.1.5 Relation Between Penalty Parameter and Beam Properties ...................................................... 27 3.2 Two-Dimensional Formulation .................................... 28 3.2.1 Cubic Spline Interpolation Functions ...................... 28 3.2.2 A General form of Cubic Spline ............................ 29 3.2.3 Two-Dimensional Plane Stress Problem .................. 33 CHAPTER4 RELIABILITY OF THE NEW INTERFACE TECHNOLOGY .............................................. 43 4.1 Two-Dimensional problems ........................................ 43 4.2 Error Analysis ......................................................... 49 CHAPTER 5 CONCLUSIONS .............................................. 57 BIBLIOGRAPHY ................................................................... 58 vi LIST OF TABLES Table 1. Summary of Uniform Axial Load for Uelgluenmnpsf ...................... 53 Table 2. Summary of Uniform Axial Load for Uelgluenmf ........................... 53 vii Figure 1a. Figure 1b. Figure 2. Figure 3. Figure 4a. Figure 4b. Figure 5a. Figure 5b. Figure 6a. Figure 6b. Figure 7a. Figure 7b. Figure 8a. Figure 8b. Figure 9. Figure 10a. LIST OF FIGURES Example of a one-dimensional interface element in a two dimensional problem ............................................................................. 10 Example of a two-dimensional interface element in a three-dimensional problem ............................................................................. 11 Beam under an axial load applied at the tip-Two Element .................. 22 Mesh and Configuration for a two-dimensional plane stress problem. ....34 Horizontal displacement using penalty hybrid interface element ........... 37 Vertical displacement using penalty hybrid interface element .............. 38 Horizontal displacement using user element uelgluenps_wo.f (NPS being the collocation points) ............................................................ 38 Vertical displacement using user element uelgluenps_wo.f (NPS being the collocation points) ................................................................. 39 Horizontal displacement using user element uelgluenm_wo.f (N and M being the collocation points) ..................................................... 39 Vertical displacement using user element uelgluenm_wo.f (N and M being the collocation points) ............................................................ 40 Horizontal displacement using user element uelgluenmnps.f (N, M and NPS being the collocation points) ............................................... 40 Vertical displacement using user element uelgluenmnps.f (N, M and NPS being the collocation points) ..................................................... 41 Horizontal displacement using Abaqus 6.3 .................................... 41 Vertical displacement using Abaqus 6.3 ...................................... 42 Two-dimensional uniform axial loading problem with an inclined interface (63.43 deg) .............................................................. 43 Horizontal displacement using user element uelgluenmnps.f ............... 43 viii Figure 10b. Figure 1 1a. Figure 1 1b. Figure 12. Figure 13a. Figure 13b. Figure 14a. Figure 14b. Figure 15a. Figure 15b. Figure 16 Figure l7a-e Figure 18a-e Vertical displacement using user element uelgluenmnps.f .................. 43 Horizontal displacement using penalty hybrid interface element ........... 44 Vertical displacement using penalty hybrid interface element .............. 44 Two-dimensional bending problem with an inclined interface (63.43 deg) ................................................................................. 45 Horizontal displacement using Abaqus 6.3 .................................... 45 Vertical displacement using Abaqus 6.3 ....................................... 45 Horizontal displacement using user element subroutine uelgluenmnps.f ..................................................................... 47 Vertical displacement using user element subroutine uelgluenmnps.f ..................................................................... 47 Horizontal displacement using penalty hybrid interface method ........... 46 Vertical displacement using penalty hybrid interface method ............... 46 Nomenclature for Legends in the following plots ............................ 48 Deformed shape of the interface for Uniform axial load ................. 50-52 Deformed shape of the interface for bending load ........................ 53-55 ix CHAPTER 1 INTRODUCTION 1.1 Preliminary Information Technology has evolved continuously in every part of the world. Exchange of technology through information passage has been a primary mode of development in the field of engineering. Proper usage of this technology with consistency in every respect has been a difficult task. Methods have been developed over the years to allow effective, reliable and consistent data exchange. A few methods have been effective in specific cases and others have demonstrated their robustness in terms of application. Increasing exchange of information over the intemet in design and analysis of mechanical systems has lead to innovative technologies in the past decade. Easy availability of high performance and stable computing syStems has resulted in effective analysis of large-scale mechanical systems. It has become a common task to perform analysis of large-scale models of vehicles, spacecraft and other structural assemblies using available finite element software. But these models are still developed using substructures, which may be created by different engineers using different software and in different geographical locations with little or no communication between the teams of engineers creating the models. The time required to analyze these systems has reduced drastically but the amount of time used to set up these models has been still a drawback. To avoid this analyst have looked for ways to independently model substructures and finally assemble them. This still involves a lot of tedious and intelligent human work. Usually there are cases where these models are not compatible at their interface, i.e. the nodes at the interface do not match. This calls for effective modeling and assembly during analysis. Methods have been developed over a decade to resolve this issue. A great amount of success has been achieved but the unending increase in complexities of these structures has prompted for computationally effective, efficient and robust methods. One such advancement was development of the hybrid interface technology [1-3]. Incompatible meshes generated by analysts can occur with global/local analysis, where part of the structure is modeled as the area of primary interest, in which detailed stress distributions are required, and part of the structure is modeled as the area of secondary interest. Incompatibility occurs, as the finite element mesh in the primary area tends to be finer than that in the secondary area. Transitioning sometimes leads to distorted meshes, which may significantly affect the quality of the results. Interface elements can be used in such cases to achieve pr0per load transfer across the interface. Such incompatibility usually occurs in cases of sub assemblies modeled using automeshers, which are generously used in highly complicated structural assemblies. Many automeshers generate tetrahedral meshes for solids, and distorted tetrahedra may be more susceptible to poor results. One more aspect of structural analysis is optimization and this involves multiple analysis of the same structure using different meshes. A mesh of different size and shape may be created due to change in shape of the structure. The change could be only in a specific region without any noticeable change in the rest of the structure as is the case with many structures. In such cases the structure can be effectively modeled using subdomains. This warrants for an effective and easy to use technology to assemble the substructures. 1.2 Literature Review Analysts and researchers have worked over past decade to develop a reliable, effective and robust interface technology to solve the problems mentioned above [Sect 1.1]. The new interface technology, presented through this thesis is another such development, which has been developed using various methods. These methods have been studied extensively in past. Discussions on developments in each of methods and the interface technology have been presented below. 1.2.1 Penalty Method The penalty method for enforcing constraints has been studied extensively in the field of finite elements and mathematics. Penalty parameters have been used to avoid Dirichlet boundary conditions and hence modifying the functional to be minimized. The modified functional now consists of the original functional and a penalty term. When the modified functional is minimized we see that the boundary condition gets satisfied in a least squared sense [14]. Babuska was the first to use the penalty-function method in conjunction with the Finite Element method. He proved the existence and uniqueness of the finite element solution to the penalty—function formulation of the Dirichlet problem for the Poisson’s equation [15]. Element formulations based on penalty parameters have been studied. One such example is the case of shear deformable plate elements. These elements can be considered as elements derived from classical plate theory by treating the zero transverse shear strains as constraints. These constraints are formulated using penalty method and reduced integration is applied, as the penalty terms obtained are shear energy terms. This has been demonstrated by Reddy [22]. The penalty method has also been applied to formulation of contact algorithms by many researchers. Belytschko and Neal have applied penalty method to pinball algorithm for contact [23]. 1.2.2 Least Squares Method The Least Squares method has been applied extensively in engineering fields like optimization and finite elements. The fundamental advantage of Least Squares method is, the error (residual in finite element) is squared and thus the functional when minimized approaches the minimum value from only one side. This is advantageous as both positive and negative errors are squared and thus an appropriate measure of error is obtained. The least squares method provides numerous theoretical and computational advantages in the algorithmic design and implementation of finite element methods that are not present in standard Galerkin discretizations, when applied to elliptic boundary value problems in linear elasticity, fluid flows and convection-diffusion. Most notably it leads to symmetric and positive definite algebraic problems [36]. Cao has applied Least squares method to interface problems [37]. 1.2.3 Collocation Method The collocation method is well known for its ease of implementation and has been utilized mainly to gain computational efficiency, although the choice of collocation points drastically affects the solution of a problem. Ioakimidis and Theocaris [38] showed that the choice of collocation points such that sum of all additional terms in the numerical approximation of a singular integral equation vanish, gives faster convergence and a good solution. Ascher [39] while working on stiff boundary value problems proved that the best Runge-Kutta schemes known, which are symmetric and algebraically stable are those equivalent to collocation at Gaussian points. He also showed that collocation at gaussian points turns out to be more advantageous than any other choices because it is computationally cheap and has better convergence. In 1982 Redekop [40] introduced Fundamental Collocation method for solving planar elasto-static problems. Here the boundary conditions were satisfied using boundary point least squares collocation technique. De Boor and Swartz have showed that choice of Gaussian points as collocation points gives better convergence as compared to Runge Kutta method and same order of accuracy [24]. 1.2.4 Interface Element Technology A lot of work was done to provide a new effective method to resolve the issue of assembly of incompatible meshes. Work was done by Farhat [45], which involved a domain decomposition approach for partitioning the spatial domain into a set of disconnected subdomains, each assigned to an individual processor. Lagrange multipliers have been introduced to enforce compatibility at the interface. Ransom [13] worked on a hybrid interface element, which is based on a constraint formulation through Lagrange multipliers. Pantano [16-17] tried to improve on the computational aspect of the interface element by using a penalty parameter based constraint formulation along the interface. The hybrid interface element in comparision to penalty hybrid interface element is more challenging to implement in a commercial code. 1.2.5 Multipoint Constraint Method Efforts made initially to resolve the problem of mesh interfacing were mainly concentrated on moving nodes or writing multi-point constraint (MPC) equations on the interfaces [4-6]. The biggest restriction of moving nodes is that both sides should have the same number and type of elements. Therefore this method is not practical for general application. Different types of elements along an interface can be connected using multi- point constraint equations. Splines used as multi-point constraint equations handle more general cases. However, MPC equations by definition provide additional relationships for the existing degrees of freedom on the interface and in this process the number of independent degrees of freedom are reduced. If no new degrees of freedom are created this could result in additional local stiffness or other non-physical effects in the model [4- 6]. 1.2.6 Hybrid Interface Method This element developed at NASA LaRC [1-7] allows connection of independently modeled substructures with incompatible discretization along the common boundary. This method has developed into a very effective method as Lagrange multipliers apply constraints accurately. However this method results in a non-positive definite system of equations. Although robust solvers exist for solving non-positive definite system of equations, they typically require row and/or column pivoting which adversely impacts the solver efficiency. It also has input requirement and implementation issuses, which present a challenge for commercial codes. Aminpour, Pageau and Shin [12] have presented an alternate method based on hybrid interface element, which takes care of the implementation issues. The variational constraint equations obtained are not the most stable and numerical ill conditioning of the resulting system of equations may not be avoidable. 1.2.7 Penalty Hybrid Interface Method The penalty-based interface element developed at Michigan State University by Pantano [16-17] has some advantages over the lagrange-multiplier based hybrid interface method. In this case the variational formulation of the interface element enforces interface constraints via the penalty method. One of the advantages of such a variational formulation is, that a positive-definite and banded stiffness matrix with reduced DOFs is obtained. This approach enhances the computational efficiency as state of art sparse solvers can be used to solve the system of equations. From accuracy point of view the penalty method enforces the constraints only approximately, depending on the values of the penalty parameters chosen. This method has been effectively used to solve numerous problems. 1.3 Objective of Present Study The present study focuses on the effective use of the collocation method instead of piecewise integration along the interface for enforcing the constraint between two incompatible meshes. The constraint is still enforced approximately using the penalty method as in the penalty hybrid interface element [16-17]. The primary difficulty is to determine the collocation point locations. As integration along the interface is avoided by using collocation, we observe that the amount of computation is reduced considerably, at the same time all the desirable characteristics of the penalty hybrid interface method are conserved. Maintaining all those desirable characteristics depends on the choice of collocation point locations. Collocation method is numerically more efficient, especially in explicit dynamic analysis. Collocation method is a more general method, which can be applied to finite element models, as well as non-finite element models (e.g. FEM-FDM). As no integration is required, the approximation functions are irrelevant to the procedure. Thus it works well for any order of element. The accuracy and stability of the result depends on the choice of collocation points. Three different options were investigated. The most reliable and effective choice is presented in detail with a discussion on why this choice is reliable and effective. 1.4 Organization of the Thesis In Chapter 2 early developments in interface element technology are presented along with a detailed formulation of a collocation-based interface element. Chapter 3 provides formulation of one-dimensional collocation based interface element for two-dimensional problems. A discussion on the selection of collocation stations is presented. In Chapter 4 a comparison is made between the different collocation schemes and Pantano’s interface element. A variety of problems have been solved and the results for specific interesting cases are presented in this chapter to demonstrate the new interface element. Chapter 5 discusses conclusions and opportunities for further study and research in this interface technology and related areas. CHAPTER 2 FORMULATION OF INTERFACE ELEMENT Let us consider two substructures that represent either a 2D or a 3D assembly, having substructures £21 and (22 as shown in Figures 1(a) and 1(b). At the interface between the two substructures we employ an interface element, which binds these substructures like glue at the common interface. The interface element is discretized with a set of nodes, which do not depend on the number of nodes at the interface in the subdomains Q] and £22. In the collocation based interface element discussed in this chapter the coupling terms associated with the interface element are arranged in the form of a stiffness matrix and then assembled along with the global stiffness matrix. Q] qs Q; ul 0V I12 V Figure 1(a) Example of a one-dimensional interface element in a two-dimensional problem. 10 \ -p‘¢..‘,w- u. .. - N. N. __..-v..-.g.w -.- . I Figure 1(b) Example of a two-dimensional interface element in a three-dimensional problem. 2.1 Hybrid Interface Method In the process of development of hybrid interface element Aminpour, Ransom and McCleary developed two other formulations [7] based on collocation method and discrete least squares method. Both these methods failed the patch test. The failure of collocation and discrete least squares formulations stems from the fact that the load transfer mechanisms of these formulations are not, in general, consistent with the load transfer mechanism of the finite elements. The original hybrid interface method [1-9] uses lagrange multipliers to apply the constraints. The variational formulation transfers the load using the same load transfer mechanism as the finite elements, used to model the connected subdomains. The 11 constraints are exactly satisfied in this case. The lagrange multiplier DOFs are variables found by solving the system of equations. The system of equations is obtained from the first variation of total potential energy (TPE) with respect to all the variables taken one at a time (DOFs of all the nodes and lagrange multiplier DOFs). As the DOFs include Lagrange multiplier DOFs the implementation in commercial codes is challenging. The system of equations is non-positive definite. Thus state of art matrix solvers for positive definite system of equations cannot be used. Efficient solvers do exist for solving non- positive definite system of equations, but they typically involve row and/or column pivoting which adversely impacts the solver efficiency. The variational formulation involves minimization of the total potential energy (TPE) of the assembly. The form of TPE is as follows: ”=49, +7502+I41(V-u1)dS+I4Q(V—“2)d5 (1) 5 S The nodal displacements of the sub-domain Q,- are identified by q? and qj. The DOFs with o as superscript are not on the interface and ones with idenote DOFs that are on the interface. The displacements are expressed in terms of unknown nodal displacements q; as: u j = qu'j (2) where N,- can be matrices containing linear lagrange interpolation functions. The displacement field V of the interface element is approximated on the entire interface surface in terms of unknown nodal displacements q, as: V = qu (3) 12 where T is a matrix of cubic spline interpolation functions. The Lagrange multipliers are assumed to be of the form xi = Ra. The interpolating functions in the matrix R are constants for linear finite elements and linear functions for quadratic finite elements. As the Lagrange multipliers also enter the DOFs, the first variation of Iris taken with respect to all the DOFs including the vectors a1 and a2: 57! . c = O Iqi’ ,qiasai 45,01,052 (4) These equations enforce displacement continuity, with 2.,- representing the interface tractions. The total traction is zero on the interface. The first variation flyields the system of equations: "K,“ K,” 0 K 1‘” {’0 0 o 0 Kg 0 O K g" 0 O 0 MT 0 0 _ 0 0 M; _ T _ r ._ where M}. ——S[_N1Rjdsj and G}. —S[T Rjdsj. ,j—1,2. I O O K; K3” 0 O O 1 GOOD 0 G,T GI M. 0 - lql" ‘ft‘ 0 0 41" f1" 0 M2 q; 2“ 0 0 1 ”=7Z'1+7[2+—2'}’lz 1 n: ’71:] 15 where n1, p and m1 are number of collocation stations for each of the interfaces. The number of nodes on left hand side is n1, that on the right hand side is m1 and p is number of nodes on the spline. The form V|n denotes that the displacement field is evaluated at collocation station n. The displacement fields are defined in the same manner as in Penalty Hybrid Interface Element. The first variation of 7: yields required equations, which form a positive definite set after applying boundary conditions. Defining coupling terms as before, we have: if P T Gi(I’K)=7j [EN] LNKII +1 (16) Gi's(”K)=7j(NI|K +TKI1) (17) 3s nl-I-ml T Gr (”0:71 El TKIITII,+2 (19) where: j = 1, 2 N, I is the linear lagrange interpolation function of the 1“ node on the interface of j’h substructure evaluated at the spline point l. TKI is the cubic spline interpolation function for the Km node on the interface element 1 evaluated at 1” node on the interface of the j’h substructure. The stiffness matrices are assembled in the same manner, as done in penalty hybrid interface element. 16 2.4 Determination of the Penalty Parameters A method for automatically determining the value of the penalty parameters can be found in references [16-17], but for the sake of completion it is presented here in condensed form. The penalty parameters are a set of predetermined constants. The penalty method imposes constraints in an approximate form, and the accuracy of the solution is dependent on the value of the penalty parameter. This penalty parameter is dependent on the material properties and geometry of the model. Further, there is a relationship between the penalty parameter and the Lagrange multiplier that enforces a given constraint. As the lagrange multiplier method imposes the constraint exactly, it defines the upper limit to the accuracy of the penalty method. The value of penalty parameter in terms of the material and geometrical properties is found using the correct solution for the problem. A different penalty parameter should be defined for each of the different degrees of freedom in case of elements with multiple-degrees of freedom. The penalty parameter could be found by considering a fairly simple model with only one kind of loading so that penalty parameters for each of the degrees of freedom can be obtained one at a time. The formulations and solution are obtained analytically using both the Lagrange multiplier method and the penalty method. The displacement solutions for each degree of freedom are compared and the ratio between them is expressed in the form as shown: la ’run 1: u A A penalty = 1 + g (20) where f = f (element geometric properties, material properties, and loads) Once this expression is identified, the penalty parameter 7 is set equal to: 17 r=flf (2D Then the ratio between the solutions becomes independent of material and geometrical properties of the element: penalty 1 “lagrange = 1+ 5 (22) Thus we see that the accuracy of solution directly depends on the value assigned to the parameter ,6. But the accuracy of the solution cannot be indefinitely increased, since round off error would rise as ,6 is increased. But once a compromise between constraint representation error and the round off error has been reached a value of ,6 can be found which is able to produce the same level of accuracy for every combination of material and geometrical properties. More complicated cases where the value of ,6 is dependent on more than one DOF are considered in [16-17]. However it has been found that the value of ,6 in case of collocation method needs to be 1 order higher than that in case of penalty hybrid element method [16-17]. 2.5 Automatic Round-Off Error Control In order for the constraints to be properly imposed we need to satisfy certain inherent properties of the interface element stiffness matrix. One such property is that the sum of all the terms in a row or column be equal to zero. This condition cannot be usually achieved because of round-off error, and the resulting inaccuracy grows with value of the penalty parameter. An important measure thus devised is to look at row imbalances in the 18 stiffness matrix. If K. is the diagonal stiffness associated to the n-th nodal DOF, it is sufficient to consider the ratio: Q. = " (23) where ER" is the unbalance in the interface element stiffness matrix row related to the DOF n. ER" = 2 K", (24) 1' When the value of Q” exceeds about 1E-7, errors in the solution may become appreciable. The discussed row imbalance is proportional to the value or the penalty parameter, ER" cc )6 It is also approximately true that: ER.°<7°<fl-K.=>Q.°<fl (25) Thus, an algorithm was developed in [16-17] to control round-off error. Its steps can be summarized as follows: 0 Stiffness terms are computed in the interface element stiffness matrix. 0 For each row in the stiffness matrix 0 The highest term is selected and stored in a variable X o The row imbalance of the stiffness matrix is stored in a variable ER 0 Q = £1?- is evaluated 0 The highest Q is found compared to a given constant C. Typically C = 1E-7 is used. IIt'W o If Q > C, the parameter 6 is reduced according to: 6 = 6 mln 19 0 The interface stiffness matrix is thus calculated using the 6 = 6 "‘W. Thus we see that this kind of approach reduces the risks of round-off error, which could affect the solution. 20 CHAPTER 3 PENALTY HYBRID INTERFACE ELEMENT USING COLLOCATION In this chapter hybrid interface element, penalty hybrid interface element and penalty hybrid interface element using collocation are compared for one-dimensional problems. A formulation of two-dimensional problems involving one-dimensional interface element is presented. Multiple problems using the user element implemented in Abaqus 5.8/6.2/6.3 are presented to show the reliability of the collocation based interface technology. In the latter part of the chapter an error analysis is performed on a set of problems to assess the accuracy of the method in finite element simulations. The first objective is to obtain a stiffness matrix using the general formulation. Once the stiffness matrix for the interface element is obtained, it is assembled along with other element stiffness matrices to form the global stiffness matrix. Then the usual procedure is followed which involves applying boundary conditions and inverting the global stiffness matrix to obtain the unknown displacements. 3.1 Axial load — Two bar elements A very suitable problem to show the differences between the lagrange multiplier method, penalty method and the collocation based interface element is presented. In this problem we have axial load applied to a bar, which is modeled as a two-element mesh with an interface element between these two elements. The problem is depicted in Figure 2. 21 1 2 v 3 4 e c o % : u1 E, A "2 V "3 E, A “4 Figure 2. Beam under an axial load applied at the tip — Two Element 3.1.1 Lagrange Multiplier Method The hybrid interface method introduces lagrange multipliers xi, and 112 in the total potential energy (TPE) of the system to satisfy the displacement continuity condition. Thus the TPE takes the form as shown: 7r=7r1+zr2 +III°(V-l12)+/IQ -(v—u3) where 7n, 7r; are, respectively, the TPE of the first and of the second bar. It: [1(ox-ex)dV+ j —1—(ox-ex)dV-Pu4+/il-(v—u2)+,i/,'(v—u3) V12 V22 1 2 21 2 -(I)EA(%‘) dx+%- {EA(%) dx—Pu4 +1, -(v—u2)+A/Zo(v-u3) Displacements are approximated linearly as Substituting the approximated displacements into TPE gives the following result. 22 (26) (27) (28) (29) ”_EA ’- dN dN, _E_A 21- dN de ”-7 (Ilia “‘ dx —I§"17x_]dx+ 2 {(Zu "' dx _I§urd_x_}x_P'u4+ (30) A-(v—ufl-I-Z/Z-(v—ufl Setting the first variation with respect to all the DOFs to zero we get: L dN- 9—” =EA-j[iN—1 Eu, 1 x=0 (31) all] 0 (ix 1' dx d_7_r__ L dN2 de Buzz {,l dx I?” dx ’1‘ ( ) 87: L le de =EA- — - — =0 33 Big ({[de§“} (3%) 12 ( ) L dN- 31:5A.,(__dN2 z“, I —P=0 (34) 8114 0 dx 1' dx a7£=xll+xilx2==0 (35) av 87: 5/1—1=v—u2=0 (36) air Ezv—u3=0 (37) After applying the boundary conditions (in: 0) these equations can be written in the same form as the classical equation involving stiffness, displacement and force. 23 95. __EA 0 0 0 0 L L -55. 5:4. 0 0 0 -1 L L o 0 E __E_A. o 0 L L 0 -54 12/: 0 0 L L 0 O O O 1 O —1 0 O l O _ O 0 —1 O l 0 Solving this system of equations we have: n=—P xil—P 12-40 u —u -s—52 2 3 EA 2PL u4—— EA 3.1.2 Penalty Method (38) (39) (40) (41) (42) (43) In the penalty method the displacement continuity constraint is imposed through two penalty parameters y. and y2. Thus the TPE of the system shown is Figure 2 changes its form. 1 l 7Z'=7'[1+7Z'2 +§yl -(v—u2)2 +EY2'(V"u3)2 where 7n and 7:2 are respectively, the TPE of the first and of the second bar. (44) 1 l l 1 7t: {—2-(0'x-8x)dV+ I 5(O'x-8x)dV—Pu4+§yl-(v—u2)2+§y2-(v—u3)2 (45) V1 V2 24 2 L L d 2 12L d 1 1 -[EA(Z:-)dn — [EA(d:]2 dx—Pu4+§y, ~(v—u2)2+-2—72-(v-u3)2 0 The displacement is approximated using linear langrange interpolation functions. 2 = Z uJN i=1 all—sluts) Substituting the approximation for displacements we have: ét’I'W—LQ) +—2-rz- EA L dN EA LLL dN, --2- ,(s, .1; 152-9) +— :12: Tie dx N . d J x—P-u4+ dx ] The first variation of It with respect to all the DOF assumes the following forms: L 35 EA- {VI—NZ. 8114 0 (ix 872' —=71‘(V—“2)+72°(V—“3)=0 av Thus the FE model is obtained as: 25 (46) (47) (48) (49) (50) (51) (52) (53) The solution to this system of equations is as follows: fi='P PL 112 - -— EA 3.1.3 Collocation based Interface Element E -123: 0 0 0 L L (0 __E_[4_ E+y -7 O O L L l l 112 0 -71 71+72 -72 0 W A 0 0 -72 _E_A_+,2 -15.. u;; L L [u4 0 0 0 -54; .55; ‘ L L L - (54) (55) (56) (57) (58) (59) From section 2.3 it is known that in this technique the displacement continuity constraint is evaluated at each of the nodes along the interface. Thus we see that the total potential energy for the problem at hand is as follows: 1 l 71:72] +7[2 +371 -(v—u2)2 +—2-}’2 ~(v—u3)2 (60) It is observed that the TPE for collocation based interface element and penalty hybrid interface element are same. Thus the formulation of stiffness matrix, which 26 follows similar procedures as in penalty based interface element, is the same as in Eq. (54). 3.1.4 Comparison of the three methods It is seen that lagrange multiplier method enforces the constraint accurately and thus is 3 upper limit for the latter two methods in terms of accuracy. The convergence for the last two methods as the penalty parameters are increased is presented. . . . . 1 L L Lzmztv= lelt —+— = — (61) 71—>°° 71->°° 7’1 EA EA L l 1 L Limit “3 = Limit [—+—+— = (_)P (62) 71.72 ->°° 71,72 ->°° EA 71 72 EA 2L 1 1 2L Limit “4 = Limit [—+—+— = (—)P (63) 71.72 —><>o r1124... EA 7’1 7’2 EA 3.1.5 Relation between Penalty Parameter and Beam Properties It is seen that the result from lagrange multiplier method is same as the theoretical result. Both penalty method and collocation based interface element converge to the theoretical value based on the penalty parameter value. It is very important as mentioned in sections 2.5-2.6 that an appropriate method to find the value of penalty parameter be developed. This value of penalty parameter can be obtained by observing the results obtained through theoretical calculation and penalty method. A method to obtain this penalty parameter is presented below [16]. 27 “exact 2PL = 1 + 2 L (64) EA EA Thus by substituting EA 71:72 ='BO(—L—l (65) for the penalty parameters we obtain “(penalty 1 = 1 + — (66) “Exact fl 3.2 Two-Dimensional Formulation In this section formulation for a one dimensional interface element for two- dimensional problems is presented with the help of an example. The interface element is formed using cubic spline interpolation functions as approximation functions for the displacement field. 3.2.1 Cubic Spline Interpolation Functions Spline functions are mathematical tools to get a better understanding of available data by passing a smooth curve through the data points. Splines exist in various orders but cubic splines are most widely used in engineering practice. It is important to note that cubic splines do not have the “wiggle” problem associated with higher order interpolating polynomials. Cubic spline functions can connect as many data points as possible. A general form of cubic spline as used in [13,17] is used. In this general form of cubic spline 28 interpolation functions the second derivatives at the end points are calculated by differentiating a cubic function passing through the first or last four points depending on the end point at which the second derivative is desired. This general form of cubic spline interpolation function is found to be computationally efficient. 3.2.2 A general form of Cubic Spline A description of the general form of cubic spline is presented for completion. A very efficient matrix method is followed which makes it easy to program. Given a series of points x,- (i=0,l,...,n) which are generally not evenly spaced, and the corresponding function values f(x,-), the cubic spline function denoted by g(x) may be written as: g(x) : g,xx(xi) (x,-+1-x)3 6 Ax,- f (x,- )[-(£‘+Xlx:—x)-] + f (xi+1)[(xA—x:’ )] -Axi(x—xi):|+ . _ .3 _Axi(xi+l -x)] + g,xx(x,+1) [0‘ x1) 6 Ax,- (67) where Ax) = x,- - x,- and g,” denotes double differentiation with respect to x. This equation provides the interpolating cubics over each interval for i=0,l,.....,n and may be given in matrix form as: g = ri‘lgxx + T2f (68) For each of the k values of x at which the spline function is to be evaluated, x, S x, S x. 1+1 ’ k=1,2,....,p, and p is the number of evaluation points. The T, and T2 matrices can be written in the form: 29 T2: where (L1)1,1 (6)1,2 ~ _ (502.1 (102.2 [(61 )p,l (fl )p,2 (L2)1,1 (52)1,2 ~ (f2)2.1 (’2)2,2 (2),, (mpg A 7 (L1)1,n+1 (f1)2.n+1 (69) (t1)p,n+l_ (12)1,n+1 (f2)2,n+l (70) (12)p,n+1_ for xk < x,- 01' xk >xi+1 O .. l (x,-+1—xk)3 ._. (11);“: = < - —Ax,-(x,-+1—xk) for x, S xk Sxi+1 andJ—1+1 (71) 6 Ax,- _l_ (xk —x,-)3 —Ax°(x —x-) for x, Sxk Sx,+1 andj=i+2 L 6 Ax, r k 1 and [3,xx(x0)‘ i8, (x1) gxx = ”i I (72) 84,0"), rf(xo)‘ f z, f(:x1) , (73) If(xn), Note that there are, at most, two non-zero coefficients in each row of the T1 and T2 matrices in equations (69) and (70) respectively. Applying additional smoothness 3O conditions (i.e., equating the first and second derivatives of adjacent interpolation cubics at x,) yields a set of simultaneous equations of the form: Ax-_ 2 i _ i |: Air—i1 ]g,xx(xi-l) + [L2%l]g,xx(x,) + [l]g.xx(xi+1) = i=1,2,...,n-1 (74) 6 f(xi+1)- f(xi) _ f(xi) " f(xi—l) (M02 (AxiXAxi—l) if the x,- are evenly separated with spacing Ax, then the above equation transforms to: [llamas—1H I4Is,..+111g,..(x..s> = élf‘xHWZLOZ“ “xi-1’ (75) - (Axi) Now we can write the equations for evenly separated and unevenly separated as: Ag,xx = Pf (76) The coefficients of matrices A and P are dependent upon the end conditions. Whether the equations are of the form of that for evenly spaced data or unevenly spaced data, there are n-l equations and n+1 unknowns g,xx(xo), g,xx(x1),..., g,x,(x,.). The two necessary additional equations are obtained by specifying conditions on g,xx(xo) and g,xx(x,.). In general spline, these second derivatives are calculated by differentiating (twice) a cubic function, which passes through the first four pseudo-nodes along the interface path and another cubic function that passes through the last pseudo-nodes along the interface path. Evaluating this cubic function: ._ 2 3 g(xo) - a0 + alx + azx + a3x (77) at the first four points gives: 31 < g(xr) g(xz) lg(xo)‘ (8063), f II pup—Hp. x0 x1 x2 x3 Solving for the coefficients yields a = N'1 g or Iao :.._ "“‘lvf:°5'\n' A... ("an "._. .r_,, I... . Figure 13a. Horizontal displacement using abaqus 6.3 NKOQAI—‘lOmwOI—hhqo N (o m I o .h- Figure 13b. Vertical displacement using abaqus 6.3 46 U, 01 Horizontal displacement using user element subroutine uelgluenmnps. f 5 0555444444444 0 0000000000000 . +..___._..._. e eeeeeeeeeeeee 4 0628888776665 3 0998543210987 fl 0876147036814 1 0258111222233 + +.______..__. Figure 14a. 132 Vertical displacement using user element subroutine uelgluenmnps.f 47 Figure 14b. Eggauflumaanr... ., . -. Inlagggllgufinazmssmwu— ' ’ ”II-Illllflliigg‘fn'idmh a. “ Figure 15a. Horizontal Displacement using penalty hybrid interface method Figure 15b. Vertical displacement using penalty hybrid interface method The tip displacement for this bending problem obtained using the beam theory is PL3 3.125E-4 (UZTIP :5). From Figures 13a through 15b it can be observed that the solution obtained using the collocation based interface element matches that of the 48 penalty hybrid interface element. This result also matches the result obtained by using Abaqus 6.3 on a uniform mesh. The displacement field is observed to be continuous across the interface. To make sure that there is no discontinuity of displacement at the interface, the interface area was observed closely. 4.2 Error Analysis In this section a set of problems have been solved using various methods viz. penalty hybrid interface method, collocation methods (two of the choices) and using Abaqus v6.3. For the collocation based interface technology two choices of collocation points have been used. One of the choices for the points is, the nodes on either sides of the interface element (N and M) and the other is, the nodes on the interface element and those on either sides of the interface (N, M and NPS). Both uniform axial load and bending load have been used for error analysis. The plots below show, how the edges on either sides of the interface look in case of different analysis schemes. The legends have been named in a particular fashion to make them more intuitive. The nomenclature used is as shown in Figure 16. /® (uéFlg uenmnps) m Interface side : a - left / \ b - right Value of N. M and NPS respectively. Default: NPS = 4 if not specified Analysis Type : uelglue - Penalty hybrid method uelgluenm - Collocation method with N and M uelgluenmnps - Collocation method with N. M and NPS abaqus - Abaqus 6.3 without interface element Note: abaqus doesn't have interface side and NM and NPS Figure 16. Nomenclature for legends in the following plots 49 10.0 ----~-,- 9.0 Interface plots (N=3. M82. NPs-4) i-e- 8361613652 " " “a; 1‘ 611515111332 W 787 77.136151011016327 .i‘...;.zsu..'gtsaasz“ i ..f.i.*i:-9!.°§¥1€"£‘0_833?..:'::PI!!'9.'."!E"E'.'B§3? .i.:i‘~':;: 1.99999? - __ 9.0 7.0 8.0 5.0 ...gg 4.0 3.0 2.0 1.0 0.0 -1.o---w~~- 90.9 Y Coor'chde of the interface I: r I r g 1 90.95 90.97 90.99 91.01 cte-- - Figure 17a. Deformed shape of the interface for Uniform axial load 9.0 Interface plots (N-4, Ill-3. NPs-4) —- -- . -1‘r.- ~ - auelgluenmnps43 - -—-- buolgiuenmnpsfl , - ~21~abequs — <> - euelgluo43 — L‘Sr- buelglue43 --/‘~- uelgluenm“ -~ ->< -- buelgluenm43 9.0 7.0 8.0 5.0 I): I -— 5+? 3* ~10 )- . [30- - 1:1). D 4.0 .1," .1)... 3.0 2.0 Y Com. of the Irrtertece 1.0 0.0 : L - .- i ) é) a. 1 _,.-- | xix. 'j‘ if" I'i\ .‘7‘. ....... ,_ .... - . k . T 1 a T I 1 -1.0 -~» - - ~ - -........._.- 90.99750 Figure 17b. 90.99850 90.99950 91 .00050 91 .00150 91 .00250 91 .00350 X Coordhate otthe Interface Deformed shape of the interface for Uniform axial load 50 Interface plots (NIB, M810. NPSI4) — <>- euelglue8-10-4 - 1-3- . buelglu08-10-4 MA euelgluenm8-10-4 --/-- buelgluenm8-10-4 - -x—-- auelgluenmnpsB-I 0-4 —---- buelgluenmnps8-10-4 ~-::~--: abaqus 100 ..._. _. -_..._-.. . ._____.W_..._-__........_-_-“_.__.____- .._.-1_....--.-L.....-.._..---.-.._ 9.0 g ......................... 3,9. . 9'0 hfi ..... f ” S 7.0 ‘ :3 _‘~. 4 6.0 a “It? 5 5.0 9’9- a 3 g I g "0 é 7’ ; {35‘ 3.0 ‘ ""ié’f1'7: 99 ILA \. ' 2.0 ' *9; . i >- ‘ ‘~ . 1.0 g * 0.0 r r 1 E I r r I '''''''''''' ;-... £1 1 -1_o_, .-_-- . - .. ...... - L”- --,--...-,- . . __---_. 90.9993 90.9995 90.9997 90.9999 91.0001 91.0003 91.0005 91.0007 91.0009 91.0011 XCOordInate ofthe Menace Figure 170. Deformed shape of the interface for Uniform axial load Interface plots (N=8, M=10, NPS=8) - <>- auelglue8-10-8 .- [3— buelglue8-10-8 ,3..- auelgluenm8-10-8 ~\ buelguenm8-10-8 --‘=1r--~ auelgluenmnp58-10-8 —-—~ buelgluenmnps8-10-8 ~-«-.'~\- abaqus 10.0 — - L--. ~ , 9.0 . . .... % . ..:'.:::.:::.:' .::‘.'...;'.'_."".f.‘. illijilffj ::: ~ u .u E g 8.0 {5% 7.} 7.0 f" 1;';';:f!' .A .8 ' i 6.0 a ‘ , 1...: 5.0 is if i 4.0 fl 0 2.0 % ..... , _ > 1.0 . _ a ...... I 9 9 . 0,0 ‘2. ------- : """" H r a . . i -1.0‘ ~ . . ' - ,-. - 3 90.99800 90.99850 90.99900 90.99950 91.00000 91.00050 91.00100 91.00150 X Coordinae of the interface Figure 17d. Deformed shape of the interface for Uniform axial load 51 Interface plots (N=2, M=7, NPS=4) — <>- auelglue27 ~ $.3- buelglue27 W1}. auelgluenm27 ---><-- buelgluenm27 ~~ {~19- . auelgluenmnps27 —- ~-~ buelgluenmnp527 . -~O~ . abaqus 9.0 Eh . ._ ‘ 3 8.0 ' r-fi " '7‘ - i g 7.0 f; ‘ , E - -; 6.0 9.9g ‘/ 5 5 0 ,. .5 ‘5 ' 99 0 4.0 -' .15 :g 3.0 1..; .- 4‘ o 2.0 #45; 7. I o 831‘; _.....->< >- 1.0 2.}; 99 0.0; 93%; 9 . . . _1 '0 .. . . . . .. .. .-_.. . - .. ... -,H-..... ..... ..-_.. ..- ._-. .._-.--.____._ 1' 90.95 91.00 91.05 91.10 91.15 91.20 91.25 X Coordinate of the interface Figure l7e. Deformed shape of the interface for Uniform axial load Table 1 and 2 contain a summary of the results presented in Figures 17a through 176. It can be clearly seen from the Figures 17 a through 17e, that the choice of collocation points suggested in Chapter 2 (p. 15-16, i.e. N, M and NPS as collocation points) matches the results obtained using penalty hybrid interface element and Abaqus 6.3. Also the other choice where only points on either sides (N and M) of the interface are used, fails to match the results obtained using penalty hybrid interface and Abaqus 6.3. 52 S. No. Mesh Type N-M-NPS Max I U1 _ U, I Tip Deflection lUmax | 1 3-2-4 0 1 2 4-3-4 0 1 3 8-10-4 0 1 4 8-10-8 0 1 5 2-7-4 0 1 Table 1. Summary of Uniform Axial Load for Uelgluenmnps. f s. No. Mesh Type N-M-NPS Max I U’ _ U' n Tip Deflection l U max I 1 3-2-4 0.1191 1.061 2 4-3-4 0.0020 1.001 3 8-10-4 0 1.001 4 8-10-8 0.0001 1.001 5 2-7-4 0.1723 1.213 Table 2. Summary of Uniform Axial load for Uelgluenm. f Figure 18a through 18e contain the error analysis for bending load, which is presented in very similar format as uniform axial load above. It can be clearly observed from the plots that the suggested choice of collocation points provides results, in accordance with those obtained using penalty based interface element and Abaqus 6.3. 53 Interface plots (N=3, M=2, NPS=4) — O- auelglue32 ~9 {3 - buelglue32 MA auelgluenm32 --/-- buelgluenm32 — -t~- auelquenmnps32 —— buelgluenmnpsSZ ~ 0 -- abaqus 8.0 /§ . ‘ /” . ' - . . . ' g /’/ 1A,.- A 6.0 297 ' 3 4 o // 0 , f/ . 3 / l a: .‘ - g 2.0 // 8 . .9 ‘ >- , ,‘fi/ 0.0 b 7- ' 11"" I I r 1 (Ir " 5? 89.875 89.925 89.975 90.025 90.075 90.125 X coordinate of interface Figure 18a. Deformed shape of the interface for bending load Interface plots (N-4, M83, ups-4) — o- aueIque43 . r} - bueIque43 -- Am weighenm43 ,9... bueIquenm43 —+-~ auelgluenmnps43 ---—-- bueIquenmnpuS 0 ~ abaqus 9.0 -- 8.0 1'77 7.0 xfififf‘ ' fil § ”/(g' —‘... 6.0 p. . '5 996%: g 5.0 .14" ‘ . (flag? g 4.0 «3‘71 . / ' . 30 /1H 1 g 2'0 /'za' ! 1.0 / >' . , h) f 0.0 V II r 1 7 fl '1 -° if" % '2.° ..L..~99-9 ..........._.._ . .9 . -..-9..._.....,,,,..m._-. 99-99.9W9-9-~9999. . W... . . . ...... .. . . ....... ..- .. . . 89.87 89.925 89.975 90.025 90.075 90.125 X Coordhate ofthe Interface Figure 18b. Deformed shape of the interface for bending load 54 Interface plots (NIB. III-10, NPs-4) — <>- auelglue8-10-4 ~~ £31 — buelglue8-10-4 WA eueIquean—1 0-4 ---><-- bueIquenm8-10-4 -»m—-- welgluenmnps8-10-4 --—-- buelquenmnp38-10-4 “+3... abaqus 90 T ....... -_......_ ......-___-___. i 8.0 r. / E 7.0 ’ L, § (*8 6.0 [.4 E ”/53 j s 5.0 XE’ 5 4.0 E? 4” i '5 ,./ ~ 1 " 3.0 . . I" g 2.0 9.7/5 a" " i O a g u 1.0 .7‘ ‘. >- //E i 0.0 [3&1 l I I 1| -1.0 Ex ! 89.875 89.925 89.975 90.025 90.075 90.125 X Coordlnde ofthe Interface Figure 18c. Deformed shape of the interface for bending load Interface plots (N=8, M810, NPS-8) — <>- auelglueO—10-8 {.14 buelglueB—IO-e ~~A~~euelgluenm8-10-8 -~-----bueIquenrn8-10-8 ~x—-- auelgluenmnpsB—IO—B —«—--buelgluenmnpsO-10-8 n abaqus 9'0 -. ..-N.......,.,..-...-..-.......m......,...... --.. ... i 8.0 r 7 o // 0 ' £537 ; ° 6.0 ’1’” 1 § {/9 E 5.0 E" i o x’ i g 4.0 4’W‘ ; g ,4“ . 1r 3-0 ’5 § 20 I,/ . 5 1.0 ‘ y/fl/g : > \ '/'/’E 0.0 Rate. . . . . 89.925 89.975 90.025 90.075 90.125 X Coordinate of the Interface Figure 18d. Deformed shape of the interface for bending load 55 Interface plots (N=2, M=7, NPS=4) - <> - auelglu027 ~ E} - buelglue27 We. euelgluonm27 —+~ auelgluenmnp527 —-—~- buelgluenmnp327 we."— abaqus ——.\— buelgluenm27 8.0 6.0 4.0 2.0 Y COOKING 0f the Intemce 0.0 89.875 Figure 186. f— , ”I." /’/ ..-r , w“ \" ,‘VA...V: '/ 4“, III. erjr f I,“ I I ‘4 I x Coordinate ofthe Interface Deformed shape of the interface for bending load 56 90.025 CHAPTER 5 CONCLUSIONS An effective and reliable method has been developed and implemented in popular commercial finite element software. This method based on the collocation method improves the computational efficiency. An effort was made to identify the collocation points. Starting with an initial guess as nodes on the interface three different choices were made. The most reliable choice was compared with penalty hybrid interface element, analytical solution and commercial finite element software. The accuracy of the results was observed to be the same as the other methods used for comparison. The resulting interface element can be used for large-scale problems in finite element and non-finite element problems (e. g FEM-FDM). This method is completely compatible with available commercial software. It can be used to model composite structures conveniently. It has all the features of penalty hybrid interface element. When it comes to computational efficiency the method proves to be much better for explicit analysis and has wide applications in non-finite element method. The collocation based interface element has been implemented in Abaqus as a User Element Subroutine (UEL) for two-dimensional problems. This made the testing of the collocation based interface element easier. This method can now be used to perform global/local analysis and optimization problems where the change occurs only for a small part of the model. The collocation method for interface technology can be applied to explicit solvers to capture its computational efficiency. 57 BIBLIOGRAPHY 58 [1]. [2]. [3]. [4]- [5]. [6]. [7]. [8]. [9]. BIBLIOGRAPHY C. J. Davila, J. B. Ransom and M. A. Aminpour, “Cross-Surface Interface Element for Coupling Built-up Structural Subdomains”, NASA, Langley Research Center, Technical Memorandum 109125, 1994. J. B. Ransom, S. L. McCleary and M. A. Aminpour. “ A New Interface Element for connecting Independently Modeled Substructures”, AIAA Paper, Number 93- 1503, 1993. J. M. Housner, M. A. Aminpour, C. G. Davila, J. E. Schiermeier, W. F. Stroud, J. B. Ransom, and R. E. Gillian. “An Interface Element for Global/Local and Substructuring Analysis”, presented at the MSC World Users’ Conference, Los Angeles, CA, May 8-12, 1995. J. E. Schierrneier, J. M. Housner, M. A. Aminpour and W. J. Stroud, “The Application of interface Elements to Dissimilar Meshes in Global/Local Analysis”, presented at the MSC World Users’ Conference, Newport Beach, California, June 3-7, 1996. J. E. Schiermeier, J. M. Housner, M. A. Aminpour and W. J. Stroud, “Interface Elements in Global/Local Analysis-Part 2: Surface Interface Elements, presented at the MSC Aerospace User’s Conference, Newport Beach, California, November 17-20, 1997. J. E. Schiermeier, R. K. Kansakar, J. B. Ransom, M. A. Aminpour and W. J. Stroud, “Interface Elements in Global/Local Analysis-Part 3: Shell-to-Solid Transition” presented at the MSC Worldwide Aerospace Conference, Long Beach, California, June 7-10, 1999. M. A. Aminpour, J. B. Ransom, and S. L. McCleary. “A Coupled Analysis Method for Structures with Independently Modeled Finite Element Subdomains”, International Journal for Numerical Methods in Engineering, Vol. 38, 1995, p. 3695-3718. J. T. Wang and J. B. Ransom, “Application of Interface Technology in Nonlinear Analysis of a Stitched/RFI Composite Wing Stub box”, AIAA paper, v 3 1997, p. 2295-2310. J. B. Ransom. “Interface Technology for Geometrically Nonlinear Analysis of Multiple Connected Subdomains”, AIAA Paper, Number 97-1298, 1997 59 [10]. [11]. [12]. [13]. [14]. [15]. [16]. [17]. [18]. [19]. M. A. Aminpour, and T. Krishnamurthy. “A Two-Dimensional Interface Element for Multi-Domain Analysis of Independently Modeled Three-Dimensional Finite Element Meshes”, AIAA Paper, Number 97-1297, p.1853-1861. 1997. M. A. Aminpour, T. Krishnamurthy and T. D. Fadale. “Coupling of Independently Modeled Three-Dimensional Finite Element Meshes with Arbitrary Shape Interface Boundaries”, AIAA Paper, Number 98-2060, p. 11853- 1861,1998. M. A. Aminpour, Stephane Pageau and Youngwon Shin, “ An Alternative Method for the Interface Modeling Technology”, AIAA Paper, Number 2000- 1352, p. 1-13, 2000. Jonathan B. Ransom, “On Multifunctional Collaborative Methods in Engineering Science”, Langley Research Center, Hampton, Virginia, NASA-2001-tm211046, 2001. R. Courant, K. Friedrichs and J. Lewy, “Uber die Partieller Differenzengleichungen der Mathmatischen Physik,” Mathematische Annalen. Vol. 100, 1928, p. 32; the English translation can be found in IBM journal, 1967, p 215. Babuska, I., “The Finite Element Method with Penalty,” Tech. Note BN-710, The Institute for Fluid Dynamics and Applied Mathematics, University of Maryland, August 1971. A. Pantano and Ronald C Averill, “A Penalty-Based Finite Element Interface Technology”, Computers and Structures, v80, n22, September 2002, p. 1725- 1748. A. Pantano, “A Penalty-Based Interface Technology for Connecting Independently modeled Substructures and for Simulating Growth of Delamination in Composite Structures”, PhD. Thesis, Department of Mechanical Engineering, Michigan State University, 2002. A. Pantano and Ronald C. Averill, “A Finite Element Interface Technology for modeling delamination growth in Composite Structures”, Submitted to the 43'd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, 2002. J. D. Whitcomb and K. Woo, “Evaluation of iterative Global/Local stress analysis”, ASME Enhancing Analysis Techniques for Composite Materials, 1991, p. 201-207. 60 [20]. [21]. [22]. [23]. [24]. [25]. [26]. [27]. [28]. [29]. [30]. [31]. [32]. [33]. [34]. [35]. T. Zhu and S. N Atluri, “A Modified Collocation Method and a penalty formulation for enforcing the essential boundary conditions in the element free Galerkin method”, Computational Mechanics, Springer Verlag, 1998. N. Kikuchi, “Finite Element Methods in Mechanics”, Cambridge University Press, 1986. J. N. Reddy, “The Penalty Function Method in Mechanics a review of recent advances”, ASME Conference, 1982. T. Belytschko and M. 0. Neal, “Contact-Impact by the Pinaball Algorithm with Penalty and Lagrangian Methods”, International Journal of Numerical Methods in Engineering, v 31, 1991, p. 547-572. Carl De Boor and Blair Swartz, “Collocation at Gaussian Points”, SIAM Journal of Numerical Analysis, Vol. 10, No. 4, September 1973. H. B. Nielsen, “Cubic Splines”, IMM, Department of Mathematical Modeling, Technical University of Denmark, 1998. Arthur Sard and Sol Weintraub. “A Book of Splines”, Wiley, 1971. Larry L. Shumaker, Spline functions: basic theory”. New York: Wiley, c 1981. P. M. Prenter. “ Splines and variational methods”, New York: Wiley, 1975. J. H. Ahlberg, E. N. Nilson. “ The theory of splines and their applications”. New York, Academic Press, 1967. H. Spath. “Two dimensional spline interpolation algorithms”, Wellesley, Mass: AK Peters, 1995. Abaqus v6.2 Manuals, Abaqus Inc. 2001. Calfem, A Finite Element Toolbox to Matlab Version 3.3, The Division of Structural Mechanics and the Department of Solid Mechanics at Lund University, February 1999. J. N. Reddy. “An Introductions to the Finite Element Method”. Highstown, NJ, McGraw-Hill, 1993. V. Rajaraman, “Computer Programming in Fortran 90 and 95”, Prentice-Hall of India, New Delhi, 1997. W. E. Mayo and M Cwiakala, “Programming with Fortran 77”, Schaum’s Outline Series, McGraw-HILL, 1995. 61 [36]. [37]. [38]. [39]. [40]. [41]. C. Chang and M. Gunzburger, “A Subdomain Galerkin/Least Squares Method for First Order Elliptic Systems in the Plane”, SIAM Journal of Numerical Analysis 27,1990, p. 1197-1211. Y. Cao and M. Gunzburger, “Least-Squares Finite Element Approximations to Solutions of Interface Problems”, SIAM Journal of Numerical Analysis 35, 1998, 393-405. I. N. Ioakimidis and P. S. Theocaris, “On the selection of Collocation Points for the numerical solution of singular integral equations with generalized kernels appearing in Elasticity problems”, Computers and Structures, v11, 11 4 April 1980, p 289-295. U. Ascher and G. Bader, “ Stability of collocation at Gaussian points”, SIAM Journal on Numerical Analysis, v23, n 2 April 1986, p412-422. D. Redekop, “ Fundamental solutions for the collocation method inplanar elastostatics”, Applied Mathematical Modelling, v 6, n 5, October 1982, p 240- 244. C. Farhat and M. Gerardin, “Using a reduced number of Lagrange Multipliers for assembling parallel solution algorithm”, Journal for Numerical Methods in Engineering, v 32, 1991, 1205-1227. 62 3 1293 2504