Citation |

- Permanent Link:
- http://ufdc.ufl.edu/UFE0021901/00001
## Material Information- Title:
- Finite Element Analysis Using Uniform B-Spline Approximation and Implicit Boundary Method
- Creator:
- Burla, Ravi Kumar
- Place of Publication:
- [Gainesville, Fla.]
- Publisher:
- University of Florida
- Publication Date:
- 2008
- Language:
- english
- Physical Description:
- 1 online resource (175 p.)
## Thesis/Dissertation Information- Degree:
- Doctorate ( Ph.D.)
- Degree Grantor:
- University of Florida
- Degree Disciplines:
- Mechanical Engineering
Mechanical and Aerospace Engineering - Committee Chair:
- Kumar, Ashok V.
- Committee Members:
- Sankar, Bhavani V.
Kim, Nam Ho Peters, Jorg - Graduation Date:
- 5/1/2008
## Subjects- Subjects / Keywords:
- Approximation ( jstor )
Boundary conditions ( jstor ) Boundary value problems ( jstor ) Contour lines ( jstor ) Mathematical extrapolation ( jstor ) Matrices ( jstor ) Modeling ( jstor ) Shear stress ( jstor ) Step functions ( jstor ) Stiffness matrix ( jstor ) Mechanical and Aerospace Engineering -- Dissertations, Academic -- UF boundary, bsplines, c1, c2, conforming, continuity, fem, grid, implicit, mesh, non, structured, xfem - Genre:
- Electronic Thesis or Dissertation
bibliography ( marcgt ) theses ( marcgt ) Mechanical Engineering thesis, Ph.D.
## Notes- Abstract:
- The finite element method (FEM) is a widely used numerical technique for solution of boundary value problems arising in engineering applications. FEM is well studied and has several advantages; however the stress and strain solutions obtained by FEM are discontinuous across element boundaries and require smoothing algorithms to make the stress/strain appear smooth in the analysis domain. Automatic mesh generation is still a difficult task in FEM for complicated geometries and often requires user intervention for generating quality meshes. Finite elements based on uniform B-spline basis functions defined on a structured grid are developed which provide at least continuous solution through out the analysis domain. Numerical techniques are presented in this work based on structured grids and implicit equations of the boundaries of analysis domain, material interface boundaries or coarse/fine grid interfaces for solution of various engineering boundary value problems. Solution structures are constructed using approximate Heaviside step functions for imposition of Dirichlet boundary conditions, treatment of material discontinuity to perform micromechanical analysis and for local refinement of grid. The use of structured grid eliminates the need for constructing a conforming mesh and results in significant savings of time in pre-processing stage of the design cycle. Numerical examples are presented to demonstrate the performance of B-spline elements. The results are compared with analytical solutions as well as traditional finite element solutions to demonstrate the ability of B-spline elements to represent continuous stress and strain through out analysis domain. Convergence studies show that B-spline elements can provide accurate solutions for many engineering problems with fewer numbers of elements and nodes as compared to traditional FEM. Solution structure for treatment of material boundary is validated by performing a convergence analysis on a problem involving circular inclusion in a square matrix and by determining effective properties of fiber reinforced composite. Solution structure for local grid refinement is validated by analyzing classical stress concentration problems. ( en )
- General Note:
- In the series University of Florida Digital Collections.
- General Note:
- Includes vita.
- Bibliography:
- Includes bibliographical references.
- Source of Description:
- Description based on online resource; title from PDF title page.
- Source of Description:
- This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
- Thesis:
- Thesis (Ph.D.)--University of Florida, 2008.
- Local:
- Adviser: Kumar, Ashok V.
- Statement of Responsibility:
- by Ravi Kumar Burla
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright by Ravi Kumar Burla. 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/11/2008
- Classification:
- LD1780 2008 ( lcc )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

In this expression, the terms involving only inclusion functions and its gradients are separated from the terms containing the shape functions and their gradients. This decomposition leads to efficient and accurate computation of the stiffness matrix. The volume integral has been decomposed into a line (surface) integral along the normal to the boundary and a surface integral along the boundary. The normal for the boundary represented by implicit equation 0 = 0 is given by V at the boundary. The component involving only inclusion function and its gradients is termed as C2 in Eq. (5-32) and is evaluated in the similar manner as explained in section 5.1.3. 5.3.4 Periodic Boundary Conditions Composite materials have periodic arrangement of cells and their behavior can be predicted by analyzing the behavior of the unit cell which is often termed as representative volume element (RVE). The variations in the micro-scale (RVE) are rapid and periodic as compared to the variations at the macro-scale. In order to analyze the behavior at the micro- scale, the RVE is analyzed with appropriate loading and boundary conditions. The RVEs tile the computational domain by translation and the neighboring cells fit into each other in both deformed and un-deformed states as shown in Figure 5-6. Hence the boundary conditions for the RVE are periodic in order to preserve the continuity of displacements, strains and stress across each RVE. The periodic boundary conditions are expressed as linear constraints and they are implemented as multipoint constraints in a typical finite element code. The periodic boundary conditions are expressed in the following form: uF -uF = const or in general uF -AuF + u= 0 (5-33a) uF -uF = const or in general uF -AuFi + p = 0 (5-33b) LIST OF TABLES Table page 5-1 Set of periodic boundary conditions for 6 different strains .................. ............... 102 6-1 Average stresses for various loading cases................................... ...............162 6-2 Comparison of effective properties................................................ 163 6-3 Average stresses for various loading cases................................... ...............164 N2 1)+,(r,s) =N(r)N (s); i = 1..4;j = 1..4 (4-10) NAT1k 1)+4(J )+ ,(r,s,t)=N,(r)NJ(s)Nk(t); i = 1..4;j = 1..4;k = 1..4 (4-11) The geometry mapping between parametric and physical space can be defined using Eq. (4-7). The lower bounding node for this element is node 6 and the upper bounding point corresponds to the node 11. Again, some of the support nodes lie outside the domain of the element and are shared with neighboring elements. The 2-D bi-cubic B-spline element has 16 basis functions (Eq. 4-10) while the 3D cubic B-spline has 64 basis functions (Eq. 4-11). In these equations, N (r) represents the cubic basis functions as defined in Eq. (4-5) *4L2 a 4a t Figure 6-22. Analysis model for a square plate with a circular hole IlMeo Figure 6-23. Contour plot of o-using FEM Q4 elements ^ ^, ^ ^ ^^ ^ ^^ ^ ^ , .111U I .mB3E LIST OF REFERENCES Adams, D.F., Crane, D.A., 1984. Finite element micromechanical analysis of a unidirectional composite including longitudinal shear loading. Computers and Structures 6 1153-1165. Amestoy P.R., Davis T.A., Duff I.S., 1996. An approximate minimum degree ordering algorithm. SIAM Journal on Matrix Analysis and Applications 17(4) 886-905. Atluri S.N., Kim H.G., Cho J.Y., 1999. Critical assessment of the truly Meshless Local Petrov- Galerkin (MLPG), and Local Boundary Integral Equation (LBIE) methods. Computational Mechanics 24(5) 348-372. Atluri S.N., Zhu T.L., 2000a. Meshless local Petrov-Galerkin (MLPG) approach for solving problems in elasto-statics. Computational Mechanics 25(2-3) 169-179. Atluri S.N., Zhu T., 2000b. New concepts in meshless methods. International Journal for Numerical Methods in Engineering 47(1) 537-556. Babuska I., Banerjee U., Osborn, J.E., 2002. On principles for the selection of shape functions for the Generalized Finite Element Method, Computer Methods in Applied Mechanics and Engineering 191(49-50) 5595-5629 Belytschko T., Krongauz Y., Fleming M., Organ D., Snm L., Wing K., 1996. Smoothing and accelerated computations in the element free Galerkin method. Journal of Computational and Applied Mathematics. 74(5) 111-126. Belytschko T., Krongauz Y., Organ D., Fleming M.,Krysl P., 1996. Meshless methods: an overview and recent developments. Computer Methods in Applied Mechanics and Engineering 139(1-4) 3-47. Belytschko T., Lu Y.Y., Gu L., 1994. Element-free Galerkin methods. International Journal for Numerical Methods in Engineering 37(2) 229-256. Belytschko T., Parimi C., Moes N., Usui S., Sukumar N., 2003. Structured extended finite element methods for solids defined by implicit surfaces. International Journal for Numerical Methods in Engineering 56(4) 609-635. Boresi A.P., Schmidt R.J., Sidebottom O.M., 1993. Advanced mechanics of materials. John Wiley and Sons Inc. New York. Cai, Y.C., Zhu, H.H., 2004. Direct imposition of essential boundary conditions and treatment of material discontinuities in the EFG method. Computational Mechanics 34(4) 330-338. Chamis, C.C., 1984. Simplified composite micromechanics equations for hygral, thermal and mechanical properties. SAMPE Quarterly 14-23. Clark B.W., Anderson D.C., 2003a. The penalty boundary method. Finite Elements in Analysis and Design 39(5-6) 387-401. u = Dug + " a / / \. / D - / \ te i e Uo eUI e5 Figure 3-1. Solution structure for imposing essential boundary conditions n ot n=T u= Figure 3-2. Elastic solid under equilibrium aHjc aH' 0 Hinc Ng2 0 Hnc 0 AHJnC AnCF Ng2 0 B22= 0 Nn. = 0 0 = H B22 (5-31) 2 2N2 0 Ng OHawnc 2 ti Using the above procedure, the strain-displacement matrix can be rewritten as B2 = HncB21 + H2B2 Similarly, the first component of strain-displacement matrix is written as B, =(1- HincB,) H B These matrices can be defined for 3D elements also in a similar manner. The last three components in stiffness matrix as shown in Eq. (5-29) contain gradients of inclusion functions which are large near the inclusion boundary when 6 is made very small and they vanish away from the boundary. Therefore the terms in these integrals will be non-zero only in a thin band along the boundary. This allows the volume integral to be converted to a surface integral. The contribution from the gradients of the shape functions can be assumed to be constant across the band but varies along the boundary and the contribution from the gradients of inclusion function is constant along the boundary but varies normal to the boundary. The band is represented in terms of two parameters t,n as shown in the Figure 5-3. The stiffness matrix as decomposed in Eq. (5-29) has four components; the second component JBfCB12dV can be E rewritten as shown below. J BTCB12dV= f BT ((i-Hn)CI2 )12dV SB (5-32) ~- 1 -r --r = B1T 1-H n,)C [-2 do B12dF f Bl1CB212dF F 110VF D, ON + N, O Oa, a, [B]3,2m 0 D D N + N (5-5) ax, 2 2 aN 8x 8x 8xD D, + N D, ND OD 2 ax2 ax ax ax The above equation shows that the strain-displacement matrix has gradients of both shape functions and Dirichlet functions. It is convenient to decompose the [B] matrix into two parts [B] = [B] + [B2] such that [B,] contains terms involving the gradient of the shape functions and [B2] contains terms with the gradients of the D function as shown below: N1 0 D, D1 0 0 0 l D 9X2 DvI B12 (5-6) [B1],, = 0 D2 =0 D,2 0 0 =[D B] (5-6) DX2 0 0 D, D2 9 0 N, 0 0 .D 0 a 2 al 0 9NAx ax, 1 xo, [B, 32 0 N, 2 0 N = 2 0 [D2 B2l (5-7) x, 1x 0 N, 3x2 2x2m [ND ]Nm = o D OD N D aD, a2, a & Dx, x^, 8x, Using the above procedure, the strain-displacement matrix can be rewritten as [B]= [DI] [B i + D2 ][ Similar matrices can be defined for three dimensional elements also. Using the decomposition of strain displacement matrix, the stiffness matrix is decomposed into the following components: solution structure demonstrated in Figure 3-1, D represents the Dirichlet function which is zero at the first node of element el and reaches unity within this element and remains unity through out the domain. The boundary value function is defined by ua which takes on the value of the prescribed boundary value of u, at the boundary and ug represents the grid variable which is interpolated on the grid consisting of elements el-e5. The approximate solution is given by u = Dug + ua and it can be observed that the approximated solution satisfies the boundary condition exactly. The solution structure for linear elastostatic and steady state heat transfer problems is first discussed in this section for a general D-function and a detail discussion on construction D-functions is provided later in this chapter. 3.1.1 Modified Weak Form for Linear Elasticity The governing differential equation for an elastic body under equilibrium as shown in Figure 3-2 is given as: Solve: V.b +b = 0 in Q Subject to: u = u0 on F, (3-2) a.n=to on rF In the above equation, Q is the analysis domain bounded by essential boundary F, and natural boundary Ft. {(} is Cauchy stress tensor and {b} is the body force and the Dirichlet boundary conditions are displacements specified as uo on essential boundary and the natural boundary conditions are specified traction to. The principle of virtual work is the weak form for the boundary value problem and is given as: fJ{} T{} dQ = Su}T {t}dF+ + Ju}T {b}dQ (3-3) Q F, Q Analytical FE Q4 FE Q9 CBS SQBS -100 Figure 6-32. 03 0 32 0 34 0 36 0 38 04 0 42 0 44 0 46 0 48 05 Distance from center of cylinder, in Comparison of radial stresses along the thickness of the cylinder Analytical FE Q4 FE Q9 CBS SQBS S900 ---- 600 800 I- 700 600 03 04 05 Distance from center of cylinder, in Figure 6-33. Comparison of tangential stresses along the thickness of the cylinder Aspect Ratio = L/t Figure 6-6. Thin cantilever beam subjected to end loading 10-- --- % error in CBS deflection -e- % error in Euler-beam deflection 15 Aspect ratio of beam Figure 6-7. Plot of error in the tip deflection with varying aspect ratio 124 --------- -S------------------- Grid Boundary Solid Boundary SD=D,=0 Consistent Nodal Loads Figure 6-15. Typical non-conforming grid using cubic B-spline elements 10-5 10 101 10 10 Log (Num Nodes) Figure 6-16. Convergence plot of L2 error for cantilever beam CBS QBS FEQ4 FEQ9 Kumar A.V., Lee J., 2006. Step function representation of solid models and application to mesh free engineering analysis. Journal of Mechanical Design. Transactions of the ASME 128(1) 46-56. Kumar A.V., Lee J.H., Burla R.K., 2005. Implicit solid modeling for mesh free analysis. In: Proceedings of International Design Engineering Technical Conferences & Computers in Engineering IDETC/CIE Long Beach, California. Leung A.Y.T., Au F.T.K., 1990. Spline finite elements for beam and plate. Computers and Structures 37(5) 717-729. Li, S., Lim, S.H., 2005. Variational principles for generalized plane strain problems and their applications. Applied Science and Manufacturing 36(3) 353-365. Liszka T.J., Duarte C.A.M., Tworzydlo W.W., 1996. hp-meshless cloud method. Computer Methods in Applied Mechanics and Engineering 139(1-4) 263-288. Liu G.R., Gu Y.T., 2001. A point interpolation method for two-dimensional solids. International Journal for Numerical Methods in Engineering 50(4) 937-951. Lucy L.B., 1977. A numerical approach to the testing of the fission hypothesis. The Astronomical Journal 82(12) 1013-24. Marrey, R.V., Sankar, B.V., 1997. A Micromechanical Model for Textile Composite Plates. Journal of Composite Materials 31(12) 1187-1213. Melenk J.M., Babuska I., 1996. The partition of unity finite element method: Basic theory and applications. Computer Methods in Applied Mechanics and Engineering 139(1-4) 289- 314. Moes N., Bechet E., Tourbier M., 2006. Imposing Dirichlet boundary conditions in the extended finite element method. International Journal for Numerical Methods in Engineering 67(12) 1641-1669. Monaghan J.J., 1987. An introduction to SPH. Computer Physics Communications 48(1) 89 96. Mortenson M.E., 1997. Geometric Modeling. Wiley, New York. Nayroles B., Touzat G., Villon P.,1992. Generalizing the finite element method: Diffuse approximation and diffuse elements. Computational Mechanics 10(5) 307-318. Osher S.J., Fedkiw R.P., 2002. Level Set Methods and Dynamic Implicit Surfaces. Springer Verlag, Philadelphia. Padmanabhan, S., 2006. Analysis using non-conforming structured grid and implicit boundary method. Doctoral Dissertation, University of Florida. Stress Y 2.9443E6 2.2016E6 1.4589E6 7.1621E5 -2.6487E4 ,-7.6918E5 1.5119E6 -2.2546E6 S-2.9973E6 SA I -3.74E6 S-4.4827E6 Figure 6-13. Contour plot of o-, using cubic B-spline elements Y B D X P A L Figure 6-14. Thick cantilever beam used for convergence analysis performed by computing the L, errors and energy errors for a series of models with varying grid densities. For comparison purposes, errors are also computed for a series of models with FE Q4 (bi-linear) and FE Q9 (bi-quadratic) elements. The plot for L, error is shown in Figure 6-25 and plot for energy error is shown in Figure 6-26. Cubic B-spline elements perform much better than the quadratic B-spline elements and FE elements. For a given number of nodes, the errors for B- splines are far smaller when compared to the traditional finite elements. Figure 6-27 shows comparison with the analytical solution along the line 0 = It can be 2 observed that cubic B-spline element computes the stress concentration factor as 3.06 which is very close to the theoretical value of 3.0. 6.1.5 Thick Cylinder with Internal Pressure A hollow thick walled cylinder of internal radius 0.3in and external radius 0.5in is modeled as shown in Figure 6-28 and for analysis purpose only a quarter of the cylinder is used. The Dirichlet boundary conditions are applied on the straight edges of the quarter such that radial displacement is allowed while tangential displacement is restricted. The inner curved edge is subjected to a pressure of 500 psi. The material is assumed to be made of steel with a modulus of elasticity as 30 x 106 psi and a Poisson's ratio of 0.3. The analytical solution is obtained from the governing differential equations for displacements, strains and stresses in a thick walled cylinder with open ends (Boresi et.al. 2003). The displacements and stresses as a function of radial position r is: (r ((1v)+a22 (6-6a) ur = 2 2 2- (6-6a) E b 2-a 2) 6-40 Comparison of temperature along line AB ................................... .......... ................. 143 6-41 Three dimensional clamp used for analysis................................ ....................... 143 6-42 Convergence plot for strain energy for clamp model ............... ................... ...........144 6-43 Contour plot of bending stresses using H8 elements in FEM for 3D clamp ...................144 6-44 Contour plot of bending stresses using H8 elements for 3D clamp ................................. 145 6-45 Contour plot of bending stresses using cubic B-spline elements for 3D clamp ..............145 6-46 Three dimensional rib model used for analysis .................................... ............... 146 6-47 Convergence plot for strain energy for rib model .................. ................................146 6-48 Contour plot of bending stresses using tetrahedral elements for 3D rib ........................147 6-49 Contour plot of bending stresses using B-spline elements for 3D rib .............................147 6-50 Contour plot of bending stresses using H8 elements for 3D rib.................................. 148 6-51 Three dimensional dial bracket used for analysis .................. ................................148 6-52 Convergence plot for strain energy for dial bracket model ..........................................149 6-53 Contour plot of bending stresses using a cubic B-spline element model for dial b rack et ........... ... ........ .......... .......... ............................................... 14 9 6-54 Contour plot of bending stresses using a tetrahedral element model for dial bracket.....150 6-55 Contour plot of bending stresses using a H8 element model for dial bracket ...............150 6-56 Square plate with a circular inclusion ..................................................... ............ 151 6-57 Contour plot of o using Q4 elements in IBFEM ......................................................151 6-58 Contour plot of o- using Q4 elements in FEM.........................................................152 6-59 Convergence plot for error in strain energy for square matrix with inclusion model......152 6-60 Comparison of ux along the line AB ................................................... ..................153 6-61 Comparison of stresses along the line AB .......................................................... 153 6-62 Com prison of strains along the line AB ........................................................ .............. 154 6-63 Plot of stresses for cases 1 and 2............................................. ............................. 155 or smaller. The Inclusion function is used to construct the solution structure for material discontinuity as shown in Eq. (3-19). u = (1 H'nc)ugl+ H'ncug2 (3-19) This solution structures combines the solution from both grids to represent the solution in the analysis domain. When the Inclusion function is unity, the solution is given by u = ug2, which is the solution from the inclusion part and when the Inclusion function is zero, the solution is given by u = ugl, which is the solution from the matrix. In the region where the inclusion function varies from unity to zero, the solution is a blend of the solutions from matrix grid and inclusion grid. This way of constructing solution structure ensures the displacement continuity of the solution through out the analysis domain. This solution structure allows the normal strain discontinuity; this can be seen from the gradients of the solution structure as shown below OBu (uI' OH( l +H Ou, g2 +OH" g2 '= (1-H ) gu + Hnc _+- u (3-20) Z ua aH1 xu8 a i In this expression, the first term: (1-H H"n) and third term H"nc- are continuous at H"inc Hnc"" the interface boundary while the second term: u1' and fourth term u2 are 8x Dx discontinuous at the interface boundary and vanish in the matrix region and are non-zero in the inclusion region. Therefore these terms provide independent slopes at both sides of the interface resulting in discontinuous normal strains. Modified Weak Form. The weak form of the elastostatic boundary value problem is the principle of virtual work which can be written in the following form as shown in Eq. (3-21). Here {6S} is the virtual strain, {0} is Cauchy stress tensor, {t} is the traction vector, {6u} is Pengcheng S., Peixiang H., 1995. Bending analysis of rectangular moderately thick plates using spline finite element method. Computers and Structures 54(6) 1023-1029. Rvachev V.L., Shieko T.I., 1995. R-functions in boundary value problems in mechanics. Applied Mechanics Reviews 48(4) 151-188. Rvachev V.L., Sheiko T.I., Shapiro V., Tsukanov I., 2001. Transfinite interpolation over implicitly defined sets. Computer-Aided Geometric Design 18(3) 195-220. Sankar, B.V., and Marrey, R.V., 1997. Analytical Method for Micromechanics of Textile Composites. Composites Science and Technology 57(6) 703-713. Shapiro V., 1998. Theory of R-functions and applications: A primer. CPA Technical Report CPA88-3, Cornell Programmable Automation. Sibley School of Mechanical Engineering. Shapiro V., Tsukanov I., 1999. Meshfree simulation of deforming domains. Computer Aided Design 31(7) 459-471. Sokolnikoff I.S., 1956. Mathematical Theory of Elasticity, McGraw-Hill, New York. Sukumar N., Moran B., Belytschko T., 1998. Natural element method in solid mechanics. International Journal for Numerical Methods in Engineering 43(5) 839-887. Sun, C.T., Chen, J.L., 1990. A micromechanical model for plastic behavior of fibrous composites. Composites Science and Technology 40 115-29. Sun, C.T., Vaidya, R.S., 1996. Prediction of composite properties from a representative volume element. Composites Science and Technology 56(2) 171-179. Surana K.S., Rajwani A., Reddy J.N., 2006. The k-version finite element method for singular boundary-value problems with application to linear fracture mechanics, International Journal for Computational Methods in Engineering Science and Mechanics 7(3), 217-239. Taliercio, A, 2005. Generalized plane strain finite element model for the analysis of elastoplastic composites. International Journal of Solids and Structures 42(8) 2361-2379. Timoshenko S.P., Goodier J.N., 1969. Theory of Elasticity, McGraw Hill, New York. Tyrus, J.M., Gosz, M., DeSantiago, E., 2007. A local finite element implementation for imposing periodic boundary conditions on composite micromechanical models. International Journal of Solids and Structures 44(9) 2972-2989. Vermeulen A.H., Heppler G.R., 1998. Predicting and avoiding shear locking in beam vibration problems using the B-spline field approximation method. Computer Methods in Applied Mechanics and Engineering 158(3-4) 311-327. Whitney, J.M., Riley, M.B., 1966. Elastic properties of fiber reinforced composite materials. A1AA Journal 4 1537-42. Df"t/udoa 9E 1 6 i E I ilE-2 21E.55 E-12 (a) ,E-I (b) Figure 6-3. Plot of Dirichlet functions (a) D, (b) D, kVT. fi = qh h,TO T=T Figure 3-3. Steady state heat transfer problem Curves represented by implicit equations 05- 0.0 Union: D1u2 = Di+D2-DiD2 11 Intersection: Din2 = D1D2 Subtract: D-2 = DiD-2 Figure 3-4. Boolean operations on Dirichlet functions interpolation/approximation using the element basis functions. The nodal values of u" can be constructed similar to the grid variable ug using the grid element basis functions. u = -Nu (3-17) In the above equation, N, are the shape functions of the grid element and u' are the nodal values of u". It is important that the function u" should belong to the space spanned by the shape functions used to represent the grid variable. Otherwise, the grid variable has difficulty compensating for u~ in the solution structure leading to poor convergence. Indeed, the solution structure will not satisfy completeness requirement, that is the ability to accurately represent constant strains, if basis functions used to construct u" does not span a sub-space of the basis functions that are used to represent the grid variable ug. If the basis functions N, in Eq. (3-17) are identical to the basis functions used for approximating the grid variable then this requirement is automatically satisfied. If solutions that have C1 or C2 continuity is desired then quadratic or cubic B-spline basis should be used in Eq. (3-17). The nodes of the boundary elements should be assigned nodal values of ua such that the resulting approximation will have the desired values at the boundary passing through the element. This is very easy to achieve when the assigned value is constant or even linearly varying along the boundary. When the specified boundary condition is a constant value of displacement along boundary, all the support nodes of the element have to be assigned value equal to the specified boundary condition. In the case of the 2D cubic B-spline element, the approximation within the element is controlled by 16 support nodes. Therefore, the value at all these nodes should be assigned to the same value to ensure that the resulting approximation has a constant value within the element. Similarly, if the specified boundary value has to vary linearly (as in hydrostatic pressure which varies linearly in the z-direction), the values 3.3.2 Modified Weak Form The solution structure for grid refinement using the step function for refinement, H'r can be written as shown in the previous section and the virtual displacements are written as: u = (1- Href)(D3u')+Href (D3u2) (3-27) The grid solutions in both grids are approximated by ug' = {N1{Xel with N1 = (1-Hr)N and ug2 = {N2Xe2} with N2 = H N, where Nrepresents the shape functions for grid variables. The boundary values are approximated u"' = {N"' } X"' with N" =(i-Href)N and ua2 = {a2) Xa2 } with Na2 =HrfN where N represents the basis functions used to approximate the boundary value function. Using the above approximations, the expressions for strain takes the following forms: a8u Eu X el X al 12 2, a In the above expression, B, and B2 are modified strain displacement matrices which contain the terms involving the gradients of H'ref, Dirichlet functions and the gradients of the shape functions used to approximate the grid solution and Bal, Ba2 are strain displacement matrices which contain the gradients of the shape functions used to construct boundary value functions. The expression for stress is given as o- = C ,, where C is the elasticity stiffness matrix, and in computation of stress. The modified weak form takes the following discrete form when the above expressions for stress and strain are incorporated into the weak form (Eq. 3-29). virtual displacement and {b} is the body force. For elastostatic problems, the essential boundary conditions are displacements {u} = { u } specified on F, and the natural boundary conditions are traction {t} = {to} specified on Ft. J {E T{a}d=Q S{u}T {t}dF + 6u}T {b}dQ (3-21) n Ft, When the solution structure shown in Eq. (3-19) is used and the displacements are interpolated as ug1 = [A1]u l }, with N 1 = (1- H') N in grid corresponding to matrix and ug2 = [g2 ]u 2} with Ng2 = H"nCN in inclusion grid, the expression for strains takes the following form Ou1 Ox, Ei =2 [Bg B 2]{U}gl (3-22) 712 2 +- In this expression, [B1 B2] is the modified strain displacement matrix and is a combination of strain displacement matrices from the two grids. The expression for stress is given as o-, = C,,c, where C, is the elasticity stiffness matrix, and in computation of stress, the inclusion (or matrix) elasticity stiffness matrix is used when the strains are evaluated in the inclusion (or matrix) region. The modified weak form takes the following discrete form when the above expressions for stress and strain are incorporated into the weak form (Eq. 3-23). AVE TE T z(X'e} J[B, B2] [C][B, B2]{Xe}d= Z{X'e}T (I N 2} {b}d2+ el a e-l (3-23) {Nxe} { F1 2} {t}dF e=1 Fp 5 NUMERICAL IMPLEMENTATIONS.....................................67 5.1 Numerical Implementations for Linear Elastic Problems ..........................................67 5.1.1 Computation of Stiffness Matrix for Internal Elements.......................................68 5.1.2 Computation of Stiffness Matrix for Boundary Elements................................. 68 5.1.3 Computation of Stiffness Matrix for Elements with Essential Boundaries............69 5.1.4 Com putation of Load V ector ....................................................... ............... 73 5.2 Numerical Implementations for Steady State Heat Transfer................ .............. ....74 5.3 Numerical Implementations for M icromechanics ................................. ................ 76 5.3.1 Computation of Stiffness Matrix in the Matrix Region.......................................79 5.3.2 Computation of Stiffness Matrix in the Inclusion Region............... .................79 5.3.3 Computation of Stiffness Matrix in the Interface Region ....................................80 5.3.4 Periodic B oundary Conditions ........................................ .......................... 82 5.3.5 Computation of Effective Properties ....................................................... 84 5.3.6 Com putation of Load V ector.............................................. .............................. 87 5.4 Numerical Implementations for Local Grid Refinement.................. .................87 5.4.1 Computation of Stiffness Matrix in Regions 1 and 2..........................................88 5.4.2 Computation of Stiffness M atrix in Region 3 ................................ ............... 89 5.4.3 Computation of Load Vector................ .... ...... ..................90 5.5 Extrapolation of Solution for Elements with Low Volume Fraction ............................90 5.5.1 Construction of Extrapolation M atrix ......................... ........... ............... .... 92 5.5.2 Construction of Extrapolation Matrix for a Boundary Element..........................94 5.5.3 Computation of Stiffness M atrix ............................................................. 94 6 A N A L Y SIS A N D R E SU L T S...................................................................... ...................103 6.1 Validation of Solution Structure for Essential Boundary Conditions ..........................103 6.1.1 P watch T est ................... ............ ... ..................... .... .... ........ 103 6.1.2 V alidation of Extrapolation Technique ............................................................. 106 6.1.3 Convergence Analysis of a Cantilever Beam ...................................................... 107 6.1.4 Square Plate with a Circular Hole ........................................................ .......... 110 6.1.5 Thick Cylinder with Internal Pressure ....................................... ................111 6.1.6 Steady State Heat Transfer Analysis ...................................... ...............112 6.1.7 A analysis of a 3D Clam p ...................................................................... 113 6.1.8 Analysis of a 3D Rib Structure ................................................................. 114 6.1.9 Analysis of a 3D Dial Bracket .................................................... ...............115 6.2 Validation of Solution Structure for Micromechanics..............................116 6.2.1 Square Plate with Circular Inclusion.................... ..... ................116 6.2.2 Determination of Effective Properties of a Fiber Reinforced Composite ..........17 6.3 Validation of Solution Structure for Grid Refinement ............................118 6.3.1 Application of Grid Refinement to Constant Stress Problem ............................1118 6.3.2 Stress Concentration Problem ..................................... ......................... .......... 119 7 CONCLUSIONS AND FUTURE WORK............................................... ............... 165 7 .1 C o n clu sio n s....................................................... ................ 16 5 7 .2 F utu re W ork ...........................................................................167 Computation of Stiffness Matrix. The stiffness matrix comprises of a volume integral and a surface integral. The volume integral is performed as explained in section 5.1.1 and 5.1.2 after subdividing the boundary elements into triangles or tetrahedra. The surface integral for load computation is evaluated as explained in section 5.1.4. For the boundary elements with essential boundaries, the [B] matrix can be decomposed as shown below: D N +N 0 [Bl] = D Ol (5-15) 0 D '2 + In the above equation, [B] matrix contains the terms involving gradients of shape functions and also the gradients ofD function. It is advantageous to decompose this matrix further to [B] = [B1] + [B2] such that contains terms involving the gradient of the shape functions and [B2] contains terms with the gradients of the D function as shown below: D 0 0 ax, Qx, 0 D-N 0 N Nx, Ox, N1D 0 8D, 8x, ~N, . [B 2]2 0 N, = N 0 =D2 B2 (5-17) a2 00 N1 10 = 2DJLo N1 Sx2 0 N . )2 1X N, N- Using the above procedure, the gradient matrix can be rewritten as [B] = D [B + D2 [B, Similar matrices can be defined for 3D elements also. Using this decomposition, the stiffness matrix is decomposed into the following components: BCB 0 (5-43) K= [B1 B2 TC[B1 B2]dV= V (5-43) Q, Q, 0 0 The computation of stiffness matrix is performed in a similar manner as explained in sections 5.1.1-5.1.3. In the Region 2, the step function of refinement is equal to one and the solution structure becomes u = Dug2 + ua2 and the strain displacement matrix will be similar to the one shown in Eq. (5-42) and hence the stiffness matrix becomes: K= [B, B2]TC[B1 B2]dV= J B CB V (5-44) Q, Q, -- 120 B CB The stiffness matrix is computed as explained in sections 5.1.1-5.1.3. 5.4.2 Computation of Stiffness Matrix in Region 3 In the Region 3, the step function of refinement is varying and the Dirichlet functions are unity. It is to be noted that there is a small region of the size 6 where both Dirichlet functions and step function of refinement may be active, since this region is very small, the contribution to the stiffness matrix in this region is negligible. The solution structure takes the form: u = (l- Href) (ug" + ual) + Href (Ug2 + ua2) and hence the strain displacement matrix takes the form B = [B B2] (1- H') 0 NO 0-Ngl 0 ax, ax, B, = B, +B,2 = 0 (1-H) x + 0 N --- (5-45) a(1-H)2 (1-H X 2 6x2 1x, 2 1i SCBS QBS -- FEQ4 FEQ9 10-4 10 103 Log (Num Nodes) Figure 6-26. Convergence plot for error in strain energy for stress concentration problem Plot of X-Stress along 0 = i/2 Analytical CBS QBS 1.5 2 2.5 3 Distance from center of circular hole, m 3.5 4 Figure 6-27. Plot of along the line 0 N ricr -----i~t- r- w "N CHAPTER 3 SOLUTION STRUCTURES FOR IMPLICIT BOUNDARY METHOD The essential boundary conditions are imposed in FEM by specifying the values at the boundary nodes on which essential boundary conditions are specified. When structured grids are used for discretizing the analysis domain, nodes are not guaranteed to be present on essential boundaries; hence special techniques need to be devised to be able to impose essential boundary conditions. In this chapter, a detail discussion on construction of solution structure for exact imposition of Dirichlet boundary conditions using the implicit equations of boundaries is presented. Approximate Heaviside step functions are used to construct a solution structure based on the technique proposed by Kantorovich and Krylov (Kantorovich and Krylov, 1958 and is termed as Implicit Boundary Method (Padmanabhan, 2006 ; Kumar et. al. 2007). Non homogenous essential boundary conditions are imposed by constructing piecewise continuous functions which satisfy the specified values at the essential boundaries and are termed as boundary value functions. This solution structure is used for modifying the weak form corresponding to linear elasticity problems and steady state heat transfer problems. Implicit Boundary Method has been extended in this dissertation to be able to use B-spline basis. A solution structure for treatment of material discontinuity based on the implicit equation of the material boundaries and approximate Heaviside step function is presented later in this chapter. Techniques for imposing periodic boundary conditions using multipoint constraints are also discussed in detail. In the final part of this chapter, a solution structure for adaptive refinement of grid based on the implicit equations of the grid element boundaries and approximate Heaviside step functions is discussed in detail. The weak form corresponding to linear elasticity is modified using this solution structure. Hillig K., Reif V., Wipper J., 2001. Weighted extended B-spline approximation of Dirichlet problems. SIAM Journal on Numerical Analysis 39(2) 442-462. Hillig K., Reif V., 2003. Nonuniform web-splines. Computer Aided Geometric Design 20(5) 277-294. Hillig K., 2003. Finite element methods with B-Splines. SIAM, Philadelphia. Kagan P., Fischer A., Bar-Yoseph P.Z., 1998. New B-spline finite element approach for geometrical design and mechanical analysis. International Journal for Numerical Methods in Engineering 41(3) 435-458. Kagan P., Fischer A., Bar-Yoseph P.Z., 2003. Mechanically based models: Adaptive refinement for B-spline finite element. International Journal for Numerical Methods in Engineering 57(8) 1145-1175. Kantorovich L.V., Krylov V.I., 1958. Approximate Methods of Higher Analysis. Interscience, New York. Kim, H.J., Swan C.C., 2003. Algorithms for automated meshing and unit cell analysis of periodic composites with hierarchical tri-quadratic tetrahedral elements. International Journal for Numerical Methods in Engineering 58(11) 1683-1711. Kohno H., Tanahashi T., 2004. Numerical analysis of moving interfaces using a level set method coupled with adaptive mesh refinement. International Journal for Numerical Methods in Fluids 45(9) 921-944. Krongauz Y., Belytschko T., 1996. Enforcement of essential boundary conditions in meshless approximations using finite elements. Computer Methods in Applied Mechanics and Engineering 131(1-2) 133-145. Krongauz, Y., Belytschko, T.,1998. EFG approximation with discontinuous derivatives. International Journal for Numerical Methods in Engineering 41(7) 1215-1233. Krongauz Y., Belytschko T., 1997. Consistent pseudo-derivatives of meshless methods. Computer Methods in Applied Mechanics and Engineering 146(3-4) 371-386. Krysl P., Grinspun E., Schroder P., 2003. Natural hierarchical refinement for finite element methods. International Journal for Numerical Methods in Engineering 56(8) 1109-1124. Krysl P., Belytschko T., 1997. Element-free Galerkin method: convergence of the continuous and discontinuous shape functions. Computer Methods in Applied Mechanics and Engineering 148(3-4) 257-277. Kumar A.V., Padmanabhan S., Burla R., 2007. Implicit boundary method for finite element analysis using non-conforming mesh or grid. International Journal for Numerical Methods in Engineering. (In press). Clark B.W., Anderson D.C., 2003b. The penalty boundary method for combining meshes and solid models in finite element analysis. Engineering Computations 20(4) 344-365. Cook, R.D., Malkus D.S., Plesha M.E., Witt, R.J., 2003. Concepts and Applications of Finite Element Analysis. John Wiley and Sons Inc., New York. Cordes, L.W., Moran, B., 1996. Treatment of material discontinuity in the element-free Galerkin method. Computer Methods in Applied Mechanics and Engineering 139(1-4) 75-89. Cottrell J.A., Reali A., Bazilevs Y., Hughes T.J.R., 2006. Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering 195 (41-43) 5257- 5296. Cottrell J.A., Hughes T.J.R., Reali A., 2007. Studies of refinement and continuity in isogeometric structural analysis. Computer Methods in Applied Mechanics and Engineering 196(41-44) 4160-4183. Dang, T.D., Sankar, B.V., 2007. Meshless local Petrov-Galerkin formulation for problems in composite micromechanics. AIAA Journal 45(4) 912-921. De S., Bathe K.J., 2000. Method of finite spheres. Computational Mechanics 25(4) 329-345. De S., Bathe K.J., 2001. The method of finite spheres with improved integration. Computers and Structures 79(22-25) 2183-2196. Dolbow J., Belytschko T., 1999. Numerical integration of the Galerkin weak form in meshfree methods. Computational Mechanics 23(3) 219-230. Duarte C.A., Oden J.T., 1996. h-p adaptive method using clouds. Computer Methods in Applied Mechanics and Engineering 139(1-4) 237-262. Duarte C.A., Babuska I., Oden J.T., 2000. Generalized finite element methods for three- dimensional structural mechanics problems. Computers and Structures 77(2) 215-232 Farin G., 2002. Curves and Surfaces for CAGD: A Practical Guide. Morgan Kaufmann Publishers, San Francisco. Gibson, R. F., 1994. Principles of Composite Material Mechanics. McGraw- Hill, New York. Gingold R.A., Monaghan J.J., 1982. Kernel estimates as a basis for general particle methods in hydrodynamics. Journal of Computational Physics 46 429-53. Hughes T.J.R., Cottrell J.A., Bazilevs Y., 2005. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering 194 (39-41) 4135-4195. Hughes T.J.R., 2000. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications, New York. Ta=30, h=15 ^=80 J=50 T=80 0000 0 80 Ta=30, h=15 (a) 1m Ta=30, h=15 T=80 T=50 .5m 0.25 m A B (b) Figure 6-34. Heat conduction in a block with series of pipes Dfunction T 1EO 9E-1 8E-1 7E-1 6E-1 2E-1 3F- 1 1E-1 OEO Figure 6-35. Plot of Dirichlet function Figure 6-24. Contour plot of o using cubic B-spline elements 105 10 103 Log (Num Nodes) Figure 6-25. Convergence plot for L2 error norm for stress concentration problem 2 766IEO 2,436E0 2.10BSEO 1.77MED 1.45B2ED 1.1211EO 70201E.1 4A29E-1 1 15WIVE-1 -10532E-1 e CBS QBS FEQ4 FEQ9 -6 x 106 Theoretical displacement --- 8=1.0-1 5 -e- 8 =5.0-2 a- 8 =1.0-2 -*-- 8 =1.0-3 S4 I 8 =1.0-5 E O. 2 -I In 1 0 0.2 0.4 0.6 0.8 Distance from left edge of plate, m Figure 6-4. Effect of range parameter on the quality of solution Figure 6-5. Contour plot of o- for one element model using cubic B-splines Stress X 1E6 9.9E5 9.8E5 9.7E5 9.6E5 9.5E5 9.4E5 9.3E5 9.2E5 9.1E5 9E5 Temperature T 8.0039E1 7.5163E1 7.0287E1 6.5411E1 6.0535E1 5.566E1 5.0784E1 4.5908E1 4.1032E1 S.6156E1 3.1281E1 Figure 6-38. Contour plot of temperature distribution using cubic B-spline elements -- --::: i--: -J I1III-f1 1UN U E1t:. ....... /" Temperature T 8.0078E1 7.5199E1 7.0319E1 6.544E1 6.0561El 5.568?FF1 5.'f0OE I 4.5923EEI 4.1044E1 I3.6165E1 3.1286E1 Figure 6-39. Contour plot of temperature distribution using quadratic B-spline elements IgLIIIIIIIIIIIIIIIIl/l = 3vy2L-x)+ 4+5v)Dx+3L-x)x2 (6-2b) P P (D2 > cx= (L-x)y ; = -y ;oY=0 (6-2c) I 21 4 A typical non-conforming grid used for analysis is shown in Figure 6-15. Even though this problem can be modeled using a conforming grid, a non-conforming grid is used in order to demonstrate the convergence behavior of the solution structure along with B-spline basis. The B- spline elements have nodes outside the elements which form an extra layer of nodes as can be seen in Figure 6-15. The consistent nodal loads are distributed on the nodes corresponding to the boundary elements on which traction is specified. A convergence analysis is performed for this model where a series of models with increasing grid density are used for analysis and the corresponding error norms are computed and plotted. The L2 error norm represents the errors in the displacements and the expression which defines error is shown in Eq. (6-3). Energy error norm represents the error in stresses and strains and the expression for defining energy error is shown in Eq (6-4). In these expressions, u, s and a represent displacement, strain and stress respectively, while the superscript e represents exact solution and h represents approximate solution. L = K(ue ue ) .(u euh)Q (6-3) E = j(ae _h(Teh)dQ2 (6-4) The convergence plot for L2 error norm and energy error norms are computed for various grid densities and shown in Figure 6-16 and Figure 6-17 respectively. [K,]= f[B1]T [k][B, ]dV + f[B,]T [k][B2]dV + f[B2 ]T [k][B, ]dV + f[B2 ]T [k][B2]dV Sv (5-18) =[K + ([K2]+[K ]T) +[K] The first component [K, is computed as described earlier by decomposing the boundary elements into triangles or tetrahedra to perform the area or volume integrals respectively. The second and third components contain gradients of Dirichlet functions which are large near the boundary when 6 is made very small and they vanish away from the boundary, hence the terms in the integrals will be non-zero in a thin band along the boundary. This allows us to convert the volume integral into a surface integral as explained in detail in section 5.1.3. The second term in Eq. (5-18) can be written as: [K2]= [B1]T [k][B2]dV = j [Bl] f[D1j[k][D2 dBB2 dF E r (5-19) = [B C2 2[B,] dF This term is evaluated as a surface integral as explained in section 5.1.3. 5.3 Numerical Implementations for Micromechanics The modified weak form for the solution structure corresponding to micromechanics was discussed in a previous chapter and was shown to be: NE T IT {8Xe} JI B2] [C][B, B2] {XedQ {Xe} T{N N2} {b}d+ e1 C el (5-20) {jXe}T IN1 N2 }{t}dF e=1 pt {Xe} = {g ug2 T (5-21) In this section, a detailed discussion on the numerical implementations for the computation of the terms involved in this weak form is presented. The left hand side of the weak form corresponds to the stiffness matrix while the right hand side corresponds to the load vector. The solution structure for micro-mechanics is such that the solution from the matrix grid inside the inclusion is zero due to the step function of inclusion and similarly the solution from inclusion grid outside of the inclusion vanishes. Therefore the nodal values of the matrix grid inside the inclusion are set to zero, provided these nodes are outside the influence ofH'"C and the nodal values of the inclusion grid outside the inclusion are set to zero. This leads to the reduction in the total number of degrees of freedom for the analysis model. The finite element solution for the grid variables are written as ug1 = [Ng' ]{ug)} and u2 = [Ng2 ]u2 }. For two dimensional case, the strain vector using the strain displacement relations is written as shown below: uu S a2 [B1 B2 Ug (5-22) 712I +- a 2 1, The components of the strain displacement matrix B1,B2 are expressed as: (1 H ) agl w 0 N 0 (1 H O Ha) )x, 1 x, Hx, xN a2 a21 OX1, O [K]= f[B]T [C][B1]dV+ [B,]T [C][B]dV+ [B2 ]T [C[B1]dV+[B2 ]T[C][B]dV Sv v (5-8) =[K1]+ ([K2]+[K2 ]T +[K3,] The first component [K, is computed as described earlier by decomposing the boundary elements into triangles or tetrahedra to perform the area or volume integrals respectively. The second and third components contain gradients of Dirichlet functions which are large near the boundary when 6 is made very small and they vanish away from the boundary, hence the terms in the integrals will be non-zero in a thin band along the boundary. This allows us to convert the volume integral into a surface integral. The contribution from the gradients of finite element shape functions can be assumed to be constant across the band but varies along the boundary and the contribution from the gradients of Dirichlet function are constant along the boundary but vary normal to the boundary. The band is represented in terms of two parameters t,n as shown in the Figure5-3 and the terms in stiffness matrix take the following form. The stiffness matrix as decomposed in Eq. (5-8) has four components; the second component [K2] can be further decomposed as shown in the following equation. [K= J[B T[c][B2]dV D [C] D2 d1B d VE (5-9) S[BJTf J[C2 [B2]dF In this expression, the terms involving only Dirichlet functions and their gradients are separated from the terms containing the basis functions and their gradients. This decomposition allows us to compute the stiffness matrix efficiently and accurately. The volume integral has been decomposed into a line integral along the normal to the boundary and a surface integral along the boundary. The normal for the boundary represented by implicit equation 0 = 0 is given The multipoint constraints in Implicit Boundary Method can be implemented exactly as is done in traditional finite element method (Cook et. al. 2003). It is to be noted that the RVE is always a rectangle/cuboid, hence when Lagrange interpolations are used, the nodes are guaranteed to be present on the boundaries on which multipoint constraints are imposed. Furthermore, since structured grids are used, identical nodes are always guaranteed to be present on the opposite faces and facilitate the imposition of periodic boundary conditions. However, when FEM is used for analysis, care needs to be taken to generate mesh such that there are identical nodes present on opposite faces. The nodes on the faces (or edges in 2D) of the RVE can be partitioned into two parts: independent or master nodes (s) and dependent or slave nodes (s). The nodal displacement vector can be similarly partitioned into u = [um u], where um are displacement of the master nodes and u, are those of slave nodes. The displacements for the slave nodes are specified in terms of the master nodes in the following form us = Giu + g1, where gi is the non-zero periodic boundary conditions (e.g., Case 1: u =u k + L). Hence the displacements are expressed as : U = [Um Us= IG [um]+ = G[um]+g (5-34) where, I denotes the identity matrix and 0 is the zero vector. The matrices G and g are constructed globally for all the master degrees of freedom present in the analysis model. Using this expression for displacement and a similar expression of the corresponding virtual displacement 6u = G [6um ] in the weak form, the global stiffness matrix is transformed into K' = GTKG and the load vector is increased by -GTKg to impose the multipoint constraints. It is to be noted that in a general case of multipoint constraints, G and g are constructed such that Penalty boundary method developed by Clark and Anderson (Clark and Anderson, 2003a, 2003b) uses a regular grid and a constrained variational formulation to include both essential and natural boundary conditions which is then solved using penalty method. Shapiro (Shapiro, 1998 ; Shapiro and Tsukanov, 1999) and Hollig (Hollig et. al. 2001; Hollig and Reif 2003 ; Hollig 2003) used R-functions to define the boundary of the problem domain and the solution structure proposed by Kantorovich and Krylov (Kantorovich and Krylov 1958) was constructed to satisfy Dirichlet boundary conditions. Shapiro used transfinite Lagrange interpolation to impose Dirichlet boundary conditions while Hollig used weight functions based on distance functions and R-functions. Hollig developed extended B-spline basis to represent the solution in the analysis domain. The method presented in this dissertation is very similar to the approach used by Hollig, however approximate Heaviside step function of essential boundaries are used in this work to construct weight functions and the behavior of solution structure for very small values of range parameter used to construct step functions is also studied. All the internal elements in our method have identical stiffness matrices unlike Hollig's method. 2.4 Meshless Methods A number of methods have been developed in the past two decades which avoid mesh generation. These methods are collectively termed as meshless or meshfree methods. Some of these methods provide continuous gradients in the analysis domain. Smoothed particle hydrodynamics (SPH) is one of the first meshless methods in literature (Lucy, 1977 ; Gingold, 1982); it was first used to model astrophysical phenomenon and was modified to include continuum mechanics problems (Monahan, 1982). This approach uses Shepard shape functions as interpolants and uses point collocation method for discretization and nodal integration for stiffness matrix computation. The drawbacks of this approach are poor accuracy and spatial CHAPTER 2 ANALYSIS TECHNIQUES FOR BOUNDARY VALUE PROBLEMS In this chapter, a detailed discussion on several analysis techniques that are used for solution of engineering boundary value problems is presented. In the first part of this chapter, a brief overview of the finite element method (FEM) is presented along with its advantages and shortcomings. Techniques for employing B-spline approximations for analysis are discussed in the next part of this chapter. A detailed discussion on structured grid methods and meshless methods is presented in the later part of this chapter Micromechanical analysis is an important engineering application where FEM is used extensively. Microstructures are modeled using a representative volume element and periodic boundary conditions are applied using multipoint constraint techniques. Meshless methods have been extended for treatment of material discontinuity and for performing micromechanical analysis. A survey on these techniques is presented in this chapter along with a survey on employing FEM for determination of effective properties. Another important feature of FEM is its capability for local refinement. A survey on the refinement techniques available in structured grid and meshless framework is presented in the last part of this chapter. 2.1 Finite Element Method Finite element method (Hughes, 2000 ; Cook et. al. 2003) is a well established numerical technique and is widely used in solving engineering problems in academia as well as in industry. For finite element analysis, the domain of analysis is subdivided into elements and the resulting finite element mesh approximates the geometry and is also used to approximate the solution by piece-wise interpolation within each element. Using Galerkin's approach, the governing differential equation in its strong form is converted into a weak form which allows the use of Table 6-1. Average stresses for various loading cases Case Unit oa,(GPa) o22 (GPa) r33 (GPa) 23 (GPa) r31 (GPa) 12 (GPa) Strain 1 e-, =1 230.84 41.14 41.14 0 0 0 2 22 =1 40.35 161.20 46.22 0 0 0 3 E33 =1 40.35 46.22 161.20 0 0 0 4 23 = 1 0 0 0 46.07 0 0 5 731 = 10 0 0 0 54.44 0 6 12 = 1 0 0 0 0 0 54.44 k the basis functions form a partition of unity N, = 1 is used to obtain (k +1)2 equations needed 1=0 to solve for all coefficients of the k +1 basis functions. Using this approach, the basis functions for quadratic and cubic B-splines can be derived as shown in Eq. (4-4) and the basis functions are plotted in Figure 4-2. Quadratic B-spline element. The shape functions for quadratic B-spline element are constructed using the procedure explained above. Quadratic B-spline element in one-dimension has three nodes and hence three shape functions. The plots of the shape functions are shown in Figure 4-2 and the corresponding equations are shown below: N, =(1-2r+r2);N2= (6-2r2);N3= l(1+2r+r2) (4-4) 8 8 8 Cubic B-spline element. The plots of basis functions for cubic B-spline are shown in Figure 4-3 where it can be seen that, unlike traditional finite element shape functions, the B- spline basis functions are not unity at the corresponding node and do not vanish at other nodes. Therefore, these basis functions do not satisfy the Kronecker-delta property and the approximation constructed using these basis functions does not interpolate nodal values. 1 1 N, = (1- 3r+3r2 -r);N2= (23-15r -3r2 +3r3); 48 48 (4-5) 1 1 N= --(23+15r- 3r2- 3r3); N4= (1+3r +3r2 +r3) 48 48 4.2 Two and Three Dimensional B-Spline Elements The basis functions for the higher dimensional B-spline elements are constructed by taking product of the basis functions for one-dimensional B-splines. A typical grid used for analysis is illustrated in Figure 4-4. This grid corresponds to cubic B-spline elements where the span corresponding to each element is controlled by 16 support nodes. The grid elements used here are regular quadrilaterals (rectangle/square) or regular hexahedra (cube/cuboid), hence the CHAPTER 7 CONCLUSIONS AND FUTURE WORK 7.1 Conclusions In this dissertation, finite elements based on uniform B-spline basis defined over a structured grid are constructed to obtain solutions with continuous gradients throughout the analysis domain. Finite elements are developed using quadratic and cubic B-spline basis for problems in both two and three dimensions. B-spline basis do not satisfy Kronecker-delta property and hence Implicit Boundary Method is extended for using B-spline basis and for applying essential boundary conditions. Solution structures are developed using the approximate Heaviside step functions based on the implicit equations of the essential boundaries to impose essential boundary conditions. The weak form for linear elasticity and steady state heat transfer was modified using this solution structure such that essential boundary conditions are automatically imposed. Methods were developed for efficiently and accurately evaluating the terms involved in the modified weak form in the limit as the approximate Heaviside step functions approach the exact Heaviside step functions. The primary motivation for the method is the desire to obtain continuous gradients in the entire analysis domain using B-spline basis on a structured grid. The B-spline elements studied in this work exhibited good convergence behavior and the comparison with traditional FEM solutions showed that very accurate solutions can be obtained with relatively fewer nodes. The performance of B-spline elements does not degrade with increasing aspect ratio of the elements and one element can be used to compute accurate solutions for beam like structures. Although B- spline basis is used in this work to achieve continuous stress and strains, the method is quite general and can be applied to other basis functions which may or may not satisfy Kronecker- they may be a fully populated matrix and vector respectively. However, when structured grids are used for analysis, identical nodes are guaranteed to be present on the opposite faces and hence with proper choice of node numbers, G matrix can be constructed to be a diagonal matrix. This simplifies the computations involved in the transformation of the stiffness matrix and hence leads to efficient implementation. The periodic boundary conditions involved in a square array of fibers are shown in Figure 3-10 and the corresponding multipoint constraint equations are shown in Table 1. 5.3.5 Computation of Effective Properties Effective properties of the composite materials with periodic distribution of microstructures are computed by applying far field strains to the RVE as shown in Figure 5-7. The displacements, strains and stresses are assumed to be periodic in the RVE. Hills lemma states that the average of the strain energy in RVE is same as the strain energy due to average stress and strains. It can be stated in the form (a : ) = () : (E) .This property is used to compute the effective properties from an RVE. Using the stress-strain relation (a) = CH(), the effective properties are represented as (a):(: ) = (CHE: ). The average stresses and strains are computed by solving the boundary value problem six times with different applied strains. These average stresses and strains are used in the relation (a) = CH(") to determine a row (column) of the homogenized stiffness matrix CH. Six different strain states are used to determine all the components of CH. The strain-stress relation for orthotropic materials is given as follows. Once the homogenized stiffness matrix is determined, the effective properties can be computed by using the strain-stress relation. 4.25 LU C 2 4.2 Cl) Log (Num Nodes) Figure 6-42. Convergence plot for strain energy for clamp model S, S11 (Avg: 75%) +1.249e+05 +1.041e+05 +8. 36e-+04 S+6.2553e+04 +-4.1 9e+04 +2.100oe+'04 +2. 189e+02 -2.057e+04 -4.135e+04 -6.214e+C'4 -8.292e+C'4 -1.037e+C05 -1.245e+C'5 Figure 6-43. Contour plot of bending stresses using H8 elements in FEM for 3D clamp HinC aNg2 + 2 H +-N 0 2, ax n Ng2 Ngl 8H nc B2 =0 H+"c -N, .. (5-23b) H +-NN2 H nc _+ N+g2 a8x2, x, x , In the above equations, the decomposed components of strain-displacement matrix have terms involving gradients of finite element shape functions Nk and the gradients of the ax inclusion function -- It will be advantageous to further split these components into a term ix containing gradients of shape function only and another term containing gradients of inclusion function only. This is accomplished as shown in the following equations. xl, &9H, (81 H= a, 0 (1 -A)gl 0 B, =B +B,12= 0 (1 H'-) .+ 0 --- gl" (5-24) a' a' 2H"" a2H" (1-H"") (1-H") 'g' ax, ax a2x, C9x H 0 n a'g2 0 H 8x, C' B2 = B21 +B22 = H"- .+ 0 N2 (5-25) Cx, ax, 8x2 C9x, x2 aHx ' It can be noticed that the modified strain displacement matrices B, B12, B21 and B22 vanish in certain regions of the computational domain and are non zero in other parts of computational domain depending on where the inclusion function and its gradients vanish. 6-16 Convergence plot of L2 error for cantilever beam .................................. ...............129 6-17 Convergence plot for error in strain energy for cantilever beam................................130 6-18 Comparison of deflection curves for 10 element model....................... ...............131 6-19 Comparison of normal stress along the line AB for 10 element model..........................132 6-20 Comparison of shear stress along the line AB for 10 element model............................ 133 6-21 Convergence plot for tip deflection in near incompressible limit.................................133 6-22 Analysis model for a square plate with a circular hole....................................................134 6-23 Contour plot of ox using FEM Q4 elements ............................................... ............. 134 6-24 Contour plot of o- using cubic B-spline elem ents......................................................... 135 6-25 Convergence plot for L2 error norm for stress concentration problem ..........................135 6-26 Convergence plot for error in strain energy for stress concentration problem ..............136 6-27 P lot of o- along the line 0 = ........................................................................... 136 2 6-28 Hollow thick cylinder subjected to uniform pressure............ ............... ............... 137 6-29 Contour plot of radial displacements using cubic B-spline element model.....................137 6-30 Contour plot of radial displacements using FEM Q4 elements.................. ............. 138 6-31 Comparison of radial displacements along the thickness of the cylinder........................138 6-32 Comparison of radial stresses along the thickness of the cylinder ...............................139 6-33 Comparison of tangential stresses along the thickness of the cylinder............................139 6-34 H eat conduction in a block with series of pipes ............................................................. 140 6-35 P lot of D irichlet function ....................................... ............ ............ ............................140 6-36 Plot of boundary value function.............................................. ............................. 141 6-37 Contour plot of temperature distribution using FE Q4 elements................................... 141 6-38 Contour plot of temperature distribution using cubic B-spline elements ........................142 6-39 Contour plot of temperature distribution using quadratic B-spline elements................142 Stress y 2.9166E6 2.1803E6 1.4439E6 7.0747E5 2.8916E4 7.6531E5 1.5017E6 2.2381E6 -2.9745E6 -3.7109E6 -4.4473E6 Figure 6-12. Contour plot of o-y for Q4 elements using extrapolation technique 6-64 Plot of stresses for cases 3 and 4............................................. ............................. 155 6-65 Plot of stresses for cases 5 and 6............................................. ............................. 155 6-66 Plot of stresses for cases 1 and 2 for 3D model .................................... ............... 156 6-67 Plot of stresses for cases 3 and 4 for 3D model .................................... ............... 156 6-68 Plot of stresses for cases 5 and 6 for 3D model .................................... ............... 156 6-69 Contour plot of u using Q4 elements and grid refinement............................................157 6-70 Contour plot of v using Q4 elements and grid refinement ..........................................157 6-71 Contour plot of u using cubic B-spline elements and grid refinement...........................158 6-72 Contour plot of v using cubic B-spline elements and grid refinement...........................158 6-73 Square plate w ith a circular hole......... .................................................. ............... 159 6-74 Contour plot of o- for local refinement using Q4 elements .......................................159 6-75 Contour plot of o- for local refinement using cubic B-spline elements.........................160 6-76 Convergence plot for L2 error in local grid refinement ..............................................160 6-77 Convergence plot for energy error in local grid refinement ................. ................161 6-78 Contour plot of o- for local refinement using Q4 and cubic B-spline elements ............161 LIST OF FIGURES Figure page 3-1 Solution structure for imposing essential boundary conditions............... ...................49 3-2 Elastic solid under equilibrium ............................................... ............................... 49 3-3 Steady state heat transfer problem ............................................. ............................ 50 3-4 Boolean operations on Dirichlet functions ............................................. ............... 50 3-5 A CSG tree representing Dirichlet function for multiple boundaries..............................51 3-6 Boundary value function using cubic B-spline basis.............................. ............... 52 3-7 Representation of grids for material discontinuity.............. ..... .................. 53 3-8 One dimensional inclusion function for treatment of material discontinuity ....................53 3-9 Two dimensional inclusion function for an elliptic inclusion in a square matrix..............54 3-10 Representation of solution structure for local grid refinement............... ................ ...54 3-11 Approximate step function for refinement................................... .......................... 55 3-12 Coarse and refined grid interface with multiple implicit equations..............................55 3-13 Step function for refinement for a 2D model.................... ..........................56 4-1 Approxim ation using ID cubic B-splines...................................... ........................ 64 4-2 Shape functions for ID quadratic B-spline element .................................. ............... 64 4-3 Shape functions for ID cubic B-spline element ..................................... .................64 4-4 Typical grid used for analysis and the mapping used for geometry and field variable.....65 4-5 Two dimensional quadratic B-spline elements in parametric space..............................66 4-6 Two dimensional cubic B-spline elements in parametric space .....................................66 5-1 Typical non-conforming grid used for analysis...................................... ............... 95 5-2 Triangulation of boundary element for volume integration........................................96 5-3 Representation of band on essential boundaries......................... ....... ........................... 96 5-4 Approximation of natural boundary for computation of load vector..............................97 In Figure 5-11, the approximation for displacements in the internal element is given by u = Nui + Nu + Nu + N4u' and the approximation in the boundary element is given by U = NIu + N2,u + Nu' + N4 . The shape functions of cubic B-splines and the their third derivative is given by 1 1 N (1 -3r+3r2 -r3);N2= (23-15r-3r2 +3r3) 48 48 (5-52) 1 1 N, -(23+15r- 3r2- 3r);N (1+ 3r + 3r2 +r3) 48 48 d3N1 -1 d3N2 3. d3N3 -3. d3N 1 (5-53) dr3 8 dr3 8' dr3 8 dr3 8 Imposing C3 continuity at the intersection of internal element and boundary element results in the following equation and the extrapolation relation is given as: 1 2 3 1 (5-54) Su' = -u' + 4u t 6u' + 4u' Using the above relation, the extrapolation matrix in two dimensions is given as: \Uo ~-1 0 4 0 -6 0 4 01 )T} { o = 00 v uI v, u, v, ui, v, (5-55) Lv, 0 -1 0 4 0 -6 0 41 "2 "' "' " Quadratic B-spline approximation. The quadratic b-spline approximations have C1 continuity. Hence C2 continuity is imposed (Figure 5-12) between the piecewise polynomial segments between internal nodes and external nodes. The extrapolation relation between the external node and internal node is computed as: uz = u 3u + 3u (5-56) The extrapolation in two dimensions is constructed as: { = [ 7 3 03 0 .; v ; v ; (5-57) v 0 1 0 -3 0 3 u 3 S, .11 (Avg: 75%) 1- S p1 Tim 1 - B ti "e i.- Figure 6-48. Contour plot of bending stresses using tetrahedral elements for 3D ribr 1.1.7073E5 *q -, L 5, St, 5.3909E 4 ncrer I Sp Time 1-1.7974E5 Figure 6-48. Contour plot of bending stresses using tetrahedral elements for 3D rib Stress xx 2.0967E5 9.285E4 A 5.3909E4 -2.3973E4 I-.'6.2914E4 -1.7974E5 Figure 6-49. Contour plot of bending stresses using B-spline elements for 3D rib 147 (a) (b) Figure 6-28. Hollow thick cylinder subjected to uniform pressure 1.0742E-5 Figure 6-29. Contour plot of radial displacements using cubic B-spline element model Figure 6-29. Contour plot of radial displacements using cubic B-spline element model The solution structure for adaptive refinement can be extended for performing analysis of 3D problems. The h-, p- and k- adaptivity can be incorporated into this method. B-spline basis provides a convenient methodology for raising the order and continuity of the approximation space, this property can be used to enable k-adaptivity and can also be used for k-version FEM. for analysis, the nodes may not lie on the material interface and hence there is challenge in imposing interface conditions for material discontinuity. In Implicit Boundary Method, this is achieved by construction of a solution structure which ensures that interface conditions are satisfied for displacements. The implicit equations of the inclusion boundaries are used to construct a step function for inclusion such that it is unity inside the inclusion and smoothly blends to zero at the inclusion boundary. Such a function can be used to patch the solution from the inclusion grid and matrix grid to form a solution structure. Figure 3-8 shows the inclusion function in one-dimension and Figure 3-9 shows the inclusion function for an elliptic inclusion in a square matrix. It can be seen that the inclusion function is unity inside the inclusion and transitions to zero at the boundary of the inclusion, the rate with which this function reaches zero is specified by the parameter 6. The expression for inclusion function is similar to the Dirichlet function described earlier in this chapter and shown below 0 0_<0 H"nc(0)= 1 1O' 0<0 1 _>8 This function inherits all the properties of Dirichlet function and is a shifted approximate step function that tends to the Heaviside step function in the limit as 5- 0 but has a value of zero on the boundary unlike in traditional approximations of step functions where it has value of 0.5 Moreover, for allowing normal strain discontinuity, H" nC() is constructed to have non-zero gradient at 0 = 0. At 0 = S this function has Ck 1 continuity. In this work, the step function approximation of HI"C with k= 2 and S 0 is used as Inclusion function with values of S w 10 2 H"f ONg2 0 ~ N g2 0 B 21 B22, = 0 H:' + 0 N .g (5-46) ax2 2 2 1 212 i 1 The strain displacement matrix takes the similar form as explained in section 5.3.3. Hence the stiffness matrix will take the following form: iTC B T C CB1 BT CB2V = [K K, K= [BB B2 C[B B2dV= 2 V K K12 (5-47) B, T BCB, B CB, K2 K,, The computation of stiffness matrix is performed in the similar procedure as explained in detail in section 5.3.3. 5.4.3 Computation of Load Vector The load vector is computed in the similar manner as explained in section 5.3.6. 5.5 Extrapolation of Solution for Elements with Low Volume Fraction In structured grid based methods, the grid elements do not conform to the geometry of the analysis domain and hence there may be few elements in the model with very small solid region. When the volume of the solid region inside the grid elements approach zero, the stiffness matrix corresponding to these elements will have components with magnitudes which are very small in comparison to the internal element stiffness matrices. This leads to ill-conditioning of the system of linear equations and hence poor convergence behavior. In order to avoid the problems associated with grid elements having small volume fractions, an extrapolation scheme is presented in this dissertation. The degrees of freedom corresponding to the elements with small volume fractions are constrained to be an extrapolation of the degrees of freedom from the internal elements. Due to this extrapolation scheme, the degrees of freedom for external nodes will not be arbitrarily large due to near singularity of the - Solid Region Approximation of Boundary Natural Boundary Figure 5-4. Approximation of natural boundary for computation of load vector Figure 5-4. Approximation of natural boundary for computation of load vector X .016E0 "$.9 LI 2.086E0 I.'"hFb l.46H9Ei 1.]10LIPI .499 -1 5.405E-1 .2-2__________ Figure 6-75. Contour plot of o, for local refinement using cubic B-spline elements 10-1 CBS SI 1 1 -- --- -I- Q4 L _ _, ,_ _, _,_ _. ,_ _, ,_ j ^ Log (Num Nodes) Figure 6-76. Convergence plot for L2 error in local grid refinement -I 1-14 - \iLt 3.1 Solution Structure for Enforcing Essential Boundary Conditions The essential boundary conditions are specified as prescribed displacements or prescribed temperatures on the essential boundaries. Let u be a trial function defined over Q c R2 or R3 that must satisfy the boundary condition u = a along F, which is the essential boundary of Q. Here, u corresponds to a displacement component in linear elasticity and temperature in steady state heat transfer. If D, (x) = 0 (x e Q and a : Q -> R ) is the implicit equation of the boundary F,, then the trial function can be defined as u(x) = (,(x)U(x) + a(x) (3-1) The solution structure for the trial function defined in Eq. (3-1) is then guaranteed to satisfy the condition u = a along the boundary defined by the implicit equation I (x) = 0 for any U: Q -> R. The variable part of the solution structure is the function U(x) which can be replaced by a finite dimensional approximate function Uh(x) defined by piece-wise interpolation/approximation within elements of a structured grid. In the context of such solution structures, the function (a (x) is often referred to as a weighting function. Since these functions assist in imposing the Dirichlet boundary conditions, we have referred to these functions as Dirichlet functions or just D-functions. D-function is constructed such that it takes a value of zero and has non-zero gradient on the essential boundary and has a value of one far from the essential boundary. Instead of directly using the implicit equation of the boundary as D- functions, approximate step functions constructed using these implicit equations have been used to construct solution structures in this work. The use of approximate step functions for construction of D-functions localizes them to the boundary elements and hence the stiffness matrix for all the internal matrices is identical as the D-function is unity inside them. In the ..- ............ -- .. .. ..... o ----------i---t-- -s-- S: u=10(1-xl)+5(1-x2) -- +---*----c- --c- *----< -r4 -- f------tr--- .----+--4-----... + --+ - e----- ----- -- -- + --- ... -. 4-4.....4 ..... 4....... + *. *- .-"... ... o -4 4 i ----+t--*t-----t-6- O ~~0* -4- 4 - o 1 (a) Non-zero boundary value specified at the edges I 1i- 5 -- - O ' C) A 0 1 5 1.5 (b) Boundary value function constructed using cubic B-spline basis Figure 3-6. Boundary value function using cubic B-spline basis u: o o 4 06 --Figure 3-13. Step function for refinement for a 2D model OB Os, 4 -02 0 02 Q04 O 88 Figure 3-13. Step function for refinement for a 2D model The components involving only the gradients of Dirichlet function is termed as [c3] and its integral across the band is given as 1 w 4 s1 n2 4 1 & Sc n i 4 C11 +C33 1 s (C12 + C33) CT3 e=% -rh 3( n 3s JoJ 3sh(5- (5-12) V(1 1^ 9(-f)}f _,q f C 4C 41 4 +C134 J(12 J,,33 s s 22 2 T s 33 J s J11J22 1 CS8 32 12 9s, 31 / 2 \ 3 The computation of the integral along the boundary F is performed numerically by Gauss quadrature. This approach leads to accurate and efficient computation of stiffness matrix for boundary elements with essential boundaries. 5.1.4 Computation of Load Vector The right hand side of the modified weak form as shown in Eq(5-1) has two terms, the first term is a volume integration for the body loads and applied non-homogenous essential boundary conditions which is done similar to the computation of stiffness matrix as explained in the previous sections. The second term is a surface integration over the natural boundaries. The boundary is approximated as a set of line segments in 2D as shown in Figure 5-4 and surface triangles in 3D and integration is performed using Gauss Quadrature as shown below: {ft2m = f{N}T {tdF = w, {N((I)}T {t(Sl)} L, 1=1 (5-13) {f}2ml =Z{f} 1=] As the grid is made finer, the approximation of the boundary approaches the exact representation and the errors associated with approximation of boundaries with lines (planes) approach zero. In traditional FEM, the natural boundary is approximated by the surface formed by the nodes on the boundary and associated shape functions. This results in a load vector which has non-zero values corresponding to the nodes which are present on the surface. However, in parameter range for the element is r e [-1,1], the parameter value for first support node is k+1 k+1 r = and for the (k + 1)t support node is r = for the parameterization to be uniform. 2 2 The approximated function shown in Fig 4-1, is a linear combination of the basis functions within each element so that the span of the element e is represented as k k k f (r) = -NJu = Y aIr' (4-1) 1=0 1=0 j=O- In the above equation, u are the support nodal values and it is assumed for convenience that the first support node for element e has index 0. B-spline approximation of order k has Ck-1 continuity at the junction between adjacent elements. This continuity requirement can be stated as f() f, m + ) 0,1m= ...k-1 (4-2) where, e and e+1 are adjacent elements whose derivatives are evaluated at the end point (r=--) and start point (r=1)respectively. m = 0 represents the continuity equation f/(1)= f, (-1) Using Eq. (4-2), the continuity equations can be expressed as k k k k g *] S(j m) i (j-)! ( jm)! (4-3) k j! )J-m 0 tJ- ak+l 1)J k+ = 0 =0 (j m)! Since this equation is valid for arbitrary values of u, i = O.k +1, the coefficients corresponding to each u,i = O..k +1 should vanish. This results in a set of k + 2 linearly independent equations. Since k+2 equations are obtained for each value of m = 0,1,..k -1, a total of k(k + 2) linearly independent equations are obtained. An additional equation stating that also demonstrates that cubic B-spline elements do not undergo shear locking in thin limit, and when the structures are made of thin beams, one element can be used to get acceptable solutions. 6.1.2 Validation of Extrapolation Technique In this example, the extrapolation technique presented in this dissertation is validated. A curved beam as shown in Figure 6-8 is used for analysis. The curved beam is a three-quarter circular beam which is fixed on one end and loaded on the other end. The material of the beam is assumed to be made of steel with a modulus of elasticity of 200GPa and a Poisson's ratio of 0.3. The geometry for this beam is such that there are several elements on the boundary which may have very small support from the nodes lying inside the geometry of the analysis domain and hence the resulting stiffness matrix will be almost singular and can lead to ill-conditioning of the system of linear equations and hence leads to poor convergence behavior. The extrapolation technique presented in this dissertation is used to extrapolate the solution for the external nodes in terms of the internal nodes. Bilinear Lagrange interpolation are used to perform a convergence analysis where a series of models with varying grid densities are analyzed and the error in strain energy is plotted as shown in Figure 6-9. It can be observed that the convergence behavior is not smooth due to near- singular stiffness matrices. Figure 6-10 shows the error in strain energy using the extrapolation scheme. It can be observed that the convergence plot is very smooth and gives an expected behavior. When B-spline elements are used for modeling this example, the quality of solution was not severely affected even when several elements are present on the boundary with small support. This may be because B-spline elements have a wide span as compared to Lagrange elements and hence the contribution to stiffness from a neighboring element leads to non- singular stiffness matrices and hence does not lead to ill-conditioning. Figure 6-9 shows the plot it is to be noted that these grids gl and g2 are conceptual grids and analysis is performed on the grid g. 3.3.1 Solution Structure The approximate step function that is used to blend the solutions from the coarse and the refined grids is constructed by using the implicit equations of the coarse/refined grid interface. This step function for refinement is defined as: 0 0<0 i(20T o0 < < 2\)) 2 Hre' (9) 2 S 20 ( 1k (3-24) 1- 12 1 #_>8 In the above equation, # is the implicit equation of the coarse/refined grid interface, 6 is the range parameter which defines the region in which the step function for refinement varies from zero to one, and k is an integer which defines the order of the step function used. The step 1 S function for refinement is such that it takes on a value of zero at # = 0, at = and one at 2 2 S= S. It has a zero slope and Ck1 continuity at both # = 0 and # = S. Hence, when B-spline basis is used to approximate the solution, an appropriate k needs to be chosen so as to maintain the continuity of the solution across the coarse/refined grid interface. The plot of Href and 1- Href is shown in Figure 3-11, and it can be seen that the above mentioned properties are satisfied by the step function of refinement. The above definition for H're is valid only when the coarse/refined grid interface is defined by a single implicit equation. However, in many practical cases the interface is defined by Xa T 8E1 7.2E1 6.4E1 5.6E1 4.8E1 4El 3.2E 1 F.4E1 1.6E1 8EO) OEO Figure 6-36. Plot of boundary value function NT11 +8.000e+01 +7.594e+01 +7.188e+01 +6.782e+01 +6.376e+01 +5.970e+01 +5.563e+01 +5.157e+01 +4.751e+01 +4.345e+01 +3.939e+01 +3.533e+01 +3.127e+01 Figure 6-37. Contour plot of temperature distribution using FE Q4 elements Therefore the computational domain is divided into three regions as shown in Figure 5-5 and they are explained in detail in the following. 5.3.1 Computation of Stiffness Matrix in the Matrix Region In the matrix region of the microstructure, the degrees of freedom for grid corresponding to matrix are active, therefore the strain-displacement matrix as shown in Eq. (5-23) becomes B = [B1 B2] = [B,1 0] with Hn" = 0 in Eqs. (5-24) and (5-25). The stiffness matrix will take the following form in matrix region: -tB T 1CmatB 0 K= [B1 B2 ]Tmat[B, B2]dV= B 0 V (5-26) Q, Q, I This stiffness matrix is computed by performing the volume integration over the elements in matrix region. The grid element in the matrix region will be an internal element and the volume integration is performed as discussed in the previous section 5.1.1. 5.3.2 Computation of Stiffness Matrix in the Inclusion Region This region corresponds to the inclusion region. In this region the degrees of freedom for grid corresponding to inclusion are active, therefore the strain-displacement matrix as shown in Eq. (5-23) becomes B = [B B2] = [0 B21] with H" = 1 in Eqs. (5-24) and (5-25). The stiffness matrix takes the following form: K= f[B B2 c]Tinc [B l B2]dV= f T V (5-27) QQ, 21, 21 The stiffness matrix is computed by performing the volume integration over the elements in inclusion region. The grid element in this region will be an internal element and the volume integration is performed as discussed in the previous sections 5.1.1. specified at all the support nodes should be on a plane corresponding to the linear distribution. This is illustrated in Figure 3-6, where a simple geometry is used to demonstrate the construction of the boundary value function. The left edge has a constant value of ua = 10 while the right edge has a bilinear variation of u, = 10(1- x) + 5(1- y). The boundary value function is constructed by identifying the elements for the left boundary and specifying the nodal value of 10 for all the support nodes, and then by identifying the boundary elements for right boundary and specifying the nodal values computed using 10(1- x) + 5(1 y) on the corresponding support nodes. The nodes which are not support nodes for the boundary elements are set to zero. The plot for the boundary value function is shown in Figure 3-6 and it can be observed that the boundary value function is continuous over the entire domain and it smoothly transitions from the specified value of 10 from the left boundary to zero in the middle of the domain and then smoothly transitioning from zero to the specified bilinear distribution along the right edge. Even though a simple geometry is used to illustrate the concept of constructing boundary value function, the procedure is quite general. However, care has to be taken to ensure that the grid has sufficient density to be able to construct the appropriate boundary value function. For more general variation of boundary values, the nodal values can be computed by least-square minimization. The rest of the nodes in the grid that are not part of any boundary element can have any arbitrary value and therefore can be set equal to zero. The boundary value function u" contributes to the load computation on the right hand side of Eq. (3-10) for all those elements in which the gradient of u" has non-zero magnitude. 3.2 Solution Structure for Micromechanics In the finite element method, a conforming mesh is created such that the element edges/faces form the material boundary. When structured grids as shown in Figure 3-7 are used CHAPTER 1 INTRODUCTION 1.1 Overview Numerical methods are necessary for solving many problems arising in engineering analysis. Even though analytical solutions are available for several engineering problems with simple geometries, for complicated domains it is difficult to obtain analytical solutions and hence numerical methods are needed to obtain accurate solutions. The Finite Element Method (FEM) is widely used in engineering analysis not only due to its applicability to solve a variety of physical problems but also due to its ability to handle complex geometry and boundary conditions. In FEM, the analysis domain is subdivided into finite elements which approximate the geometry of the analysis domain. The field variable is approximated on these finite elements using piecewise continuous interpolations. Even though FEM is a well established technique, it has certain limitations. Approximation schemes used in FEM are typically Co continuous and hence there is a discontinuity of gradients across the elements. The lack of gradient continuity in FEM necessitates the smoothing of results. There are a number of approximation schemes like B-splines, moving least squares and Hermite interpolations which provide at least C1 continuity in the domain of interest. B-spline basis functions have compact support like traditional finite element shape functions and lead to banded stiffness matrices. They form partition of unity and hence can exactly represent the state of constant field in the domain of interest. The B-spline approximations have good reproducing properties and when a cubic B-spline is used, it can represent all polynomials up to third order exactly. B-spline approximation provides high order of continuity and is capable of providing accurate solutions with continuous gradients through out the domain. The uniform B-spline basis functions are polynomials and accurate integration for stiffness matrices and load vectors can be The weak form takes the following form when the above approximations are incorporated into Eq. (3-13) J[B] k[B] k {TdQ l + I[N]h[N] {T }dF =- JB k[B]{Ta}dQ, Fa rC a, (3-14) + [N]Tfd + f I [N]Tq,od- f [N] h(Ta- T")dF ,e Th Fr In the above expression, all the known quantities are moved to the right hand side and the unknown quantities are placed in the left hand side. The first term on the right hand side is the contribution to the load due to applied essential boundary conditions, the second term is the load due to the heat generation, the third term is load due to applied heat flux and the last term is load due to convection boundary conditions. The first term on the left hand side represents the stiffness matrix while the second term represents the stiffness due to convection boundary condition. 3.1.3 Dirichlet Function The Dirichlet function D is a function which vanishes on all essential boundaries where ith component of displacement (or temperature) is prescribed. In addition, DZ should have a non- vanishing gradient on the essential boundaries to ensure that the gradients of displacements are not constrained. Furthermore, it is necessary to ensure that D, is non-zero inside the analysis domain so that the solution is not constrained anywhere within the domain of analysis. Implicit equation of the curve or surface representing the boundary can be used for this purpose. However if multiple curves are involved then it is necessary to combine them into a single equation. Some authors have proposed the use of R-functions to construct the Boolean combination of implicit functions (Shapiro, 1998 ; Rvachev et. al., 2001). Hollig (2003) recommends constructing a weight function using distance functions as: w(x) = 1 max(0,1 dist(x, OD) / )' This function Figure 6-46. Three dimensional rib model used for analysis 0.245 0.24 0.235 2) a) W 0.23 0.225 Abaqus Linear Tetra IBFEM H8 IBFEMCBS 0.21502 102 103 Log (Num Nodes) Figure 6-47. Convergence plot for strain energy for rib model the implicit boundary method, nodes are not guaranteed to be present on the surface, therefore the load vector for a boundary element with natural boundary will have non-zero values for all the nodes corresponding to this element. The computed load vector is a work-equivalent load similar to the consistent load vector used in traditional FEM. When the solid region inside the element is very small, its contribution towards the stiffness matrix will be negligible and may lead to singular stiffness matrices. To overcome this problem, the elements with low volume fraction, i.e. r <5%, = vol(Solid) are identified and the vol(Element) nodal displacements in these elements are constrained to be an extrapolation of the displacements from the neighboring internal elements. The details of this extrapolation approach are presented in section 5.5. 5.2 Numerical Implementations for Steady State Heat Transfer The modified weak form for the solution structure for steady state heat transfer was discussed in a previous chapter and was shown to be: f I[B]T k [B] {T }d, + f J[N] h[N] {J}dr = F r (5-14) -Z R] k[B]{T'd + [Nrde [N]Tfdo, + [N] q dF i[N] h(Ta T.")dF ,e Fe Th rc In this equation, the left hand side corresponds to stiffness matrix and the right hand side corresponds to load vector. The left hand side has two terms the first term is stiffness matrix and the second term is a contribution to the stiffness matrix due to convection boundary condition. The right hand side has four terms, the first term is load due to applied temperature on essential boundary condition, the second term is load due to heat generation in the body, the third term is due to applied heat flux on natural boundary and the last term is load due to convection on convection boundary. NE NE J{8X'j}J[B, BT]2[C][B, B2]{X'jdQ {Xe}T {Fe}+ el NE e-1 (3-29a) NBE I X' Tf {N N2I {t} dl e=l pt {Fe = N N2)T {b} d2- {Bar Ba2 T {a d (3-29b) Q, Q, In this expression {Xe = {ug ug2 and the summation is over the number of elements (NE) and the number of boundary elements (NBE) in the grid used for analysis. As in the finite element method these element matrices can be assembled to form global equations, where the left hand side terms form a stiffness matrix and the right hand side a force vector. The Dirichlet function and the boundary value function are constructed as explained in sections 3.1.3 and 3.1.4. In this chapter, modified weak forms for the solution structures corresponding to enforcing essential boundary conditions, treatment of material discontinuity and local refinement of grid are presented in detail. These modified weak forms are used in finite element implementation for solving several engineering boundary value problems. In the next chapter, a detail discussion on the construction ofB-spline based finite elements is presented. 5-5 Division of computational domain into three regions............................................ 98 5-6 Undeformed and deformed configuration of a set of 4 RVEs ........................ ..........98 5-7 Representative volume element used for computation of effective properties.................99 5-8 Division of computation domain into various regions...........................................99 5-9 Typical analysis grid with elements having low volume fraction .................................100 5-10 Extrapolation of external node for linear Lagrange element ...................................100 5-11 Extrapolation of external nodes for cubic B-spline element...................................100 5-12 Extrapolation of external nodes for quadratic B-spline element .................................... 101 5-13 Extapolation for boundary elem ent...................................................................... ...... 101 6-1 Plate subjected to uniaxial loading ....................................................................... .... 121 6-2 Non-conforming grid using one cubic B-spline element................................................121 6-3 P lot of D irichlet functions............................ .................................................. .......... 122 6-4 Effect of range parameter on the quality of solution .............. ....... ............... 123 6-5 Contour plot of o-x for one element model using cubic B-splines ...............................123 6-6 Thin cantilever beam subjected to end loading............................................................124 6-7 Plot of error in the tip deflection with varying aspect ratio............................................124 6-8 Curved beam subjected to uniform loading ............................... ................................. 125 6-9 Plot of error in strain energy for curved beam (without extrapolation) ...........................125 6-10 Plot of error in strain energy for curved beam using extrapolation technique...............126 6-11 Contour plot of o- for Q4 elements (without extrapolation)........................... ........126 6-12 Contour plot of o- for Q4 elements using extrapolation technique .............................127 6-13 Contour plot of o-, using cubic B-spline elements .................................. ............... 128 6-14 Thick cantilever beam used for convergence analysis.........................................128 6-15 Typical non-conforming grid using cubic B-spline elements............... ..................129 14l.Eill L il3Ell 2L.2iTII 1I21-SE11l .L5!1I1 L084ELL 1.7253Eli I : 40JIH ]:l .(3TELO 16I !5~SB Figure 6-66. Plot of stresses for cases 1 and 2 for 3D model (A) Plot of 1,, normal stress for Case 1 and (B) Plot of 72, normal stress for Case 2 Sim. Hi 7IEll l.086EII11 1.72T3EII I'rIaI' I i--E II 5i l l .l l S.255jE10 |.412UE10 S., 6 FLI I.' 4 6-11E 10 4*Ij!E10 4 n1(EIP 3.5163E10 9i 4EID Figure 6-67. Plot of stresses for cases 3 and 4 for 3D model (A) Plot of 33 normal stress for Case 1 (B) Plot of a23 shear stress for Case 4 . 1 IEL 1 1s.i I[, 2.4234EI .6949El6 I~^~ 6.0657E10 I11 l"[1 14234E10 I9'4FIO B Figure 6-68. Plot of stresses for cases 5 and 6 for 3D model (A) Plot of -31 shear stress for Case 5 and (B) Plot of z12 shear stress for Case 6 156 refinement is obtained when the order of continuity of the basis functions is increased. In meshless methods, refinement is achieved by inserting nodes in the regions of stress concentrations. Hierarchical mesh refinement is studied by Krysl et. al. (2003) in the finite element method framework. In this method, the nodal basis functions on the coarse mesh are split into the basis functions of same order in the finer mesh. They have provided algorithms for refinement as well as un-refinement for finite element meshes. Hollig (2003) also presented a hierarchical refinement of B-spline basis functions, where the B-spline basis are split into the B-spline basis of same order but on a finer grid with half the size of the coarse grid. The hierarchical refinement provides a framework for getting a multi-level resolution, where the basis functions are split in a recursive manner. In this dissertation, a technique for local grid refinement is presented based on construction of solution structure using approximate step functions. Stress xx 1.248E5 9.966E4 2",E .I-_F I 4.937E4 2.6".1- 4 -7.635E4 1-1.015E5 -1.266E5 Figure 6-44. Contour plot of bending stresses using H8 elements for 3D clamp Stress xx 1.233E5 9.846E4 k6 -- 4.888E4 .L ..1 LF S-7.506E4 5 9.985E4 -1.246E5 Figure 6-45. Contour plot of bending stresses using cubic B-spline elements for 3D clamp expected tip deflection. This example demonstrates that the B-spline elements do not undergo volumetric shear locking. 6.1.4 Square Plate with a Circular Hole The performance of quadratic and cubic B-spline elements in regions of stress concentration is studied in this example. An infinite plate with a central circular hole as shown in Figure 6-22 is considered for this analysis. The plate is subjected to uniform unit traction (o = 1) along the x axis at infinity. The exact solution for this problem is shown in Eq. (6-5) (Sokolnikoff, 1956). Symmetry of the problem is exploited to model only a quadrant of the analysis domain. The radius of the hole is chosen to be 1 (a = 1) and length and width of the plate is chosen to be 4 (1 = w = 4a). The exact solution is used to compute the tractions on the traction boundaries (right and top boundaries) and symmetry boundary conditions are applied on the bottom and left edges while the inner edge is traction free. Ia2 3 3a4 x =0 1-- -cocos20+cos40 + 4 cos40 (6-5a) r 2 ) 2r a 21 3a4 c L =r -- ,sin20+i4O sin440 (6-5b) r 2 2 2r4 S= cos20- cos40 cos40 (6-5c) y r 2 2r4 -- /-l v lFv 1 2 cosO+--- cos 3- -cos30 l+v 1 2 a 1 a2 1 a4 u=- a -- rcos+- -cos+--cos30---cos30 (6-5d) E 1+v +v r 2 r 2 r3 l+v -v -v a2 a2 a4 (6-e) v=-- cL rsinO-- -sinO+--sin3O-- sin30 (6-5e) E 1+v l+v r 2 r 2 r3 The typical analysis models used for convergence analysis along with the contour plot of stresses are shown in Figure 6-23 and Figure 6-24. The error analysis for this problem is Polynomial segment(span) C2Continuity l Approximal Nodal Values 1 E2 E1 3 E2 4 E3 5 Nodes Elements Figure 4-1. Approximation using ID cubic B-splines I N2 I I I I 0.5 I NI NI N3 I I I I 0 ---- ---'- ----- 0 1 2 mi 3 Element Figure 4-2. Shape functions for ID quadratic B-spline element 10.5 1I 3 I I I I I I 0.5i 1 2 N1 3 o 0 r-1 r-1 Element Figure 4-3. Shape functions for ID cubic B-spline element Function 6 6 4 The comparison of ux along the line is shown in Figure 6-60. It can be observed that the displacement solution obtained by IBFEM matches exactly with the FEM solution. A slope discontinuity can be observed at the inclusion boundary. Figure 6-61 shows the line plots o- and o- along the line AB and Figure 6-62 shows the comparison of Ex and c It can be observed that the solution for both stresses and strains match exactly with the FEM solution. The discontinuity in stress and strain can be observed at the inclusion interface. 6.2.2 Determination of Effective Properties of a Fiber Reinforced Composite 2D model. In this example a two-dimensional model is used for computation of effective properties of a Boron-Aluminum fiber reinforced composite with volume fraction of 0.47. For determining the effective properties in the transverse directions, a plane strain assumption is used. In order to determine the homogenized material properties in the longitudinal direction, a generalized plane strain model (Li and Lim, 2005) is used. The homogenized shear properties in the longitudinal direction are determined by using a longitudinal shear loading model with one degree of freedom in the longitudinal direction is used. This assumption results in non-zero shear strains in z-y and z-x directions only while all other strain and stress components are identically zero. The results obtained by IBFEM are compared with the results available in literature (Sun and Vaidya, 1996, Sun and Chen, 1990, Chamis, 1984 and Whitney, 1966). The periodic boundary conditions shown in Table 5-1 are used for analysis. The average stresses for each case are tabulated in Table 6-1. Figures 6-63, 6-64 and 6.65 show the stress distribution for various cases with the deformed shape. The effective properties are computed by using the average stresses and strains and are compared in Table 6-2. It can be observed that the effective strains in the longitudinal directions are specified. A constant strain in the x3 is assumed and the strains and stresses are expressed as follows 8, 1 C11 ou E2 an2 2 C21 2 2 and I E3 03 C31 L712 u, 3'12 0 ax, ax, Ox2 1Xl Using the constant strain in 3-direction C2 12 22 0 C2 + 23 T12 0 0 C66 12 and-3 [C13 C23 0] e2 +C33 3 12/2 0 Ez C66 E3 C66 /12 i, the stresses are written as follows: 3 = [C] {} + Ck3 The principle of virtual work as shown in Eq. (3.3) is modified by choosing the virtual strain vector as 6E = {(E, 8E2 712 } to get the following equation (Eq. 5-39). In this expression, the last term is the load due to strain in the third direction. The rest of the formulation is exactly similar to plane strain formulation. {&6} T{( } dQ n fu}T {t}dF + f{6u}T {b} d- {6} CkE3 (5-39) n F, Q Longitudinal shear loading. The homogenized shear properties in the longitudinal direction are computed by using a longitudinal shear loading, where there is only one degree of freedom which is the displacement in the 3rd-direction. This assumption results in non-zero shear strains in 3-2 or 3-1 directions only while all other strain and stress components are identically zero. The equations for shear strains and shear stress are given by the following expression: (5-37) (5-38) Matrix Inclusion Grid 1 Figure 3-7. Representation of grids for material discontinuity Zero Slope 1 HInc I-H HInc Non t t t t 1 Matrix Boundary( = o) Inclusion Boundary (0 = 0) v Figure 3-8. One dimensional inclusion function for treatment of material discontinuity Grid 2 latrix stiffness matrix of partial elements. An added advantage of this procedure is that the total number of degrees of freedom is reduced and hence less computational time for solution of linear equations. The procedure for this extrapolation strategy is explained with reference to Figure 5-9.The external nodes in the grid are denoted by n and the internal nodes are represented by n'. The displacement of an external node is associated with the displacement of a set of internal nodes and the number of inner nodes that are used will depend on the approximation scheme used for the grid variable. For example, for linear Lagrange elements, two internal nodes are associated for a given external node and for cubic B-spline elements four internal nodes are associated with a given external node. The extrapolation (constraining of external nodes using internal nodes) can be written as follows Un = I(u,"...,Un, ) (5-48) In the above expression, I is a function which maps the degrees of freedom of a set of internal nodes to the displacement of an external node. This function can be expressed as an extrapolation matrix and the transformation can be represented as a matrix multiplication: U dxl = [C]ndxm {U, } (5-49) In the above equation m is the number of internal nodes used to constrain the displacements of the outer nodes and the dimension of the problem as nd and the size of [C] is nd m. L IST O F R E F E R E N C E S .............................................................................. ..........................169 B IO G R A PH IC A L SK E T C H ............................................................................... ............... ..... 175 Here {6S} is the virtual strain and {6u} is virtual displacement. A trial solution structure for the displacement field can be constructed such that the Dirichlet boundary conditions are satisfied u}= [D]{ug}+ {u) =u)+ {u} (3-4) In this solution structure, {ug} is a grid variable represented by piece-wise approximation over the grid, and {u0} is the boundary value function which is a vector field constructed such that it takes on the prescribed values of displacements at boundaries. The variable part of the solution structure is {u"} = [D] {ug which satisfies the homogeneous boundary conditions. [D] = diag(D,.., D ) is a diagonal matrix whose components DZ are Dirichlet functions that vanish on boundaries on which the ith component of displacement is specified, nd is the dimension of the problem. The same D-functions used to construct the trial solution are also used to construct the test functions. {(u} = [D]S6ug (3-5) Using these solution structures, the strains and stresses can be expressed as the sum of a homogeneous part (E s, u ) and a boundary value part ( e, o ). S (Usj +us, )+ (u +ua = Es+ (3-6) 2 \ 1 2 J 1 o, = Csk kl = Uk = ( + a o-y (3-7) Substituting ({c =( }+ {a in the weak form, we can rewrite it as J~S)} {sdQ= { Su}T {t}dF + JSu}T {b}dQ- J{6} {a ) d (3-8) Q F, Q Q imposed by placing nodes on the boundary so that the shape functions corresponding to these nodes emulate Kronecker-delta behavior. Shepard functions are used as partition of unity functions and quartic splines are used for weight functions. Method of HP Clouds is a partition of unity based meshless method developed by Duarte and Oden (Duarte and Oden 1996 ; Duarte et. al. 1996 ; Liszka et. al. 1996). HP-Cloud functions are created by multiplying the partition of unity functions with local basis functions. Shepard functions or MLS functions are used as partition of unity. Anisotropic approximation can be incorporated by using the basis function which has higher order in one direction and lower order in other direction. P-refinement is done by increasing the size of the local domain; hence more number of points is included in the sub-domain and also by increasing the polynomial order in the basis function. Boundary conditions are applied by modified variational principle. Natural element method developed by Sukumar et. al.(1998) uses the natural neighbor interpolants which are based on Voronoi tessellation of the nodes. These interpolants are smooth over the whole domain excepting at the nodes where they are continuous. The essential boundary conditions are applied on nodes for convex domains. 2.5 Techniques for Determination of Effective Properties of Composite Microstructures The finite element method (FEM) has been used extensively in micromechanical analysis of composite microstructure to determine effective properties (Sun and Vaidya; 1996 ; Taliercio, 2005 ; Xia et. al. ; 2003 ; Marrey and Sankar, 1997 ; Zhu et. al. 1998). The effective elastic properties are determined by performing six different analyses on the representative volume element (RVE) with unit non-zero strain component in each analysis and by imposing appropriate periodic boundary conditions on the boundaries of the RVE. The mesh is generated in such a way that the element boundaries conform to the matrix-inclusion interfaces. However, for complicated microstructures, automatic mesh generation becomes difficult and special Solid Region D() = D(,l2)n3 = (D(q) + -D(() D()D(q2))(D(03)) Figure 3-5. A CSG tree representing Dirichlet function for multiple boundaries Numerical examples are presented to demonstrate the performance of B-spline elements. The results are compared with analytical solutions as well as traditional finite element solutions to demonstrate the ability of B-spline elements to represent continuous stress and strain through out analysis domain. Convergence studies show that B-spline elements can provide accurate solutions for many engineering problems with fewer numbers of elements and nodes as compared to traditional FEM. Solution structure for treatment of material boundary is validated by performing a convergence analysis on a problem involving circular inclusion in a square matrix and by determining effective properties of fiber reinforced composite. Solution structure for local grid refinement is validated by analyzing classical stress concentration problems. -A- ---15 ---- -B Tx=1MPa Figure 6-1. Plate subjected to uniaxial loading Nodes Solid Boundary -Element Boundary Figure 6-2. Non-conforming grid using one cubic B-spline element 121 1.0 m Z "11 1 v21 31 0 0 0 E1 E2 E3 12 V32 0 0 0 1' E1 E2 E3 011 _22 V13 V23 1 0 0 0 o 22 E33 E1 E2 E3 22 [ t (5-35) 723 0 0 0 1 0 0 r23 731 G23 131 7122 0 0 0 0 1 0 1"12 G31 0 0 0 0 0 G12 In two-dimensional problems, the effective properties in the transverse directions are determined by analyzing the model in plane strain conditions. The longitudinal modulus is computed by analyzing the model in generalized plane strain conditions. The longitudinal shear moduli are computed by analyzing the model in longitudinal shear loading conditions. Each of these cases are explained in brief in the following sub-sections. Plane strain. In plane strain formulation, the strains in the x3 direction (longitudinal direction) are assumed to be zero and the stress-strain relation is constructed as follows: r; 1-v V 0 2 (+v)(l 2v) 1 v 0 2 (5-36) 1"12 0 0 1-2v 712 2 This stress-strain relation is used in the modified weak form shown in Eq. (3.3) to solve the boundary value problem with specified strain loading and periodic boundary conditions to compute the transverse material properties. Generalized plane strain. The homogenized material properties in the longitudinal direction are computed by using a generalized plane strain model (Li and Lim, 2005), where the Figure 6-8. Curved beam subjected to uniform loading IBFEM Q4 IBFEMCBS -10 10 Log (Num Nodes) Figure 6-9. Plot of error in strain energy for curved beam (without extrapolation) instability (tensile instability). The behavior of solution depends on the proper distribution of nodes. So there is always an effort needed to place the nodes properly in order to avoid instability. Element Free Galerkin Method was developed by Belytschko and his coworkers (Belytschko et. al., 1994 ; Belytschko et. al., 1996 ; Belytschko et. al.; 1996, Dolbow and Belytschko, 1999 ; Krongauz and Belytschko, 1996 ; Krongauz and Belytschko; 1997 ; Krysl and Belytschko 1997). This method is based on moving least square approximation and uses global Galerkin method for discretization. Quadrature was performed on the background cells. Boundary conditions were applied using Lagrange multipliers or by modification of variational principle, where Lagrange multipliers were replaced by traction forces. The former method leads to an increased number of algebraic equations and the stiffness matrix is not positive definite. The latter method reduces the accuracy of the results and is not as convenient as imposing essential boundary conditions directly. A new method for imposing essential boundary conditions by using a string of Finite Elements around the boundary of the domain was later incorporated into this method. Meshless Local Petrov-Galerkin (MLPG) was introduced as a truly meshless method by Atluri and his coworkers (Atluri et al. 1999 ; Atluri and Zhu 2000a, 2000b). This method uses moving least square approximation as the interpolation functions. It uses the Local Petrov- Galerkin weak form which is developed in the local sub-domain. The trial and test spaces are different; the test functions vanish on the boundary of the sub-domain, as long as these boundaries are internal to the problem domain. If the sub-domain intersects the boundary of global domain, then the boundary condition is applied on this local boundary by penalty method. The stiffness matrix obtained is banded but not symmetric. The boundary conditions are applied coordinate system (4, q) with origin at node 11 for Q1. In order to use a parametric coordinate (r, s) with origin at the center of the element, the (4, q) coordinate system needs to be translated for each quadrant as Q1: 4=r-1;u=s-1 Q2: =r+l;q=s-1 Q3: r=r-1; r=s+l Q4: =r+l;r/=s+1 ( The basis functions at any given point (r, s) are computed by first determining which quadrant the point lies in and then by determining the corresponding support nodes and basis functions. As an example, the basis functions for a point in Q1 are ND (r, s) = N(r )NU (s 1); N2D (r, s) = N (r 1)N (s 1) N2D (r, s)= N, (r 1)N (s 1); N (r, s) = N (r 1)N (s 1) N2D(r,s)= N,(r 1)N (s -1);N (r,s)= N3(r 1)N(s 1) (4-9) N2 (r,s)= Nl (r 1)N3(s -1); N12(r,s)= N (r 1)N, (s -1) ND(r,s)=N3(r-1)N3(s- 1);ND (r,s)= ND(r,s)= 0 ND (r,s)= ND (r,s)= N2 (r,s)= N, (r,s)= N (r,s) = 0 N,2D(r,s) refers to bi-quadratic basis function corresponding to ith node of the element and N, (r) represents the one-dimensional quadratic basis function as defined in Eq. (4-4). Three dimensional quadratic B-spline element can be constructed in a similar manner and has eight octants which are assembled to construct the 64-node element. The basis functions are constructed using the product of ID basis functions and in any given octant, only 27 basis functions are active. The mapping for the geometry is done as shown in Eq. (4-7). Cubic B-Spline elements. The two dimensional cubic B-spline elements are constructed by product of one-dimensional cubic B-splines. In this case, there are 16 support nodes as shown in Figure 4-6 and there is a node at each vertex of the element ([-1,1] x [-1,1] in parametric space). In this section, a discussion on the numerical implementations for the computation of the terms involved in this weak form is presented. The left hand side of the weak form corresponds to the stiffness matrix while the right hand side corresponds to the load vector. The solution structure for local grid refinement is such that the solution from the coarse grid inside the refined grid is zero due to the step function for refinement. Therefore the nodal values of the coarse grid inside the refined are set to zero, provided these nodes are outside the influence ofHrf This leads to the reduction in the total number of degrees of freedom for the analysis model. The computational domain can be divided into three regions as shown in Figure 5-8. It can be observed that the solution is uI outside the refined grid and transitions from uI to u2 inside the refined grid. It can also be observed that in the regions where step function for refinement is active, the Dirichlet functions are inactive and vice versa. This is because of the fact the step functions for refinement are active (varying from zero to one) only on the interface of the edges which are shared between the refined and unrefined grid elements and the Dirichlet functions are active only on essential boundaries and the essential boundaries cannot form the interface edges. The computation of stiffness matrix is explained in detail in the next section. 5.4.1 Computation of Stiffness Matrix in Regions 1 and 2 The step function of refinement is not active (H"r = 0) in Region 1 and the solution structure becomes u = Dugl + ual and hence the strain displacement matrix is given as follows: dN dD ON, OD cx 'xax, D + 0 ;[BX]32m 0 0 =[B]320 0 0 (5-42) 20 0 D + N, ,- D + N, a2i2 s,1 1x This matrix is similar to the one explained in section 5.1.3 and the stiffness matrix is given as: Varying Hinc Region ,.. Inclusion Region j Matrix Region Figure 5-5. Division of computational domain into three regions F2 j F3 I F11 F31 Fl| I Igure 5-6. Undeformed and deformed configuration of a set of 4 RVEs I _F4 I Figure 5-6. Undeformed and deformed configuration of a set of 4 RVEs Str 3.8702EXII 3.2872Ell .Z154GI1 .953El1' 2. E 1 .I6.'E I I 1.5114r11 A i_2-,it 3ISIEII _11 4 lll 7.21~111 ga B Figure 6-63. Plot of stresses for cases 1 and 2 (A) Plot of a,, normal stress Case 1 and (B) Plot of ,22 normal stress for Case 2 7A781E0LI 6.0 18t10 2.-199E0L S1.4402EO1 .1044EE n S%h.- 5.1312E10Q ."S 1110 I|.B5i6EllI 3.326E10 I :49E:: B Figure 6-64. Plot of stresses for cases 3 and 4 (A) Plot of ,33 normal stress for Case 1 and (B) Plot of ,23 shear stress for Case 4 Sm.. ' hIT E |WIEI 7,47SiE30 d 411F 1.1 I " I.O IFIIrl 13599EIO 2.M5PE1B I L640ZEIg S.1044E9 87E 10 7.4781E10 6.4flss3l 2.018'f0l 2 I6lE I0 S140. 10E1O 9.1044E9 Figure 6-65. Plot of stresses for cases 5 and 6 (A) Plot of z31 shear stress for Case 5 and (B) Plot of ,12 shear stress for Case 6 Uno ONNOMMISM 110500: :00 H : sINNUIP I.M.: a.: ::.a:: REMO, I -I a, Moan MM II.MOVIN "I' I U.M a. .00. AM.- I- a NE, manHIR, i: Masson -7-lummommmm Mmfflm:am'V Congo :::: MMMMIII'' 'Immonar, IIIIINNAMEW now, : IS 00 ma In, no HHU, 12 :::::am MEN Now. .... USMIL MWOUNIN: is :00 off so no s no Is no 16:111 MUNK-Momma-HOUNN M11"OFIff HHIff HUH :::: Momommommor AU Mir IIIHMMMHUPO- 1 :28 T man Hir IMOV 11 L El .0 MEE 2:0:111-! Mae L ARMS! Room .RNE 011006 UNIMEMEMEMILL AMMER N ::: 601110:21111111"IN: "MORE: mosomm Immmommom HOWNUMMISIC ComomHOM: asommossm! I'MUMNON MUNFI, -10 :901151. 11 IMSEE MMNNIAUI ANIU.N.86. 'Emma- -;Wmmm' aa a 'a"HIPIL ::EMmM"6=e'IMRMMMM HHHO"002 NAMMEME "I'll `7 'I"' 13 p on: go 11 11 2. an I IF EMENS6 IUH: :::::::a.- -.11MMIN HHUMN: 20HUNN: MOUNMEMEMN; 4mmommommom using penalty method. Another advantage of this procedure over FEM is that smoothing techniques are not required for computation of derivative. There is no need of special integration schemes, and there is flexibility in choosing the size and shape of the local domain, which can lead to convenient methods for integration in non-linear problems. The drawbacks of this method are that the stiffness matrices generated by this approach are not symmetric positive definite and penalty method imposes boundary conditions approximately. Another drawback is that we have to make sure that the union of all the local domains should cover the global boundary. Partition of unity based methods (PUFEM) form a class of meshless methods which use partition of unity functions to blend the local approximation functions which represent the local behavior (Melenk and Babuska, 1996). The main advantage of PUFEM is that we can use the knowledge of the local behavior of solutions efficiently to get good numerical approximation. This method also yields good results in the cases where FEM gives poor results or when FEM is computationally very expensive. One drawback of this approach is that the stiffness matrix obtained from this procedure can be singular and the basis function needs to be modified appropriately or special techniques are needed to solve the singular system. Generalized Finite Element Method was developed as a hybrid of Finite Element Method and Partition of Unity (Babuska et. al. 2002 ; Duarte et. al. 2000). The desirable features like inter-element continuity, local approximability and Kronecker-delta property from FEM are retained in GFEM. Method of finite spheres (De and Bathe 2000 ; De and Bathe 2001) was developed as meshless method which is based on partition of unity; this approach does not use background cells for numerical quadrature. It uses compactly supported partition of unity functions and Bubnov-Galerkin method as weighted residual scheme. Integration of weak form is performed on spheres (circles in 2D), using specialized cubature rules. Essential boundary conditions are rU3 732 a 32=3 G32 OJ{732{ (5-40) 73 J 31 0 CG31 jy31 The principle of virtual work as shown in Eq. 3-3 is modified by choosing the virtual strain vector as 6b = {Y32 731 } This formulation is used to specify unit shear strains and hence the corresponding shear moduli are computed. 5.3.6 Computation of Load Vector The right hand side of the modified weak form Eq. (5-20) comprises of load terms. The first term is load due to body load and the volume integration is performed as explained in section 5.1.1 and 5.1.2. The second term is the load term due to applied traction. The shape function vector will be either (N1 0) or (0 N2} depending on which region the natural boundary condition is applied. The computation of load vector is performed in a similar manner as explained in section 5.1.4. 5.4 Numerical Implementations for Local Grid Refinement The modified weak form for the solution structure corresponding to local grid refinement was discussed in a previous chapter and was shown to be: NE NE {(X'} JI[B, B2] [C][B, B22] X'ed S6X'e} {FT el a e=l (5-41a) NBE X {e}'T IN1 N2} {t}dF e=l pt {F'e} N1 N 2} b} dQ- (Ba, Ba2 T a d (5-41b) X u u2 (5-41c) {Xe = Inu1 ug2}T (5-41c) |

Full Text |

PAGE 1 FINITE ELEMENT ANALYSIS USING UNIFORM B-SPLINE APPROXIMATION AND IMPLICIT BOUNDARY METHOD By RAVI KUMAR BURLA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2008 1 PAGE 2 2008 Ravi Kumar Burla 2 PAGE 3 To my parents 3 PAGE 4 ACKNOWLEDGMENTS I would like to express my sincere gratitude to my advisor and chairman of my supervisory committee, Dr. Ashok V. Kumar, for his guidance, encouragement, enthusiasm and constant support throughout my doctoral research. I would like to thank him for the numerous insights he provided during every stage of the research. Without his assistance it would not have been possible to complete this dissertation. I would also like to thank the members of my advisory committee, Dr. Bhavani V. Sankar, Dr. Nam-Ho Kim and Dr. Jorg Peters. I am grateful for their willingness to serve on my committee, for providing help whenever required, for involvement and valuable suggestions in my oral qualifying examination, and for reviewing this dissertation. My colleagues at Design and Rapid Prototyping Laboratory at the University of Florida deserve special thanks for their help and many fruitful discussions. I would especially thank Sanjeev Padmanabhan, Greg Freeman and Sung Uk-Zhang. I would like to thank my parents and my brothers for their constant love and support without which this would not have been possible. 4 PAGE 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ...............................................................................................................4 LIST OF TABLES ...........................................................................................................................8 LIST OF FIGURES .........................................................................................................................9 ABSTRACT ...................................................................................................................................14 CHAPTER 1 INTRODUCTION..................................................................................................................16 1.1 Overview...........................................................................................................................16 1.2 Goals and Objectives........................................................................................................19 1.3 Outline..............................................................................................................................19 2 ANALYSIS TECHNIQUES FOR BOUNDARY VALUE PROBLMS................................21 2.1 Finite Element Method.....................................................................................................21 2.2 B-Spline Approximation Based Methods.........................................................................22 2.3 Structured Grid Methods..................................................................................................23 2.4 Meshless Methods............................................................................................................24 2.5 Techniques for Determination of Effective Properties of Composite Microstructures....27 2.6 Techniques for Local Refinement....................................................................................28 3 SOLUTION STRUCTURES FOR IMPLICIT BOUNDARY METHOD.............................30 3.1 Solution Structure for Enforcing Essential Boundary Conditions....................................31 3.1.1 Modified Weak Form for Linear Elasticity............................................................32 3.1.2 Modified Weak Form for Steady State Heat Transfer............................................34 3.1.3 Dirichlet Function...................................................................................................36 3.1.4 Boundary Value Function.......................................................................................38 3.2 Solution Structure for Micromechanics............................................................................40 3.3 Solution Structure for Grid Refinement...........................................................................44 3.3.1 Solution Structure...................................................................................................45 3.3.2 Modified Weak Form.............................................................................................47 4 B-SPLINE FINITE ELEMENT FORMULATION...............................................................57 4.1 One Dimensional B-Spline Elements...............................................................................58 4.2 Two and Three Dimensional B-Spline Elements.............................................................60 5 PAGE 6 5 NUMERICAL IMPLEMENTATIONS..................................................................................67 5.1 Numerical Implementations for Linear Elastic Problems................................................67 5.1.1 Computation of Stiffness Matrix for Internal Elements.........................................68 5.1.2 Computation of Stiffness Matrix for Boundary Elements......................................68 5.1.3 Computation of Stiffness Matrix for Elements with Essential Boundaries............69 5.1.4 Computation of Load Vector..................................................................................73 5.2 Numerical Implementations for Steady State Heat Transfer............................................74 5.3 Numerical Implementations for Micromechanics............................................................76 5.3.1 Computation of Stiffness Matrix in the Matrix Region..........................................79 5.3.2 Computation of Stiffness Matrix in the Inclusion Region......................................79 5.3.3 Computation of Stiffness Matrix in the Interface Region......................................80 5.3.4 Periodic Boundary Conditions...............................................................................82 5.3.5 Computation of Effective Properties......................................................................84 5.3.6 Computation of Load Vector..................................................................................87 5.4 Numerical Implementations for Local Grid Refinement..................................................87 5.4.1 Computation of Stiffness Matrix in Regions 1 and 2.............................................88 5.4.2 Computation of Stiffness Matrix in Region 3........................................................89 5.4.3 Computation of Load Vector..................................................................................90 5.5 Extrapolation of Solution for Elements with Low Volume Fraction...............................90 5.5.1 Construction of Extrapolation Matrix....................................................................92 5.5.2 Construction of Extrapolation Matrix for a Boundary Element.............................94 5.5.3 Computation of Stiffness Matrix............................................................................94 6 ANALYSIS AND RESULTS...............................................................................................103 6.1 Validation of Solution Structure for Essential Boundary Conditions............................103 6.1.1 Patch Test.............................................................................................................103 6.1.2 Validation of Extrapolation Technique................................................................106 6.1.3 Convergence Analysis of a Cantilever Beam.......................................................107 6.1.4 Square Plate with a Circular Hole........................................................................110 6.1.5 Thick Cylinder with Internal Pressure..................................................................111 6.1.6 Steady State Heat Transfer Analysis....................................................................112 6.1.7 Analysis of a 3D Clamp.......................................................................................113 6.1.8 Analysis of a 3D Rib Structure.............................................................................114 6.1.9 Analysis of a 3D Dial Bracket..............................................................................115 6.2 Validation of Solution Structure for Micromechanics....................................................116 6.2.1 Square Plate with Circular Inclusion....................................................................116 6.2.2 Determination of Effective Properties of a Fiber Reinforced Composite............117 6.3 Validation of Solution Structure for Grid Refinement...................................................118 6.3.1 Application of Grid Refinement to Constant Stress Problem..............................118 6.3.2 Stress Concentration Problem..............................................................................119 7 CONCLUSIONS AND FUTURE WORK...........................................................................165 7.1 Conclusions.....................................................................................................................165 7.2 Future Work....................................................................................................................167 6 PAGE 7 LIST OF REFERENCES.............................................................................................................169 BIOGRAPHICAL SKETCH.......................................................................................................175 7 PAGE 8 LIST OF TABLES Table page 5-1 Set of periodic boundary conditions for 6 different strains.............................................102 6-1 Average stresses for various loading cases......................................................................162 6-2 Comparison of effective properties..................................................................................163 6-3 Average stresses for various loading cases......................................................................164 8 PAGE 9 LIST OF FIGURES Figure page 3-1 Solution structure for imposing essential boundary conditions.........................................49 3-2 Elastic solid under equilibrium..........................................................................................49 3-3 Steady state heat transfer problem.....................................................................................50 3-4 Boolean operations on Dirichlet functions........................................................................50 3-5 A CSG tree representing Dirichlet function for multiple boundaries................................51 3-6 Boundary value function using cubic B-spline basis.........................................................52 3-7 Representation of grids for material discontinuity.............................................................53 3-8 One dimensional inclusion function for treatment of material discontinuity....................53 3-9 Two dimensional inclusion function for an elliptic inclusion in a square matrix..............54 3-10 Representation of solution structure for local grid refinement..........................................54 3-11 Approximate step function for refinement.........................................................................55 3-12 Coarse and refined grid interface with multiple implicit equations...................................55 3-13 Step function for refinement for a 2D model.....................................................................56 4-1 Approximation using 1D cubic B-splines..........................................................................64 4-2 Shape functions for 1D quadratic B-spline element..........................................................64 4-3 Shape functions for 1D cubic B-spline element................................................................64 4-4 Typical grid used for analysis and the mapping used for geometry and field variable.....65 4-5 Two dimensional quadratic B-spline elements in parametric space..................................66 4-6 Two dimensional cubic B-spline elements in parametric space........................................66 5-1 Typical non-conforming grid used for analysis.................................................................95 5-2 Triangulation of boundary element for volume integration...............................................96 5-3 Representation of band on essential boundaries................................................................96 5-4 Approximation of natural boundary for computation of load vector.................................97 9 PAGE 10 5-5 Division of computational domain into three regions........................................................98 5-6 Undeformed and deformed configuration of a set of 4 RVEs...........................................98 5-7 Representative volume element used for computation of effective properties..................99 5-8 Division of computation domain into various regions.......................................................99 5-9 Typical analysis grid with elements having low volume fraction...................................100 5-10 Extrapolation of external node for linear Lagrange element...........................................100 5-11 Extrapolation of external nodes for cubic B-spline element............................................100 5-12 Extrapolation of external nodes for quadratic B-spline element.....................................101 5-13 Extapolation for boundary element..................................................................................101 6-1 Plate subjected to uniaxial loading..................................................................................121 6-2 Non-conforming grid using one cubic B-spline element.................................................121 6-3 Plot of Dirichlet functions................................................................................................122 6-4 Effect of range parameter on the quality of solution.......................................................123 6-5 Contour plot of x for one element model using cubic B-splines..................................123 6-6 Thin cantilever beam subjected to end loading................................................................124 6-7 Plot of error in the tip deflection with varying aspect ratio.............................................124 6-8 Curved beam subjected to uniform loading.....................................................................125 6-9 Plot of error in strain energy for curved beam (without extrapolation)...........................125 6-10 Plot of error in strain energy for curved beam using extrapolation technique.................126 6-11 Contour plot of y for Q4 elements (without extrapolation)...........................................126 6-12 Contour plot of y for Q4 elements using extrapolation technique................................127 6-13 Contour plot of y using cubic B-spline elements..........................................................128 6-14 Thick cantilever beam used for convergence analysis.....................................................128 6-15 Typical non-conforming grid using cubic B-spline elements..........................................129 10 PAGE 11 6-16 Convergence plot of 2 L error for cantilever beam..........................................................129 6-17 Convergence plot for error in strain energy for cantilever beam.....................................130 6-18 Comparison of deflection curves for 10 element model..................................................131 6-19 Comparison of normal stress along the line AB for 10 element model...........................132 6-20 Comparison of shear stress along the line AB for 10 element model..............................133 6-21 Convergence plot for tip deflection in near incompressible limit....................................133 6-22 Analysis model for a square plate with a circular hole....................................................134 6-23 Contour plot of x using FEM Q4 elements....................................................................134 6-24 Contour plot of x using cubic B-spline elements..........................................................135 6-25 Convergence plot for 2 L error norm for stress concentration problem...........................135 6-26 Convergence plot for error in strain energy for stress concentration problem................136 6-27 Plot of x along the line 2 .......................................................................................136 6-28 Hollow thick cylinder subjected to uniform pressure......................................................137 6-29 Contour plot of radial displacements using cubic B-spline element model.....................137 6-30 Contour plot of radial displacements using FEM Q4 elements.......................................138 6-31 Comparison of radial displacements along the thickness of the cylinder........................138 6-32 Comparison of radial stresses along the thickness of the cylinder..................................139 6-33 Comparison of tangential stresses along the thickness of the cylinder............................139 6-34 Heat conduction in a block with series of pipes..............................................................140 6-35 Plot of Dirichlet function.................................................................................................140 6-36 Plot of boundary value function.......................................................................................141 6-37 Contour plot of temperature distribution using FE Q4 elements.....................................141 6-38 Contour plot of temperature distribution using cubic B-spline elements........................142 6-39 Contour plot of temperature distribution using quadratic B-spline elements..................142 11 PAGE 12 6-40 Comparison of temperature along line AB......................................................................143 6-41 Three dimensional clamp used for analysis.....................................................................143 6-42 Convergence plot for strain energy for clamp model......................................................144 6-43 Contour plot of bending stresses using H8 elements in FEM for 3D clamp...................144 6-44 Contour plot of bending stresses using H8 elements for 3D clamp.................................145 6-45 Contour plot of bending stresses using cubic B-spline elements for 3D clamp..............145 6-46 Three dimensional rib model used for analysis...............................................................146 6-47 Convergence plot for strain energy for rib model............................................................146 6-48 Contour plot of bending stresses using tetrahedral elements for 3D rib..........................147 6-49 Contour plot of bending stresses using B-spline elements for 3D rib.............................147 6-50 Contour plot of bending stresses using H8 elements for 3D rib......................................148 6-51 Three dimensional dial bracket used for analysis............................................................148 6-52 Convergence plot for strain energy for dial bracket model.............................................149 6-53 Contour plot of bending stresses using a cubic B-spline element model for dial bracket..............................................................................................................................149 6-54 Contour plot of bending stresses using a tetrahedral element model for dial bracket.....150 6-55 Contour plot of bending stresses using a H8 element model for dial bracket.................150 6-56 Square plate with a circular inclusion..............................................................................151 6-57 Contour plot of x using Q4 elements in IBFEM............................................................151 6-58 Contour plot of x using Q4 elements in FEM...............................................................152 6-59 Convergence plot for error in strain energy for square matrix with inclusion model......152 6-60 Comparison of along the line AB................................................................................153 xu 6-61 Comparison of stresses along the line AB.......................................................................153 6-62 Comparison of strains along the line AB.........................................................................154 6-63 Plot of stresses for cases 1 and 2......................................................................................155 12 PAGE 13 6-64 Plot of stresses for cases 3 and 4......................................................................................155 6-65 Plot of stresses for cases 5 and 6......................................................................................155 6-66 Plot of stresses for cases 1 and 2 for 3D model...............................................................156 6-67 Plot of stresses for cases 3 and 4 for 3D model...............................................................156 6-68 Plot of stresses for cases 5 and 6 for 3D model...............................................................156 6-69 Contour plot of using Q4 elements and grid refinement.............................................157 u 6-70 Contour plot of using Q4 elements and grid refinement.............................................157 v 6-71 Contour plot of using cubic B-spline elements and grid refinement...........................158 u 6-72 Contour plot of using cubic B-spline elements and grid refinement...........................158 v 6-73 Square plate with a circular hole......................................................................................159 6-74 Contour plot of x for local refinement using Q4 elements...........................................159 6-75 Contour plot of x for local refinement using cubic B-spline elements.........................160 6-76 Convergence plot for 2 L error in local grid refinement..................................................160 6-77 Convergence plot for energy error in local grid refinement............................................161 6-78 Contour plot of x for local refinement using Q4 and cubic B-spline elements............161 13 PAGE 14 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy FINITE ELEMENT ANALYSIS USING UNIFORM B-SPLINE APPROXIMATION AND IMPLICIT BOUNDARY METHOD By Ravi Kumar Burla May 2008 Chair: Ashok V. Kumar Major: Mechanical Engineering The finite element method (FEM) is a widely used numerical technique for solution of boundary value problems arising in engineering applications. FEM is well studied and has several advantages; however the stress and strain solutions obtained by FEM are discontinuous across element boundaries and require smoothing algorithms to make the stress/strain appear smooth in the analysis domain. Automatic mesh generation is still a difficult task in FEM for complicated geometries and often requires user intervention for generating quality meshes. Finite elements based on uniform B-spline basis functions defined on a structured grid are developed which provide at least continuous solution through out the analysis domain. Numerical techniques are presented in this work based on structured grids and implicit equations of the boundaries of analysis domain, material interface boundaries or coarse/fine grid interfaces for solution of various engineering boundary value problems. Solution structures are constructed using approximate Heaviside step functions for imposition of Dirichlet boundary conditions, treatment of material discontinuity to perform micromechanical analysis and for local refinement of grid. The use of structured grid eliminates the need for constructing a conforming mesh and results in significant savings of time in pre-processing stage of the design cycle. 1C 14 PAGE 15 Numerical examples are presented to demonstrate the performance of B-spline elements. The results are compared with analytical solutions as well as traditional finite element solutions to demonstrate the ability of B-spline elements to represent continuous stress and strain through out analysis domain. Convergence studies show that B-spline elements can provide accurate solutions for many engineering problems with fewer numbers of elements and nodes as compared to traditional FEM. Solution structure for treatment of material boundary is validated by performing a convergence analysis on a problem involving circular inclusion in a square matrix and by determining effective properties of fiber reinforced composite. Solution structure for local grid refinement is validated by analyzing classical stress concentration problems. 15 PAGE 16 CHAPTER 1 INTRODUCTION 1.1 Overview Numerical methods are necessary for solving many problems arising in engineering analysis. Even though analytical solutions are available for several engineering problems with simple geometries, for complicated domains it is difficult to obtain analytical solutions and hence numerical methods are needed to obtain accurate solutions. The Finite Element Method (FEM) is widely used in engineering analysis not only due to its applicability to solve a variety of physical problems but also due to its ability to handle complex geometry and boundary conditions. In FEM, the analysis domain is subdivided into finite elements which approximate the geometry of the analysis domain. The field variable is approximated on these finite elements using piecewise continuous interpolations. Even though FEM is a well established technique, it has certain limitations. Approximation schemes used in FEM are typicallycontinuous and hence there is a discontinuity of gradients across the elements. The lack of gradient continuity in FEM necessitates the smoothing of results. There are a number of approximation schemes like B-splines, moving least squares and Hermite interpolations which provide at least continuity in the domain of interest. B-spline basis functions have compact support like traditional finite element shape functions and lead to banded stiffness matrices. They form partition of unity and hence can exactly represent the state of constant field in the domain of interest. The B-spline approximations have good reproducing properties and when a cubic B-spline is used, it can represent all polynomials up to third order exactly. B-spline approximation provides high order of continuity and is capable of providing accurate solutions with continuous gradients through out the domain. The uniform B-spline basis functions are polynomials and accurate integration for stiffness matrices and load vectors can be 0C 1C 16 PAGE 17 performed by using well known techniques like Gauss Quadrature. These properties of B-splines make them attractive for use in analysis. B-splines have not been used extensively for analysis because they do not satisfy Kronecker-delta property and hence essential boundary conditions cannot be imposed directly by specifying nodal values. B-splines are defined on a grid which limits the geometry of the analysis domain to a parallelepiped. In order to be able to use uniform B-spline basis for approximation of finite element solutions, structured grids need to be used since the uniform B-spline basis are defined only on a uniform structured grid. A structured grid is much easier to generate than a finite element mesh and all the elements in the grid have regular geometry (rectangles/cuboids). The use of structured grids for analysis eliminates the need to generate a conforming mesh which is a significant advantage over traditional FEM. Automatic mesh generation algorithms for FEM are well developed 2D problems but can be unreliable for 3D geometries often resulting in poor or distorted elements in some regions that can lead to large errors in solution. Significant amount of user intervention is needed to correct such problems making mesh generation the most time consuming step in the design process. These problems associated with mesh generation can be avoided with the use of a structured grid. When structured grids are used to approximate solution in the analysis domain, the boundaries of the domain need to be defined independent of the structured grid and this can lead to several elements on the boundary which have partial volume inside the grid elements. Special techniques are therefore needed to compute the stiffness matrices and load vectors in these partial elements. B-spline approximations do not satisfy Kronecker-delta property and hence boundary conditions cannot be applied by specifying nodal values. Furthermore, the nodes in a structured grid are not guaranteed to lie on the boundary, therefore the traditional methods used 17 PAGE 18 in FEM for applying essential boundary conditions cannot be used and new techniques are needed for applying essential boundary conditions. In order to enforce essential boundary conditions when approximations are used on structured grids, Implicit Boundary Method is used in this work. In Implicit Boundary Method, solution structures are constructed by using the approximate Heaviside step functions based on the implicit equations of the boundaries for exact imposition of essential boundary conditions. 1C Implicit Boundary Method is extended in this work by developing solution structure for treatment of material discontinuity using approximate Heaviside step functions based on the implicit equations of matrix/inclusion boundaries. This technique is further used for determination of effective properties of fiber reinforced composite. In many practical problems local refinement of mesh/grid is necessary to accurately compute the stresses in the regions of stress concentrations; this problem is addressed in this work by developing solution structure for local grid refinement using approximate Heaviside step functions based on the implicit equations of the coarse/refined grid boundaries. The performance and capabilities of B-spline based finite elements are studied by solving several problems in structural and thermal analysis in both 2D and 3D. Convergence analysis is also performed using several examples to study the behavior of B-spline elements and it has been demonstrated that B-spline element provide accurate solutions with fewer number of elements when compared to Lagrange elements used in traditional FEM. Micromechanical analysis for determining effective properties of fiber reinforced composites is performed by employing the solution structure for treatment of material discontinuity using Lagrange finite elements. Solution structure for local refinement is used to analyze stress concentration problems where local 18 PAGE 19 resolution of solution near regions of stress concentrations is necessary for performing design analyses. 1.2 Goals and Objectives The goal of this research is to develop finite elements based on B-spline approximations and to extend Implicit Boundary Method by developing solution structures for analysis of engineering boundary value problems involving material discontinuity and local grid refinement. The main objectives of this dissertation are To extend the Implicit Boundary Method for use of B-spline basis functions. To develop finite elements based on B-spline approximation which provides at least continuity in the analysis domain. 1C To develop solution structures for modeling material discontinuity and imposition of periodic boundary conditions involved in the micromechanical analysis of composite microstructures and to determine effective elastic properties of fiber reinforced composites. To develop solution structures for local refinement of grid to accurately capture stress concentrations and study the use of the B-spline approximations and Lagrange interpolations in local refinement. To develop efficient numerical techniques for computation of various terms involved in the Galerkin weak form. To study the convergence behavior of the solution structures and B-spline finite elements using several numerical examples and to study the advantages and disadvantages of using the B-spline based finite element. 1.3 Outline The rest of the dissertation is organized as follows: Chapter 2 discusses several analysis techniques that are reported in literature which use B-splines to approximate solution. A survey on structured grid methods, meshless methods, treatment of material discontinuity and adaptive grid/mesh refinement is also presented. Chapter 3 discusses the construction of solution structure for exact imposition of essential boundary conditions for structural and thermal analyses. Solution structures are developed for 19 PAGE 20 modeling material discontinuity and for imposition of periodic boundary conditions which are used for micromechanical analysis of composites. Finally, solution structure is developed for local refinement of grid to improve the local resolution of the solution. In chapter 4, finite elements based on B-spline approximation are developed which provide and continuity in one, two and three dimensions. 1C 2C Chapter 5 discusses efficient numerical implementations for computation of stiffness matrices and load vectors involved in the modified weak form corresponding to the solution structures presented in chapter 3. These techniques are implemented in JAVA and several problems in structural mechanics, steady state heat transfer, composite micromechanics and localized stress concentrations are studied in detail in chapter 6. Convergence analyses are also performed for these problems to study the behavior B-spline elements and solution structures in terms of accuracy and efficiency. Finally, in chapter 7, several conclusions are drawn from the numerical experiments performed in chapter 6. The advantages and shortcomings of proposed solution structures and B-spline finite elements are discussed in detail. The scope for future work is presented towards the end of the dissertation. 20 PAGE 21 CHAPTER 2 ANALYSIS TECHNIQUES FOR BOUNDARY VALUE PROBLMS In this chapter, a detailed discussion on several analysis techniques that are used for solution of engineering boundary value problems is presented. In the first part of this chapter, a brief overview of the finite element method (FEM) is presented along with its advantages and shortcomings. Techniques for employing B-spline approximations for analysis are discussed in the next part of this chapter. A detailed discussion on structured grid methods and meshless methods is presented in the later part of this chapter Micromechanical analysis is an important engineering application where FEM is used extensively. Microstructures are modeled using a representative volume element and periodic boundary conditions are applied using multipoint constraint techniques. Meshless methods have been extended for treatment of material discontinuity and for performing micromechanical analysis. A survey on these techniques is presented in this chapter along with a survey on employing FEM for determination of effective properties. Another important feature of FEM is its capability for local refinement. A survey on the refinement techniques available in structured grid and meshless framework is presented in the last part of this chapter. 2.1 Finite Element Method Finite element method (Hughes, 2000 ; Cook et. al. 2003) is a well established numerical technique and is widely used in solving engineering problems in academia as well as in industry. For finite element analysis, the domain of analysis is subdivided into elements and the resulting finite element mesh approximates the geometry and is also used to approximate the solution by piece-wise interpolation within each element. Using Galerkins approach, the governing differential equation in its strong form is converted into a weak form which allows the use of 21 PAGE 22 interpolation/approximation schemes with lower continuity requirements than the corresponding strong form. This leads to a set of linear equations which are then solved using Gauss elimination, Cholesky decomposition or iterative schemes like Gauss-Siedel. FEM has been applied to a wide variety of physical problems in linear/non-linear static and dynamic elasticity, steady state heat transfer, electrostatics and magnetostatics. The approximation space used in FEM is typically Lagrange interpolation that yields continuous solutions. However, continuous solutions have discontinuous stresses and strains (gradients) across the element boundaries. The lack of stress and strain continuity necessitates the use of smoothing schemes to plot these quantities such that they appear to be continuous. Solutions can be represented with higher degree of continuity using several other basis functions like B-splines, moving least squares (MLS) and Hermite functions which provide at least continuity in the domain of interest. In this dissertation, B-spline approximations are studied in detail which provide and continuous solutions in the analysis domain. 0C 0C 1C 1C 2C Another limitation of the FEM is the need to generate a conforming mesh, which may need user intervention for correcting the problems associated with distorted meshes. Automatic mesh generation algorithms work satisfactorily for 2D but are still unreliable for 3D. In order to avoid the problems associated with mesh generation, several methods have been proposed which fall into two categories: meshless methods and structured grid methods. In this chapter, a detailed survey on various structured grid and meshless methods is presented. 2.2 B-Spline Approximation Based Methods A number of researchers have studied the use of B-spline approximations to solve specific problems where the geometry of the analysis domain is simple (e.g. beams or rectangular plates). The boundary conditions were applied by using multiple knots at the boundary of the analysis 22 PAGE 23 domain which emulates Kronecker-delta behavior on the boundaries. Leung and Au (Leung and Au, 1990) introduced the boundary transformation method by constructing an appropriate transformation matrix based on the boundary configuration to impose Dirichlet boundary conditions. Kagan et. al.(1998, 2003) have provided a scheme to integrate B-spline finite element method with a mechanically based computer aided design and the Dirichlet boundary conditions were enforced by imposing multiple knots on the essential boundaries. They also developed methods for adaptive refinement of B-spline elements. In summary, these methods have limitations on the geometry of analysis domain. Hughes and his coworkers (Hughes et. al. 2005 ; Cottrell et. al. 2006 ; Cottrell et. al. 2007) have developed a more general approach for using B-splines called the isogeometric analysis method, where the geometry is represented by Non-Uniform Rational B-Spline (NURBS) basis functions to get an exact geometric representation. The same basis functions can be used for the analysis or their order can be raised without changing the geometry or parameterization. This enables the traditional h-, and padaptive refinement as well as k-refinement where the degree of continuity of the solution can be raised. Several alternatives for imposing Dirichlet boundary conditions were proposed including direct application by specifying control variables and use of constraint equations, but this still remains a topic of research. 2.3 Structured Grid Methods A number of techniques have been developed based on structured grids to approximate the solution. Extended Finite Element Method or X-FEM was developed by Belytschko and co-workers (Belytschko et. al. 2003). This method uses implicit equations for definition of the geometry of the analysis domain and the Dirichlet boundary conditions are applied using Lagrange multipliers. 23 PAGE 24 Penalty boundary method developed by Clark and Anderson (Clark and Anderson, 2003a, 2003b) uses a regular grid and a constrained variational formulation to include both essential and natural boundary conditions which is then solved using penalty method. Shapiro (Shapiro, 1998 ; Shapiro and Tsukanov, 1999) and Hollig (Hollig et. al. 2001; Hollig and Reif 2003 ; Hollig 2003) used R-functions to define the boundary of the problem domain and the solution structure proposed by Kantorovich and Krylov (Kantorovich and Krylov 1958) was constructed to satisfy Dirichlet boundary conditions. Shapiro used transfinite Lagrange interpolation to impose Dirichlet boundary conditions while Hollig used weight functions based on distance functions and R-functions. Hollig developed extended B-spline basis to represent the solution in the analysis domain. The method presented in this dissertation is very similar to the approach used by Hollig, however approximate Heaviside step function of essential boundaries are used in this work to construct weight functions and the behavior of solution structure for very small values of range parameter used to construct step functions is also studied. All the internal elements in our method have identical stiffness matrices unlike Holligs method. 2.4 Meshless Methods A number of methods have been developed in the past two decades which avoid mesh generation. These methods are collectively termed as meshless or meshfree methods. Some of these methods provide continuous gradients in the analysis domain. Smoothed particle hydrodynamics (SPH) is one of the first meshless methods in literature (Lucy, 1977 ; Gingold, 1982); it was first used to model astrophysical phenomenon and was modified to include continuum mechanics problems (Monahan, 1982). This approach uses Shepard shape functions as interpolants and uses point collocation method for discretization and nodal integration for stiffness matrix computation. The drawbacks of this approach are poor accuracy and spatial 24 PAGE 25 instability (tensile instability). The behavior of solution depends on the proper distribution of nodes. So there is always an effort needed to place the nodes properly in order to avoid instability. Element Free Galerkin Method was developed by Belytschko and his coworkers (Belytschko et. al., 1994 ; Belytschko et. al., 1996 ; Belytschko et. al. ; 1996, Dolbow and Belytschko, 1999 ; Krongauz and Belytschko, 1996 ; Krongauz and Belytschko ; 1997 ; Krysl and Belytschko 1997). This method is based on moving least square approximation and uses global Galerkin method for discretization. Quadrature was performed on the background cells. Boundary conditions were applied using Lagrange multipliers or by modification of variational principle, where Lagrange multipliers were replaced by traction forces. The former method leads to an increased number of algebraic equations and the stiffness matrix is not positive definite. The latter method reduces the accuracy of the results and is not as convenient as imposing essential boundary conditions directly. A new method for imposing essential boundary conditions by using a string of Finite Elements around the boundary of the domain was later incorporated into this method. Meshless Local Petrov-Galerkin (MLPG) was introduced as a truly meshless method by Atluri and his coworkers (Atluri et al. 1999 ; Atluri and Zhu 2000a, 2000b). This method uses moving least square approximation as the interpolation functions. It uses the Local Petrov-Galerkin weak form which is developed in the local sub-domain. The trial and test spaces are different; the test functions vanish on the boundary of the sub-domain, as long as these boundaries are internal to the problem domain. If the sub-domain intersects the boundary of global domain, then the boundary condition is applied on this local boundary by penalty method. The stiffness matrix obtained is banded but not symmetric. The boundary conditions are applied 25 PAGE 26 using penalty method. Another advantage of this procedure over FEM is that smoothing techniques are not required for computation of derivative. There is no need of special integration schemes, and there is flexibility in choosing the size and shape of the local domain, which can lead to convenient methods for integration in non-linear problems. The drawbacks of this method are that the stiffness matrices generated by this approach are not symmetric positive definite and penalty method imposes boundary conditions approximately. Another drawback is that we have to make sure that the union of all the local domains should cover the global boundary. Partition of unity based methods (PUFEM) form a class of meshless methods which use partition of unity functions to blend the local approximation functions which represent the local behavior (Melenk and Babuska, 1996). The main advantage of PUFEM is that we can use the knowledge of the local behavior of solutions efficiently to get good numerical approximation. This method also yields good results in the cases where FEM gives poor results or when FEM is computationally very expensive. One drawback of this approach is that the stiffness matrix obtained from this procedure can be singular and the basis function needs to be modified appropriately or special techniques are needed to solve the singular system. Generalized Finite Element Method was developed as a hybrid of Finite Element Method and Partition of Unity (Babuska et. al. 2002 ; Duarte et. al. 2000). The desirable features like inter-element continuity, local approximability and Kronecker-delta property from FEM are retained in GFEM. Method of finite spheres (De and Bathe 2000 ; De and Bathe 2001) was developed as meshless method which is based on partition of unity; this approach does not use background cells for numerical quadrature. It uses compactly supported partition of unity functions and Bubnov-Galerkin method as weighted residual scheme. Integration of weak form is performed on spheres (circles in 2D), using specialized cubature rules. Essential boundary conditions are 26 PAGE 27 imposed by placing nodes on the boundary so that the shape functions corresponding to these nodes emulate Kronecker-delta behavior. Shepard functions are used as partition of unity functions and quartic splines are used for weight functions. Method of HP Clouds is a partition of unity based meshless method developed by Duarte and Oden (Duarte and Oden 1996 ; Duarte et. al. 1996 ; Liszka et. al. 1996). HP-Cloud functions are created by multiplying the partition of unity functions with local basis functions. Shepard functions or MLS functions are used as partition of unity. Anisotropic approximation can be incorporated by using the basis function which has higher order in one direction and lower order in other direction. P-refinement is done by increasing the size of the local domain; hence more number of points is included in the sub-domain and also by increasing the polynomial order in the basis function. Boundary conditions are applied by modified variational principle. Natural element method developed by Sukumar et. al.(1998) uses the natural neighbor interpolants which are based on Voronoi tessellation of the nodes. These interpolants are smooth over the whole domain excepting at the nodes where they are continuous. The essential boundary conditions are applied on nodes for convex domains. 2.5 Techniques for Determination of Effective Properties of Composite Microstructures The finite element method (FEM) has been used extensively in micromechanical analysis of composite microstructure to determine effective properties (Sun and Vaidya ; 1996 ; Taliercio, 2005 ; Xia et. al. ; 2003 ; Marrey and Sankar, 1997 ; Zhu et. al. 1998). The effective elastic properties are determined by performing six different analyses on the representative volume element (RVE) with unit non-zero strain component in each analysis and by imposing appropriate periodic boundary conditions on the boundaries of the RVE. The mesh is generated in such a way that the element boundaries conform to the matrix-inclusion interfaces. However, for complicated microstructures, automatic mesh generation becomes difficult and special 27 PAGE 28 techniques and algorithms need to be developed (Kim and Swan, 2003) to generate a quality mesh with quadrilateral or hexahedral elements. For example, in the analysis of composite microstructures involved in textile composites, the automatic mesh generation often becomes difficult and one of the main problems is generating a quality mesh in the interfacial region (Kim and Swan, 2003). This region is multiply connected and complications may arise for generation of mesh with identical nodes on the opposite faces of the RVE in order to impose periodic boundary conditions using multipoint constraints. Meshless methods have been extended for treatment of material discontinuity. Element Free Galerkin Method has been extended to incorporate discontinuous derivatives (Cai and Zhu, 2004 ; Cordes and Moran, 1996 ; Krongauz and Belytschko, 1998) for application to material discontinuity problems and Meshless Local Petrov-Galerkin Method has been extended to incorporate material discontinuity and for performing micromechanical analysis (Dang and Sankar, 2007). 2.6 Techniques for Local Refinement In many practical problems, there may be regions with stress concentrations and hence a uniform mesh/grid for the entire model can become very expensive computationally. In order to reduce the computational cost local mesh/grid refinement is necessary. In this section, a brief discussion on the local refinement techniques is presented. In the finite element analysis, local refinement is typically obtained by h-refinement or p-refinement. In h-refinement, nodes are inserted in the regions of stress concentration which increases the number of elements and reduces the element size in the region of interest while p-refinement is obtained by increasing the order of polynomial basis used for interpolation. Another type of refinement is termed as r-refinement, where the total number of nodes and elements in the refined model is same as the coarse model but the nodes are rearranged so that the number of nodes and elements in the regions of stress concentration is increased and k28 PAGE 29 refinement is obtained when the order of continuity of the basis functions is increased. In meshless methods, refinement is achieved by inserting nodes in the regions of stress concentrations. Hierarchical mesh refinement is studied by Krysl et. al. (2003) in the finite element method framework. In this method, the nodal basis functions on the coarse mesh are split into the basis functions of same order in the finer mesh. They have provided algorithms for refinement as well as un-refinement for finite element meshes. Hollig (2003) also presented a hierarchical refinement of B-spline basis functions, where the B-spline basis are split into the B-spline basis of same order but on a finer grid with half the size of the coarse grid. The hierarchical refinement provides a framework for getting a multi-level resolution, where the basis functions are split in a recursive manner. In this dissertation, a technique for local grid refinement is presented based on construction of solution structure using approximate step functions. 29 PAGE 30 CHAPTER 3 SOLUTION STRUCTURES FOR IMPLICIT BOUNDARY METHOD The essential boundary conditions are imposed in FEM by specifying the values at the boundary nodes on which essential boundary conditions are specified. When structured grids are used for discretizing the analysis domain, nodes are not guaranteed to be present on essential boundaries; hence special techniques need to be devised to be able to impose essential boundary conditions. In this chapter, a detail discussion on construction of solution structure for exact imposition of Dirichlet boundary conditions using the implicit equations of boundaries is presented. Approximate Heaviside step functions are used to construct a solution structure based on the technique proposed by Kantorovich and Krylov (Kantorovich and Krylov, 1958 and is termed as Implicit Boundary Method (Padmanabhan, 2006 ; Kumar et. al. 2007). Non homogenous essential boundary conditions are imposed by constructing piecewise continuous functions which satisfy the specified values at the essential boundaries and are termed as boundary value functions. This solution structure is used for modifying the weak form corresponding to linear elasticity problems and steady state heat transfer problems. Implicit Boundary Method has been extended in this dissertation to be able to use B-spline basis. A solution structure for treatment of material discontinuity based on the implicit equation of the material boundaries and approximate Heaviside step function is presented later in this chapter. Techniques for imposing periodic boundary conditions using multipoint constraints are also discussed in detail. In the final part of this chapter, a solution structure for adaptive refinement of grid based on the implicit equations of the grid element boundaries and approximate Heaviside step functions is discussed in detail. The weak form corresponding to linear elasticity is modified using this solution structure. 30 PAGE 31 3.1 Solution Structure for Enforcing Essential Boundary Conditions The essential boundary conditions are specified as prescribed displacements or prescribed temperatures on the essential boundaries. Let u be a trial function defined over that must satisfy the boundary condition 23or RR u a along u which is the essential boundary of Here, ucorresponds to a displacement component in linear elasticity and temperature in steady state heat transfer. If ( and ()0ax x :a R ) is the implicit equation of the boundary then the trial function can be defined as u ()()()()auUxxx ax (3-1) The solution structure for the trial function defined in Eq. (3-1) is then guaranteed to satisfy the condition u along the boundary defined by the implicit equation for any. The variable part of the solution structure is the function which can be replaced by a finite dimensional approximate function defined by piece-wise interpolation/approximation within elements of a structured grid. In the context of such solution structures, the function is often referred to as a weighting function. Since these functions assist in imposing the Dirichlet boundary conditions, we have referred to these functions as Dirichlet functions or just D-functions. D-function is constructed such that it takes a value of zero and has non-zero gradient on the essential boundary and has a value of one far from the essential boundary. Instead of directly using the implicit equation of the boundary as D-functions, approximate step functions constructed using these implicit equations have been used to construct solution structures in this work. The use of approximate step functions for construction of D-functions localizes them to the boundary elements and hence the stiffness matrix for all the internal matrices is identical as the D-function is unity inside them. In the a ()0ax :UR ()Ux ()hUx ()ax 31 PAGE 32 solution structure demonstrated in Figure 3-1, D represents the Dirichlet function which is zero at the first node of element e1 and reaches unity within this element and remains unity through out the domain. The boundary value function is defined by which takes on the value of the prescribed boundary value of at the boundary and au 0u g u represents the grid variable which is interpolated on the grid consisting of elements e1-e5. The approximate solution is given by g auDuu and it can be observed that the approximated solution satisfies the boundary condition exactly. The solution structure for linear elastostatic and steady state heat transfer problems is first discussed in this section for a general D-function and a detail discussion on construction D-functions is provided later in this chapter. 3.1.1 Modified Weak Form for Linear Elasticity The governing differential equation for an elastic body under equilibrium as shown in Figure 3-2 is given as: Solve: in .+b=0 Subject to: on 0u=u u (3-2) on 0.n=t t In the above equation, is the analysis domain bounded by essential boundary and natural boundary { is Cauchy stress tensor and { is the body force and the Dirichlet boundary conditions are displacements specified as on essential boundary and the natural boundary conditions are specified traction The principle of virtual work is the weak form for the boundary value problem and is given as: u t } }b 0u 0t tddTTTutub d (3-3) 32 PAGE 33 Here {is the virtual strain and { is virtual displacement. A trial solution structure for the displacement field can be constructed such that the Dirichlet boundary conditions are satisfied } }u g asuDu+uu+u a (3-4) In this solution structure, { } g u is a grid variable represented by piece-wise approximation over the grid, and { is the boundary value function which is a vector field constructed such that it takes on the prescribed values of displacements at boundaries. The variable part of the solution structure is {} }au []{} s guDu which satisfies the homogeneous boundary conditions. is a diagonal matrix whose components are Dirichlet functions that vanish on boundaries on which the component of displacement is specified, is the dimension of the problem. The same D-functions used to construct the trial solution are also used to construct the test functions. [](,..,)dindiagDDD iD thi dn g uDu (3-5) Using these solution structures, the strains and stresses can be expressed as the sum of a homogeneous part ( s ij s ij ) and a boundary value part ( aij aij ). ,,,,1122 s saasijijjiijjiijijuuuu a (3-6) () s asijijklklijklklklijijCC a (3-7) Substituting s=+ a in the weak form, we can rewrite it as tdddTTTTsautub d (3-8) 33 PAGE 34 The grid variable is approximated within each element as g gijuNu ij where are the shape functions and jN g iju are the nodal values of the grid variable for component of displacement. The virtual strain and stress within each element can then be expressed as and {}, where { and { represents the nodal values of the grid variable for the element e and [ is the strain-displacement matrix. Using these notations, the weak form can be expressed in a discrete form very similar to the finite element equations. thi {}[]{}eBX [][]{}se=CBX }T eX }eX ]B 111etNENENBEeeeeeeeeddTTTTXBCBXXFXNt (3-9a) iiedTTaFNbB d (3-9b) The summation is over the number of elements (NE) and the number of boundary elements (NBE) in the grid used for analysis. As in the finite element method these element matrices can be assembled to form global equations, where the left hand side terms form a stiffness matrix and the right hand side a force vector. Note that unlike in the finite element equations, there is a forcing term that arises from the boundary value functions. Furthermore, the nodal variables to be solved for are grid variables and not displacements. Since no assumptions have been made restricting the choice of shape functions used to represent the grid variables, this approach can be used with a variety of basis functions including Lagrange interpolation, B-spline basis functions or meshless approximations. 3.1.2 Modified Weak Form for Steady State Heat Transfer The governing differential equation for general steady state heat transfer problem shown in Figure 3-3 is given as 34 PAGE 35 Solve: T0inkf Subject to: TT on g (3-10) Tonhhkqn The corresponding weak form is given as: TTdTdTdTTTdhchkfqh (3-11) In the above equations, g is the portion of the boundary with prescribed temperature while is the portion with heat flux conditions and h c represents the portion with convection conditions and T is the ambient temperature. is the convection coefficient and is the thermal conductivity of the material which is assumed to be isotropic for simplicity and is the heat source. A trial solution structure for the temperature field is constructed such that the essential boundary conditions are satisfied h k f s agTTTDTT+ a (3-12) Proceeding in the similar manner as for the linear elastic problems discussed in the previous section and choosing the test function gTDT such that it vanishes on g the weak form modifies to the following: 0chcsssgsassasTkTdThTdTkTdfTdqTdhTTTd (3-13) The nodal degrees of freedom for temperature are approximated as gT g iiNT and the homogenous part of solution becomes sTD g iiNT where is the essential boundary function. The gradients of temperature which are used in the weak form are constructed as D s giijTBT j 35 PAGE 36 The weak form takes the following form when the above approximations are incorporated into Eq. (3-13) a0ddddTeceehcTTTggeeTTTeBkBTNhNTBkBTNfNqNh dTda (3-14) In the above expression, all the known quantities are moved to the right hand side and the unknown quantities are placed in the left hand side. The first term on the right hand side is the contribution to the load due to applied essential boundary conditions, the second term is the load due to the heat generation, the third term is load due to applied heat flux and the last term is load due to convection boundary conditions. The first term on the left hand side represents the stiffness matrix while the second term represents the stiffness due to convection boundary condition. 3.1.3 Dirichlet Function The Dirichlet function is a function which vanishes on all essential boundaries where component of displacement (or temperature) is prescribed. In addition, should have a non-vanishing gradient on the essential boundaries to ensure that the gradients of displacements are not constrained. Furthermore, it is necessary to ensure that is non-zero inside the analysis domain so that the solution is not constrained anywhere within the domain of analysis. Implicit equation of the curve or surface representing the boundary can be used for this purpose. However if multiple curves are involved then it is necessary to combine them into a single equation. Some authors have proposed the use of R-functions to construct the Boolean combination of implicit functions (Shapiro, 1998 ; Rvachev et. al., 2001). Hollig (2003) recommends constructing a weight function using distance functions as: iD thi iD iD ()1max(0,1(,)/)wxdistxD This function 36 PAGE 37 varies between 0 and 1 over a narrow strip of width near the boundary. Evaluating this function requires numerically computing the distance function. In this work, approximate step functions similar to the weighting function used by Hollig (Hollig, 2003) are proposed that can be constructed using any type of implicit function of the boundary. If the boundary condition is specified on a boundary with implicit equations then a step function at any given point x is defined as shown below 00()1101kD (3-15) This function is a shifted approximate step function that tends to the Heaviside step function in the limit as 0 but has a value of zero on the boundary instead of 0.5 as in traditional approximations of step functions (Osher, 2002). Moreover, for use as a Dirichlet function it is constructed to have non-zero gradient at 0 At this function has 1kC continuity. In this paper, the step function approximation of D with 0 is used as the Dirichlet function with values of 510 or smaller. When a single boundary passes through a boundary element the D function can be computed as in Eq. (3-15), but when multiple boundaries exist that have Dirichlet boundary conditions on them then it is necessary to construct D functions that are Boolean combinations of the individual step functions. The intersection, subtraction and union between two step functions 1D and 2D denoted as 12D 12D and 12 D respectively, can be defined as shown in Eq. (3-16). 37 PAGE 38 121DDD 2 (3-16a) 121DDD 2 (3-16b) 12121DDDDD 2 (3-16c) Figure 3-4 shows plots of D-functions created as a Boolean combination of two D-functions that represent circles. Note that the D-functions created by the Booleans inherit important properties including zero value at boundaries, non-zero slopes at boundaries and they plateau to a value of 1 for Since the step function defined in Eq. (3-15) has a value equal to zero at the boundary (and not 0.5 as in traditional approximations of step functions), the Boolean operations defined in Eq. (3-16) do not round off corners defined by the intersection of the implicit equations. If more than two boundaries exist within a boundary element, then the Dirichlet function is constructed using the Boolean tree similar to CSG tree as shown in Figure 3-5. 3.1.4 Boundary Value Function The boundary value function (displacement component for linear elasticity and temperature for steady state heat transfer) must be defined such that they have value equal to the imposed essential boundary conditions on the appropriate boundaries. If the essential boundary conditions are specified at multiple boundaries and the assigned values are not the same for all boundaries then it is necessary to construct a continuous function that takes on appropriate values at respective boundaries while transitioning smoothly in between. While there is no unique method for constructing boundary value functions, it is beneficial to ensure that the trial function is a polynomial of the same order as the shape functions used for interpolating the grid variable. This can be accomplished by constructing the boundary value function by piecewise au 38 PAGE 39 interpolation/approximation using the element basis functions. The nodal values of can be constructed similar to the grid variable au g u using the grid element basis functions. aiiiuN au (3-17) In the above equation, are the shape functions of the grid element and are the nodal values of It is important that the function should belong to the space spanned by the shape functions used to represent the grid variable. Otherwise, the grid variable has difficulty compensating for in the solution structure leading to poor convergence. Indeed, the solution structure will not satisfy completeness requirement, that is the ability to accurately represent constant strains, if basis functions used to construct does not span a sub-space of the basis functions that are used to represent the grid variable iN aiu au au au au g u If the basis functions in Eq. (3-17) are identical to the basis functions used for approximating the grid variable then this requirement is automatically satisfied. If solutions that have or continuity is desired then quadratic or cubic B-spline basis should be used in Eq. (3-17). The nodes of the boundary elements should be assigned nodal values of such that the resulting approximation will have the desired values at the boundary passing through the element. This is very easy to achieve when the assigned value is constant or even linearly varying along the boundary. When the specified boundary condition is a constant value of displacement along boundary, all the support nodes of the element have to be assigned value equal to the specified boundary condition. In the case of the 2D cubic B-spline element, the approximation within the element is controlled by 16 support nodes. Therefore, the value at all these nodes should be assigned to the same value to ensure that the resulting approximation has a constant value within the element. Similarly, if the specified boundary value has to vary linearly (as in hydrostatic pressure which varies linearly in the z-direction), the values iN 1C 2C au 39 PAGE 40 specified at all the support nodes should be on a plane corresponding to the linear distribution. This is illustrated in Figure 3-6, where a simple geometry is used to demonstrate the construction of the boundary value function. The left edge has a constant value of 10au while the right edge has a bilinear variation of The boundary value function is constructed by identifying the elements for the left boundary and specifying the nodal value of for all the support nodes, and then by identifying the boundary elements for right boundary and specifying the nodal values computed using 1 10(1)5(1)aux y 10 0(1)5(1) x y on the corresponding support nodes. The nodes which are not support nodes for the boundary elements are set to zero. The plot for the boundary value function is shown in Figure 3-6 and it can be observed that the boundary value function is continuous over the entire domain and it smoothly transitions from the specified value of 10 from the left boundary to zero in the middle of the domain and then smoothly transitioning from zero to the specified bilinear distribution along the right edge. Even though a simple geometry is used to illustrate the concept of constructing boundary value function, the procedure is quite general. However, care has to be taken to ensure that the grid has sufficient density to be able to construct the appropriate boundary value function. For more general variation of boundary values, the nodal values can be computed by least-square minimization. The rest of the nodes in the grid that are not part of any boundary element can have any arbitrary value and therefore can be set equal to zero. The boundary value function contributes to the load computation on the right hand side of Eq. (3-10) for all those elements in which the gradient of has non-zero magnitude. au au 3.2 Solution Structure for Micromechanics In the finite element method, a conforming mesh is created such that the element edges/faces form the material boundary. When structured grids as shown in Figure 3-7 are used 40 PAGE 41 for analysis, the nodes may not lie on the material interface and hence there is challenge in imposing interface conditions for material discontinuity. In Implicit Boundary Method, this is achieved by construction of a solution structure which ensures that interface conditions are satisfied for displacements. The implicit equations of the inclusion boundaries are used to construct a step function for inclusion such that it is unity inside the inclusion and smoothly blends to zero at the inclusion boundary. Such a function can be used to patch the solution from the inclusion grid and matrix grid to form a solution structure. Figure 3-8 shows the inclusion function in one-dimension and Figure 3-9 shows the inclusion function for an elliptic inclusion in a square matrix. It can be seen that the inclusion function is unity inside the inclusion and transitions to zero at the boundary of the inclusion, the rate with which this function reaches zero is specified by the parameter The expression for inclusion function is similar to the Dirichlet function described earlier in this chapter and shown below 00()1101kincH (3-18) This function inherits all the properties of Dirichlet function and is a shifted approximate step function that tends to the Heaviside step function in the limit as 0 but has a value of zero on the boundary unlike in traditional approximations of step functions where it has value of 0.5 Moreover, for allowing normal strain discontinuity, ()incH is constructed to have non-zero gradient at 0 At this function has 1kC continuity. In this work, the step function approximation of with and incH 2k 0 is used as Inclusion function with values of 210 41 PAGE 42 or smaller. The Inclusion function is used to construct the solution structure for material discontinuity as shown in Eq. (3-19). (1)incincHHg1g2uu u (3-19) This solution structures combines the solution from both grids to represent the solution in the analysis domain. When the Inclusion function is unity, the solution is given by which is the solution from the inclusion part and when the Inclusion function is zero, the solution is given by which is the solution from the matrix. In the region where the inclusion function varies from unity to zero, the solution is a blend of the solutions from matrix grid and inclusion grid. This way of constructing solution structure ensures the displacement continuity of the solution through out the analysis domain. This solution structure allows the normal strain discontinuity; this can be seen from the gradients of the solution structure as shown below g2u=u g1u=u 121(1)ggincincincgincgiiiijjjjuuuHHuHxxxxx 2ijHu (3-20) In this expression, the first term: 1(1) g incijuH x and third term 2 g incijuH x are continuous at the interface boundary while the second term: 1inc g ijHux and fourth term 2inc g ijHux are discontinuous at the interface boundary and vanish in the matrix region and are non-zero in the inclusion region. Therefore these terms provide independent slopes at both sides of the interface resulting in discontinuous normal strains. Modified Weak Form. The weak form of the elastostatic boundary value problem is the principle of virtual work which can be written in the following form as shown in Eq. (3-21). Here {is the virtual strain, { is Cauchy stress tensor, {} is the traction vector, { is } } t }u 42 PAGE 43 virtual displacement and { is the body force. For elastostatic problems, the essential boundary conditions are displacements {} specified on }b {}0u=u u and the natural boundary conditions are traction {} specified on {}0t=t t tddTTTutub d (3-21) When the solution structure shown in Eq. (3-19) is used and the displacements are interpolated as 11[]{ggiiNug1u= } with 1incHg1N N in grid corresponding to matrix and 22[]{ggiiNug2u= } with incHg2N N in inclusion grid, the expression for strains takes the following form 111222121221uxuxuuxx g 112 g 2uBBu (3-22) In this expression, 21BB is the modified strain displacement matrix and is a combination of strain displacement matrices from the two grids. The expression for stress is given as iijC j where is the elasticity stiffness matrix, and in computation of stress, the inclusion (or matrix) elasticity stiffness matrix is used when the strains are evaluated in the inclusion (or matrix) region. The modified weak form takes the following discrete form when the above expressions for stress and strain are incorporated into the weak form (Eq. 3-23). ijC 22111etNENEeeeeeNBEeeddd TTT2T111TT21XBBCBBXXNNbXNNt (3-23) 43 PAGE 44 In this expression Teg1g2X=uu and the summation is over the number of elements (NE) and the number of boundary elements (NBE) in the grid used for analysis. As in the finite element method these element matrices can be assembled to form global equations, where the left hand side terms form a stiffness matrix and the right hand side a force vector. The boundary conditions arising in micromechanical analysis are periodic boundary conditions which can be imposed in this method similar to FEM as the RVE is always a rectangle/cuboid and the nodes are always guaranteed to be on the boundaries of the unit cell. The effective properties of the composite microstructure are computed by performing six different analyses on the RVE subjected to different unit strains in each case. The average values of stresses and strains in the RVE are computed for each case and are then used to determine the effective properties. A detailed discussion on the periodic boundary conditions and computation of effective properties is presented in chapter 5. 3.3 Solution Structure for Grid Refinement In many practical problems local refinement is needed to capture the localized behavior of the solution (e.g.: stress concentrations). A uniform grid for the entire model is not suitable for this purpose. In order to be able to refine a grid at required locations, a solution structure based on approximate step function is presented in this chapter. The approximate step function blends the solutions from the coarse grid and refined grid such that the continuity of the solution across the refined grid and coarse grid interface is maintained. When higher order continuous approximations are used for approximating the field variable, a suitable step function needs to be used to maintain the continuity of the solution across the refine grid and course grid interface. In Figure 3-10, the coarse grid is represented by g1 and the locally refined grid is represented by g2, 44 PAGE 45 it is to be noted that these grids g1 and g2 are conceptual grids and analysis is performed on the grid g. 3.3.1 Solution Structure The approximate step function that is used to blend the solutions from the coarse and the refined grids is constructed by using the implicit equations of the coarse/refined grid interface. This step function for refinement is defined as: 0012022()1212221krefkH (3-24) In the above equation, is the implicit equation of the coarse/refined grid interface, is the range parameter which defines the region in which the step function for refinement varies from zero to one, and is an integer which defines the order of the step function used. The step function for refinement is such that it takes on a value of zero at k 0 12 at 2 and one at It has a zero slope and continuity at both 1kC 0 and Hence, when B-spline basis is used to approximate the solution, an appropriate needs to be chosen so as to maintain the continuity of the solution across the coarse/refined grid interface. The plot of and is shown in Figure 3-11, and it can be seen that the above mentioned properties are satisfied by the step function of refinement. k refH 1refH The above definition for is valid only when the coarse/refined grid interface is defined by a single implicit equation. However, in many practical cases the interface is defined by refH 45 PAGE 46 multiple implicit equations as shown in Figure 3-12. In order to define ref H for multiple implicit equations, a Boolean tree as shown in Figure 3-12b is used to represent the step function for refinement and the Boolean operations used to combine two step functions 1refH and 2refH denoted as and respectively are be defined as 12refH 12refH 12refH 121refrefrefHHH 2 (3-25a) 121refrefrefHHH 2 (3-25b) 12121refrefrefrefrefHHHH 2H (3-25c) The step function for refinement for the model shown in Figure 3-12 is shown in Figure 3-13. This step function is constructed by using the Boolean operations as defined in Eq. (3-32). The solution structure for the displacements is represented in the following form: (1)refrefHH1uu 2ua2a2 (3-26) In the above expression, is the step function for the coarse and refined grid interface and it vanishes outside the refined grid and transitions from zero to unity inside the refined grid. The transition from zero to unity takes place only across the boundary of the coarse/refined grid interface. In Figure 3-10, the refined step functionvaries from zero to unity across the edges represented by AB and BC as these edges are shared between the refined elements of grid and unrefined elements of grid .The refined step function is unity and hence does not vary across the edges AF, FE, ED and DC as these edges are not shared between the refined and unrefined grids. The expressions for and have solution structure as proposed for imposition of essential boundary conditions and can be written as and. refH refH 2g 1g 1u 2u 1g1u=Du+u 2g2u=Du+u 46 PAGE 47 3.3.2 Modified Weak Form The solution structure for grid refinement using the step function for refinement, ref H can be written as shown in the previous section and the virtual displacements are written as: 1refrefHH g 1 g 2uDuD u (3-27) The grid solutions in both grids are approximated by g 11eu=NX 1 with 1refH1NN and g 22eu=NX 2 with refH2N N where represents the shape functions for grid variables. The boundary values are approximated N a1a1a1u=NX with 1refHa1NN and a2a2a2u=NX with refHa2NN where represents the basis functions used to approximate the boundary value function. Using the above approximations, the expressions for strain takes the following forms: N 1112222121221uxuxuuxx e1a11a1a2e2a2XBBBBXX Xj (3-28) In the above expression, and are modified strain displacement matrices which contain the terms involving the gradients of Dirichlet functions and the gradients of the shape functions used to approximate the grid solution and,are strain displacement matrices which contain the gradients of the shape functions used to construct boundary value functions. The expression for stress is given as 1B 2B refH a1B a2B iijC where is the elasticity stiffness matrix, and in computation of stress. The modified weak form takes the following discrete form when the above expressions for stress and strain are incorporated into the weak form (Eq. 3-29). ijC 47 PAGE 48 22111etNENEeeeeNBEeedd ee TTT11TT12XBBCBBXXFXNNt (3-29a) iiedTT12aa1a2FNNbBB d (3-29b) In this expression Teg1g2X=uu and the summation is over the number of elements (NE) and the number of boundary elements (NBE) in the grid used for analysis. As in the finite element method these element matrices can be assembled to form global equations, where the left hand side terms form a stiffness matrix and the right hand side a force vector. The Dirichlet function and the boundary value function are constructed as explained in sections 3.1.3 and 3.1.4. In this chapter, modified weak forms for the solution structures corresponding to enforcing essential boundary conditions, treatment of material discontinuity and local refinement of grid are presented in detail. These modified weak forms are used in finite element implementation for solving several engineering boundary value problems. In the next chapter, a detail discussion on the construction of B-spline based finite elements is presented. 48 PAGE 49 Figure 3-1. Solution structure for imposing essential boundary conditions Figure 3-2. Elastic solid under equilibrium 49 PAGE 50 Figure 3-3. Steady state heat transfer problem Figure 3-4. Boolean operations on Dirichlet functions 50 PAGE 51 Figure 3-5. A CSG tree representing Dirichlet function for multiple boundaries 51 PAGE 52 Figure 3-6. Boundary value function using cubic B-spline basis 52 PAGE 53 Figure 3-7. Representation of grids for material discontinuity Figure 3-8. One dimensional inclusion function for treatment of material discontinuity 53 PAGE 54 Figure 3-9. Two dimensional inclusion function for an elliptic inclusion in a square matrix Figure 3-10. Representation of solution structure for local grid refinement 54 PAGE 55 Figure 3-11. Approximate step function for refinement Figure 3-12. Coarse and refined grid interface with multiple implicit equations 55 PAGE 56 Figure 3-13. Step function for refinement for a 2D model 56 PAGE 57 CHAPTER 4 B-SPLINE FINITE ELEMENT FORMULATION B-splines provide a methodology to approximate a given set of points with smooth polynomial functions. The basis functions are represented by piecewise polynomials and the approximation provides a high degree of continuity depending on the basis used. For example, continuity is obtained when quadratic basis is used and continuity is obtained when cubic basis is used and so on. It is to be noted that the Lagrange interpolations which are used in traditional FEM are continuous irrespective of the basis used. B-spline basis functions are classified into uniform, non-uniform and NURBS. Uniform B-spline basis functions have equally spaced nodes in the parameter space and the basis functions are periodic. Adjacent nodes in non-uniform B-splines are not equidistant in the parameter space. NURBS are generalization of non-uniform B-splines, where the basis functions are rational functions and are capable of representing several analytical curves exactly. In this work, uniform B-splines are used for approximating the solution space. The discussion on B-splines in the remaining of this chapter is limited to uniform B-splines. 1C 2C 0C The B-spline basis functions have compact support and lead to banded stiffness matrices. B-spline approximation provides high order of continuity and is capable of providing accurate solutions with continuous gradients through out the domain. This is a major advantage over standard finite element shape functions. The B-spline basis functions form a partition of unity and are able to accurately represent the state of constant field in the domain of interest. This property is very important for convergence of the approximate solutions provided by B-splines. The basis functions used are polynomials and accurate integration can be performed by using Gauss quadrature. The B-spline approximation has good reproducing properties and when a 57 PAGE 58 cubic B-splines are used for approximating the solution, all the polynomial solutions up to third order can be represented exactly. B-spline shape (basis) functions are traditionally constructed using recursive definition (Farin, 2002). The parameter space is partitioned into elements using a knot vector (equivalent to a collection of nodes). Very general methods are available to insert knots and elevate the order (or degree) of the polynomial basis functions. This framework therefore can be used to implement h-pand k-refinement and has been studied in the context of isogeometric analysis by Hughes et al (Hughes, 2005). Instead of using recursive definitions as in geometric modeling applications, polynomial definitions for B-spline basis functions are used here. Furthermore, it is convenient to define an independent parameter space for each element in the grid instead of a continuous parameter space (or knot vector) for all elements. Therefore it is assumed here that for each element in the grid, parameter values vary within [1as is done for traditional isoparametric finite elements. This is illustrated for the one-dimensional case in Figure 4-1 where a cubic B-spline approximation is shown. In Figure 4-1, element is defined between vertex nodes 2 and 3, but the approximation defined over (referred to as its span) is controlled by four nodes (1-4) which can be refer to as its support nodes. ,1] 4.1 One Dimensional B-Spline Elements The polynomial expressions for the basis functions can be derived using the continuity requirements between adjacent elements. A order B-spline has thk 1k support nodes with some of the support nodes lying outside the element. The span corresponding to the element is defined by the basis functions of all the support nodes. The basis functions are polynomials of the form whose coefficients are determined using continuity requirements. Since the 0kjiijjNa r ija 58 PAGE 59 parameter range for the element is [1,1]r the parameter value for first support node is 12kr and for the ( support node is 1)thk 12kr for the parameterization to be uniform. The approximated function shown in Fig 4-1, is a linear combination of the basis functions within each element so that the span of the element e is represented as 000()kkkjeiiijiij i f rNuar u (4-1) In the above equation, are the support nodal values and it is assumed for convenience that the first support node for element has index 0. B-spline approximation of order has iu e k 1kC continuity at the junction between adjacent elements. This continuity requirement can be stated as 1(1)(1)mmeemmffrr 0,1...1k ,m (4-2) where, e and e+1 are adjacent elements whose derivatives are evaluated at the end point (r=-1) and start point (r=1)respectively. 0m represents the continuity equation Using Eq. (4-2), the continuity equations can be expressed as 1(1)(1)eeff 000100110!!!(1)()!()!()!!(1)0()!kkkkjmjijijijjkjmkjkjjjjauaaujmjmjmjaujm ji1 (4-3) Since this equation is valid for arbitrary values of ,0..iuik the coefficients corresponding to each should vanish. This results in a set of linearly independent equations. Since equations are obtained for each value of a total of linearly independent equations are obtained. An additional equation stating that ,0..iuik 1) 2k 2k 0,1,..1mk (2kk 59 PAGE 60 the basis functions form a partition of unity 01kiiN is used to obtain 2(1)k equations needed to solve for all coefficients of the basis functions. Using this approach, the basis functions for quadratic and cubic B-splines can be derived as shown in Eq. (4-4) and the basis functions are plotted in Figure 4-2. 1k Quadratic B-spline element. The shape functions for quadratic B-spline element are constructed using the procedure explained above. Quadratic B-spline element in one-dimension has three nodes and hence three shape functions. The plots of the shape functions are shown in Figure 4-2 and the corresponding equations are shown below: 22123111(12);(62);(12)888NrrNrNr 2r (4-4) Cubic B-spline element. The plots of basis functions for cubic B-spline are shown in Figure 4-3 where it can be seen that, unlike traditional finite element shape functions, the B-spline basis functions are not unity at the corresponding node and do not vanish at other nodes. Therefore, these basis functions do not satisfy the Kronecker-delta property and the approximation constructed using these basis functions does not interpolate nodal values. 23231223233411(133);(231533);484811(231533);(133)4848NrrrNrrNrrrNrr rr (4-5) 4.2 Two and Three Dimensional B-Spline Elements The basis functions for the higher dimensional B-spline elements are constructed by taking product of the basis functions for one-dimensional B-splines. A typical grid used for analysis is illustrated in Figure 4-4. This grid corresponds to cubic B-spline elements where the span corresponding to each element is controlled by 16 support nodes. The grid elements used here are regular quadrilaterals (rectangle/square) or regular hexahedra (cube/cuboid), hence the 60 PAGE 61 mapping for geometry of elements from parameter space to the physical space is linear. Note that even if an iso-parametric formulation is used the mapping will degenerate to a linear map because the element geometry in physical space is regular and can be obtained by scaling the element geometry in parametric space. In the following sub-sections, detailed discussion on higher dimensional quadratic and cubic B-spline elements is presented. Quadratic B-Spline elements. The basis functions for the bi-quadratic B-spline element shown in Figure 4-5 can be constructed using the product of the one-dimensional quadratic basis functions (expressed in terms of parametric coordinates (,) ). 23(1)(,)()();1..3;1..3DjiijNNNij (4-6) The span is defined by 9 support nodes, as shown in Figure 4-5(a), ,but none of these are at the vertices of the element defined in the parameter space as the square [1 ,1][1,1] It is advantageous to have nodes at the vertices of the element. In order to accomplish this, the quadratic B-spline finite element was defined as shown in Figure 4-5(b) as an assembly of four quadrants (Q1-Q4). This representation facilitates implementation because each element will have nodes at its four vertices. The mapping of the geometry of element from the parameter space to the physical space is linear, as explained earlier, and can be represented as 1122liiirr uii x x x (4-7) Equation (4-7) is written using index notation, where li x and ui x are the lower and upper bounds respectively for the nodal coordinates of the element. For each quadrant, only 9-nodes (and hence 9 basis functions) are active. For example, the active nodes of the first quadrant Q1 are based on the node numbering shown in Figure 4-5(b). The basis functions obtained by product of one-dimensional basis functions will have a local 6,7,8,10,11,12,14,15 and 16 61 PAGE 62 coordinate system (,) with origin at node 11 for Q1. In order to use a parametric coordinate with origin at the center of the element, the (, (,)rs ) coordinate system needs to be translated for each quadrant as 1:1;12:1;13:1;14:1;1QrsQrsQrsQrs (4-8) The basis functions at any given point are computed by first determining which quadrant the point lies in and then by determining the corresponding support nodes and basis functions. As an example, the basis functions for a point in Q1 are (,)rs 2261172122831101222112212322214131523(,)(1)(1);(,)(1)(1)(,)(1)(1);(,)(1)(1)(,)(1)(1);(,)(1)(1)(,)(1)(1);(,)(1)(1DDDDDDDDNrsNrNsNrsNrNsNrsNrNsNrsNrNsNrsNrNsNrsNrNsNrsNrNsNrsNrNs22216331222222345913)(,)(1)(1);(,)(,)0(,)(,)(,)(,)(,)0DDDDDDDDNrsNrNsNrsNrsNrsNrsNrsNrsNrs (4-9) 2(,)DiNrs refers to bi-quadratic basis function corresponding to node of the element and represents the one-dimensional quadratic basis function as defined in Eq. (4-4). thi ()jNr Three dimensional quadratic B-spline element can be constructed in a similar manner and has eight octants which are assembled to construct the 64-node element. The basis functions are constructed using the product of 1D basis functions and in any given octant, only 27 basis functions are active. The mapping for the geometry is done as shown in Eq. (4-7). Cubic B-Spline elements. The two dimensional cubic B-spline elements are constructed by product of one-dimensional cubic B-splines. In this case, there are 16 support nodes as shown in Figure 4-6 and there is a node at each vertex of the element ([1 ,1][1,1] in parametric space). 62 PAGE 63 24(1)(,)()();1..4;1..4DjiijNrsNrNsij (4-10) 316(1)4(1)(,,)()()();1..4;1..4;1..4DkjiijkNrstNrNsNtijk (4-11) The geometry mapping between parametric and physical space can be defined using Eq. (4-7). The lower bounding node for this element is node 6 and the upper bounding point corresponds to the node 11. Again, some of the support nodes lie outside the domain of the element and are shared with neighboring elements. The 2-D bi-cubic B-spline element has 16 basis functions (Eq. 4-10) while the 3D cubic B-spline has 64 basis functions (Eq. 4-11). In these equations, represents the cubic basis functions as defined in Eq. (4-5) ()jNr 63 PAGE 64 Figure 4-1. Approximation using 1D cubic B-splines Figure 4-2. Shape functions for 1D quadratic B-spline element Figure 4-3. Shape functions for 1D cubic B-spline element 64 PAGE 65 Figure 4-4. Typical grid used for analysis and the mapping used for geometry and field variable 65 PAGE 66 Figure 4-5. Two dimensional quadratic B-spline elements in parametric space Figure 4-6. Two dimensional cubic B-spline elements in parametric space 66 PAGE 67 CHAPTER 5 NUMERICAL IMPLEMENTATIONS The weak form was modified for various solution structures and was discussed in a previous chapter. In this chapter, a detailed discussion on the efficient and accurate computation of the terms involved in the modified weak form is presented. In the first part of this chapter, a discussion on the computation of stiffness matrix and load vector corresponding to the solution structure for imposing essential boundary conditions is presented. In the next part of this chapter, a discussion on the construction of stiffness matrices involved in the solution structure for micromechanical analysis is presented. Finally, a discussion on the computation of various terms involved in the solution structure for adaptive refinement is presented. 5.1 Numerical Implementations for Linear Elastic Problems The modified weak form in discrete form for linear elasticity was presented in a previous chapter and was shown to be: 111etNENENBEeeeeeeeeddTTTTXBCBXXFXNt T (5-1a) iiedTTaFNbB d (5-1b) The stiffness matrix is given by the expression ed TBCB which needs to be computed for all the elements present in the structured grid. A typical structured grid used for analysis is shown in Figure 5-1. The grid elements can be classified into three types and the computation of stiffness matrix depends on the type of element that is being integrated upon. In the following, a discussion on the computation of stiffness matrix in each type of element is presented. 67 PAGE 68 5.1.1 Computation of Stiffness Matrix for Internal Elements The internal elements lie completely inside the analysis domain and the strain displacement relation for two dimensional problems can be written as follows: B 11132211210..0..mNxNxNNxxB .. (5-2) In the above equation, are the shape functions used to approximate the displacement solution and is the number of nodes for the element. The stiffness matrix is computed using Gauss Quadrature as shown below: iN m 2211(,)(,)yxnnijijijmmjiww TK=BCB (5-3) In the above equation, are the number of Gauss points in x and y directions respectively, is the Gaussian weight associated to the Gauss point ,xynn iw i at which the integrand is evaluated. It is to be noted that the stiffness matrix for all the internal elements is identical and stiffness matrix is computed only once and is used for all the internal elements. This leads to a significant savings in computation cost associated with computation of stiffness matrices when compared to FEM, where stiffness matrix is computed for all the elements. 5.1.2 Computation of Stiffness Matrix for Boundary Elements The boundary elements lie partly in the analysis domain; this leads to partial elements and hence needs special techniques for computation of stiffness matrix. The strain-displacement matrix for these elements is identical to the strain-displacement matrix for internal elements as 68 PAGE 69 shown in Eq. (5-2). To compute the stiffness matrix, the partial solid region is divided into triangles as shown in Figure 5-2. The stiffness matrix is computed by adding the contribution from all the triangles generated by triangulation as shown below: 221tnmmiiK=K (5-4) In this equation, is the number of total number of triangles and tn iK is the stiffness matrix computed using Gauss Quadrature for triangle. Several algorithms are available for automatic triangulation of a given set of points on a plane. In this work, Delaunay triangulation algorithm is used for triangulation. It is to be noted that curved edges are approximated by straight lines which can leads to error in the computation of stiffness matrix. However, this error is limited to boundary elements and the errors incurred due to this approximation are similar to the errors involved in FEM where the boundaries are approximated with conforming meshes. Furthermore, when the grid is made finer, the errors will approach zero. In three dimensions, the stiffness matrix is computed by tetrahedralizing the solid region and the errors involved in this approach are same as the ones involved in triangulation in two dimensions. thi 5.1.3 Computation of Stiffness Matrix for Elements with Essential Boundaries For the boundary elements containing essential boundaries, the Dirichlet functions are active and hence the corresponding strain-displacement matrix has gradients of Dirichlet function. Since Dirichlet functions are defined on a narrow band with 0 their gradients can have large values in the small band. In order to accurately compute the contribution due to Dirichlet function and its gradient, special techniques are required. In this section, one such technique is presented and discussed in detail. The strain-displacement matrix can be written as: 69 PAGE 70 111111123221221112112122110.[]0....mNDDNxxNDDNxxNDNDDNDNxxxxB (5-5) The above equation shows that the strain-displacement matrix has gradients of both shape functions and Dirichlet functions. It is convenient to decompose the B matrix into two parts such that contains terms involving the gradient of the shape functions and contains terms with the gradients of the D function as shown below: 1[][][]BBB 2] 1[]B 2[B 1111111211113222344212122111212110..0..0..000[]0..0000..00..0..mmNxNDNxDxNDDNxDDxNNDDNxxxBD B (5-6) 1111112222321322212212121121210..00..[]0..00....mmDDNxxNDDNNxxDDDDNNxxxx2BD B (5-7) Using the above procedure, the strain-displacement matrix can be rewritten as 1122BDBDB Similar matrices can be defined for three dimensional elements also. Using the decomposition of strain displacement matrix, the stiffness matrix is decomposed into the following components: 70 PAGE 71 eeeeeVVVVdVdVdVdVTTTT111221221223TKBCBBCBBCBBCBK+K+K+K (5-8) The first component 1K is computed as described earlier by decomposing the boundary elements into triangles or tetrahedra to perform the area or volume integrals respectively. The second and third components contain gradients of Dirichlet functions which are large near the boundary when is made very small and they vanish away from the boundary, hence the terms in the integrals will be non-zero in a thin band along the boundary. This allows us to convert the volume integral into a surface integral. The contribution from the gradients of finite element shape functions can be assumed to be constant across the band but varies along the boundary and the contribution from the gradients of Dirichlet function are constant along the boundary but vary normal to the boundary. The band is represented in terms of two parameters as shown in the Figure5-3 and the terms in stiffness matrix take the following form. ,tn The stiffness matrix as decomposed in Eq. (5-8) has four components; the second component 2K can be further decomposed as shown in the following equation. 2221ETTTVTdVdddd1121212KBCBBDCDBBCB (5-9) In this expression, the terms involving only Dirichlet functions and their gradients are separated from the terms containing the basis functions and their gradients. This decomposition allows us to compute the stiffness matrix efficiently and accurately. The volume integral has been decomposed into a line integral along the normal to the boundary and a surface integral along the boundary. The normal for the boundary represented by implicit equation 0 is given 71 PAGE 72 by at the boundary. The component involving only Dirichlet function and its gradients is termed as 2C in Eq. (5-9) and this needs to be integrated only across the band (along the normal to boundary), which can be done analytically as shown in Eq. (5-10). In computation of the integral across the band, it is assumed that the gradients of the implicit function with respect to the parametric coordinates of the element are constant; this is a valid assumption as the integration is performed in a very thin band 510 and hence the gradients remain approximately constant. 1211111211121211122212212222122212111221233133133222211223323321DD11DCDCCCxx22DD11DCDCCCxx2211DD1CCDCDC2xxDDDCDCxx 2 J sJ s J sJddJs s 2C 3311133332221111211CC22 J s J sJ s (5-10) In the above equation, the terms ii J denote the components of the Jacobian matrix of the transformation described in Eq. (4.7) which transforms the element from parametric space to the physical space 123(,,)sss 123(,,) x xx It is to be noted that the Jacobian matrix is diagonal, as the grid elements are regular quadrilaterals or hexahedra. Proceeding in the similar manner as explained above, the third component of the decomposed stiffness matrix in Eq. (5-8) can be written as shown in the following equation 31ETTTVTdVdddd222232222KBCBBDCDBBCB (5-11) 72 PAGE 73 The components involving only the gradients of Dirichlet function is termed as 3C and its integral across the band is given as 221133123322121122121122322123322332211221212112214141433311414133CCCCssJJssJJdCCCCJJssssJJ 43 C (5-12) The computation of the integrals along the boundary is performed numerically by Gauss quadrature. This approach leads to accurate and efficient computation of stiffness matrix for boundary elements with essential boundaries. 5.1.4 Computation of Load Vector The right hand side of the modified weak form as shown in Eq(5-1) has two terms, the first term is a volume integration for the body loads and applied non-homogenous essential boundary conditions which is done similar to the computation of stiffness matrix as explained in the previous sections. The second term is a surface integration over the natural boundaries. The boundary is approximated as a set of line segments in 2D as shown in Figure 5-4 and surface triangles in 3D and integration is performed using Gauss Quadrature as shown below: ii211211()()LilniimiLnimidwTTfNtNtf=f (5-13) As the grid is made finer, the approximation of the boundary approaches the exact representation and the errors associated with approximation of boundaries with lines (planes) approach zero. In traditional FEM, the natural boundary is approximated by the surface formed by the nodes on the boundary and associated shape functions. This results in a load vector which has non-zero values corresponding to the nodes which are present on the surface. However, in 73 PAGE 74 the implicit boundary method, nodes are not guaranteed to be present on the surface, therefore the load vector for a boundary element with natural boundary will have non-zero values for all the nodes corresponding to this element. The computed load vector is a work-equivalent load similar to the consistent load vector used in traditional FEM. When the solid region inside the element is very small, its contribution towards the stiffness matrix will be negligible and may lead to singular stiffness matrices. To overcome this problem, the elements with low volume fraction, i.e. 5% ()(volSolidvolElement ) are identified and the nodal displacements in these elements are constrained to be an extrapolation of the displacements from the neighboring internal elements. The details of this extrapolation approach are presented in section 5.5. 5.2 Numerical Implementations for Steady State Heat Transfer The modified weak form for the solution structure for steady state heat transfer was discussed in a previous chapter and was shown to be: a0dddddTeceehcTTggeTTTaTeeBkBTNhNTBkBTNfNqNh Td (5-14) In this equation, the left hand side corresponds to stiffness matrix and the right hand side corresponds to load vector. The left hand side has two terms the first term is stiffness matrix and the second term is a contribution to the stiffness matrix due to convection boundary condition. The right hand side has four terms, the first term is load due to applied temperature on essential boundary condition, the second term is load due to heat generation in the body, the third term is due to applied heat flux on natural boundary and the last term is load due to convection on convection boundary. 74 PAGE 75 Computation of Stiffness Matrix. The stiffness matrix comprises of a volume integral and a surface integral. The volume integral is performed as explained in section 5.1.1 and 5.1.2 after subdividing the boundary elements into triangles or tetrahedra. The surface integral for load computation is evaluated as explained in section 5.1.4. For the boundary elements with essential boundaries, the [B matrix can be decomposed as shown below: ] 1111211220.[]0.mNDDNxxNDDNxxB ..2 (5-15) In the above equation, [Bmatrix contains the terms involving gradients of shape functions and also the gradients of D function. It is advantageous to decompose this matrix further to such that contains terms involving the gradient of the shape functions and contains terms with the gradients of the D function as shown below: ] 1[][][]BBB 2[]B 111112211220..0..[]0..0..mmNNDxxDNNDxx D BB (5-16) 11111222212221221211210..00..[]0..0..0..mmDNDxxNDNNDxxDDNNxx 2BD B (5-17) Using the above procedure, the gradient matrix can be rewritten as B 1DB 22DB Similar matrices can be defined for 3D elements also. Using this decomposition, the stiffness matrix is decomposed into the following components: 75 PAGE 76 eeeeeVVVVkdVkdVkdVkdTTTT111221221223TKBBBBBBBBK+K+K+K V (5-18) The first component 1K is computed as described earlier by decomposing the boundary elements into triangles or tetrahedra to perform the area or volume integrals respectively. The second and third components contain gradients of Dirichlet functions which are large near the boundary when is made very small and they vanish away from the boundary, hence the terms in the integrals will be non-zero in a thin band along the boundary. This allows us to convert the volume integral into a surface integral as explained in detail in section 5.1.3. The second term in Eq. (5-18) can be written as: 2221ETTTVTkdVkdddd1121212KBBBDDBBCB (5-19) This term is evaluated as a surface integral as explained in section 5.1.3. 5.3 Numerical Implementations for Micromechanics The modified weak form for the solution structure corresponding to micromechanics was discussed in a previous chapter and was shown to be: 22111etNENEeeeeeNBEeeddd TTT2T111TT21XBBCBBXXNNbXNNt (5-20) Teg1g2X=uu (5-21) In this section, a detailed discussion on the numerical implementations for the computation of the terms involved in this weak form is presented. The left hand side of the weak form 76 PAGE 77 corresponds to the stiffness matrix while the right hand side corresponds to the load vector. The solution structure for micro-mechanics is such that the solution from the matrix grid inside the inclusion is zero due to the step function of inclusion and similarly the solution from inclusion grid outside of the inclusion vanishes. Therefore the nodal values of the matrix grid inside the inclusion are set to zero, provided these nodes are outside the influence of inc H and the nodal values of the inclusion grid outside the inclusion are set to zero. This leads to the reduction in the total number of degrees of freedom for the analysis model. The finite element solution for the grid variables are written as and For two dimensional case, the strain vector using the strain displacement relations is written as shown below: 11[]{ggiiNug1u= }} 22[]{ggiiNug2u= 111222121221uxuxuuxx 2 g11g2uBBu (5-22) The components of the strain displacement matrix are expressed as: 12B,B 1111112211112211(1)0..0(1)(1)(1)..gincincgiigincincgiiggincincincgincgiiiiNHHNxxNHHNxxNNHHHNHNxxxx1B .. (5-23a) 77 PAGE 78 22111222212222110.0..gincincgiigincincgiiggincincincgincgiiiiNHHNxxNHHNxxNNHHHNHNxxxx2B ... (5-23b) In the above equations, the decomposed components of strain-displacement matrix have terms involving gradients of finite element shape functions g kijN x and the gradients of the inclusion function incjH x It will be advantageous to further split these components into a term containing gradients of shape function only and another term containing gradients of inclusion function only. This is accomplished as shown in the following equations. 1111112211112121(1)0..0.0(1)..0...(1)(1)..gincincgiigincincgiiggincincggincinciiiiNHHNxxNHHNxxNNHHNNHHxxxx .. 11112B=B+B= (5-24) 22111222212221210..0.0..0....gincincgiigincincgiiggincincggincinciiiiNHHNxxNHHxxNNHHNNHHxxxx ...N 22122B=B+B= (5-25) It can be noticed that the modified strain displacement matrices ,,and vanish in certain regions of the computational domain and are non zero in other parts of computational domain depending on where the inclusion function and its gradients vanish. 11B 12B 21B 22B 78 PAGE 79 Therefore the computational domain is divided into three regions as shown in Figure 5-5 and they are explained in detail in the following. 5.3.1 Computation of Stiffness Matrix in the Matrix Region In the matrix region of the microstructure, the degrees of freedom for grid corresponding to matrix are active, therefore the strain-displacement matrix as shown in Eq. (5-23) becomes 2111B=BB=B0 with 0incH in Eqs. (5-24) and (5-25). The stiffness matrix will take the following form in matrix region: 22eedV TmatTmat111111BCB0K=BBCBBdV=00 (5-26) This stiffness matrix is computed by performing the volume integration over the elements in matrix region. The grid element in the matrix region will be an internal element and the volume integration is performed as discussed in the previous section 5.1.1. 5.3.2 Computation of Stiffness Matrix in the Inclusion Region This region corresponds to the inclusion region. In this region the degrees of freedom for grid corresponding to inclusion are active, therefore the strain-displacement matrix as shown in Eq. (5-23) becomes 212 1 B= with BB=0B1incH in Eqs. (5-24) and (5-25). The stiffness matrix takes the following form: 22eedV Tinc11Tinc212100K=BBCBBdV=0BCB (5-27) The stiffness matrix is computed by performing the volume integration over the elements in inclusion region. The grid element in this region will be an internal element and the volume integration is performed as discussed in the previous sections 5.1.1. 79 PAGE 80 5.3.3 Computation of Stiffness Matrix in the Interface Region In this region, the step function for the inclusion is active. This means the step function is varying in this region and strain-displacement matrix as shown in Eq. (5-23) becomes 21B=BB The stiffness matrix is expressed as: 22eedV TTT1112111211TT21222122KKBCBBCBK=BBCBBdV=KKBCBBCB (5-28) In this expression depending on which region is being integrated upon. Each of the terms in this matrix can be decomposed as a sum of volume integral and a surface integral. The procedure for this decomposition is described in detail for the first term in the stiffness matrix. orincmatC=CC 11K 2222eeedVdVTT111111121112TTTT111111111111K=BCBdV=B+BCB+B=BCB+BCB+BCB+BCB (5-29) In this expression, the first term is integrated as a volume integral after the partial elements are subdivided into triangles or tetrahedral as explained in section 5.1.2. The remaining three terms are converted into surface integrals, as these terms are non-zero in a thin band in which the inclusion function is non-zero. The strain-displacement matrices shown in Eq. (5-25) are rewritten in the following form: 11111122111122220..0..0..0......gginciiggincincinciiggggincinciiiiNNHxxNNHHxxNNNNHHxxxx2121BB H (5-30) 80 PAGE 81 2112222222221210..00..0..00....incincgigincincgiigiincincincincggiiHHNxxNHHNxxNHHHHNNxxxx2222BH B (5-31) Using the above procedure, the strain-displacement matrix can be rewritten as 2incH21222BBH B Similarly, the first component of strain-displacement matrix is written as 11incH11111BB HB These matrices can be defined for 3D elements also in a similar manner. The last three components in stiffness matrix as shown in Eq. (5-29) contain gradients of inclusion functions which are large near the inclusion boundary when is made very small and they vanish away from the boundary. Therefore the terms in these integrals will be non-zero only in a thin band along the boundary. This allows the volume integral to be converted to a surface integral. The contribution from the gradients of the shape functions can be assumed to be constant across the band but varies along the boundary and the contribution from the gradients of inclusion function is constant along the boundary but varies normal to the boundary. The band is represented in terms of two parameters as shown in the Figure 5-3. The stiffness matrix as decomposed in Eq. (5-29) has four components; the second component can be rewritten as shown below. ,tn EdVT1112BCB 12122111eeincincdVHdVHddTT1112112TT1121112BCBBCHBBCHBBC dB (5-32) 81 PAGE 82 In this expression, the terms involving only inclusion functions and its gradients are separated from the terms containing the shape functions and their gradients. This decomposition leads to efficient and accurate computation of the stiffness matrix. The volume integral has been decomposed into a line (surface) integral along the normal to the boundary and a surface integral along the boundary. The normal for the boundary represented by implicit equation 0 is given by at the boundary. The component involving only inclusion function and its gradients is termed as 2C in Eq. (5-32) and is evaluated in the similar manner as explained in section 5.1.3. 5.3.4 Periodic Boundary Conditions Composite materials have periodic arrangement of cells and their behavior can be predicted by analyzing the behavior of the unit cell which is often termed as representative volume element (RVE). The variations in the micro-scale (RVE) are rapid and periodic as compared to the variations at the macro-scale. In order to analyze the behavior at the micro-scale, the RVE is analyzed with appropriate loading and boundary conditions. The RVEs tile the computational domain by translation and the neighboring cells fit into each other in both deformed and un-deformed states as shown in Figure 5-6. Hence the boundary conditions for the RVE are periodic in order to preserve the continuity of displacements, strains and stress across each RVE. The periodic boundary conditions are expressed as linear constraints and they are implemented as multipoint constraints in a typical finite element code. The periodic boundary conditions are expressed in the following form: 13FFuuconst or in general 130FFuu (5-33a) 24FFuuconst or in general 240FFuu (5-33b) 82 PAGE 83 The multipoint constraints in Implicit Boundary Method can be implemented exactly as is done in traditional finite element method (Cook et. al. 2003). It is to be noted that the RVE is always a rectangle/cuboid, hence when Lagrange interpolations are used, the nodes are guaranteed to be present on the boundaries on which multipoint constraints are imposed. Furthermore, since structured grids are used, identical nodes are always guaranteed to be present on the opposite faces and facilitate the imposition of periodic boundary conditions. However, when FEM is used for analysis, care needs to be taken to generate mesh such that there are identical nodes present on opposite faces. The nodes on the faces (or edges in 2D) of the RVE can be partitioned into two parts: independent or master nodes (s) and dependent or slave nodes (s). The nodal displacement vector can be similarly partitioned into where are displacement of the master nodes and are those of slave nodes. The displacements for the slave nodes are specified in terms of the master nodes in the following form whereis the non-zero periodic boundary conditions (e.g., Case 1: msu=[uu] mu su s1mu=Gu+g 1L 1g 11fbkuu ). Hence the displacements are expressed as : msmm11I0u=[uu]=u=GugGg (5-34) where, I denotes the identity matrix and is the zero vector. The matrices Gand are constructed globally for all the master degrees of freedom present in the analysis model. Using this expression for displacement and a similar expression of the corresponding virtual displacement 0 g m u=Gu in the weak form, the global stiffness matrix is transformed into and the load vector is increased by to impose the multipoint constraints. It is to be noted that in a general case of multipoint constraints, and gare constructed such that 'TK=GKG TGKg G 83 PAGE 84 they may be a fully populated matrix and vector respectively. However, when structured grids are used for analysis, identical nodes are guaranteed to be present on the opposite faces and hence with proper choice of node numbers, matrix can be constructed to be a diagonal matrix. This simplifies the computations involved in the transformation of the stiffness matrix and hence leads to efficient implementation. The periodic boundary conditions involved in a square array of fibers are shown in Figure 3-10 and the corresponding multipoint constraint equations are shown in Table 1. G 5.3.5 Computation of Effective Properties Effective properties of the composite materials with periodic distribution of microstructures are computed by applying far field strains to the RVE as shown in Figure 5-7. The displacements, strains and stresses are assumed to be periodic in the RVE. Hills lemma states that the average of the strain energy in RVE is same as the strain energy due to average stress and strains. It can be stated in the form :=: .This property is used to compute the effective properties from an RVE. Using the stress-strain relation H=C the effective properties are represented as The average stresses and strains are computed by solving the boundary value problem six times with different applied strains. These average stresses and strains are used in the relation H:=C: H=C to determine a row (column) of the homogenized stiffness matrix Six different strain states are used to determine all the components of The strain-stress relation for orthotropic materials is given as follows. Once the homogenized stiffness matrix is determined, the effective properties can be computed by using the strain-stress relation. HC HC 84 PAGE 85 312112332121111123222213233312322232323313112123112100010001000100000100000100000EEEEEEEEEGGG-1C (5-35) In two-dimensional problems, the effective properties in the transverse directions are determined by analyzing the model in plane strain conditions. The longitudinal modulus is computed by analyzing the model in generalized plane strain conditions. The longitudinal shear moduli are computed by analyzing the model in longitudinal shear loading conditions. Each of these cases are explained in brief in the following sub-sections. Plane strain. In plane strain formulation, the strains in the 3 x direction (longitudinal direction) are assumed to be zero and the stress-strain relation is constructed as follows: 11212121010(1)(12)12002E 2 (5-36) This stress-strain relation is used in the modified weak form shown in Eq. (3.3) to solve the boundary value problem with specified strain loading and periodic boundary conditions to compute the transverse material properties. Generalized plane strain. The homogenized material properties in the longitudinal direction are computed by using a generalized plane strain model (Li and Lim, 2005), where the 85 PAGE 86 strains in the longitudinal directions are specified. A constant strain in the 3 x is assumed and the strains and stresses are expressed as follows 1111121311221222322231323333366121212122100and0000uxCCCuCCCxCCCCuuxx 123 (5-37) Using the constant strain in 3-direction, the stresses are written as follows: 1111211321222223312661213132323331200000and 0kCCCCCCCCCCCC 3 (5-38) The principle of virtual work as shown in Eq. (3.3) is modified by choosing the virtual strain vector as to get the following equation (Eq. 5-39). In this expression, the last term is the load due to strain in the third direction. The rest of the formulation is exactly similar to plane strain formulation. 1212T 3tTkddd C TTTutub (5-39) Longitudinal shear loading. The homogenized shear properties in the longitudinal direction are computed by using a longitudinal shear loading, where there is only one degree of freedom which is the displacement in the 3 rd -direction. This assumption results in non-zero shear strains in 3-2 or 3-1 directions only while all other strain and stress components are identically zero. The equations for shear strains and shear stress are given by the following expression: 86 PAGE 87 323232323231313131310;0uxGGux (5-40) The principle of virtual work as shown in Eq. 3-3 is modified by choosing the virtual strain vector as This formulation is used to specify unit shear strains and hence the corresponding shear moduli are computed. 3231T 5.3.6 Computation of Load Vector The right hand side of the modified weak form Eq. (5-20) comprises of load terms. The first term is load due to body load and the volume integration is performed as explained in section 5.1.1 and 5.1.2. The second term is the load term due to applied traction. The shape function vector will be either 01N or 20N depending on which region the natural boundary condition is applied. The computation of load vector is performed in a similar manner as explained in section 5.1.4. 5.4 Numerical Implementations for Local Grid Refinement The modified weak form for the solution structure corresponding to local grid refinement was discussed in a previous chapter and was shown to be: 111etNENEeeeeNBEeedd ee TTT1212TT12XBBCBBXXFXNNt (5-41a) iiedTT12aa1a2FNNbBB d (5-41b) Teg1g2X=uu (5-41c) 87 PAGE 88 In this section, a discussion on the numerical implementations for the computation of the terms involved in this weak form is presented. The left hand side of the weak form corresponds to the stiffness matrix while the right hand side corresponds to the load vector. The solution structure for local grid refinement is such that the solution from the coarse grid inside the refined grid is zero due to the step function for refinement. Therefore the nodal values of the coarse grid inside the refined are set to zero, provided these nodes are outside the influence of ref H This leads to the reduction in the total number of degrees of freedom for the analysis model. The computational domain can be divided into three regions as shown in Figure 5-8. It can be observed that the solution is outside the refined grid and transitions from to inside the refined grid. It can also be observed that in the regions where step function for refinement is active, the Dirichlet functions are inactive and vice versa. This is because of the fact the step functions for refinement are active (varying from zero to one) only on the interface of the edges which are shared between the refined and unrefined grid elements and the Dirichlet functions are active only on essential boundaries and the essential boundaries cannot form the interface edges. The computation of stiffness matrix is explained in detail in the next section. 1u 1u 2u 5.4.1 Computation of Stiffness Matrix in Regions 1 and 2 The step function of refinement is not active ( 0refH ) in Region 1 and the solution structure becomes and hence the strain displacement matrix is given as follows: 1g1au=Du+u 1111111213221232221112112122110..00..[]0..;[]00..00....mmNDDNxxNDDNxxNDNDDNDNxxxx BB (5-42) This matrix is similar to the one explained in section 5.1.3 and the stiffness matrix is given as: 88 PAGE 89 eedV TT111212BCB0K=BBCBBdV=00 (5-43) The computation of stiffness matrix is performed in a similar manner as explained in sections 5.1.1-5.1.3. In the Region 2, the step function of refinement is equal to one and the solution structure becomes and the strain displacement matrix will be similar to the one shown in Eq. (5-42) and hence the stiffness matrix becomes: g2a2u=Du+u 22eeTdV T121200K=BBCBBdV=0BCB (5-44) The stiffness matrix is computed as explained in sections 5.1.1-5.1.3. 5.4.2 Computation of Stiffness Matrix in Region 3 In the Region 3, the step function of refinement is varying and the Dirichlet functions are unity. It is to be noted that there is a small region of the size where both Dirichlet functions and step function of refinement may be active, since this region is very small, the contribution to the stiffness matrix in this region is negligible. The solution structure takes the form: and hence the strain displacement matrix takes the form 221refrefHHg1a1gau=u+u+u+u 12BBB 1111112211112121(1)0..0.0(1)..0...(1)(1)..grefrefgiigrefrefgiiggrefrefggrefrefiiiiNHHNxxNHHNxNNHHNNHHxxxx ..x 11112B=BB= (5-45) 89 PAGE 90 22111222212221210..0.0..0....grefrefgiigrefrefgiiggrefrefggrefrefiiiiNHHNxxNHHxxNNHHNNHHxxxx ...N 22122B=B+B= (5-46) The strain displacement matrix takes the similar form as explained in section 5.3.3. Hence the stiffness matrix will take the following form: 22eedV TTT1112111211TT21222122KKBCBBCBK=BBCBBdV=KKBCBBCB (5-47) The computation of stiffness matrix is performed in the similar procedure as explained in detail in section 5.3.3. 5.4.3 Computation of Load Vector The load vector is computed in the similar manner as explained in section 5.3.6. 5.5 Extrapolation of Solution for Elements with Low Volume Fraction In structured grid based methods, the grid elements do not conform to the geometry of the analysis domain and hence there may be few elements in the model with very small solid region. When the volume of the solid region inside the grid elements approach zero, the stiffness matrix corresponding to these elements will have components with magnitudes which are very small in comparison to the internal element stiffness matrices. This leads to ill-conditioning of the system of linear equations and hence poor convergence behavior. In order to avoid the problems associated with grid elements having small volume fractions, an extrapolation scheme is presented in this dissertation. The degrees of freedom corresponding to the elements with small volume fractions are constrained to be an extrapolation of the degrees of freedom from the internal elements. Due to this extrapolation scheme, the degrees of freedom for external nodes will not be arbitrarily large due to near singularity of the 90 PAGE 91 stiffness matrix of partial elements. An added advantage of this procedure is that the total number of degrees of freedom is reduced and hence less computational time for solution of linear equations. The procedure for this extrapolation strategy is explained with reference to Figure 5-9.The external nodes in the grid are denoted by and the internal nodes are represented by The displacement of an external node is associated with the displacement of a set of internal nodes and the number of inner nodes that are used will depend on the approximation scheme used for the grid variable. For example, for linear Lagrange elements, two internal nodes are associated for a given external node and for cubic B-spline elements four internal nodes are associated with a given external node. on in The extrapolation (constraining of external nodes using internal nodes) can be written as follows oi1mnnnu=I(u,...,u) i (5-48) In the above expression, is a function which maps the degrees of freedom of a set of internal nodes to the displacement of an external node. This function can be expressed as an extrapolation matrix and the transformation can be represented as a matrix multiplication: I 1[]{}iondndmnu=Cu (5-49) In the above equation is the number of internal nodes used to constrain the displacements of the outer nodes and the dimension of the problem as and the size of [C is m dn ] dnm 91 PAGE 92 5.5.1 Construction of Extrapolation Matrix In this section, construction of extrapolation matrix Cis described in detail. The extrapolation matrix is constructed by imposing the continuity requirements at the boundaries of the elements. Q4/H8 interpolation. The extrapolation matrix is constructed using the continuity conditions at the nodes and is explained for a 1D case and the extrapolation matrix for higher dimensions is constructed from the 1D extrapolation. When a linear interpolation for grid variable is used, the extrapolation matrix is constructed such that the displacement along the line joining the internal and outer nodes is linear, which means that a continuity requirement is imposed for the piecewise polynomial at the intersection of internal and boundary element. In Figure 5-10, the interpolation between nodes is given by and the interpolation between nodes is given by Imposing continuity at node leads to the following equation (Eq. 5-50). Solving this equation we get and 2D the extrapolation will be as shown in Eq. 5-51 1C 1andin 2iniono1iu 1122iNuNu 2andin 122iNuNu 1C 2in 22-oiuu 121212211iiiorrdNdNdNdNuuuudrdrdrdr (5-50) 112220100201iiooiinnnnnnuvuvuv (5-51) Cubic B-spline approximation. The cubic B-spline approximations have continuity. Hence continuity is imposed between the piecewise polynomial segments between internal nodes and external nodes. 2C 3C 92 PAGE 93 In Figure 5-11, the approximation for displacements in the internal element is given by and the approximation in the boundary element is given by 11223344iiiuNuNuNuNu io 12233441iiiuNuNuNuNu The shape functions of cubic B-splines and the their third derivative is given by 23231223233411(133);(231533)484811(231533);(133)4848NrrrNrrNrrrNrr rr (5-52) 33333123333133;;;888dNdNdNdNdrdrdrdr 4181o3i (5-53) Imposing continuity at the intersection of internal element and boundary element results in the following equation and the extrapolation relation is given as: 3C 1234234112343333464iiiiiiioiiiiuuuuuuuuuuuuu (5-54) Using the above relation, the extrapolation matrix in two dimensions is given as: 112233441040604001040604oiiiiiiiioTnnnnnnnnnnuuvuvuvuvv (5-55) Quadratic B-spline approximation. The quadratic b-spline approximations have continuity. Hence continuity is imposed (Figure 5-12) between the piecewise polynomial segments between internal nodes and external nodes. The extrapolation relation between the external node and internal node is computed as: 1C 2C 11233oiiuuuu (5-56) The extrapolation in two dimensions is constructed as: 112233103030010303ooTniiiiiinuuvuvuvv (5-57) 93 PAGE 94 5.5.2 Construction of Extrapolation Matrix for a Boundary Element In the above discussion, the extrapolation matrix for a single external node is constructed. For a boundary element, there may be more than one such external node. In this section, the construction of extrapolation matrix for a boundary element is described using the example shown in Figure 5-13 In Figure 5-13, there are 3 external nodes for the element of interest and one internal node. An extrapolation matrix needs to be constructed which relates the displacements of the outer nodes using the inner nodes as shown in the following equation: 104,774112,55282,552()2()32()2uIuuuuuIuuuuuIuuuu (5-58) Using the above relations, the extrapolation matrix for the boundary element shown in Figure 5-13 is constructed as: 82774105111002010002101003uuuuuuuu i (5-59) In the above discussion, construction of extrapolation matrix is presented for a specific example and this approach can be extended for construction of extrapolation matrix for Lagrange interpolations as well as B-spline approximations in both two and three dimensions. 5.5.3 Computation of Stiffness Matrix Once the extrapolation matrix for a boundary element is constructed, it can be used to compute the stiffness matrix in terms of the inner degrees of freedom (nodes). Assuming that the external nodes are related to internal nodes with the relation and the left hand side contains the degrees of freedom for the boundary element, the strain displacement matrix is ou=Cu 94 PAGE 95 given as Hence the stiffness matrix is given by Hence if is the stiffness matrix in terms of the degrees of freedom of boundary element, the modified stiffness matrix in terms of the degrees of freedom of internal nodes is given by By using this approach, the modified stiffness matrix for each boundary element is computed and later used for assembly. Once the finite element equations are solved for the internal degrees of freedom, the C matrix can be used to compute the displacements of the external nodes. o=Bu=BCu iC dVTTKCBDBC TdVTCBDB oK iToKCKC Figure 5-1. Typical non-conforming grid used for analysis 95 PAGE 96 Figure 5-2. Triangulation of boundary element for volume integration Figure 5-3. Representation of band on essential boundaries 96 PAGE 97 Figure 5-4. Approximation of natural boundary for computation of load vector 97 PAGE 98 Figure 5-5. Division of computational domain into three regions Figure 5-6. Undeformed and deformed configuration of a set of 4 RVEs 98 PAGE 99 Figure 5-7. Representative volume element used for computation of effective properties Figure 5-8. Division of computation domain into various regions 99 PAGE 100 Figure 5-9. Typical analysis grid with elements having low volume fraction Figure 5-10. Extrapolation of external node for linear Lagrange element Figure 5-11. Extrapolation of external nodes for cubic B-spline element 100 PAGE 101 Figure 5-12. Extrapolation of external nodes for quadratic B-spline element Figure 5-13. Extapolation for boundary element 101 PAGE 102 Table 5-1. Set of periodic boundary conditions for 6 different strains Case Unit Strain Front and Back Face Left and Right Face Top and Bottom Face 1 11 112233;; f bkfbkfbkuuLuuuu 112233;;rlrlrluuuuuu 11223;;tbtbtuuuuuu 3b 2 22 112233;; f bkfbkfbkuuuuuu 112233;;rlrlrluuuuLuu 11223;;tbtbtuuuuuu 3b3 33 112233;; f bkfbkfbkuuuuuu 112233;;rlrlrluuuuuu 112233;;tbtbtbuuuuuuL 4 23 112233;; f bkfbkfbkuuuuuu 112233;;rlrlrluuuuuuL 11223;;tbtbtuuuuuu 3b5 13 112233;;fbkfbkfbkuuuuuuL 112233;;rlrlrluuuuuu 11223;;tbtbtuuuuuu 3b6 12 112233;; f bkfbkfbkuuuuuu 112233;;rlrlrluuuuuu 11223;;tbtbtuuLuuuu 3b 102 PAGE 103 CHAPTER 6 ANALYSIS AND RESULTS The solution structures developed in this dissertation have been implemented by modifying a finite element code. The first part of this chapter presents numerical examples for validating the solution structure corresponding to enforcing essential boundary conditions exactly for both structural and steady state heat transfer problems. It has been be demonstrated that B-spline elements provide good convergence capabilities with continuous stresses (heat flux) in the entire analysis domain. The solutions obtained from B-spline elements are compared with the FE results obtained using bi-linear (4-noded quadrilateral) and bi-quadratic (9-noded quadrilateral) finite elements for two dimensional problems and linear and quadratic tetrahedral elements for three dimensional problems. In the second part of this chapter, the solution structure for treatment of material discontinuity and micromechanical analysis is validated. Numerical examples are presented to study the convergence capability of this solution structure. Fiber reinforced composite is used for determination of effective properties using 2D models and 3D models. In the final part of this chapter, numerical examples are presented for validation of solution structure for local refinement of grid. The classical stress concentration example is presented to study the convergence behavior using Lagrange and B-spline approximations. 6.1 Validation of Solution Structure for Essential Boundary Conditions 6.1.1 Patch Test This example consists of two parts. The first part demonstrates that a linear solution can be exactly represented using B-spline elements using the solution structure presented in this dissertation and the second part demonstrates that a cubic solution in displacement is represented accurately with only one B-spline element using a thin cantilever beam as an example. 103 PAGE 104 Uniaxial tension with constant stress. The first part of this example demonstrates the capability of B-spline elements to represent a state of constant stress. A rectangular plate of dimensions by 1.0 is subjected to uniaxial loading of 1.0 2.0m m M Pa Only a quarter of the plate was modeled because of symmetry as shown in Figure 6-1. The material of the plate is assumed to be steel with modulus of elasticity equal to and Poissons ratio of The analytical solution is given in Eq. (6-1) 200GPa 0.3 666(5.010);(1.510)10;0;0xyxyuxv y (6-1) A non-conforming grid used for this analysis is shown in Figure 6-2. It can be seen that several nodes lie outside the analysis domain as the B-spline elements are being used for analysis. The Dirichlet functions for this problem are shown in Figure 6-3. The D-function for x-component of displacement is zero on the left edge and rapidly reaches to one in the interior of the domain as can be seen in Figure 6-3a. The D-function corresponding to the y-displacements is zero on bottom edge and reaches to one inside the interior region as shown in Figure 6.3b. Since the applied value of Dirichlet boundary condition is zero, the boundary value function is zero in the entire domain. In order to study the effect of range parameter on the quality of solution, a parameter study is performed on and the solution of x-displacements along the line AB is plotted in Fig 6-4 for varying values of it can be observed that for values of 210 one element gives an exact solution. This study demonstrates that a state of constant stress can be exactly represented using Implicit Boundary Method with B-spline approximation. The contour plot for x is shown in Figure 6-5. This example serves as a patch test for implicit boundary method. 104 PAGE 105 Deflection of thin cantilever beam. The second part of this example demonstrates the ability of cubic B-spline element to represent a cubic solution accurately. A thin cantilever beam as shown in Figure 6-6 is considered for this purpose, and according to Euler-Bernoullis beam theory, the deflection curve of for beam is a cubic polynomial and tip deflection is given by 33tipPL E I where is the end load, P L is the length of beam, E is the modulus of elasticity and I is the moment of inertia of beam cross section. The tip deflection according to Timoshenko beam theory is given as 2314538tipPtLLEI (Timoshenko and Goodier, 1969). When cubic B-spline elements are used to model a cantilever beam, we expect to get accurate tip deflection as predicted by Timoshenko beam theory when only one element is used. However, we may not get exact tip deflection as predicted by Euler-Bernoullis beam theory as shear deformations are neglected in this theory, while the plane stress elements account for shear deformations. As the aspect ratio of the beam is increased, the solution from Timoshenko beam theory and Euler-Bernoulli beam theory converge to the same solution and we expect that the solution obtained by one element cubic B-splines converges to the same solution. In order to test this, a cantilever beam with varying aspect ratio is modeled using one-element cubic B-splines. The error in tip deflection for Euler-Bernoulli theory and one element model with respect to the aspect ratio of the beam is plotted in Figure 6-7. It can be observed that the solution obtained by one element model is accurate even for smaller aspect ratio of the beam and the quality of solution does not degrade with increasing aspect ratio of the beam. It can also be observed that the solution from one element model converges to the solution from Euler-Bernoulli theory and Timoshenko beam theory at large aspect ratios of the beam. This example 105 PAGE 106 also demonstrates that cubic B-spline elements do not undergo shear locking in thin limit, and when the structures are made of thin beams, one element can be used to get acceptable solutions. 6.1.2 Validation of Extrapolation Technique In this example, the extrapolation technique presented in this dissertation is validated. A curved beam as shown in Figure 6-8 is used for analysis. The curved beam is a three-quarter circular beam which is fixed on one end and loaded on the other end. The material of the beam is assumed to be made of steel with a modulus of elasticity of and a Poissons ratio of 0.3. 200GPa The geometry for this beam is such that there are several elements on the boundary which may have very small support from the nodes lying inside the geometry of the analysis domain and hence the resulting stiffness matrix will be almost singular and can lead to ill-conditioning of the system of linear equations and hence leads to poor convergence behavior. The extrapolation technique presented in this dissertation is used to extrapolate the solution for the external nodes in terms of the internal nodes. Bilinear Lagrange interpolation are used to perform a convergence analysis where a series of models with varying grid densities are analyzed and the error in strain energy is plotted as shown in Figure 6-9. It can be observed that the convergence behavior is not smooth due to near-singular stiffness matrices. Figure 6-10 shows the error in strain energy using the extrapolation scheme. It can be observed that the convergence plot is very smooth and gives an expected behavior. When B-spline elements are used for modeling this example, the quality of solution was not severely affected even when several elements are present on the boundary with small support. This may be because B-spline elements have a wide span as compared to Lagrange elements and hence the contribution to stiffness from a neighboring element leads to non-singular stiffness matrices and hence does not lead to ill-conditioning. Figure 6-9 shows the plot 106 PAGE 107 of error in strain energy when extrapolation technique was not used for cubic B-spline elements. Figure 6-10 shows the plot of error in strain energy when extrapolation was used. It can be observed that there is no appreciable difference in the convergence plots. It is to be noted that the x-axis of the convergence plot when extrapolation technique is used shows the number of active nodes. The contour plot of y without extrapolation is shown in Figure 6-11 and the contour plot of normal stress is shown in Figure 6-12 using extrapolation technique. In both these plots, bi-linear Lagrange elements were used for analysis and it can be observed that the stress solution shows expected behavior only when extrapolation is used and hence it is necessary to use extrapolation to avoid singular stiffness matrices and to obtain accurate solutions. Figure 6-13 shows the contour plot for y using cubic B-spline elements without the use of extrapolation technique and it can be observed that accurate solution is obtained for B-spline elements even when extrapolation technique was not used. However, in all the remaining examples of this section, extrapolation technique is used for Lagrange as well and B-spline elements. 6.1.3 Convergence Analysis of a Cantilever Beam This example involves the convergence study of quadratic and cubic B-spline elements and its comparison with four-node (Lagrange bi-linear element) and nine-node quadrilateral (quadratic Lagrange element) elements used in traditional FEM. A cantilever beam of length and 0. thickness is fixed at the one end and a uniform shear load of 1.0m 2m 0.1 M Pa is applied at the other end as shown in Figure 6-14. The exact solution for this problem is shown in Eq. (6-2) (Timoshenko and Goodier, 1969). 22163264PyuLxxyEI D (6-2a) 107 PAGE 108 221345364PvyLxDxLxxEI 2 (6-2b) 22;24xxyPPDLxyyII ;0y (6-2c) A typical non-conforming grid used for analysis is shown in Figure 6-15. Even though this problem can be modeled using a conforming grid, a non-conforming grid is used in order to demonstrate the convergence behavior of the solution structure along with B-spline basis. The B-spline elements have nodes outside the elements which form an extra layer of nodes as can be seen in Figure 6-15. The consistent nodal loads are distributed on the nodes corresponding to the boundary elements on which traction is specified. A convergence analysis is performed for this model where a series of models with increasing grid density are used for analysis and the corresponding error norms are computed and plotted. The 2 L error norm represents the errors in the displacements and the expression which defines error is shown in Eq. (6-3). Energy error norm represents the error in stresses and strains and the expression for defining energy error is shown in Eq (6-4). In these expressions, u and represent displacement, strain and stress respectively, while the superscript represents exact solution and represents approximate solution. e h 122LTehehu-u.u-u d (6-3) 1212ETeheh-d (6-4) The convergence plot for error norm and energy error norms are computed for various grid densities and shown in Figure 6-16 and Figure 6-17 respectively. 2L 108 PAGE 109 The plots are shown in a log-log scale, where the y-axis represents the error norm and the x-axis represents the number of nodes in the model. B-spline elements exhibit good convergence rates and for a given number of nodes, the errors are much smaller as compared to the bi-linear and quadratic Lagrange elements. In order to study the errors involved in deflection curve and stress distribution for a specific model, a 10 element model with 5 elements along the length of the beam and 2 elements through the thickness is used. Corresponding FE models with bilinear and quadratic Lagrange elements are also used for comparison purposes. Analytical solution as shown in Eq (6-2) is used a benchmark solution. The comparison of deflection curve is shown in Figure 6-18. It can be observed that the deflection obtained by quadratic and cubic B-spline elements as well as quadratic Lagrange elements used in FEM match closely with the analytical solution. The variation of normal and shear stresses along the line AB are shown in Figure 6-19 and Figure 6-20. It can be observed that B-spline elements as well as quadratic Lagrange elements used in traditional FEM give accurate solutions and match closely with the analytical solutions. However, B-spline elements provide continuous stress in the analysis domain and hence do not need any smoothing algorithms. In order to study the performance of B-spline elements in near incompressible limit, the cantilever beam shown in Figure 6-14 is analyzed under plane strain conditions with a Poissons ratio as 0.4999999. Convergence analysis is performed by analyzing a series of models with increasing grid densities and the tip deflection is plotted in Figure 6-21 for bi-linear Lagrange elements in FEM and quadratic and cubic B-spline elements. It can be observed that bi-linear Lagrange elements undergo volumetric locking while B-spline elements converge to the 109 PAGE 110 expected tip deflection. This example demonstrates that the B-spline elements do not undergo volumetric shear locking. 6.1.4 Square Plate with a Circular Hole The performance of quadratic and cubic B-spline elements in regions of stress concentration is studied in this example. An infinite plate with a central circular hole as shown in Figure 6-22 is considered for this analysis. The plate is subjected to uniform unit traction ( 1 ) along the x axis at infinity. The exact solution for this problem is shown in Eq. (6-5) (Sokolnikoff, 1956). Symmetry of the problem is exploited to model only a quadrant of the analysis domain. The radius of the hole is chosen to be 1 ( 1a ) and length and width of the plate is chosen to be 4 (). The exact solution is used to compute the tractions on the traction boundaries (right and top boundaries) and symmetry boundary conditions are applied on the bottom and left edges while the inner edge is traction free. 4lwa 242331cos2cos4cos422xarr 4a (6-5a) 24213sin2sin4sin422xyarr 4a (6-5b) 24213cos2cos4cos422yarr 4a (6-5c) 224311211coscoscos3cos31122aaaurErr r (6-5d) 22431111sinsinsin3sin31122aaavrErr r (6-5e) The typical analysis models used for convergence analysis along with the contour plot of stresses are shown in Figure 6-23 and Figure 6-24. The error analysis for this problem is 110 PAGE 111 performed by computing the 2 L errors and energy errors for a series of models with varying grid densities. For comparison purposes, errors are also computed for a series of models with FE Q4 (bi-linear) and FE Q9 (bi-quadratic) elements. The plot for 2 L error is shown in Figure 6-25 and plot for energy error is shown in Figure 6-26. Cubic B-spline elements perform much better than the quadratic B-spline elements and FE elements. For a given number of nodes, the errors for B-splines are far smaller when compared to the traditional finite elements. Figure 6-27 shows comparison with the analytical solution along the line 2 It can be observed that cubic B-spline element computes the stress concentration factor as 3.06 which is very close to the theoretical value of 3.0. 6.1.5 Thick Cylinder with Internal Pressure A hollow thick walled cylinder of internal radius 0.3and external radius 0.5 is modeled as shown in Figure 6-28 and for analysis purpose only a quarter of the cylinder is used. The Dirichlet boundary conditions are applied on the straight edges of the quarter such that radial displacement is allowed while tangential displacement is restricted. The inner curved edge is subjected to a pressure of The material is assumed to be made of steel with a modulus of elasticity as and a Poissons ratio of 0.3. The analytical solution is obtained from the governing differential equations for displacements, strains and stresses in a thick walled cylinder with open ends (Boresi et.al. 2003). The displacements and stresses as a function of radial position ris: in in 500psi 63010psi 22222211r p abruparEba (6-6a) 111 PAGE 112 222222rrarbprba (6-6b) 222222arbprba (6-6c) This example demonstrates the capability of the Implicit Boundary Method to impose zero essential boundary conditions in a non-conforming structured grid. The Dirichlet functions are constructed such that they vanish on the boundaries where essential boundary conditions are specified. The displacement solution for the radial displacement is shown in Figure 6-29 and Figure 6-30. In order to verify the accuracy of solution along the line AB, comparison of analytical solution and the computed solution for radial displacements is shown in Figure 6-31. The comparison of radial and tangential stresses is shown in Figure 6-32 and Figure 6-33. It can be observed that the solutions obtained by both quadratic and cubic B-spline elements match closely with analytical solution for displacement, radial and tangential stresses. Furthermore, the stresses computed by B-spline elements are continuous across the element boundaries. This example demonstrates that B-spline elements provide smooth and accurate solutions for stresses as well as displacements with B-spline elements. 6.1.6 Steady State Heat Transfer Analysis A series of pipes carrying hot fluids are arranged in a block to enable heat transfer as shown in Figure 6-34. Each set of alternate pipes carry fluid with temperatures 80 and 50. The surfaces of the block are exposed to air with ambient temperature of 30.The heat transfer takes place due to conduction between the pipes and the block and due to convection from the oC oC oC 112 PAGE 113 ambient air and the surfaces of the block. The conduction coefficient of the block is assumed to be 0.3and the natural convection coefficient is taken to be /oWmC 215/oWmC This problem is analyzed as a steady state heat transfer along with convection. Symmetry of the problem is exploited to model only a quarter as shown in Figure 6-34 (b). The Dirichlet functions for this problem are constructed such that they vanish at the curved edges and are unity inside the domain. The contour plot for Dirichlet function is shown in Figure 6-35. This problem demonstrates the imposition of non-zero Dirichlet boundary conditions and the boundary value function is constructed such that it takes the value of imposed boundary values at the curved edges and the contour plot of the boundary value function is plotted in Figure 6-36. The boundary value function takes on the imposed values at the curved edges and smoothly vanishes to zero in the interior of the domain. For validating the solution obtained by Implicit Boundary Method using B-spline elements, the solution is compared with the solution obtained by FE using ABAQUS package. The contour plot of the temperature distribution in the entire domain is shown in Figure 6-37, Figure 6-38 and Figure 6-39. It can be observed that the solution obtained using cubic B-spline elements has similar patterns when compared to FEM solution using bi-linear Lagrange elements. The temperature distribution along the line AB is plotted in Figure 6-40. The comparison of plots show that the temperature distribution matches exactly with the solution obtained by FEM. 6.1.7 Analysis of a 3D Clamp The clamp shown in Figure 6-41 is modeled using a three-dimensional model with 8-node hexahedral elements and B-spline elements to demonstrate some advantages of using a structured grid with undistorted elements. The structure is clamped on one end while the other end is 113 PAGE 114 subjected to traction in the Y direction. Traction of 10000 Pa is applied on the flat face as shown in Figure 6-41. The Youngs modulus is assumed to be 200 GPa while the Poissons ratio is taken as 0.3. This problem is analyzed using both Implicit Boundary Method (IBFEM) and using traditional finite element software (ABAQUS). Convergence analysis for this problem was performed by computing strain energy for a series of models with varying mesh densities as plotted in Figure 6-42. The plot is shown on semi-log scale, where the x-axis represents the number of nodes and the y-axis represents the strain energy. As the number of nodes increases both methods converge to the same solution. However, in this example, the Implicit Boundary Method provides a better estimate of the strain energy using fewer nodes and elements when cubic B-spline elements are used. For small number of elements, it is harder to accurately represent the geometry without having elements that are distorted in the finite element mesh, whereas the structured grid has regular shaped elements regardless of the grid density. The contour plot for the bending stresses is plotted in Figures 6-43, 6-44 and 6-45 using FEM, Implicit Boundary Method with linear Lagrange elements and cubic B-spline elements respectively. The results obtained using Implicit Boundary Method and relatively coarse grid is able to give an answer that is close to the answer obtained by a higher density FEM mesh. 6.1.8 Analysis of a 3D Rib Structure A 3D rib structure shown in Figure 6-46 is used to study the performance of B-spline elements and tri-linear Lagrange hexahedral elements. A constant pressure loading of 1000 is applied on the top surface and one end of the structure is clamped as shown in Figure 6-46. The material of the structure is assumed to be steel with modulus of elasticity and Poissons ratio Pa 200EG Pa 0.3 The geometry of the structure is quite complex and meshing with regular 114 PAGE 115 hexahedral elements is difficult, so finite element models with linear tetrahedral meshes are used as benchmark solutions. A sequence of models with varying grid densities is modeled and strain energy is computed for each model and the convergence in strain energy is plotted for FEM tetrahedral elements, Implicit Boundary Method linear hexahedral elements and cubic B-spline elements is plotted in Figure 6-47. It can be observed that the strain energy obtained by the sequence of models converge to the same strain energy value. In order to validate the computed stress distribution obtained by linear hexahedral elements and B-spline elements, the contour plots of the normal stress distributions is plotted in Figures 6-48, 6-49 and 6-50. It can be observed that the stress distribution matches closely with the one obtained by FEM. 6.1.9 Analysis of a 3D Dial Bracket A 3D dial bracket structure shown in Figure 6-51 is used to study the performance of B-spline elements and linear hexahedral elements. A constant pressure loading of 1000 is applied on the top flat surface and the bottom flat surface of the structure is constrained. The material of the structure is assumed to be steel with modulus of elasticity and Poissons ratio Pa 200EG Pa 0.3 The geometry of the structure is quite complex and meshing with regular hexahedral elements is difficult, so finite element models with quadratic and linear tetrahedral meshes are used as benchmark solutions. Convergence analysis is performed in terms of strain energy and the results are plotted in Figure 6-52. It can be seen that the solutions obtained by B-spline elements converge to the solution obtained from very dense FE model and the number of nodes needed to reach a convergent solution for B-spline elements is significantly smaller as compared 115 PAGE 116 to the traditional tetrahedral (linear, quadratic)elements. The contour plots of bending stresses are shown in Figures 6-53, 6-54 and 6-55. It can be observed that stress distribution obtained by the Implicit Boundary Method shows similar pattern as with the finite element solutions. This example demonstrates the ease of constructing regular grid for complicated 3D geometries without compromising the quality of solution. 6.2 Validation of Solution Structure for Micromechanics 6.2.1 Square Plate with Circular Inclusion This example demonstrates the validity of the solution structure for treatment of material discontinuity. A circular inclusion of volume fraction 0.47fv as shown in Figure6-56 is used for analysis in this example. The matrix has material properties 68.3m E GPa and 0.3m and the fiber has material properties 379.3fEG P a and 0.1f and plane stress idealization is used. The solution obtained by Implicit Boundary Method (IBFEM) is compared with the solution obtained by FEM using ABAQUS software. The contour plot of x is shown in Figures 6-57 and 6-58 and it can be observed that the stress pattern as computed by FEM matches closely with the stress pattern as computed by IBFEM. The maximum and minimum values for the stresses also match very closely. Convergence analysis is performed by analyzing a sequence of models and the error in strain energy is plotted in Figure 6-59. The solution from a highly dense FE model with 12001 quadratic quadrilateral elements is used to as an exact solution for convergence analysis. The strain energy from the dense FE model was computed to be and the error in strain energy for a given grid density is computed as the difference in computed strain energy and this value. It can be observed that the solution converges as the error approaches to zero when the grid density is increased. -64.27302094610 116 PAGE 117 The comparison of along the line is shown in Figure 6-60. It can be observed that the displacement solution obtained by IBFEM matches exactly with the FEM solution. A slope discontinuity can be observed at the inclusion boundary. Figure 6-61 shows the line plots xu x and y along the line AB and Figure 6-62 shows the comparison of x and y It can be observed that the solution for both stresses and strains match exactly with the FEM solution. The discontinuity in stress and strain can be observed at the inclusion interface. 6.2.2 Determination of Effective Properties of a Fiber Reinforced Composite 2D model. In this example a two-dimensional model is used for computation of effective properties of a Boron-Aluminum fiber reinforced composite with volume fraction of 0.47. For determining the effective properties in the transverse directions, a plane strain assumption is used. In order to determine the homogenized material properties in the longitudinal direction, a generalized plane strain model (Li and Lim, 2005) is used. The homogenized shear properties in the longitudinal direction are determined by using a longitudinal shear loading model with one degree of freedom in the longitudinal direction is used. This assumption results in non-zero shear strains in z-y and z-x directions only while all other strain and stress components are identically zero. The results obtained by IBFEM are compared with the results available in literature (Sun and Vaidya, 1996, Sun and Chen, 1990, Chamis, 1984 and Whitney, 1966). The periodic boundary conditions shown in Table 5-1 are used for analysis. The average stresses for each case are tabulated in Table 6-1. Figures 6-63, 6-64 and 6.65 show the stress distribution for various cases with the deformed shape. The effective properties are computed by using the average stresses and strains and are compared in Table 6-2. It can be observed that the effective 117 PAGE 118 properties determined by the current method match very closely with the effective elastic constants as reported in literature. 3D model. In this example a three-dimensional model is used for computation of effective properties of the fiber reinforced composite. The results obtained by three-dimensional analysis are compared with the results available in literature. The periodic boundary conditions shown in Table 5-1 are used for analysis. The average values of stresses and strains are tabulated in Table 6-3, and the effective elastic constants that are determined by the 3D model are tabulated in Table. 6-2. It can be observed that the effective properties determined by 3D model matches very closely with the ones determined by 2D model and with the ones available in literature. Figures 6-66, 6-67 and 6-68 show the stress distribution for various cases. 6.3 Validation of Solution Structure for Grid Refinement The local grid refinement technique presented in chapter 3 is validated in this section. The first example demonstrates the validity of solution structure for grid refinement by solving the problem with a linear displacements and constant stress and strain as exact solution. In order to demonstrate the performance of the refinement technique in the regions of stress concentrations, the classic problem of square plate with a circular hole is solved and the convergence of the solution structure is demonstrated in the second example. 6.3.1 Application of Grid Refinement to Constant Stress Problem In this example, the validity of solution structure for grid refinement is demonstrated. A rectangular plate with a uniform traction on one end as shown in Figure 6-1 is analyzed in this example. Only a quarter of the plate is modeled because of symmetry, and the exact solution for this problem is shown in Eq. (6-1). The solution structure presented in chapter 3 for local grid refinement is used in a recursive fashion to get multiple levels of refinement. Figure 6-69 shows the contour plot of X118 PAGE 119 displacements using bi-linear Lagrange elements and Figure 6-70 shows the contour plot of Y-displacements. It can be observed that the refinement is done for three levels and the exact solution is obtained without degrading the quality of the solution. It can be observed that at the first level and second levels of refinement, the central element is subdivided into 3x3 elements and at the third level of refinement, the central element is subdivided into 2x2 elements. It is to be noted that in hierarchical refinement, each element can be subdivided only into 2x2 elements, while the method presented in this dissertation does not have this limitation. The contour plots for X and Y displacements are shown in Figures 6-71 and 6-72 and it can be observed that for three levels of refinement, the solution is exactly represented. 6.3.2 Stress Concentration Problem This example involves the study of the performance of solution structure for local grid refinement in the regions of stress concentration. An infinite plate with a central circular hole as shown in Figure 6-73 is considered for this analysis. The plate is subjected to uniform unit traction ( 1x ) along the x-axis at infinity. Symmetry of the problem is exploited to model only a quadrant of the analysis domain. The radius of the hole is chosen to be 1 () and length and width of the plate is chosen to be 5 ( 1a 5lwa ). The exact solution as shown in Eq. (6-5) is used to compute the tractions on the traction boundaries (right and top boundaries) and symmetry boundary conditions are applied on the bottom and left edges while the inner edge is traction free. Bi-linear Lagrange elements and cubic B-spline elements are used in this study and convergence analysis is performed to demonstrate the convergence capabilities of the solution structure. Figure 6-74 shows the contour plot of the X-stress with 4 levels of refinement using a 119 PAGE 120 bi-linear Lagrange elements and it can be observed that the stress concentration factor is computed as 2.916 which is close to the exact value of 3.0 but not exact. Figure 6-75 shows the contour plot of x using cubic B-spline elements for 2 levels of refinement. It can be observed that very accurate solution is obtained by cubic B-spline approximations with relatively less number of elements. Convergence analysis is performed for this problem and Figure 6-76 shows the plot of errors in 2 L norm and Figure 6-77 shows the plot of energy error norm. It can be observed that in both the plots that the solution obtained with the local grid refinement solution structure converges and the error approaches to zero when the levels of refinement are increased. The B-spline elements converge faster as compared to Lagrange elements and hence very accurate solutions are obtained by B-spline approximations with relatively less number of elements. The proposed solution structure is capable of using a mix of Lagrange and B-spline elements in a given model. In order to test this, the stress concentration problem is modeled using bi-linear Lagrange elements for the coarsest mesh and for the first level of refinement and the second level of refinement is done using cubic B-spline elements. The contour plot for x is shown in Figure 6-78. It can be observed that very accurate solution is obtained by using a mix of Lagrange and B-spline elements and the stress concentration factor is computed to be 3.009 which is very close to the analytical solution of 3.0. 120 PAGE 121 Figure 6-1. Plate subjected to uniaxial loading Figure 6-2. Non-conforming grid using one cubic B-spline element 121 PAGE 122 Figure 6-3. Plot of Dirichlet functions (a) (b) uD vD 122 PAGE 123 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6x 10-6 Distance from left edge of plate, m X displacement, m Theoretical displacement = 1.0-1 = 5.0-2 = 1.0-2 = 1.0-3 = 1.0-5 Figure 6-4. Effect of range parameter on the quality of solution Figure 6-5. Contour plot of x for one element model using cubic B-splines 123 PAGE 124 Figure 6-6. Thin cantilever beam subjected to end loading 0 5 10 15 20 25 30 -10 -8 -6 -4 -2 0 2 4 6 8 10 Aspect ratio of beam% error % error in CBS deflection% error in Euler-beam deflection Figure 6-7. Plot of error in the tip deflection with varying aspect ratio 124 PAGE 125 Figure 6-8. Curved beam subjected to uniform loading 101 102 103 104 -10 0 10 20 30 40 50 60 70 80 90 Log (Num Nodes)% Error in Strain Energy IBFEM Q4IBFEM CBS Figure 6-9. Plot of error in strain energy for curved beam (without extrapolation) 125 PAGE 126 101 102 103 104 -5 0 5 10 15 20 25 30 Log (Num Nodes)% Error in Strain Energy IBFEM Q4IBFEM CBS Figure 6-10. Plot of error in strain energy for curved beam using extrapolation technique Figure 6-11. Contour plot of y for Q4 elements (without extrapolation) 126 PAGE 127 Figure 6-12. Contour plot of y for Q4 elements using extrapolation technique 127 PAGE 128 Figure 6-13. Contour plot of y using cubic B-spline elements Figure 6-14. Thick cantilever beam used for convergence analysis 128 PAGE 129 Figure 6-15. Typical non-conforming grid using cubic B-spline elements 101 102 103 104 105 10-8 10-7 10-6 10-5 Log (Num Nodes)L2 Error CBSQBSFEQ4 FEQ9 Figure 6-16. Convergence plot of 2 L error for cantilever beam 129 PAGE 130 101 102 103 104 105 10-2 10-1 100 Log (Num Nodes) Energy Error CBSQBSFEQ4 FEQ9 Figure 6-17. Convergence plot for error in strain energy for cantilever beam 130 PAGE 131 0 0.2 0.4 0.6 0.8 1 -6 -5 -4 -3 -2 -1 0x 10-5 Distance along the beam lengthDeflection along neutral axis AnalyticalCBSQBSFEQ4 FEQ9 Figure 6-18. Comparison of deflection curves for 10 element model 131 PAGE 132 -0.1 -0.05 0 0.05 0.1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2x 106 Distance along the beam thicknessNormal Stress AnalyticalCBSQBSFEQ4 FEQ9 Figure 6-19. Comparison of normal stress along the line AB for 10 element model 132 PAGE 133 -0.1 -0.05 0 0.05 0.1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2x 105 Distance along the beam thicknessShear Stress AnalyticalCBSQBSFEQ4 FEQ9 Figure 6-20. Comparison of shear stress along the line AB for 10 element model 101 102 103 104 105 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0x 10-5 Log (Num Nodes)Tip Deflection CBSQBSFEQ4Exact Figure 6-21. Convergence plot for tip deflection in near incompressible limit 133 PAGE 134 Figure 6-22. Analysis model for a square plate with a circular hole Figure 6-23. Contour plot of x using FEM Q4 elements 134 PAGE 135 Figure 6-24. Contour plot of x using cubic B-spline elements 101 102 103 104 105 10-5 10-4 10-3 10-2 Log (Num Nodes)L2 Error CBSQBSFEQ4 FEQ9 Figure 6-25. Convergence plot for 2 L error norm for stress concentration problem 135 PAGE 136 101 102 103 104 105 10-4 10-3 10-2 10-1 Log (Num Nodes) Energy Error CBSQBSFEQ4 FEQ9 Figure 6-26. Convergence plot for error in strain energy for stress concentration problem 1 1.5 2 2.5 3 3.5 4 0.5 1 1.5 2 2.5 3 3.5 Plot of X-Stress along = /2Distance from center of circular hole, mX-Stress AnalyticalCBSQBS Figure 6-27. Plot of x along the line 2 136 PAGE 137 Figure 6-28. Hollow thick cylinder subjected to uniform pressure Figure 6-29. Contour plot of radial displacements using cubic B-spline element model 137 PAGE 138 Figure 6-30. Contour plot of radial displacements using FEM Q4 elements 0.3 0.4 0.5 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3x 10-5 Distance from center of cylinder, inRadial displacement ur AnalyticalFE Q4FE Q9CBSQBS Figure 6-31. Comparison of radial displacements along the thickness of the cylinder 138 PAGE 139 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 -500 -400 -300 -200 -100 0 100 Distance from center of cylinder, inRadial stress r AnalyticalFE Q4FE Q9CBSQBS Figure 6-32. Comparison of radial stresses along the thickness of the cylinder 0.3 0.4 0.5 500 600 700 800 900 1000 1100 Distance from center of cylinder, inTangential stress t AnalyticalFE Q4FE Q9CBSQBS Figure 6-33. Comparison of tangential stresses along the thickness of the cylinder 139 PAGE 140 Figure 6-34. Heat conduction in a block with series of pipes Figure 6-35. Plot of Dirichlet function 140 PAGE 141 Figure 6-36. Plot of boundary value function Figure 6-37. Contour plot of temperature distribution using FE Q4 elements 141 PAGE 142 Figure 6-38. Contour plot of temperature distribution using cubic B-spline elements Figure 6-39. Contour plot of temperature distribution using quadratic B-spline elements 142 PAGE 143 0 0.1 0.2 0.3 0.4 0.5 45 50 55 60 65 70 75 80 85 Distance from point A along AB, mTemperature FE Q4CBSQBS Figure 6-40. Comparison of temperature along line AB Figure 6-41. Three dimensional clamp used for analysis 143 PAGE 144 102 103 104 105 4 4.05 4.1 4.15 4.2 4.25 4.3 4.35 4.4 4.45x 10-3 Log (Num Nodes)Strain Energy Abaqus H8IBFEM H8IBFEM CBS Figure 6-42. Convergence plot for strain energy for clamp model Figure 6-43. Contour plot of bending stresses using H8 elements in FEM for 3D clamp 144 PAGE 145 Figure 6-44. Contour plot of bending stresses using H8 elements for 3D clamp Figure 6-45. Contour plot of bending stresses using cubic B-spline elements for 3D clamp 145 PAGE 146 Figure 6-46. Three dimensional rib model used for analysis 102 103 104 0.215 0.22 0.225 0.23 0.235 0.24 0.245 Log (Num Nodes)Strain Energy Abaqus Linear TetraIBFEM H8IBFEM CBS Figure 6-47. Convergence plot for strain energy for rib model 146 PAGE 147 Figure 6-48. Contour plot of bending stresses using tetrahedral elements for 3D rib Figure 6-49. Contour plot of bending stresses using B-spline elements for 3D rib 147 PAGE 148 Figure 6-50. Contour plot of bending stresses using H8 elements for 3D rib Figure 6-51. Three dimensional dial bracket used for analysis 148 PAGE 149 102 103 104 105 106 6 7 8 9 10 11 12 13x 10-4 Log (Num Nodes)Strain Energy Abaqus 10-Node TetrahedronAbaqus 4-Node TetrahedronIBFEM H8IBFEM CBS Figure 6-52. Convergence plot for strain energy for dial bracket model Figure 6-53. Contour plot of bending stresses using a cubic B-spline element model for dial bracket 149 PAGE 150 Figure 6-54. Contour plot of bending stresses using a tetrahedral element model for dial bracket Figure 6-55. Contour plot of bending stresses using a H8 element model for dial bracket 150 PAGE 151 Figure 6-56. Square plate with a circular inclusion Figure 6-57. Contour plot of x using Q4 elements in IBFEM 151 PAGE 152 Figure 6-58. Contour plot of x using Q4 elements in FEM 102 103 104 10-10 10-9 10-8 Log (Num Nodes)Error in Strain Energy Figure 6-59. Convergence plot for error in strain energy for square matrix with inclusion model 152 PAGE 153 0 0.5 1 1.5 0 0.5 1x 10-8 Distance along line ABX-Displacement along line AB Abaqus uxIBFEM ux Figure 6-60. Comparison of along the line AB xu 0 0.2 0.4 0.6 0.8 1 1.2 1.4 -200 0 200 400 600 800 1000 1200 1400 Distance along line ABStress components along line AB Abaqus xxAbaqus yyIBFEM xxIBFEM yy Figure 6-61. Comparison of stresses along the line AB 153 PAGE 154 0 0.2 0.4 0.6 0.8 1 1.2 1.4 -1 -0.5 0 0.5 1 1.5 2 2.5x 10-8 Distance along line ABStrain components along line AB Abaqus xxAbaqus yyIBFEM xxIBFEM yy Figure 6-62. Comparison of strains along the line AB 154 PAGE 155 A B Figure 6-63. Plot of stresses for cases 1 and 2 (A) Plot of 11 normal stress Case 1 and (B) Plot of normal stress for Case 2 22 A B Figure 6-64. Plot of stresses for cases 3 and 4 (A) Plot of 33 normal stress for Case 1 and (B) Plot of shear stress for Case 4 23 A B Figure 6-65. Plot of stresses for cases 5 and 6 (A) Plot of 31 shear stress for Case 5 and (B) Plot of shear stress for Case 6 12 155 PAGE 156 A B Figure 6-66. Plot of stresses for cases 1 and 2 for 3D model (A) Plot of 11 normal stress for Case 1 and (B) Plot of normal stress for Case 2 22 A B Figure 6-67. Plot of stresses for cases 3 and 4 for 3D model (A) Plot of 33 normal stress for Case 1 (B) Plot of shear stress for Case 4 23 A B Figure 6-68. Plot of stresses for cases 5 and 6 for 3D model (A) Plot of 31 shear stress for Case 5 and (B) Plot of shear stress for Case 6 12 156 PAGE 157 Figure 6-69. Contour plot of u using Q4 elements and grid refinement Figure 6-70. Contour plot of using Q4 elements and grid refinement v 157 PAGE 158 Figure 6-71. Contour plot of u using cubic B-spline elements and grid refinement Figure 6-72. Contour plot of using cubic B-spline elements and grid refinement v 158 PAGE 159 Figure 6-73. Square plate with a circular hole Figure 6-74. Contour plot of x for local refinement using Q4 elements 159 PAGE 160 Figure 6-75. Contour plot of x for local refinement using cubic B-spline elements 100 101 102 103 10-4 10-3 10-2 10-1 Log (Num Nodes)L2 Error CBSQ4 Figure 6-76. Convergence plot for 2 L error in local grid refinement 160 PAGE 161 100 101 102 103 10-3 10-2 10-1 Log (Num Nodes)Energy Error CBSQ4 Figure 6-77. Convergence plot for energy error in local grid refinement Figure 6-78. Contour plot of x for local refinement using Q4 and cubic B-spline elements 161 PAGE 162 Table 6-1. Average stresses for various loading cases Case Unit Strain 11()GPa 22()GPa 33()GPa 23()GPa 31()GPa 12()GPa 1 111 230.84 41.14 41.14 0 0 0 2 221 40.35 161.20 46.22 0 0 0 3 331 40.35 46.22 161.20 0 0 0 4 231 0 0 0 46.07 0 0 5 311 0 0 0 0 54.44 0 6 121 0 0 0 0 0 54.44 162 PAGE 163 Table 6-2. Comparison of effective properties Property IBFEM 2D IBFEM 3D Ref:Sun and Vaidya, 1996 Ref: Sun and Chen, 1990 Ref: Chamis, 1984 Ref:Whitney, 1966 1() E GPa 214.83 214.85 215 214 214 215 2() E GPa 144.07 143.85 144 135 156 123 3() E GPa 144.07 143.85 144 135 156 123 23()GGPa 46.07 45.97 45.9 43.6 13()GGPa 54.44 54.35 57.2 51.1 62.6 53.9 12()GGPa 54.44 54.35 57.2 51.1 62.6 53.9 23 0.25 0.25 0.29 0.31 13 0.20 0.19 12 0.20 0.19 0.19 0.19 0.20 0.19 163 PAGE 164 Table 6-3. Average stresses for various loading cases Case Unit Strain 11()GPa 22()GPa 33()GPa 23()GPa 13()GPa 12()GPa 1 111 230.57 40.34 40.34 0 0 0 2 221 40.35 160.93 46.25 0 0 0 3 331 40.35 46.25 160.93 0 0 0 4 231 0 0 0 45.98 0 0 5 131 0 0 0 0 54.37 0 6 121 0 0 0 0 0 54.37 164 PAGE 165 CHAPTER 7 CONCLUSIONS AND FUTURE WORK 7.1 Conclusions In this dissertation, finite elements based on uniform B-spline basis defined over a structured grid are constructed to obtain solutions with continuous gradients throughout the analysis domain. Finite elements are developed using quadratic and cubic B-spline basis for problems in both two and three dimensions. B-spline basis do not satisfy Kronecker-delta property and hence Implicit Boundary Method is extended for using B-spline basis and for applying essential boundary conditions. Solution structures are developed using the approximate Heaviside step functions based on the implicit equations of the essential boundaries to impose essential boundary conditions. The weak form for linear elasticity and steady state heat transfer was modified using this solution structure such that essential boundary conditions are automatically imposed. Methods were developed for efficiently and accurately evaluating the terms involved in the modified weak form in the limit as the approximate Heaviside step functions approach the exact Heaviside step functions. The primary motivation for the method is the desire to obtain continuous gradients in the entire analysis domain using B-spline basis on a structured grid. The B-spline elements studied in this work exhibited good convergence behavior and the comparison with traditional FEM solutions showed that very accurate solutions can be obtained with relatively fewer nodes. The performance of B-spline elements does not degrade with increasing aspect ratio of the elements and one element can be used to compute accurate solutions for beam like structures. Although B-spline basis is used in this work to achieve continuous stress and strains, the method is quite general and can be applied to other basis functions which may or may not satisfy Kronecker165 PAGE 166 delta property such as Lagrange interpolations and Moving Least Squares approximation used in many meshless methods. An added advantage of the current method is the elimination of the mesh generation process during the preprocessing stage of the design cycle. Generating a uniform structured grid that encloses the geometry is easier than to generate a conforming mesh, furthermore all the elements used in the analysis are regular and undistorted unlike traditional FEM and hence problems associated with distorted elements are eliminated. As shown through several examples in the dissertation, models with fewer elements can be used for getting accurate solutions. In addition, the method would enable design engineers to perform analysis directly using solid models instead of first generating a mesh and then applying boundary conditions on the mesh, hence resulting in significant savings in pre-processing time. Another advantage of this method is that all internal elements have the same stiffness matrix resulting in significant reduction in computation required to determine the global stiffness matrix. The Implicit Boundary Method is extended for analysis of micromechanical analysis by developing a solution structure for treatment of material discontinuity in a structured grid framework. A technique for imposing multipoint constraints based on the technique for imposition of multipoint constraints in traditional FEM is presented. This approach is used for modeling an inclusion in a matrix and convergence analysis is performed to validate the solution structure. Lagrange interpolations are used for performing micromechanical analysis to determine effective properties of fiber reinforced composites and it has been shown that the solution structure presented in this work has good convergence properties and is capable of determining the effective elastic constants with very good accuracy in both two and three 166 PAGE 167 dimensions. This method holds the potential for determining effective properties of complicated microstructures using structured grids for each inclusion and a structured grid for matrix. The Implicit Boundary Method is extended for local grid refinement by constructing a solution structure based on approximate Heaviside step functions of the coarse/refined grid boundary. This solution structure is validated by analyzing the classical stress concentration problem. Convergence analysis shows that the solution structure has good convergence capabilities. In traditional adaptive techniques, multiple types of approximations (Eg. combination of Lagrange and B-spline approximations) is difficult to use, but with the proposed solution structure, it has been demonstrated that multiple types of approximations can be used. Using higher order elements (B-spline elements) in the regions of stress concentration and lower order elements (Lagrange elements) far from the regions of stress concentrations makes the analysis computationally efficient and provides accurate solutions as demonstrated in the numerical examples. 7.2 Future Work The B-spline elements developed in this dissertation are used for analysis of linear elastic and steady state heat transfer problems only. B-spline based elements can be extended for analysis of Kirchhoff plates, where the weak form requires continuous approximation space. The Implicit Boundary Method with B-spline approximations can be extending for other linear analysis types like electro-static and magneto-static problems, non-linear analysis like large elastic deformations and elasto-plastic deformations, fracture mechanics and contact mechanics. 1C The solution structure for material discontinuity can further be used to compute effective properties of complicated microstructures and for performing damage analysis of composites. The method can further be extended to incorporate B-spline approximations. 167 PAGE 168 The solution structure for adaptive refinement can be extended for performing analysis of 3D problems. The h-, pand kadaptivity can be incorporated into this method. B-spline basis provides a convenient methodology for raising the order and continuity of the approximation space, this property can be used to enable k-adaptivity and can also be used for k-version FEM. 168 PAGE 169 LIST OF REFERENCES Adams, D.F., Crane, D.A., 1984. Finite element micromechanical analysis of a unidirectional composite including longitudinal shear loading. Computers and Structures 6 1153-1165. Amestoy P.R., Davis T.A., Duff I.S., 1996. An approximate minimum degree ordering algorithm. SIAM Journal on Matrix Analysis and Applications 17(4) 886-905. Atluri S.N., Kim H.G., Cho J.Y., 1999. Critical assessment of the truly Meshless Local Petrov-Galerkin (MLPG), and Local Boundary Integral Equation (LBIE) methods. Computational Mechanics 24(5) 348-372. Atluri S.N., Zhu T.L., 2000a. Meshless local Petrov-Galerkin (MLPG) approach for solving problems in elasto-statics. Computational Mechanics 25(2-3) 169-179. Atluri S.N., Zhu T., 2000b. New concepts in meshless methods. International Journal for Numerical Methods in Engineering 47(1) 537-556. Babuska I., Banerjee U., Osborn, J.E., 2002. On principles for the selection of shape functions for the Generalized Finite Element Method, Computer Methods in Applied Mechanics and Engineering 191(49-50) 5595-5629 Belytschko T., Krongauz Y., Fleming M., Organ D., Snm L., Wing K., 1996. Smoothing and accelerated computations in the element free Galerkin method. Journal of Computational and Applied Mathematics. 74(5) 111-126. Belytschko T., Krongauz Y., Organ D., Fleming M.,Krysl P., 1996. Meshless methods: an overview and recent developments. Computer Methods in Applied Mechanics and Engineering 139(1-4) 3-47. Belytschko T., Lu Y.Y., Gu L., 1994. Element-free Galerkin methods. International Journal for Numerical Methods in Engineering 37(2) 229-256. Belytschko T., Parimi C., Moes N., Usui S., Sukumar N., 2003. Structured extended finite element methods for solids defined by implicit surfaces. International Journal for Numerical Methods in Engineering 56(4) 609-635. Boresi A.P., Schmidt R.J., Sidebottom O.M., 1993. Advanced mechanics of materials. John Wiley and Sons Inc. New York. Cai, Y.C., Zhu, H.H., 2004. Direct imposition of essential boundary conditions and treatment of material discontinuities in the EFG method. Computational Mechanics 34(4) 330-338. Chamis, C.C., 1984. Simplified composite micromechanics equations for hygral, thermal and mechanical properties. SAMPE Quarterly 14-23. Clark B.W., Anderson D.C., 2003a. The penalty boundary method. Finite Elements in Analysis and Design 39(5-6) 387-401. 169 PAGE 170 Clark B.W., Anderson D.C., 2003b. The penalty boundary method for combining meshes and solid models in finite element analysis. Engineering Computations 20(4) 344-365. Cook, R.D., Malkus D.S., Plesha M.E., Witt, R.J., 2003. Concepts and Applications of Finite Element Analysis. John Wiley and Sons Inc., New York. Cordes, L.W., Moran, B., 1996. Treatment of material discontinuity in the element-free Galerkin method. Computer Methods in Applied Mechanics and Engineering 139(1-4) 75-89. Cottrell J.A., Reali A., Bazilevs Y., Hughes T.J.R., 2006. Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering 195 (41-43) 5257-5296. Cottrell J.A., Hughes T.J.R., Reali A., 2007. Studies of refinement and continuity in isogeometric structural analysis. Computer Methods in Applied Mechanics and Engineering 196(41-44) 4160-4183. Dang, T.D., Sankar, B.V., 2007. Meshless local Petrov-Galerkin formulation for problems in composite micromechanics. AIAA Journal 45(4) 912-921. De S., Bathe K.J., 2000. Method of finite spheres. Computational Mechanics 25(4) 329-345. De S., Bathe K.J., 2001. The method of finite spheres with improved integration. Computers and Structures 79(22-25) 2183-2196. Dolbow J., Belytschko T., 1999. Numerical integration of the Galerkin weak form in meshfree methods. Computational Mechanics 23(3) 219-230. Duarte C.A., Oden J.T., 1996. h-p adaptive method using clouds. Computer Methods in Applied Mechanics and Engineering 139(1-4) 237-262. Duarte C.A., Babuska I., Oden J.T., 2000. Generalized finite element methods for three-dimensional structural mechanics problems. Computers and Structures 77(2) 215-232 Farin G., 2002. Curves and Surfaces for CAGD: A Practical Guide. Morgan Kaufmann Publishers, San Francisco. Gibson, R. F., 1994. Principles of Composite Material Mechanics. McGrawHill, New York. Gingold R.A., Monaghan J.J., 1982. Kernel estimates as a basis for general particle methods in hydrodynamics. Journal of Computational Physics 46 429-53. Hughes T.J.R., Cottrell J.A., Bazilevs Y., 2005. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering 194 (39-41) 4135-4195. Hughes T.J.R., 2000. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications, New York. 170 PAGE 171 Hllig K., Reif V., Wipper J., 2001. Weighted extended B-spline approximation of Dirichlet problems. SIAM Journal on Numerical Analysis 39(2) 442-462. Hllig K., Reif V., 2003. Nonuniform web-splines. Computer Aided Geometric Design 20(5) 277-294. Hllig K., 2003. Finite element methods with B-Splines. SIAM, Philadelphia. Kagan P., Fischer A., Bar-Yoseph P.Z., 1998. New B-spline finite element approach for geometrical design and mechanical analysis. International Journal for Numerical Methods in Engineering 41(3) 435-458. Kagan P., Fischer A., Bar-Yoseph P.Z., 2003. Mechanically based models: Adaptive refinement for B-spline finite element. International Journal for Numerical Methods in Engineering 57(8) 1145-1175. Kantorovich L.V., Krylov V.I., 1958. Approximate Methods of Higher Analysis. Interscience, New York. Kim, H.J., Swan C.C., 2003. Algorithms for automated meshing and unit cell analysis of periodic composites with hierarchical tri-quadratic tetrahedral elements. International Journal for Numerical Methods in Engineering 58(11) 1683-1711. Kohno H., Tanahashi T., 2004. Numerical analysis of moving interfaces using a level set method coupled with adaptive mesh refinement. International Journal for Numerical Methods in Fluids 45(9) 921-944. Krongauz Y., Belytschko T., 1996. Enforcement of essential boundary conditions in meshless approximations using finite elements. Computer Methods in Applied Mechanics and Engineering 131(1-2) 133-145. Krongauz, Y., Belytschko, T.,1998. EFG approximation with discontinuous derivatives. International Journal for Numerical Methods in Engineering 41(7) 1215-1233. Krongauz Y., Belytschko T., 1997. Consistent pseudo-derivatives of meshless methods. Computer Methods in Applied Mechanics and Engineering 146(3-4) 371-386. Krysl P., Grinspun E., Schroder P., 2003. Natural hierarchical refinement for finite element methods. International Journal for Numerical Methods in Engineering 56(8) 1109-1124. Krysl P., Belytschko T., 1997. Element-free Galerkin method: convergence of the continuous and discontinuous shape functions. Computer Methods in Applied Mechanics and Engineering 148(3-4) 257-277. Kumar A.V., Padmanabhan S., Burla R., 2007. Implicit boundary method for finite element analysis using non-conforming mesh or grid. International Journal for Numerical Methods in Engineering. (In press). 171 PAGE 172 Kumar A.V., Lee J., 2006. Step function representation of solid models and application to mesh free engineering analysis. Journal of Mechanical Design. Transactions of the ASME 128(1) 46-56. Kumar A.V., Lee J.H., Burla R.K., 2005. Implicit solid modeling for mesh free analysis. In: Proceedings of International Design Engineering Technical Conferences & Computers in Engineering IDETC/CIE Long Beach, California. Leung A.Y.T., Au F.T.K., 1990. Spline finite elements for beam and plate. Computers and Structures 37(5) 717-729. Li, S., Lim, S.H., 2005. Variational principles for generalized plane strain problems and their applications. Applied Science and Manufacturing 36(3) 353-365. Liszka T.J., Duarte C.A.M., Tworzydlo W.W., 1996. hp-meshless cloud method. Computer Methods in Applied Mechanics and Engineering 139(1-4) 263-288. Liu G.R., Gu Y.T., 2001. A point interpolation method for two-dimensional solids. International Journal for Numerical Methods in Engineering 50(4) 937-951. Lucy L.B.,1977. A numerical approach to the testing of the fission hypothesis. The Astronomical Journal 82(12) 1013-24. Marrey, R.V., Sankar, B.V., 1997. A Micromechanical Model for Textile Composite Plates. Journal of Composite Materials 31(12) 1187-1213. Melenk J.M., Babuska I., 1996. The partition of unity finite element method: Basic theory and applications. Computer Methods in Applied Mechanics and Engineering 139(1-4) 289-314. Moes N., Bchet E., Tourbier M., 2006. Imposing Dirichlet boundary conditions in the extended finite element method. International Journal for Numerical Methods in Engineering 67(12) 1641-1669. Monaghan J.J., 1987. An introduction to SPH. Computer Physics Communications 48(1) 89 96. Mortenson M.E., 1997. Geometric Modeling. Wiley, New York. Nayroles B., Touzat G., Villon P.,1992. Generalizing the finite element method: Diffuse approximation and diffuse elements. Computational Mechanics 10(5) 307-318. Osher S.J., Fedkiw R.P., 2002. Level Set Methods and Dynamic Implicit Surfaces. Springer Verlag, Philadelphia. Padmanabhan, S., 2006. Analysis using non-conforming structured grid and implicit boundary method. Doctoral Dissertation, University of Florida. 172 PAGE 173 Pengcheng S., Peixiang H., 1995. Bending analysis of rectangular moderately thick plates using spline finite element method. Computers and Structures 54(6) 1023-1029. Rvachev V.L., Shieko T.I., 1995. R-functions in boundary value problems in mechanics. Applied Mechanics Reviews 48(4) 151-188. Rvachev V.L., Sheiko T.I., Shapiro V., Tsukanov I., 2001. Transfinite interpolation over implicitly defined sets. Computer-Aided Geometric Design 18(3) 195-220. Sankar, B.V., and Marrey, R.V., 1997. Analytical Method for Micromechanics of Textile Composites. Composites Science and Technology 57(6) 703-713. Shapiro V., 1998. Theory of R-functions and applications: A primer. CPA Technical Report CPA88-3, Cornell Programmable Automation. Sibley School of Mechanical Engineering. Shapiro V., Tsukanov I., 1999. Meshfree simulation of deforming domains. Computer Aided Design 31(7) 459-471. Sokolnikoff I.S., 1956. Mathematical Theory of Elasticity, McGraw-Hill, New York. Sukumar N., Moran B., Belytschko T., 1998. Natural element method in solid mechanics. International Journal for Numerical Methods in Engineering 43(5) 839-887. Sun, C.T., Chen, J.L., 1990. A micromechanical model for plastic behavior of fibrous composites. Composites Science and Technology 40 115-29. Sun, C.T., Vaidya, R.S., 1996. Prediction of composite properties from a representative volume element. Composites Science and Technology 56(2) 171-179. Surana K.S., Rajwani A., Reddy J.N., 2006. The k-version finite element method for singular boundary-value problems with application to linear fracture mechanics, International Journal for Computational Methods in Engineering Science and Mechanics 7(3), 217-239. Taliercio, A, 2005. Generalized plane strain finite element model for the analysis of elastoplastic composites. International Journal of Solids and Structures 42(8) 2361-2379. Timoshenko S.P., Goodier J.N., 1969. Theory of Elasticity, McGraw Hill, New York. Tyrus, J.M., Gosz, M., DeSantiago, E., 2007. A local finite element implementation for imposing periodic boundary conditions on composite micromechanical models. International Journal of Solids and Structures 44(9) 2972-2989. Vermeulen A.H., Heppler G.R., 1998. Predicting and avoiding shear locking in beam vibration problems using the B-spline field approximation method. Computer Methods in Applied Mechanics and Engineering 158(3-4) 311-327. Whitney, J.M., Riley, M.B., 1966. Elastic properties of fiber reinforced composite materials. AlAA Journal 4 1537-42. 173 PAGE 174 Xia, Z., Zhang, Y., Ellyin, F., 2003.A unified periodical boundary conditions for representative volume elements of composites and applications. International Journal of Solids and Structures 40(8) 1907-1921. Yu, W., Tang, T., 2007. A variational asymptotic micromechanics model for predicting thermoelastic properties of heterogeneous materials. International Journal of Solids and Structures 44(22-23) 7510-7525. Zhu, H., Sankar, B.V., Marrey, R.V., 1998. Evaluation of Failure Criteria for Fiber Composites Using Finite Element Micromechanics. Journal of Composite Materials 32(8) 766-782. 174 PAGE 175 BIOGRAPHICAL SKETCH Ravi Kumar Burla was born and brought up in Hyderabad, India. He graduated with a Bachelor of Technology degree in mechanical engineering and a Master of Technology degree in computer-aided design and automation from Indian Institute of Technology, Bombay, India in June 2001. After graduation, he joined Geometric Software Solutions Limited, Pune, India as a Software Engineer where he worked for 2 years. He then enrolled in the doctoral program in Mechanical and Aerospace Engineering at the University of Florida in Fall 2003. He was awarded an Alumni Fellowship. His areas of interest include application of B-splines for finite element analysis, numerical methods and computer-aided geometric design. 175 |