Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0000705/00001
## Material Information- Title:
- Implicit Solid Element for 2d Plane Stress Analysis
- Creator:
- YOUNG NAM HWANG (
*Author, Primary*) - Copyright Date:
- 2008
## Subjects- Subjects / Keywords:
- Boundary conditions ( jstor )
Circles ( jstor ) Coordinate systems ( jstor ) Curved beams ( jstor ) Engineering ( jstor ) Geometry ( jstor ) Quadrants ( jstor ) Stress distribution ( jstor ) Triangles ( jstor ) Vertices ( jstor )
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Nam Hwang Young. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Embargo Date:
- 7/1/2003
- Resource Identifier:
- 434616470 ( OCLC )
## UFDC Membership |

Downloads |

## This item is only available as the following downloads: |

Full Text |

PAGE 1 IMPLICIT SOLID ELEMENT METHOD FOR 2D PLANE STRESS ANALYSIS By YOUNG NAM HWANG A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2003 PAGE 2 Copyright 2003 by Young Nam Hwang PAGE 3 To my lovely family PAGE 4 ACKNOWLEDGMENTS I express my sincere gratitude to my advisor and the chairman of my thesis committee, Dr. Ashok V. Kumar, for his guidance, enthusiasm and continual support throughout my study. I thank him for his tireless effort in reading several drafts of this thesis. I am also grateful to my committee members, Dr. John K. Schueller and Dr. Nagaraj K. Arakere, for their advice and patience in reviewing this thesis. I also thank my group fellows, Mr. Jongho Lee and Mr. Siva Kumar, for their kind advice and help regarding programming. Finally, I would like to thank my family for their constant support and encouragement. Without them, I would not have been able to muster the courage to come over to America, let alone complete this thesis. iv PAGE 5 TABLE OF CONTENTS Page ACKNOWLEDGMENTS.................................................................................................iv LIST OF TABLES...........................................................................................................viii LIST OF FIGURES...........................................................................................................ix ABSTRACT......................................................................................................................xii CHAPTER 1 INTRODUCTION........................................................................................................1 1.1 Goals and Objectives..............................................................................................2 1.2 Outline....................................................................................................................3 2 NUMERICAL METHODS FOR ENGINEERING ANALYSIS.................................4 2.1 Finite Element Method...........................................................................................4 2.2 Mesh Free Method..................................................................................................7 3 GEOMETRY REPRESENTATION USING IMPLICIT SOLID ELEMENTS........10 3.1 Introduction...........................................................................................................10 3.2 Solid Modeling Technique...................................................................................11 3.3 Implicit Solid Elements........................................................................................13 3.4 Simple Primitives..................................................................................................16 3.5 Defining and Editing Primitives...........................................................................19 4 INTEGRATION ALGORITHM................................................................................21 4.1 Generating Mesh...................................................................................................21 4.2 Gaussian Quadrature.............................................................................................22 4.3 Quadtree Representation.......................................................................................24 4.4 Integration Using Quadtree Subdivision..............................................................25 4.4.1 Integration within an Internal Cell.............................................................27 4.4.2 Integration within a Partial Cell.................................................................29 4.4.2.1 One vertex of the partial cell is inside the geometry........................30 4.4.2.2 Two vertices of the partial cell are inside the geometry..................30 4.4.2.3 Three vertices of the partial cell are inside the geometry................31 v PAGE 6 4.4.3 Computing the Stiffness Matrix.................................................................31 4.5 Analytical Solution of Integration........................................................................32 5 GRAPHICAL DISPLAY ALGORITHM...................................................................40 5.1 Creating Geometry................................................................................................40 5.2 Creating TriangleArray.........................................................................................40 5.2.1 TriangleArray for an Internal Cell..............................................................40 5.2.2 TriangleArray for a Partial Cell..................................................................41 5.2.3 Example of TriangleArray..........................................................................42 5.3 Color Distribution.................................................................................................42 5.3.1 Smoothness.................................................................................................42 5.3.2 Compute the Value to the Intersection Point and the Mid Point................44 5.4 Deformed Shape...................................................................................................45 6 ANALYSIS AND RESULTS.....................................................................................47 6.1 Circular Disk with Concentrated Load.................................................................47 6.1.1 Problem Statement......................................................................................47 6.1.2 Creating Geometry and Mesh.....................................................................48 6.1.3 Analysis and Results...................................................................................48 6.1.4 Comparing the Results between ISEM and I-DEAS..................................51 6.2 Plate with a Circular Hole Under Uniaxial Tension.............................................52 6.2.1 Problem Statement......................................................................................52 6.2.2 Creating Geometry and Mesh.....................................................................52 6.2.3 Analytical Solution.....................................................................................53 6.2.4 Analysis and Results...................................................................................56 6.2.5 Comparing Results.....................................................................................57 6.3 Curved Beam 1.....................................................................................................57 6.3.1 Problem Statement......................................................................................57 6.3.2 Creating Geometry and Mesh.....................................................................58 6.3.3 Analytical Solution.....................................................................................59 6.3.4 Analysis and Results...................................................................................62 6.3.5 Comparing Results.....................................................................................63 6.4 Curved Beam 2.....................................................................................................64 6.4.1 Problem Statement......................................................................................64 6.4.2 Creating Geometry and Mesh.....................................................................65 6.4.3 Analytical Solution.....................................................................................66 6.4.4 Analysis and Results...................................................................................68 6.4.5 Comparing Results.....................................................................................69 7 CONCLUSIONS AND FUTURE WORK.................................................................72 7.1 Summary and Discussion.....................................................................................72 7.2 Future Work..........................................................................................................73 vi PAGE 7 APPENDIX A INPUT FORMAT FOR JAVAFEMVIEWER SOFTWARE.....................................74 B SCENEGRAPH FOR JAVAFEMVIEWER SOFTWARE........................................77 LIST OF REFERENCES...................................................................................................78 BIOGRAPHICAL SKETCH.............................................................................................80 vii PAGE 8 LIST OF TABLES Table page 4-1. Weight and location....................................................................................................24 4-2. Stiffness data...............................................................................................................39 6-1. Analysis results of circular disk..................................................................................52 6-2. Analysis results of plate with hole..............................................................................57 6-3. Analysis results of curved beam 1..............................................................................64 6-4. Analysis results of curved beam 2..............................................................................69 viii PAGE 9 LIST OF FIGURES Figure page 2-1. Typical mesh from I-DEAS. A) A 2D mesh example. B) A 3D mesh example..........5 2-2. Smoothed boundary. A) In FEM. B) In the Mesh Free method...................................9 3-1. Quadratic 9-node elements. A) Parametric space(r, s coordinates). B) Real space(x, y global coordinates)...................................................................................................14 3-2. Density distribution within a 9 node quadrilateral element. A) Solid view. B) Density fringes. C) Density contours.....................................................................................18 3-3. Editing the face represented using 2D element..........................................................20 4-1. Mesh generation. A) FEA. B) ISEM..........................................................................21 4-2. Quadtree representation..............................................................................................25 4-3. Quadtree representation for integration. A) Circle geometry with 4 elements. B) Quadtree subdivision of element 4. C) Quadtree representation of subdivision......26 4-4. Quadratic 4-node element in two different coordinate systems.................................27 4-5. Triangular element. A) Case 1. B) Case 2. C) Case3.................................................30 4-6. A quarter of circle geometry with global coordinate..................................................32 5-1. TriangleArray for an internal cell...............................................................................41 5-2. TriangleArray for a partial cell. A) One vertex and two intersection points. B) Two vertices and two intersection points. C) Three vertices and two intersection points.41 5-3. TriangleArray created for circle geometry.................................................................42 5-4. Von Mises stress before applying smoothing algorithm............................................43 5-5. Von Mises stress after applying smoothing algorithm...............................................44 5-6. Value of the point within the element........................................................................45 5-7. Deformed shape with circle geometry........................................................................46 ix PAGE 10 6-1. Circular disk...............................................................................................................47 6-2. Circular geometry with mesh and boundary conditions.............................................48 6-3. Displacement Magnitude by ISEM............................................................................49 6-4. Displacement Magnitude by I-DEAS.........................................................................49 6-5. Von Mises Stress by ISEM.........................................................................................50 6-6. Von Mises Stress by IDEAS......................................................................................50 6-7. Maximum Principal Stress by ISEM..........................................................................51 6-8. Maximum Principle Stress by IDEAS........................................................................51 6-9. Plate with a hole.........................................................................................................52 6-10. Plate with mesh and boundary condition..................................................................53 6-11. Circular hole in an infinite plane under uniaxial tension. A) Infinite plate with a small hole. B) distribution for = 2/, 3/2......................................................54 6-12. Stress X for 144 elements by ISEM.........................................................................56 6-13. Stress X for 140 elements by I-DEAS......................................................................56 6-14. Curved beam 1..........................................................................................................58 6-15. Curved beam1 with mesh and boundary conditions.................................................58 6-16. Analytical expression for A, R and A m ......................................................................60 6-17. Displacement Y by ISEM.........................................................................................62 6-18. Displacement Y by I-DEAS.....................................................................................62 6-19. Stress Y by ISEM.....................................................................................................63 6-20. Stress Y by I-DEAS..................................................................................................63 6-21. Curved beam 2..........................................................................................................65 6-22. Curved beam 2 with mesh and boundary conditions................................................66 6-23. Displacement Y by ISEM. A) Q4Node. B) Q9Node...............................................68 6-24. Displacement Y by I-DEAS. A) Q4Node. B) Q8Node............................................68 6-25. Stress Y by ISEM. A) Q4Node. B) Q9Node............................................................69 x PAGE 11 6-26. Stress Y by I-DEAS. A) Q4Node. B) Q8Node........................................................69 6-27. Displacement at point A...........................................................................................70 6-28. Stress at point B........................................................................................................70 6-29. Stress at point C........................................................................................................71 B-1. Scenegraph for the JavaFemViewer software............................................................77 xi PAGE 12 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science IMPLICIT SOLID ELEMENT METHOD FOR 2D PLANE STRESS ANALYSIS By Young Nam Hwang May 2003 Chair: Ashok V. Kumar Major Department: Mechanical and Aerospace Engineering The traditional method for generating mesh in Finite Element Analysis (FEA) requires the mesh to conform to the boundaries of the geometry. This requires complex meshing algorithms for automated mesh generation. The main goal of this project is to develop techniques for engineering analysis using a solid model representation based on implicit solid elements. In this representation, the boundaries of the solid are represented using implicit equations instead of the parametric equations used in traditional solid modeling techniques based on Boundary Representation. To perform the engineering analysis a regular mesh that does not conform to the geometry is generated. Integration algorithms were developed to compute the system equations using such a nonconforming mesh. This approach avoids the need for complex meshing algorithms and reduces the time required to generate finite element models. Graphic algorithms for displaying results on a nonconforming mesh were also developed. These new algorithms for nonconforming mesh were implemented in software that uses implicit solid elements to represent geometry. These algorithms were implemented by modifying previously developed-object oriented software for finite element analysis. xii PAGE 13 Using this software, 2D plane stress analysis was performed to verify the software, and the results were compared with commercial FEA software, I-DEAS, and analytical solutions. xiii PAGE 14 CHAPTER 1 INTRODUCTION There are various numerical methods for solving engineering problems such as heat transfer, electrostatic potential, fluid mechanics, and vibration analysis. Finite element analysis is a widely used numerical method and mesh free method is a newly developed numerical method. In finite element analysis, the spatial domain where the partial differential governing equations are defined is discretized into elements. A mesh is defined as any of the open spaces or interstices between the strands of net that is formed by connecting nodes in a predefined manner (Liu 2002). By using FEM with a properly predefined mesh, complex differential or partial differential governing equations can be approximated by a set of algebraic equations for the mesh. The system of algebraic equations for the whole problem domain can be formed by assembling sets of algebraic equations for all the elements. The mesh free method is used to establish a system of algebraic equations for the whole problem domain without the use of predefined mesh. Mesh free method use a set of nodes scattered within the problem domain as well as sets of nodes scattered on the boundaries of the domain to represent the problem domain and its boundaries. These sets of scattered nodes do not form a mesh, which means that no information on the relationship between the nodes is required. These numerical methods are discussed in detail in Chapter 2. A numerical method, referred to here as implicit solid element method (ISEM), is presented here that establish a system of algebraic equations using non-conforming mesh or a regular grid using implicit solid elements. Traditional methods for generating mesh 1 PAGE 15 2 in finite element analysis require the mesh to conform to the boundaries of the geometry. This requires complex meshing algorithms for automated mesh generation. To avoid complex meshing algorithm, we use a regular grid that is not required to conform to the geometry. The geometry used here for engineering analysis is defined using implicit solid element representation. This solid model representation scheme uses implicit curves or surface to represent the boundaries of the solid (Kumar and Lee 2002). The implicit curve/surface functions are defined by piece-wise interpolation within each implicit solid element. Within each element the density is defined as a parametric function. The solid is defined as the region within the element where the density is greater than a threshold or â€œlevel setâ€ value. The boundary of the solid is therefore a contour of the density function along which the density is equal to the level set value. The detailed discussion about solid modeling scheme using implicit solid element are presented in Chapter 3. 1.1 Goals and Objectives The goal of this project is to develop algorithm for non-conforming mesh that can perform the engineering analysis and display the results with the geometry that is created using implicit solid element. To achieve the above goals, the following objectives were identified. Perform literature survey for the numerical methods for engineering analysis to develop and implement a suitable algorithm for the ISEM process. Develop integration algorithms for non-conforming mesh and display algorithm for geometry created using implicit solid method. Implement these algorithms in a software by extending previous implemented FEM software. PAGE 16 3 1.2 Outline The topics covered in the rest of this thesis are organized as follows: Chapter 2 discusses numerical methods for engineering analysis such as Finite Element Method and Mesh Free Method. The concept and advantages of both methods are briefly presented. Chapter 3 summarizes the geometry representation using implicit solid elements and discusses other solid modeling techniques. Chapter 4 discusses the integration algorithm for the non-conforming mesh using implicit solid element analysis. The result of the integration using this new integration algorithm is compared with the analytical solution in this chapter. Chapter 5 discusses the graphical display algorithms for the implicit solid element analysis. The algorithm for creating geometry array, smoothing colors, and deforming shape are presented. Chapter 6 presents the examples of analysis and discusses the results. And finally Chapter 7 draws conclusions based on the previous chapters and suggests future work in this area. PAGE 17 CHAPTER 2 NUMERICAL METHODS FOR ENGINEERING ANALYSIS There are various numerical methods, such as finite element method (FEM), finite difference method (FDM), mesh free method and so on, for engineering analysis. In this chapter, we will discuss FEM from which the implicit solid element analysis is derived, and also briefly discuss mesh free method. 2.1 Finite Element Method The finite element method is a powerful numerical method for solving engineering problems and can analyze complex, structural, mechanical and electrical system. The finite element method is used to analyze both linear and nonlinear system and to solve problems in static and dynamic systems. Because of its various applications such as heat transfer, electrostatic potential, fluid mechanics, vibration analysis and so on, the finite element method has become very popular. For many physical problems, it is very difficult to find an analytical solution for complex geometry with complex boundary conditions. In the finite element method, a complex region defining a continuum is discretized into simple geometric shapes called finite elements that are connected at specified node points as shown in Figure 2-1. The shapes of the elements are intentionally made as simple as possible, such as triangles and quadrilaterals in two-dimensional domains, and tetrahedra, pentahedra, and hexahedra in three dimensions. The entire mosaic-like pattern of elements is called a mesh. 4 PAGE 18 5 A B Figure 2-1. Typical mesh from I-DEAS. A) A 2D mesh example. B) A 3D mesh example. The material properties and the governing relationships are considered over these elements and expressed in terms of unknown values at element corners. These equations describe the physical problem that is to be analyzed and generally can be expressed using the weak form of a variational principle. VSVdVdSdV}{}{}{}{}]{][[][}{bXfXXBDBXTTTT (2-1) Eq. 2-1 represents the weak form applied to an element. V is the volume of the element, S is the surface area of the element that lies on the boundary of the object, {f} is the transaction acting on these boundaries and {b} is the body force acting on the element. The vector {X} contains the nodal variables while the {X} vector contains the corresponding virtual variables. The [D] matrix contains material properties associated with the analysis and the [B] matrix denotes an operator that, when applied to the nodal variable vector, yields some physical quantity associated with the analysis. From the Eq. 2-1, the element stiffness can be written as PAGE 19 6 VdV]][[][][BDBKTe (2-2) For different kinds of problem or analysis type, the matrix [D] is different and depends on the constitutive equation and the material properties. For two-dimensional plane-stress analysis, the [D] matrix is 210001011][2ED (2-3) where E and are Youngâ€™s Modulus and Poissonâ€™s ratio respectively. The matrix [B] also varies for different analysis types. For example, the [B] matrix for two-dimensional plane-stress analysis is xNyNxNyNxNyNyNyNyNxNxNxNBnnnn...0...000...00][22112121 (2-4) where n is the number of nodes in the element and N i is the ith interpolation (or shape) function. In each element the governing equations (usually in the integral form shown in Eq. 2-1) are transformed into algebraic equations that are an approximation of governing equations. In an assembly process, the algebraic equations for all elements are assembled to achieve a system of equations for the model as a whole. Solving these equations gives us the approximate behavior of the continuum. With a sufficiently large number of elements, the assembly of the individual element solution can represent a reasonably accurate solution over the complete region. PAGE 20 7 2.2 Mesh Free Method As we mentioned above, complex partial differential equations are largely solved using numerical method, such as the finite element method (FEA). A mesh must be predefined to provide a certain relationship between the nodes, which is the basic formulation of this numerical method. Creation of a mesh for the problem domain is a prerequisite in using FEM packages and the analyst spends the majority of their time in creating the mesh. The stresses obtained using FEM packages are discontinuous and less accurate if the mesh density and quality is poor. When handling large deformation, considerable accuracy is lost because of the element distortion. The root of these difficulties is the need to create a mesh and the idea of eliminating the mesh has evolved. The mesh free method is used to establish a system of algebraic equations for the whole problem domain without using a predefined mesh. Mesh free methods use a set of nodes scattered within the problem domain as well as node scattered on the boundaries of the domain to represent the problem domain and its boundaries. These sets of scattered nodes do not form a mesh, which means that no information on the relationship between the nodes is required, at least for field variable interpolation. In Mesh Free method (Liu 2002), adaptive schemes can be easily developed, as there is no mesh, no connectivity involved. This provides flexibility in adding or deleting points/ nodes whenever and wherever needed. The analyst can save the time spend on conventional mesh generation because there is no mesh, and the nodes can be created by a computer in a fully automated manner. This can translate into major cost and time savings in modeling and simulation projects. PAGE 21 8 The fundamental difference between FEM and Mesh Free method is the construction of the shape functions. In FEM, the shape functions are constructed using elements. These shape functions are therefore predetermined for different types of elements before the finite element analysis starts. In Mesh Free method, the shape functions constructed are usually only for a particular point of interest, and the shape function changes as the location of the point of interest changes. The construction of the element free shape function is performed during the analysis. Modeling the geometry is also a difference between these two methods. In FEM, curved parts of the geometry and its boundary can be modeled using curves and curved surfaces using high-order elements. If leaner elements are used, these curves and surfaces are straight lines or flat surfaces. Figure 2-2A shows an example of smooth boundary approximated in the finite element model by the straight line edges of triangular elements. The accuracy of representation of the curved parts is controlled by the number of elements and the order of the elements used. In Mesh Free methods, the boundary is represented by nodes as shown in Figure 2-2B. At any point between two nodes on the boundary, one can interpolate using Mesh Free shape functions. Because the Mesh Free shape functions are created using nodes in a moving local domain, the curved boundary can be approximated very accurately even if linear polynomial bases are used. It is common in Mesh Free methods to use higher-order polynomials. PAGE 22 9 A B Figure 2-2. Smoothed boundary. A) In FEM. B) In the Mesh Free method. There are a number of Mesh Free methods, such as the element free Galerkin (EFG) method (Belytschko et al. 1994), the meshless local Petrov-Galerkin (MLPG) method (Atluri and Zhu 1998), the finite point method (Onate et al. 1996), and so forth. PAGE 23 CHAPTER 3 GEOMETRY REPRESENTATION USING IMPLICIT SOLID ELEMENTS The geometry used here for engineering analysis is obtained using implicit solid element representation. This chapter summarizes this geometry representation scheme (Kumar and Lee 2002). 3.1 Introduction The boundaries of the solid are represented using implicit curves or surfaces in the implicit solid element scheme. The implicit curve/surface functions, referred to here as the shape density functions (or density functions), are defined by piece-wise interpolation within each implicit solid element. These elements can be triangular or quadrilateral in 2D and tetrahedral elements in 3D. Within each element the density is defined as a parametric function. The solid is defined as the region within the element where the density is greater than a threshold or â€œlevel setâ€ value. The boundary of the solid is therefore a contour of the density function along which the density is equal to the level set value. Simple shapes or primitives can often be defined using single implicit solid element while more complex geometry can be represented either using a grid/ mesh of elements or by Boolean combination of simple primitives. The geometry of the primitives can be edited or modified by changing the density values at the nodes of the elements, by moving the nodes as well as by changing the level set value. Since the boundaries of the solid are implicit surface equations, it is easier to perform set membership classification and to evaluate Boolean operation, thus making this representation particularly suited for representing primitives in a CSG tree. 10 PAGE 24 11 3.2 Solid Modeling Technique Solid modeling techniques used by contemporary Computer Aided Design (CAD) software evolved in the early 1970â€™s (Hoffmann 1987, Mantyla 1988, Mortenson 1997). There are two popular approaches to solid modeling. One approach, known as Constructive Solid Geometry (CSG) uses a procedural representation of solid construction. The solid model is constructed by the Regularized Boolean combination of half-spaces defined using implicit equations. The construction of the solid is represented using a tree data structure (CSG tree) where the leaf nodes are primitives and the branch nodes represent the Boolean operation. The other approach is Boundary Representation (or B-Rep) where parametric representation of curves and surfaces is used to define the boundaries of the solid and the connectivity between them (or topological/combinatorial structure) is explicitly defined as relations between vertices, edges, half-edges, loops, faces and shells (Lee 1999). Parametric representation of curves and surfaces enables the representation of arbitrarily complex geometry using B-Rep. However, creating B-Rep models directly by defining each vertex, edge and surface is extremely cumbersome. Therefore, hybrid systems were developed where the solid is represented procedurally using a CSG tree but the primitives are represented using B-Rep. The B-Rep model of the solid represented by the CSG tree is evaluated automatically using â€œBoundary Evaluatorâ€ algorithms (Requicha et al. 1985). For each regularized Boolean operation, the boundary evaluator algorithm, detects intersection between participating solid, computes intersection geometries (intersection points and curves as well as sub-divided faces) and classifies them using set membership classification algorithm to determine whether the geometry is in, on or outside the final solid defined by the Boolean operation. While set membership classification is easy when the solid is represented as a combination of half PAGE 25 12 spaces (where the surfaces are implicit equations, f(x,y,z) = 0), it is more difficult to do for B-Reps. Furthermore, constructing the topology of the resultant solid automatically is also difficult. As a result it is very expensive to build robust and reliable software that can handle every special case. However, over the years many commercial systems have been developed that have very reliable boundary evaluator algorithms that are robust and can handle almost all special and degenerate cases including non-manifold topologies (ACIS 1995). To cope with the high cost of implementing these algorithms, many companies in the CAD industry buy â€œgeometric modeling kernelsâ€ from other companies that have already implemented them. We have used an alternate representation for solid primitives using implicit solid elements. Traditionally, implicit curves and surfaces are represented as f(p) = c where p R 2 for planar curves and p R 3 for surfaces. They are sometimes referred to as level sets or isocurve/surface since the curve or surface corresponds to a constant value of the function f(p). If c = 0 then the curve or surface is called the zero set of f(p). An implicit surface divides space into regions where f(p) > 0 and f(p) < 0. Therefore, if the boundaries of the solid are represented implicitly as half-spaces one can use the sign of f(p) to perform set membership classification that is to determine whether a given point is inside, outside or on the boundary of the solid. Despite this advantage, implicit representation for curves and surfaces has traditionally been considered inferior to parametric representation. The reasons suggested for this include axis dependence of implicit curves, difficulty in tracing or generating graphics as well as difficulty in fitting and manipulating freeform shapes. Considerable progress has been made to solve some PAGE 26 13 of these problems (Bloomenthal 1997) due to which implicit curves and surfaces now have numerous application in graphics and animation. To enable graphical display a variety of polygonization, tessellation or tracing algorithms have been developed (Lorensen and Cline 1987). Ray tracing algorithms have also been used for visualization of implicit surfaces (Glassener 1989). When the implicit function, f(p) is a polynomial the surface is called algebraic surface. Most common primitive shapes such as sphere, ellipse, cone and cylinders can be expressed as algebraic surface using quadratic polynomials. For more complex shapes, implicit algebraic surface patches (or A-splines) have been developed (Bajaj et al. 1995). These patches can be used for C 1 and C 2 interpolations or approximations and also for interactive freeform modeling schemes. Implicit curves and surfaces have also been used for visualization problem such as shape reconstruction from unorganized data sets (Zhao et al. 2000). Solid modeling using implicit solid elements is explained below where the boundary of the solid models are implicitly defined curves / surfaces that are defined mainly using quadrilateral and hexahedral implicit solid elements. Within each element, the implicit curve or surface function (or the shape density functions) are defined by interpolating the density values defined at the nodes of the element. 3.3 Implicit Solid Elements The implicit curves and surfaces are represented here as the level set of a implicit surface function bsr ),( for planar faces or btsr ),,( for solids. The implicit surface function, , has been referred to as the shape density function (or simply density function). The density function has a value that greater than a threshold or level set value, b , inside the solid and less than outside. Unlike traditional approach, where the implicit PAGE 27 14 surface function is defined directly in terms of the Cartesian coordinates, the density function used here is a parametric function whose arguments are the parametric coordinates of an implicit solid element. This is illustrated using a two dimensional example in Figure 3-1 where the element has nine nodes and its geometry in the parametric space is shown in Figure 3-1A while its real geometry is shown in Figure 3-1B. In the parametric space the element is a square whereas in the real space the element can be distorted quadrilateral. The mapping between the parametric coordinates and the real coordinates x, y, and z is defined as, ),,()( and ),,()( ,),,()(111tsrMzr,s,tztsrMyr,s,tytsrMxr,s,txniiiniiiniii (3-1) In Eq. 3-1 are the coordinates of the node i and are mapping basis functions that define the mapping between parametric space ( and the real space . ),,(iiizyx ),,(tsrMi),,tsr ),,(zyx b s 5 1 5 1 2 8 2 9 r 8 6 9 y 4 6 7 3 3 4 7 x A B Figure 3-1. Quadratic 9-node elements. A) Parametric space(r, s coordinates). B) Real space(x, y global coordinates) The density function is defined within the elements by specifying the value of the function at the nodes of these elements. The value of the function within each element is PAGE 28 15 obtained by interpolating the values at the nodes using appropriate basis functions or interpolation functions. A parametric interpolation scheme is used as shown below where parameters vary between and 1. ),,(tsr niiisrNsr1),(),( for 2D elements (3-2) niiitsrNtsr1),,(),,( for 3D elements (3-3) Where, i , is density at the node i, N i are the interpolation basis functions used to interpolate the density and n is the total number of nodes in the element. Each element is a cube (or square in 2D) of side length equal to 2 in parametric space. The parameters r, s, and t vary from to 1 within each element. If the mapping basis functions M are identical to the density interpolation basis functions then the element can be referred to as â€œisoparametric elementâ€. However, it is often beneficial to use simpler basis function for mapping. ),,(tsri ),,(tsrNi In this formulation, even though the boundary of the solid is represented in the implicit form bzyx ),,( , the density function itself is represented in a parametric form, ),,(tsr with a mapping defined between the parameters ( and the global coordinates . Therefore, many of advantages traditionally associated with parametric curves and surfaces are also applicable to the implicit solid element approach including axes independence and the ease in tracing. The parametric form of the density function enables the development of simple algorithm to polygonize the surface for graphical display. Based on the above description, the solid may now be defined as the set of points S such that ),,tsr ),,zyx ( PAGE 29 16 S= 11,11,11),,(,),,(11tsrwheretsrMXXtsrNXniiinibii (3-4) Where X R 3 for 3D elements. The boundary of the solid represented by an element can be defined as the set of points B such that B= 111,),,(,),,(11,11,11),,(,),,(1111torsorrwheretsrMXXtsrNXANDtsrwheretsrMXXtsrNXniiinibiiniiinibii (3-5) A variety of elements and basis functions can be defined and used to represent the objects in 2D and 3D. Higher degrees polynomials can be used as basis function for elements with larger number of nodes per element. If the geometry to be represented is simple (e.g. rectangle, circle or ellipse) then it is often possible to define the required shape density function using just one element. Primitives thus defined using one or more elements can be combined using various Boolean operations to create a procedural definition of more complex solids, expressed as a CSG tree. In the next section, a few simple primitives defined using implicit solid elements are described. 3.4 Simple Primitives Any quadrilateral or hexahedral primitives can be represented directly using quadrilateral or hexahedral elements respectively by setting the nodal values of density to be greater than or equal to the level set value. While it is possible to create a variety of implicit solid elements, we have primarily used the 9-node quadrilateral element for two PAGE 30 17 dimensional primitives such as quadrilaterals, circles, and ellipses. The quadrilateral element shown in Figure 3-1A has the following biquadratic basis functions. )1)(1(411srrsN , )1)(1(412srrsN , )1)(1(413srrsN )1)(1(414srrsN , )1)(1(2125srsN , )1)(1(2126srrN )1)(1(2127srsN , )1)(1(2128srrN , (3-6) )1)(1(229tsN The same basis functions could be used to define the mapping between parametric and real space defined by Eq. 3-1. This would enable the element in real space to be distorted as shown in the Figure 3-1. However, we have used the following bilinear 4 node interpolation to simplify the mapping with the assumption that the edges of the element will remain straight lines in the real space. srM11411 , srM11412 srM11413 , srM11414 (3-7) Figure 3-2 shows the example of circle, which is defined using a single 9-node quadrilateral element. The boundary is defined as . Simplifying the expression for 915.0),(ibiiNsr ),(sr and setting the coefficients of r, s, rs, r 2 s, rs 2 and r 2 s 2 equal to zero, we get the following conditions for the equation 5.0),(sr to be a circle. 1 = 2 = 3 = 4 , 5 = 6 = 7 = 8 , 9 =2 5 1 , 9 5 (3-8) A possible set of nodal density values that satisfy these conditions is: 1 = 2 = 3 = 4 = 0, 5 = 6 = 7 = 8 = 0.5, and 9 =1. The density distribution within the element PAGE 31 18 for these values of nodal densities is shown in Figure 3-2. Figure 3-2A shows the solid view of circle and Figure 3-2B shows the density distribution in a gray scale (density range 0 to 1 corresponds to gray scale white to black). Figure 3-2C shows contour lines corresponding to constant values of shape density. The boundary of the solid is one of the contours, whose value is equal to the boundary value b . Assuming b = 0.5, the boundary of the solid is the circle highlighted in Figure 3-2B. This results in a circle of unit radius within the element in the parametric space whose equation is r. The geometry in the real space will also be a circle if the shape of the element is square and diameter of the circle will be equal to the size of the square. If the element shape is a rectangle in real space, then the circle in parametric space maps into an ellipse in real space with major and minor diameters equal to the length and height of the rectangle. 122 s 0 0.2 0.4 0.7.80.80.99 A B -1 -0.8 -0.6 -0.4 -0.2 C 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 rs0.40.40.50.50.50.50.50.60.60.60.60.60.60.60.60.70.70.70.70.700.80.80.80.90.0.9 Figure 3-2. Density distribution within a 9 node quadrilateral element. A) Solid view. B) Density fringes. C) Density contours. By providing a variety of element types such as quadratic 4-node, 9-node and triangular 3-node a very large variety of primitives can be created. Rectangular primitives can be created from both the quadratic 4-node and 9-node quadrilateral elements. However, circle primitives can be created only from 9-node elements. Three-dimensional primitives can be defined using hexahedral elements such as hexahedral 8-node or 18-node elements. PAGE 32 19 3.5 Defining and Editing Primitives The geometry defined within each element can be edited by modifying the density values at each node, the nodal coordinates of the element or level set value of density. This diversity in editing methods enables the creation of a variety of primitives using a single element. Figure 3-3 shows how the solid geometry is changed in response to the change in the density values. In Figure 3-3A, the density value at the nodes of the 9-node element is set such that the density at all nodes is set to b = 0.5, except the central node where density is set equal to 1. The geometry is therefore identical to that of the element since every point within has density greater than b , The geometry in Figure 3-3B is obtained by changing the density values at nodes to 1 = 2 = 3 = 4 = 1, 5 = 6 = 7 = 8 = 0.7 and 9 = 0.4. The geometry can also be modified by changing the nodal coordinates as shown in Figure 3-3C where the element nodes are moved to stretch the element into a rectangle. This method of editing the geometry enables a representation of a variety of quadrilateral primitives, such as rectangles, trapezoids and parallelograms, to be constructed by deforming the element in Figure 3-3A. Therefore, the nodal points of the element can be considered to be similar to control points in the sense that moving their position can change the geometry. PAGE 33 20 A B C D Figure 3-3. Editing the face represented using 2D element Yet another way to modify the geometry is by changing the level set, b , so that a different level set or contour of the density function becomes the boundary of solid. This is illustrated in the Figure, where the solid in Figure 3-3D is obtained from the solid in Figure 3-3C by changing the level set value from 0.5 to 0.45. PAGE 34 CHAPTER 4 INTEGRATION ALGORITHM 4.1 Generating Mesh Mesh generated for FEA is shown in Figure 4-1A. The mesh coincides with geometry, and we can integrate function over each element to get element stiffness matrix (equation number) as we discussed in Chapter 2. Generating mesh in ISEM is different from generating mesh in FEA. In Figure 4-1B, the geometry is defined using a quad 9-node solid element as we discussed in Chapter 3. Then we divide the solid element to the proper number of quad 4-node or 9-node finite elements to create nonconforming mesh as shown in Figure 4-1B. Unlike mesh in FEA, the mesh in ISEM does not coincide with geometry, so we need new algorithm for proper integration at each finite element. Mesh Geometry A B Figure 4-1. Mesh generation. A) FEA. B) ISEM. 21 PAGE 35 22 4.2 Gaussian Quadrature In traditional FEA, functions are integrated over elements using Gaussian quadrature (Tirupathi and Ashok 1997). Consider the problem of numerically evaluating a one-dimensional integration of the form I = (4-1) 11)(drrf The Gaussian quadrature approach for evaluating I is given below. This method has proved most useful in finite element analysis. Consider the n-point approximation I = (4-2) 111)()(niiirfwdrrf Where w 1 , w 2 , , w n are the weight and r 1 , r 2 , , r n are the Gauss points. The idea behind Gaussian quadrature is to select the n Gauss points and n weight such that Eq. 4-2 provides an exact answer for polynomials f(r) of as large a degree as possible. In other words, the idea is that if the n-point integration formula is exact for all polynomials up to as high a degree as possible, then the formula will work well even if f is not a polynomial. Consider the formula with n=1 as I = (4-3) 1111)()(rfwdrrf Since there are two parameters, w 1 and r 1 , we consider requiring the formula in Eq. 4-3 to be exact when f(r) is polynomial of order 1. Thus if , then we require rCCrf10)( I = 011210111110221)(CrCrCdrrCCdrrf I G = 11101110111)()(rCwCwrCCwrfw I-I G =0 2111010 rCwCwC (4-4) PAGE 36 23 From Eq. 4-4, we see that the error is zeroed if w 1 = 2 , r 1 = 0 For any linear function f, then, we have I = 11)0(2)(fdrrf Consider the formula with n = 2 as )()()(I221111rfwrfwdrrf (4-5) We have four parameters to choose: w 1 , w 2 , r 1 , and r 2 . We can therefore expect the formula in Eq. 4-5 to be exact for a cubic polynomial. Thus, choosing yields. 332210)(rCrCrCCrf 2011114332210322432)(ICCrCrCrCrCdrrf )()(I2211Grfwrfw )()(32322221023132121101rCrCrCCwrCrCrCCw 3322311222221112211021)()()()(CrwrwCrwrwCrwrwCww 0IIG (4-6) From Eq. 4-6, we see that the error is zeroed if 121ww , 5773.03121rr For any cubic function f, then, we have )5773.0()5773.0()(I11ffdrrf PAGE 37 24 From the above discussion, we can conclude that n-point Gaussian quadrature will provide an exact answer if f is a polynomial of order (2n-1) or less. Table 4-1. Weight and location. Number of points, n Weight, w i Location, r i 1 2 0.0 2 1 0.57735 3 0.55555 0.88888 0.77459 0.0 4 0.34785 0.65214 0.86113 0.33998 Gaussian quadrature can be extended to two-dimensional integration of the form I = njniiijisrfwwdrdssrf111111),(),( Consider the formula with n=1 as )0,0(*2*2)0,0(),(111111ffwwdrdssrf Consider the formula with n=2 as ),(),(),(),(),(22221212212111111111srfwwsrfwwsrfwwsrfwwdrdssrf )577.0,577.0()577.0,577.0()577.0,577.0()577.0,577.0( ffff 4.3 Quadtree Representation Quadtree representation concept has been used to create approximated representation of geometry (Lee 1999). It represents a two-dimensional object, such as the one shown in Figure 4-2A, as a set of squares of different sizes by recursively subdividing the root square enclosing the object. Figure 4-2B shows the process of subdividing object. PAGE 38 25 A B Figure 4-2. Quadtree representation A root square enclosing the object is created and subdivided into four quadrants by halving its sides. Then each quadrant is classified according to its relative location with respect to the object. If the quadrant is neither completely inside nor completely outside the object, it is subdivided into four quadrants by halving its sides. Then each quadrant is classified according to its relative location with respect to the object. If the quadrant is neither completely inside nor completely outside object, it is subdivided again. This subdivision process is repeated until the mesh density requirement is satisfied and the quadrants that are either â€œcompletely insideâ€ or â€œoverlappingâ€ the object are collected. When the â€œoverlappingâ€ quadrants are collected, they are modifed to include only the inside of the object. Thus the object is represented by the collection of the â€œcompletely insideâ€ quadrants and the modifed â€œoverlappingâ€ quadrants. 4.4 Integration Using Quadtree Subdivision Integration using quadtree subdivision is a similar concept to a quadtree representation. Instead of generating meshes, Quadtree subdivision concept is used for the proper integration within each finite element. PAGE 39 26 (3) (2) (4) (1) Geometry Mesh A 1 29 30 20 21 22 23 24 26 25 32 33 27 31 34 28 14 15 16 17 18 19 8 9 12 11 10 13 2 3 4 5 6 7 B N ode representing a partial cell N ode representing an internal cell N ode representing an external cell N ode with descendants 1 34 33 29 30 24 22 16 15 13 14 12 11 10 21 20 19 117 8 28 31 32 27 26 25 8 9 7 6 2 3 4 5 23 C Figure 4-3. Quadtree representation for integration. A) Circle geometry with 4 elements. B) Quadtree subdivision of element 4. C) Quadtree representation of subdivision. An example of a circular geometry represented using 4 finite elements is shown in Figure 4-3A. Figure 4-3B shows the process of recursively subdividing the element, and Figure 4-3C shows the quadtree representation of this subdivision. If element is completely outside the object, there is no need to integrate within this element. If the element is completely inside object, integration is carried out by gauss quadrature. If the element is neither completely inside nor completely outside the object, it is subdivided PAGE 40 27 into four quadrants by halving its sides. Then each quadrant is classified according to its relative location with respect to the object. If the quadrant is completely outside, it will be referred to as an external cell, and there is no need to integrate within this cell. If the quadrant is completely inside, it is an internal cell, and integration with this cell is accomplished by gauss quadrature. If the quadrant is neither internal cell or external cell, it is a partial cell and that is subdivided again. This subdivision process is repeated until a termination criteria is satisfied and the quadrants either â€œcompletely insideâ€ or â€œoverlappingâ€ the object are collected. A quadrant which is overlapping the object is called as a partial cell. Integration within partial cells, it is accomplished by partial cell integration method discribed later. The total integration area is therefore represented by a collection of the internal cells and the partial cells. 4.4.1 Integration within an Internal Cell u v r s (bi, di) (ai, ci) i t h cell x T2 T4 T1 3 ( -1 , -1 ) 2 (-1,1) 4 ( 1 , -1 ) 1 (1,1) T3 y Real S p ace Parametric S p ace Figure 4-4. Quadratic 4-node element in two different coordinate systems. PAGE 41 28 We use the isoparametric coordinate for integration as shown in Figure 4-4. The integration within parametric space can be represented as the sum of the integration within cells. a i , b i are the r coordinate of the i th cell and, c i, d i are the s coordinate of the i th cell. celliiiinirsdcbarsAdrdssrfdrdssrfdAyxf11111det),(det),(),(][JJKe sysxryrxJ , rysxsyrxdrdsdArs Jdet To integrate within the i th cell, we need a local coordinate system (u, v) for the cell as shown Figure 4-4. Integration within the cell can be represented as dudvvufdrdssrfuvrsrsdcbaiiiiJJJKcelldetdet),(det),(][1111 vyvxuyuxJ , uyvxvyuxdudvdrdsrs Jdet u,v is subdivided coordinate system. Considering related formula of r, s coordinate to u, v coordinate is iiiauabr)1(2)( 2)(iiabur 0 vr iiicvcds)1(2)( 0 ur 2)(iicdvs 2)(2)(detiiiiuvcdabuyvxvyux J Integral equation for each cell by quadtree subdivision drdssrfrsdcbaiiiiJKcelldet),(][ PAGE 42 29 dudvcvcdauabfuvrsiiiiiiJJdet]det[)1(2)(,)1(2)(1111 For example, stiffness matrix of the quadrant number 16 and 22 in Figure 4-4 is [ drdssrfrsJKe16det),(]75.05.025.05.0 dudvvufrs1111321)1(81,21)1(8141detJ drdssrfrsJKe22det),(][5.005.01 dudvvufrs11112)1(41,1)1(4141detJ 4.4.2 Integration within a Partial Cell Each partial cell is divided into triangular elements based on the number of vertices inside the geometry in the partial cell as shown in Figure 4-5. Then the area of the triangle is computed and the area ratio between the area of the triangle and cell is calculated. The integral of the function within the partial cell is approximated as area ratio. There are three cases for creating triangular element in partial cell as shown in Figure 4-5. The shaded area is inside the geometry. The intersection point between the cell and the geometry is represented as circle, , and the vertex inside the geometry is represented as square, . drdssrfrsdcbaiiiiJdet),( PAGE 43 30 Geometry (bi, di) A (bi, di) C (bi, di) B B B A A C D C (ai, ci) A B (ai, ci) C (ai, ci) Figure 4-5. Triangular element. A) Case 1. B) Case 2. C) Case3. 4.4.2.1 One vertex of the partial cell is inside the geometry One vertex, point C, is inside the geometry as shown in Figure 4-5A. Compute the intersection points of the geometry and the cell (point A, B). Make a triangle with intersection points (point A, B) and the vertex inside the geometry (point C). Compute the area of a triangle (triangle ABC). Compute the area ratio (area ratio=triangle area/cell area). [K partial cell ]=area ratio[K cell ]=area ratio drdssrfrsdcbaiiiiJdet),( Cell number 11, 14 and 25 in Figure 4-3B are the examples of this case. 4.4.2.2 Two vertices of the partial cell are inside the geometry Two vertices, point C and D, are inside the geometry as shown in Figure 4-5B. Compute the intersection points of the geometry and cell (point A, B) Make two triangles with intersection points (point A, B) and the vertices inside the geometry (point C, D) Compute the area of two triangles (triangle ABC and ACD) PAGE 44 31 Compute the area ratio (area ratio=(triangle area1+triangle area2)/cell area) [K partial cell ]=area ratio[K cell ]=area ratio drdssrfrsdcbaiiiiJdet),( Cell number 4, 12, 15, 18, 28 and 31 in Figure 4-3B are the examples of this case. 4.4.2.3 Three vertices of the partial cell are inside the geometry Three vertices, point D, E and F, are inside the geometry as shown in Figure 4-5C. Compute the intersection points of the geometry and cell (point A, B) Make a triangle with intersection points (point A, B) and vertex outside the geometry (point C) Compute the area of a triangle (triangle ABC) Compute the area ratio (area ratio=(mesh area-triangle area)/cell area) [K partial cell ]=area ratio[K cell ]=area ratio drdssrfrsdcbaiiiiJdet),( Cell number 6, 17, 26 and 34 in Figure 4-3B are the examples of this case. 4.4.3 Computing the Stiffness Matrix Both the internal cells and the partial cells are integrated and added. Stiffness matrix for element (4) in Figure 4-3B is computed below. There are two material properties for isotropic materials. We have used the values for Youngâ€™s Modulus E = 2.0710 11 and Poissonâ€™s ratio = 0.3. Stiffness matrix of shaded area in Figure 4.-3A is celliiiinirsdcbaAdrdssrfdAyxf1det),(),(][JKe ][...][][][][...][][][cell34cell14cell12cell11cell9cell5cell4eKKKKKKKK [K e ] = PAGE 45 32 1056.9835877773.06870314-090.59533997940.51486518-74.17890830-32.7885435183.40001943-350.7950248273.06870314-08.71248308550.00292685-96.05985836-62.6555515674.17890830-620.4160784361.52628359090.59533997550.00292685-69.6434267023.5864951996.05985836-940.51486518-74.17890830-73.06870314-940.51486518-96.05985836-23.5864951969.64342670550.00292685-090.5953399773.06870314-74.17890830-74.17890830-62.6555515696.05985836-550.00292685-08.7124830873.06870314-61.52628359210.4607843632.7885435174.17890830-940.51486518-090.5953399773.06870314-56.98358777350.7950248283.40001943-83.40001943-620.4160784374.17890830-73.06870314-61.52628359350.7950248296.0526441471.85759988350.7950248261.5262835973.06870314-74.17890830-620.4160784383.40001943-71.8575998896.05264414e 4.5 Analytical Solution of Integration To verify stiffness matrix computed above, we have computed the stiffness matrix of element (4) in Figure 4-3B. We take a quarter of the circle and integrate [B] T [C][B]t over the area as shown in Figure 4-6. 4 y (0.5, 0.5) 122yx r r(s) (0, 0.5) x (0, 0) (0.5, 0) Figure 4-6. A quarter of circle geometry with global coordinate. Stiffness matrix of circle geometry as shown Figure 4-6 is drdstsrJBCBKTedet11)(1 For isotropic materials, the two materials properties are Youngâ€™s modulus E and Poissonâ€™s ratio . PAGE 46 33 [C]= 000202 21E , 12E Where Youngâ€™s Modulus E is 2.0710 11 and Poissonâ€™s ratio is 0.3, then = 2.2747e11, 2 =6.8242e10, and = 7.9615e10. Shape functions of four-node isoparametric quadrilateral element are srN11411 , srN11412 srN11413 , srN11414 In isoparametric formulation, we use the same shape functions N i to also express the coordinates of a point within the element in terms of nodal coordinates. Thus, )1(41)1)(1(415.0)1)(1(415.0114441rsrsrNxNxNxxiii )1(41)1)(1(415.0)1)(1(415.0221141ssrsrNyNyNyyiii (4-7) Expressing the derivatives of a function in x, y coordinates in terms of its derivatives in r, s coordinates as sNrNyNxNiiii1J , 410041sysxryrxJ , det[J] =1/16 (4-8) PAGE 47 34 sNrNsNrNyNxNiiiiii441004116 srN1411 , srN1412 , srN 1413 , srN1414 rsN1411 , rsN1412 , rsN 1413 , rsN1414 sxN11 , sxN12 , sxN 13 , sxN14 ryN11 , ryN12 , ryN13 , ryN14 [B] T [C][B] = xNyNxNyNxNyNxNyNyNyNyNyNxNxNxNxNxNyNyNxNxNyNyNxNxNyNyNxNxNyNyNxN443322114321432144443333222211110000000000020200000000 For easy calculating the matrix multiplication, is substituted by a, is substituted by b, and is substituted by c. 2 PAGE 48 35 = rsrrsscabbarssrrs111...0101010000.....101110101 Using maple program to get the result of above matrix multiplication > B:=linalg[matrix](3,8,[(1+s),0,-(1+s),0,-(1-s),0,(1-s),0,0,(1+r),0,(1-r),0,-(1-r),0,-(1+r),(1+r),(1+s),(1-r),-(1+s),-(1-r),-(1-s),-(1+r),(1-s)]); := B1s0 1s0 1s0 1s001r01r01r01r1r1s1r1s1r1s1r1s > C:=linalg[matrix](3,3,[a,b,0,b,a,0,0,0,c]); := C ab0ba000c > BT:=transpose(B); := BT ()transpose B BCB:=evalm(BT &* C &* B); PAGE 49 36 B CB := ()1s2a()1r2c()1sb()1r()1rc()1s[,,()1sa()1s()1rc()1r()1sb()1r()1rc()1s,()1sa()1s()1rc()1r()1sb()1r()1rc()1s,()1sa()1s()1rc()1r()1sb()1r()1rc()1s,()1sb()1r()1rc() ,,] 1s () 1r2a() 1s2c[,,()1rb()1s()1sc()1r()1ra()1r()1sc()1s,()1rb()1s()1sc()1r()1ra()1r()1sc()1s,()1rb()1s()1sc()1r()1ra()1r()1sc()1s,()1sa()1s()1rc()1r()1rb()1s()1sc()1r[,()1s2a()1r2c ,,], () 1sb() 1r() 1rc() 1s,,()1sa()1s()1rc()1r()1sb()1r()1rc()1s,()1sa()1s()1rc()1r()1sb()1r()1rc()1s,()1sb()1r()1rc()1s()1ra()1r()1sc()1s[,()1sb()1r()1rc() ,], 1s () 1r2a() 1s2c,,()1rb()1s()1sc()1r()1ra()1r()1sc()1s,()1rb()1s()1sc()1r()1ra()1r()1sc()1s,()1sa()1s()1rc()1r()1rb()1s()1sc()1r[,()1sa()1s()1rc()1r()1rb()1s()1sc()1r,()1s2a()1r2c ,],, () 1sb() 1r() 1rc() 1s,()1sa()1s()1rc()1r()1sb()1r()1rc()1s,()1sb()1r()1rc()1s()1ra()1r()1sc()1s[,()1sb()1r()1rc()1s()1ra()1r()1sc()1s,()1sb()1r() ,],, 1rc() 1s () 1r2a() 1s2c,,()1rb()1s()1sc()1r()1ra()1r()1sc()1s,()1sa()1s()1rc()1r()1rb()1s()1sc()1r[,()1sa()1s()1rc()1r()1rb()1s()1sc()1r,()1sa()1s()1rc()1r()1rb()1s()1sc()1r,()1s2a()1r2c ],,, () 1sb() 1r() 1rc() 1s,]()1sb()1r()1rc()1s()1ra()1r()1sc()1s[,()1sb()1r()1rc()1s()1ra()1r()1sc()1s,()1sb()1r()1rc()1s()1ra()1r()1sc()1s,()1sb()1r() ,,, 1rc() 1s () 1r2a() 1s2c,] PAGE 50 37 Arranging the result as 8 by 8 matrix form is = 222222)1()1(..........)1()1()1)(1()1)(1(...111111scrascrarscsrbsrcsrbrcsa Considering circle geometry as 4122yx Where )1(41rx , )1(41sy from Eq. 4-7. 41)1(161)1(1612222sryx 4)1()1(22sr 1)1(4)(2ssr Stiffness matrix of the first quarter of circle is drdstsrJBCBKTedet11)(1 Assuming thickness t = 1, and computing det[J] = 1/16 from Eq. 4-8. The stiffness of each matrix component can be derived as follows. drdsrcsaksr11)(1221111161 drdssrcsrbksr11)(1121111161 . PAGE 51 38 drdsscraksr11)(12288)1()1(161 Integrate above equations using maple program > int(int((1+s)^2*a+(1+r)^2*c,r=-1..sqrt(4-(1+s)^2)-1),s=-1..1)/16; .60295998761011 > int(int((1+s)*b*(1+r)+(1+r)*c*(1+s),r=-1..sqrt(4-(1+s)^2)-1),s=-1..1)/16; .18482125001011 .. > int(int((-1-s)^2*a+(1-r)^2*c,r=-1..sqrt(4-(1+s)^2)-1),s=-1..1)/16; .69748806881011 Assembling stiffness matrix with above results is [K e ] = 10974880688.6080354167.36074526348.05187599812.0206566792.4793493314.2375766543.38056208331.0080354167.3730385231.80046266475.0076551900.6658459981.2206566792.44265208333.0552733455.16074526348.00046266475.0675666050.9603740794.3076551900.65187599812.0206566792.4080354167.35187599812.0076551900.6603740794.3675666050.90046266475.06074526348.0080354167.3206566792.4206566792.4658459981.2076551900.60046266475.0730385231.8080354167.3552733455.14265208333.0793493314.2206566792.45187599812.06074526348.0080354167.3974880688.68056208331.0375766543.3375766543.34265208333.0206566792.4080354167.3552733455.18056208331.0029599877.6848212500.18056208331.0552733455.1080354167.3206566792.44265208333.0375766543.3848212500.1029599876.6e Compared stiffness data between ISEM result and Analytical solution PAGE 52 39 Table 4-2. Stiffness data FEM Result Analytical Solution Error k 11 6.052644149 6.029599876 0.381% k 12 1.857599887 1.848212500 0.505% k 13 -3.400019438 -3.375766543 0.713% k 14 0.416078436 0.426520833 2.510% k 15 -4.178908307 -4.206566792 0.662% k 16 -3.068703147 -3.080354167 0.380% k 17 1.526283596 1.552733455 1.733% k 18 0.795024823 0.805620833 1.333% PAGE 53 CHAPTER 5 GRAPHICAL DISPLAY ALGORITHM 5.1 Creating Geometry In 3D computer graphics, everything is modeled and rendered using line of polygon arrays with vertex-based data. In our implementation, Java3D is used to display graphics. Geometry is an abstract superclass in Java3D, so the referenced object is an instance of a subclass of Geometry. The Geometry subclasses may be used to specify points, lines, and filled polygons (triangles and quadrilaterals). Arrange of these vertex-based primitives can be stored in subclasses of the GeometryArray abstract class, which has arrays that maintain data per vertex. The subclasses of GeometryArray are LineArray, PointArray, QuadArray, and TriangleArray. TriangleArray is used to display the result of analysis in ISEM program. 5.2 Creating TriangleArray The triangle array is created using the concept of quadtree representation (Lee 1999). Instead of generating mesh or integrating, the quadtree method is used to create triangle array in each cell to display the geometry. As in the previous chapter, triangle array is created in both internal cells and partial cells. 5.2.1 TriangleArray for an Internal Cell In case of an internal cell, we calculate mid point of the cell, and create four triangle arrays as shown in Figure 5-1. 40 PAGE 54 41 Figure 5-1. TriangleArray for an internal cell. 5.2.2 TriangleArray for a Partial Cell Each partial cell is divided into the triangle elements based on the number of vertices inside of geometry in the quadrant. We calculate the intersection points of geometry and cell. When one vertex of the cell is inside the geometry, we calculate the mid point of the triangle defined by the intersection points and the vertex inside the geometry. Then we create three triangle arrays with those points shown in Figure 5-2A. In the same way, when two and three vertices of cell are inside the geometry, we create four and five triangles as shown in Figure 5-2B and 5-2C. A B C Figure 5-2. TriangleArray for a partial cell. A) One vertex and two intersection points. B) Two vertices and two intersection points. C) Three vertices and two intersection points. PAGE 55 42 5.2.3 Example of TriangleArray TriangleArray is created in the internal cell and the partial cell. An example of TriangleArray created for a circle is shown Figure 5-3. Figure 5-3. TriangleArray created for circle geometry. 5.3 Color Distribution 5.3.1 Smoothness After creating TriangleArray, we need to assign the proper color value to the vertices of the TriangleArray to display color. In case of plane-stress, we can display displacements, strains, and stresses as a color distribution. The displacements are continuous across the element boundaries. The element stresses or strains are calculated using derivatives of displacements which are discontinuous across the element boundaries. Values obtained at the element edge when calculated in adjacent elements may differ substantially if a coarse finite element mesh is used. The stress differences at PAGE 56 43 the element boundaries decrease as the finite element mesh is refined, and the rate at which this decrease occurs is determined by the order of the elements in the discretization. That the element stresses and strains are, in general, not continuous across element boundaries as shown in Figure 5-4. To solve this discontinuity problem, we have a smoothing algorithm to get average value of stress/strain at each node. We loop through all the elements and get the values at each node of the element. If there is more than one value at a node, we add all values and divide the sum of values by the number of values. Figure 5-5 shows that color distribution after using average value at each node. Figure 5-4. Von Mises stress before applying smoothing algorithm. PAGE 57 44 Figure 5-5. Von Mises stress after applying smoothing algorithm. 5.3.2 Compute the Value to the Intersection Point and the Mid Point We can calculate the value of displacements/stresses/strains at the intersection point and mid point. We use the shape functions of 4-node or 9-node isoparametric quadrilateral elements. As we discussed in Chapter 4, shape functions of 4-node quadrilateral elements are srN11411 , srN11412 srN11413 , srN11414 (5-1) In isoparametric formulation, we use the same shape functions N i to also express the value of a point within the element in terms of nodal values. Thus, 41iiiNvv (5-2) For example, we can calculate the value at point A in Figure 5-6 as follows PAGE 58 45 41iiAiANvv = AAAAAAAAsrvsrvsrvsrv11411141114111414321 (5-3) Also we can calculate the value at point B and C using Eq. 5-3. v1 v2 A C B v4 v3 Figure 5-6. Value of the point within the element. 5.4 Deformed Shape We have already created TrianlgeArray and assigned the proper value at vertices in each TriangleArray. To get the deformed shape, we need the displacement at each vertex. We can calculate the displacement using followed eqation. 41iiiNdd (5-4) After getting the displacement at each vertex, we add the displacements to original vertices and create the new triangle array with deformed vertices. Figure 5-7 shows the deformed shape of a circular disk due to applied concentrated load. PAGE 59 46 Figure 5-7. Deformed shape with circle geometry. PAGE 60 CHAPTER 6 ANALYSIS AND RESULTS In this chapter, some models are analyzed using both ISEM and FEA (I-DEAS, version 9). The result of each case is presented and compared including the boundary conditions. Also analytical solution is presented for some examples to compare with the results. 6.1 Circular Disk with Concentrated Load 6.1.1 Problem Statement A simple circular disk of constant thickness is considered here with loading as shown in Figure 6-1. The purpose of this simple problem is to compute displacements and stress distribution and plot the results in both ISEM and I-DEAS. For this reason, it is not required to generate high density meshes, and hence low density meshes are created in both ISEM and I-DEAS with the same number of elements. 300N Youngâ€™s Modulus E: 200 GPa Poissonâ€™s ratio : 0.29 Diameter: 1 m Thickness: 0.1 m Figure 6-1. Circular disk. 47 PAGE 61 48 6.1.2 Creating Geometry and Mesh The circular geometry is created using implicit solid element as we discussed in Chapter 3 (Kumar and Lee 2002). Figure 6-2A shows the circle geometry whose diameter is 1m created using the 9-node quadrilateral solid element for two-dimensional primitives. Figure 6-2B shows the mesh generated on the circle geometry using the 4-node quadrilateral finite elements by ISEM program. After the mesh is created, the type of analysis is specified and boundary conditions are applied on the mesh. In this example, the analysis type is plane-stress and the boundary conditions are shown in Figure 6-2C. Displacements are fixed in the all directions for the nodes at the bottom and the value of the concentrated load applied at the top is 300 N. A B C Figure 6-2. Circular geometry with mesh and boundary conditions. 6.1.3 Analysis and Results Figure 6-3, 6-5, 6-7 show the results of ISEM analysis using the new integration algorithm and graphical display algorithm presented in Chapter 4 and Chapter 5. The results of FEA analysis are shown in Fig 6-4, 6-6 6-8. PAGE 62 49 Figure 6-3. Displacement Magnitude by ISEM. Figure 6-4. Displacement Magnitude by I-DEAS. PAGE 63 50 Figure 6-5. Von Mises Stress by ISEM. Figure 6-6. Von Mises Stress by IDEAS. PAGE 64 51 Figure 6-7. Maximum Principal Stress by ISEM. Figure 6-8. Maximum Principle Stress by IDEAS. 6.1.4 Comparing the Results between ISEM and I-DEAS The analysis result in both ISEM and I-DEAS is shown in Table 6-2. NE represents the â€œNumber of Elementsâ€ that affects the analysis. The total number elements are not same as the number of elements affecting analysis in ISEM. There are the elements outside the geometry and these elements do not affect the analysis. Three finite elements at the four corners of the solid element are outside the geometry as shown in Figure 6-7. Therefore, even though the total number of elements is 144, the number of elements PAGE 65 52 affecting the analysis in ISEM is 132, which are the same as the number of elements used for analysis in I-DEAS. Table 6-1. Analysis results of circular disk. NE Displacement Mag Stress Von Mises Sress Max Prin ISEM 132 5.34E-08 4.84E+04 5.74E+04 I-DEAS 132 5.57E-08 3.58E+04 4.97E+04 As you can see from the above analysis results, the results of displacement magnitude by both ISEM and FEA are similar (less than 5% difference), but the results for stress are somewhat different because of stress concentration. However, the pattern of color distribution for both analyses is very similar as shown in Figure 6-5 ~ 6-8. 6.2 Plate with a Circular Hole Under Uniaxial Tension 6.2.1 Problem Statement Figure 6-9 shows a square plate (11m) with circular hole in the center. The thickness of the plate is 0.1 m and the diameter of the circular hole is 0.2m. The plate will be required to support a tensile load of 800N applied to each end, indicated by the arrow. Figure 6-9. Plate with a hole. 6.2.2 Creating Geometry and Mesh Figure 6-10A shows the plate with a circular hole created by two 9-node quadrilateral solid elements. A square is created, then a circle (radius is 0.2) is created. The circle is subtracted from the square. Figure 6-10B shows the mesh generated on the geometry using the 9-node quadrilateral finite elements. In this example, analysis type is PAGE 66 53 plane-stress and boundary conditions are shown in Figure 6-10C. Displacements are fixed in the x direction and free in the y direction on the left edge of the mesh, and a uniform pressure of 800 Pa is applied on right edge. A B C Figure 6-10. Plate with mesh and boundary condition. 6.2.3 Analytical Solution The concept of a stress concentration factor is used to verify the result in this section (Boresi et al. 1993). In the tension test of bar of an isotropic homogeneous material of constant cross-sectional area A, the stress is assumed to be uniformly distributed over the cross section, provided the section is far removed from the ends of the bar, where the load may be applied in a nonuniform manner. Nonuiformity of stress will occur when there are holes or notches in the cross section. This nonuniformity in stress distribution may result in a maximum stress max at a section near the hole that is considerably larger than the average stress A/ Pn , where P is the total tension load. ncS max PAGE 67 54 The ratio is called the stress concentration factor for the section. If cS max is the calculated value c of the localized stress as found from the theory of elasticity, or experimental method, becomes , called the calculated concentration factor. cS ccS Figure 6-11A shows the case of infinite plate with small circular hole of radius a under uniaxial tension . With respect to polar coordinates ( ), r , the plane stress components at any point P are given by the formulas (Boresi and Chong 1987) 2cos311212222222rarararr 2cos312122222rara 2sin31122222rarar (6-1) P r y a a x 3 A B Figure 6-11. Circular hole in an infinite plane under uniaxial tension. A) Infinite plate with a small hole. B) distribution for = 2/, 3/2 PAGE 68 55 The stress state given by Eq. 6-1 satisfies the boundary conditions at a r (0 rrr for all ) and r ( rrxx , 0 rxy for ,0 and 0rryy , 0 rxy for 2/3,2 / ). For a r , )2cos21( (6-2) attains its maximum values of 3(max) for 2/ , 2/3 , and a compressive value for ,0 . Fig 6-5B shows attains a maximum tensile value of three times the uniformly distributed stress , at the hole a r for 2/ , 2/3 . Hence the stress concentration factor = 3. However when the diameter of the hole is comparable to the width of the plate, Eq. 6-2 is considerably in error. The problem of a plate strip with a circular hole has been studied theoretically and experimentally. The resultant formula is ccS 3.013maxkkSncc (6-3) where is the ratio (width of strip/diameter of hole) and k n is the average stress over the cross-sectional area of the plate remaining at the section containing the hole. Width of strip is 1m and diameter of hole is 0.2m, so ratio,, is 5 (width of strip/diameter of hole). The weakened cross-sectional area is k 08.01.08.0 WA m 2 1000008.0800WnAP N/ m 2 From Eq. 6-3, is ccS 6415.23.5143.013maxkkSncc PAGE 69 56 Pa106415.2100006415.24maxnccS attains its maximum values of for Pa106415.24(max) /2 , 2/3 as shown in Fig 6-11A. 6.2.4 Analysis and Results The model is analyzed by ISEM and I-DEAS for stress in x-direction. The results of ISEM analysis are shown in Figure 6-12, and the results of FEA using I-DEAS program are shown in Figure 6-13. Figure 6-12. Stress X for 144 elements by ISEM. Figure 6-13. Stress X for 140 elements by I-DEAS. PAGE 70 57 6.2.5 Comparing Results The analysis result for the two different element sizes is shown in Table 6-2. Error represents difference between the results of ISEM and I-DEAS analysis and the analytical solution that is . Pa106415.24(max) Table 6-2. Analysis results of plate with hole. ISEM(Q9Node) I-DEAS(Q8Node) NEAE NE MSX Error NE MSX Error 10 100 2.52E+04 4.60% 116 2.52E+04 4.55% 12 144 2.50E+04 5.49% 140 2.50E+04 5.30% NEAE represents the â€œNumber of Elements Along the Edgeâ€. NE represents the â€œNumber of Elementsâ€ affecting the analysis. MSX represents the â€œMaximum value of Stress in X-directionâ€. As we can see the Table 6-2, the results of stress by both ISEM and FEA are very similar and the error between this results and analytical solution is around 5%. 6.3 Curved Beam 1 6.3.1 Problem Statement The curved beam in Figure 6-14 has a square cross section (0.2 m 0.2 m) and radius of curvature R = 0.4m. 6kN is loaded in its plane. This beams is made of a steel for which E=200 GPa and =0.29. PAGE 71 58 R R P = 6 kNABC Figure 6-14. Curved beam 1. 6.3.2 Creating Geometry and Mesh Figure 6-15 shows how to create the curved beam. The curved beam is created by three 9-node quadrilateral solid elements is shown in Figure 6-15A. The big circle (radius is 0.5) is created, then the small circle (radius is 0.3) is created. The small circle is subtracted from the big circle, so annular shape is created. Finally to get the curved beam as shown in Figure 6-15, we created a rectangle and subtracted it from the annular. The mesh is generated on the curved beam using 4-node quadrilateral finite elements by ISEM program. After the mesh is created, the type of analysis is specified and the boundary conditions are applied on the mesh as shown in Figure 6-15B. A B Figure 6-15. Curved beam1 with mesh and boundary conditions. PAGE 72 59 6.3.3 Analytical Solution To get the analytical solution of deflection and stresses at given point A, B, and C, Castiglianoâ€™s theorem and the equation of the circumferential stress distribution for the curved beam are used (Boresi et al. 1993). The strain energy U for a structure is equal to the sum of the strain energies of its members. The loading for the jth member of the structure is assumed to be such that the strain energy U for that member is j jTjSjMjNjUUUUU where , , , are the strain energy for axial loading, bending, shear loading, and a torsion, respectively. jNU jMU jSU jTU dzEANUN22 , dzEIMxxM22 U, dzGAkVyS22 U, dzGJTT22 U E, A, I, G, J, k represent the Youngâ€™s modulus, area, moment of inertia, shear modulus, polar moment of inertia, and the correction coefficient for the shear strain energy, respectively. With total strain energy U of the structure known, the deflection of of the structure at the location of a concentrated force in the direction of is iq iF iF iiFUq mjijjjjijjjjijjjjjijjjjFTJGTdzFMIEMdzFVAGVkdzFNAEN1 (6-4) The circumferential stress distribution for the curved beam is )()(ARAArrAAMAPmmx (6-5) PAGE 73 60 where the analytical expression for A, R and A m is shown in Figure 6-16. acbAcaRacbAmln2)( b a c Figure 6-16. Analytical expression for A, R and A m . For a cross section of the curved beam located at angle from the section on which P is applied, cosPN , cosPN sinPV , sinPV )cos1( PRM , )cos1(RPM 2/302/30)(sinsin)(coscos RdGAkPRdEAPqP 2/30)cos1()cos1( RdREIPR 2/302/30)2cos2121(2.1)2cos2121(dGAPRdEAPRqP 2/303)2cos2121cos21(dEIPR Therefore, q, the deflection at point A in Figure 6-14 is P EIPRGAPREAPRqP32494)3(2.143 PAGE 74 61 433223)200)(10200()12()400)(6000(249)200)(77500(4)400)(6000)(3(2.1)200)(10200(4)400)(6000(3 130588.0002189.0100686.74 mm m 13348.0 41033481. The circumferential stress distribution for the curved beam is computed by using Eq. 6-5 and equations in Figure 6-16. 04.02.02.0A m 2 102165.03.05.0ln2.0lniomrrbA m 4.023.05.02iorrR m B is the stress at point B in Figure 6-14. )()(2ARAArArAPRAPmBmBB )04.0102165.04.0)(3.0)(04.0()102165.03.004.0)(24.06000(04.06000 44689364318936150000 Pa Pa 6104689.4 C is the stress at point C in Figure 6-14. )()(2ARAArArAPRAPmCmCC )04.0102165.04.0)(5.0)(04.0()102165.05.004.0)(24.06000(04.06000 29213633071363150000 Pa Pa 6109214.2 PAGE 75 62 6.3.4 Analysis and Results The model is analyzed by ISEM for displacement and stress in y-direction. The same model is analyzed by I-DEAS for displacement and stress in y-direction such that the total number of elements is as close to the number of elements that is affected by the analysis in ISEM. The results of ISEM analysis are shown in Figure 6-17, 6-19 and the results of FEA using I-DEAS program are shown in Figure 6-18, 6-20. Figure 6-17. Displacement Y by ISEM. Figure 6-18. Displacement Y by I-DEAS. PAGE 76 63 Figure 6-19. Stress Y by ISEM. Figure 6-20. Stress Y by I-DEAS. 6.3.5 Comparing Results Points A, B, and C are the points at which we check the values of displacement and stresses as shown in Figure 6-14. The analytical solution is calculated in the previous section, and Error represents difference between the analytical solutions and the results of ISEM and I-DEAS analysis. PAGE 77 64 Table 6-3. Analysis results of curved beam 1. Displacement Stress NE A Error B Error C Error 192 1.23E-04 7.93% -4.15E+06 7.25% 2.92E+06 0.09% ISEM 408 1.28E-04 4.23% -4.44E+06 0.56% 2.98E+06 2.07% 206 1.30E-04 2.46% -4.33E+06 3.11% 2.92E+06 0.09% I-DEAS 395 1.30E-04 2.68% -4.49E+06 0.47% 2.91E+06 0.39% Analytical Solution 1.33E-04 -4.47E+06 2.92E+06 As you can see Table 6-3, in case of low density meshes, the ISEM results are different from the I-DEAS results and analytical solution for displacement at point A and stress at point B. Other values for both analyses are similar to each other and analytical solution, but it is difficult for us to see the stability of the analysis results from the above. To check the stability, the number of element will be gradually increased in both ISEM and I-DEAS, and the results will be compared with the analytical solution. However we need to create new example to check stability because we can not give proper boundary conditions with increasing the number of elements. Currently we only can give boundary conditions to geometry that conforms to the mesh. An appropriate example will be presented in the next section. 6.4 Curved Beam 2 6.4.1 Problem Statement The curved beam in Figure 6-21 has a square cross section (0.3m 0.3m) and the inside diameter of curved beam is 0.4m. Load P is 20kN. This beams is made of a steel for which E=200 GPa and =0.29. PAGE 78 65 P = 20kNBCA0.4 m Figure 6-21. Curved beam 2. 6.4.2 Creating Geometry and Mesh The curved beam is created by three 9-node quadrilateral solid elements as shown in Figure 6-22A. The big circle (radius is 0.5) is created, then the small circle (radius is 0.2) is created. The small circle is subtracted from the big circle, so annular shape is created. Finally to get the curved beam as shown in Figure 6-21, we created a rectangle and subtracted the rectangle from the annular. The mesh is generated on the curved beam using both the 4-node and 9-node quadrilateral finite element by ISEM program. After mesh is created, the type of analysis is specified and the boundary conditions are applied on the mesh as shown in Figure 6-22B. PAGE 79 66 A B Figure 6-22. Curved beam 2 with mesh and boundary conditions. 6.4.3 Analytical Solution For a cross section of the curved beam located at angle from the section on which P is applied, sinPN , sinPN cosPV , cosPV sinPRM , sinRPM 000)(sinsin)(coscos)(sinsinRdREIPRRdGAkPRdEAPqP 00)2cos2121(2.1)2cos2121(dGAPRdEAPRqP 03)2cos2121(dEIPR Therefore, q, the deflection at point A in Figure 6-21 is P EIPRGAPREAPRqP22)(2.123 PAGE 80 67 433223)300)(10200(2)12()350)(20000()300)(77500(2)350)(20000)((2.1)300)(10200(2)350)(20000( 009977.0001892.0101086.64 mm =1 m 01248.0 510248. The circumferential stress distribution for the curved beam is 09.03.03.0bhA m 2 2749.02.05.0ln3.0lniomrrbA m 35.022.05.02iorrR m B is the stress at point B in Figure 6-21. )()(ARAArArAPRAPmBmBB )09.02749.035.0)(2.0)(09.0()2749.02.009.0)(35.020000(09.020000 24135152191293222222 Pa Pa 6104136.2 C is the stress at point C in Figure 6.12. )()(ARAArArAPRAPmCmCC )09.02749.035.0)(5.0)(09.0()2749.05.009.0)(35.020000(09.020000 9654061187628222222 Pa Pa 5106541.9 PAGE 81 68 6.4.4 Analysis and Results The model is analyzed by both ISEM and I-DEAS for displacement and stress in y-direction with similar number of elements. Also the same model is analyzed using two different element types with the same number of elements. Figure 6-23 shows the displacement results obtained quadrilateral 4-node (Q4Node) and 9-node (Q9Node) finite elements in ISEM and the number of elements that affects analysis is 284. Figure 6-24 shows the displacement for quadrilateral 4-node and 8-node (Q8Node) finite elements in I-DEAS and the number of elements that affects analysis is 290. Figure 6-25 shows the stress computed by ISEM and Figure 6-26 shows the stress computed by I-DEAS. A B Figure 6-23. Displacement Y by ISEM. A) Q4Node. B) Q9Node. A B Figure 6-24. Displacement Y by I-DEAS. A) Q4Node. B) Q8Node. PAGE 82 69 A B Figure 6-25. Stress Y by ISEM. A) Q4Node. B) Q9Node. A B Figure 6-26. Stress Y by I-DEAS. A) Q4Node. B) Q8Node. 6.4.5 Comparing Results In Table 6-4, NE represents the â€œNumber of Elementâ€ that affects the analysis. Points A, B, and C are the points at which we check the values of displacement and stresses as shown in Figure 6-21. Error represents difference between the analytical solutions and the results of ISEM and I-DEAS analysis. Table 6-4. Analysis results of curved beam 2. Displacement Stress Element Type NE A Error B Error C Error Q4 284 -1.13E-05 9.72% -2.29E+06 4.92% 9.88E+05 2.32% ISEM Q9 284 -1.20E-05 3.89% -2.38E+06 1.21% 1.08E+06 9.73% Q4 291 -1.12E-05 10.58% -2.32E+06 3.92% 9.58E+05 0.78% I-DEAS Q9 291 -1.16E-05 7.13% -2.43E+06 0.60% 9.92E+05 2.72% Analytical Solution -1.25E-05 -2.41E+06 9.65E+05 PAGE 83 70 As you can see in Table 6-4, in case of the displacement result, the analysis result using quad 9-node and quad 8-node elements have less error than the analysis results using quad 4-node elements. And the analysis result obtained by ISEM has less error than that obtained by I-DEAS for the displacements. The results of the stress at point B for both analyses are very similar to each other. However the results of the stress at point C have large error for the analysis results obtained by quad 9-node elements. So the number of element is gradually increased from 284 to 1028 in case of using quad 4-node elements and from 284 to 600 in case of using quad 9-node elements in ISEM analysis. Also in I-DEAS analysis, the similar number of elements is gradually increased and the results are compared below. Q4Node-1.30E-05-1.25E-05-1.20E-05-1.15E-05-1.10E-05-1.05E-05-1.00E-051234567891011 AS ISEM IDEAS Q9Node(ISEM), Q8Node(I-DEAS)-1.26E-05-1.24E-05-1.22E-05-1.20E-05-1.18E-05-1.16E-05-1.14E-05-1.12E-05-1.10E-05123456 AS ISEM IDEAS Figure 6-27. Displacement at point A. Q4Node-2.55E+06-2.50E+06-2.45E+06-2.40E+06-2.35E+06-2.30E+06-2.25E+06-2.20E+06-2.15E+061234567891011 AS ISEM IDEAS Q9Node(ISEM), Q8Node(I-DEAS)-2.60E+06-2.55E+06-2.50E+06-2.45E+06-2.40E+06-2.35E+06-2.30E+06-2.25E+06-2.20E+06-2.15E+06123456 AS ISEM IDEAS Figure 6-28. Stress at point B. PAGE 84 71 Q4Node9.30E+059.40E+059.50E+059.60E+059.70E+059.80E+059.90E+051.00E+061.01E+061.02E+061234567891011 AS ISEM IDEAS Q9Node(ISEM), Q8Node(I-DEAS)9.00E+059.20E+059.40E+059.60E+059.80E+051.00E+061.02E+061.04E+061.06E+061.08E+06123456 AS ISEM IDEAS Figure 6-29. Stress at point C. Figure 6-27 shows the errors of analysis result for the displacement gradually decrease with increasing the number of elements for both element types. And the analysis results obtained by ISEM have less error than that obtained by I-DEAS. Figure 6-28 shows the stress values at point B fluctuate with increasing number of elements. At point C as you can in Figure 6-29, the stress values are stable with increasing number of elements for both analysis types and element types and analysis result using quad 9-node and quad 8-node elements have more error than the analysis results using quad 4-node elements. As we mentioned before, error represents difference between the analytical solutions and the results of ISEM and I-DEAS analysis. The equation of the circumferential stress distribution for the curved beam is used to compute stress for analytical solution. However this solution is based on the assumption that the radial stress and shear stress are sufficiently small so that the state of stress is essentially one-dimensional. Therefore analytical solution itself has an error based on this assumption. Because of this reason, even though the analysis result using quad 9-node and quad 8-node elements are not closer to analytical results than results from quad 4-node elements, we can not say the result using quad 9-node is worse than the result using quad 4-node. PAGE 85 CHAPTER 7 CONCLUSIONS AND FUTURE WORK The current work mostly focuses on the development of algorithms for 2D plane stress analysis using implicit solid element method. The following section deals with the conclusions made from the previous chapters. 7.1 Summary and Discussion In this thesis, integration and graphical display algorithms for nonconforming mesh using implicit solid element are presented. Unlike traditional Finite Element Method, a regular grid that is not required to conform to the geometry is used for analysis in Implicit Solid Element Method. The analysts spend the majority of their time in creating the mesh in FEM. The advantage using ISEM for engineering analysis is that the time for generating mesh can be considerably reduced. The analysts do not need to spend the majority of their time in creating the mesh using ISEM for engineering analysis because the mesh is a regular grid placed over the geometry. In the current implementation, there are limitations on the types of boundary conditions that can be applied because the mesh is non-conforming in ISEM. For general application of ISEM to engineering analysis, it is necessary to develop and implement the algorithms for applying these boundary conditions. Literature survey was carried out for the integration algorithms and the graphical display algorithms. Algorithms suitable for the ISEM were implemented as described in Chapter 4 and 5. These algorithms were implemented and added to an object-oriented FEM software (JavaFemViewer). Basic Pre processing capabilities were also added to the 72 PAGE 86 73 software. This software was tested using the various verification problems for 2D plane-stress analysis and the results obtained were compared with the results by I-DEAS (FEA) and the analytical solutions in Chapter 6. The following section suggests some future work for improvement of ISEM process. 7.2 Future Work 1. The mesh does not conform to the boundaries of the geometry in ISEM. Due to this reason, there are limitations on applying boundary conditions. Penalty method can be employed to apply boundary conditions on a regular non-conforming mesh that completely overlaps the geometry boundaries (Clark and Anderson Proceedings of DETC). 2. Implicit Solid Element formulation was implemented only for 2D plane stress analysis. This can be extended to 3D analysis and can be implemented in the ISEM software. Traditional methods for generating 3D mesh in finite element analysis significantly increase the time needed for engineering analysis because the mesh is required to conform to the CAD geometry. 3D ISEM can reduce the time to create complex 3D mesh. Instead of the quadtree representation used in 2D ISEM, octree representation scheme can be used in 3D ISEM (Lee 1999). PAGE 87 APPENDIX A INPUT FORMAT FOR JAVAFEMVIEWER SOFTWARE Total number of nodes Node number x y z . Number of element group Interpolation type total number of element Element number connectivity mpID gcID . Number of solution Name of solution type Solver Data SkylineSolver : no data NewMarksSolver : line1-mass density, damping coefficient, alpha, delta, time step, total time line2state, value, value (initial condition for displacement) line3state, value, value (initial condition for velocity) line4state, value, value (initial condition for acceleration) (if state is 0, then has constant value state has number, then the number represents the previous solver-no value ) EigenSolver : required number of eigen vector 74 PAGE 88 75 Number of A_Type Name of A_Type Number of associated Egroup EGroupID Number of material property set Material property set Axisymmetric : Young's Modulus, Poisson's Ratio Plane_stress : Young's Modulus, Poisson's Ratio ThreeD_stress : Young's Modulus, Poisson's Ratio, Thermal Coefficient Beam : Young's Modulus Htransf : Thermal Conductivity Number of geometric constant set Geometric constant set Axisymmetric : Thickness Plane_stress : Thickness ThreeD_stress : Thickness Beam : Moment of Inertia, C/s Area of Element, C/s height Htransf : Thickness Number of load type Load type number{1:Concentrated load, 2:Uniform load, 3:Linear load, 4:Thermal load} 1. Concentrated load Number of loaded node Loaded node number load(ndof) . PAGE 89 76 2. Uniform load Number of loaded element Loaded element number element type s, t value load(x,y,z) moment(x,y,z) . 3. Linear load Number of loaded element Loaded element number element type s, t value load(x,y,z) moment(x,y,z) load(x,y,z) moment(x,y,z) . 4. Thermal load Reference temperature Present temperature Number of constraint node Node number constraint(ndof) PAGE 90 APPENDIX B SCENEGRAPH FOR JAVAFEMVIEWER SOFTWARE Figure B-1. Scenegraph for the JavaFemViewer software This figure shows the Scenegraph that is implemented in the ViewCanvas class. As shown in figure, the meshShape Shape3D contains the Geometry node for which the finite element mesh is passed. 77 PAGE 91 LIST OF REFERENCES Spatial Technology Inc., 1995, Geometry Modeler Application Guide, Spatial Technology Inc., Boulder, CO. Atluri SN and Zhu T, 1998, A New Meshless Local Petrov-Galerkin (MLPG) Approach in Computational Mechanics, Computational Mechanics, 22, pp.117-127. Bajaj C, Chandrajit L, Chen J and Xu G, 1995, Modeling with Cubic A-Patches, ACM Transactions on Graphics, 14(2), pp.103-133. Belytschko T, Lu YY, and Gu L, 1994, Element-Free Galerkin Methods, Int. J. Numer. Methods Eng., 37, pp. 229-256. Bloomenthal J, Bajaj C, Blinn J, Cani-Gascuel MP, Rockwood A, Wyvill B and Wywill G, 1997, Introduction to Implicit Surfaces, Morgan Kaufmann Publisher Inc, San Francisco, CA. Boresi AP and Chong KP, 1987, Elasticity in Engineering Mechanics, Elsevier, NY. Boresi AP, Schmidt RJ and Sidebottom OM, 1993, Advanced Mechanics of Materials Fifth Edition, John Wiley and Sons, NY. Clark BW and Anderson DC, Proceedings of DETC, Finite Element Analysis in 3D using the Penalty Boundary Method. Glassener A, 1989, An Introduction to Ray Tracing, Academic Press, London. Hoffmann C, 1987, Geometric and Solid Modeling-An Introduction, Morgan Kaufmann Publishers, San Mateo, CA. Kumar AV and Lee J, 2002, Solid Modeling using Implicit Solid Elements, Submitted to Journal of Computing and Information Science in Engineering Transactions of the ASME, University of Florida. Lee K, 1999, Principles of CAD/CAM/CAE System, Addison Wesley, Readings, MA. Liu GR, 2002, Mesh Free Methods, CRC Press, NY. Lorensen W and Cline H, 1987, Marching Cubes-A High Resolution 3D Surface Construction Algorithm, Computer Graphics, 21(4), pp.163-69. 78 PAGE 92 79 Mantyla M, 1988, An Introduction of Solid Modeling, Computer Science Press, Rockville, MD. Mortenson ME, 1997, Geometric Modeling Second Edition, Wiley Computer Publishing, John Wiley and Sons, NY. Onate E, Idelsohn OC, Aienkiewicz RL and Taylor, 1996, A Finite Point Method in Computational Mechanics Application to Convective Transport and Fluid Flow, Int. J. Numer. Methods Eng., 39, pp.3839-3866. Requicha, AAG and Voelcker HB, 1985, Boolean Operations in Solid Modeling: Boundary Evaluation and Merging Algorithms, Proceedings of the IEEE, 73(2), pp. 30-40. Tirupathi RC and Ashok DB, 1997, Introduction of Finite Elements in Engineering Second Edtion, Prentice-Hall, Upper Saddle River, NJ. Zhao HK, Osher S, Merriman B and Kang M, 2000, Implicit and NonParametric Shape Reconstruction from Unorganized Data using a Variational Level Set Method, Computer Vision and Image Understanding, 80, pp.295-314. PAGE 93 BIOGRAPHICAL SKETCH The author was born on June 10, 1971, in Gwangju, Korea. In 1998, he graduated with a Bachelor of Science degree in mechanical engineering from Chosun University, Gwangju, Korea. He entered the Master of Science program at the University of Florida in fall 2000. 80 |