Ronald D. Kriz, Associate Professor
Engineering Science and Mechanics
Virginia Polytechnic Institute and State University
Blacksburg, Virginia 24061
In previous sections we studied stresses near laminate stress free edges and ply cracks that exist at an interface between two dissimilar anisotropic materials. These results were predicted by a Finite Element Model (FEM) that did not include singularities. Instead stresses gradients near the free edge and ply crack tips were approximated with a high density FEM mesh near the singurality. In this section we introduce an analytic model that predicts these singularites: 7.1) at a bimaterial interface intersecting a stress free edge, see Figure 38, and, 7.2) at a ply crack tip intersecting an interface seperating two dissimilar anisotropic media, see Figure 40. Both models are solved using Stroh's method, Ref.[24]. These singularites can be introduced into enriched FEMs that included a variety of different boundary conditions (BCs) without using high density FEM meshes. With these enriched FEMs, stress intensity factors can be solved for as unknowns along with the nodal displacements for each unique geometry and BCs. Section 8 will outline how to create a finite element model with elements "enriched" with these singularities, Benzley Ref.[29] and Heyliger & Kriz Ref.[30] .
For both models, shown in Figs. 38 and 40, it is important to note that the singularity coincides with the origin of the coordinates. Hence each model has a unique coordinate system. Since the objective here is to reproduce published results of these two models the same coordinate system established by Ting et al. will be used. This choice will require rederivation of an arbitray rotation in the 1-2 plane for orthorhombic symmetry outlined in section 4.1. The derivations below will follow the notation and procedures outlined by Ting et al, Refs. [25] and [26]. The notes below will expand on the solution introduced by Ting et al. by developing the numerical approach to solve for these singularities, These singularities were also calculated by using Green's functions at dissimilar anisotropic material interfaces terminating at a stress free edge, Tewary & Kriz Refs. [27] and [28]. .
T.C.T. Ting and S.C. Chou, Ref.[25], used Stroh's method to solve for the singularity where the interface seperates two dissimilar anisotropic regions and terminates at a stress free edge, see Figure 38. Here we only expand on the development of the numerical solution and use the same notation in Ref.[25] except the exponent on the radius, r, is changed from, δ, to -k similar to the exponent in Ref.[26].
The solution is organized into three parts: 1) transformation of the fourth order stiffness tensor, Cijkl, for arbitrary rotation in the 1-3 plane, see eqn 35 Ref.[25], 2) calcuate eigenvalues, pk, and eigen vectors, νi and τij, from equations (9) and (10) Ref.[25], and 3) calcuate exponent (singularity) on radius, r, from the determinate of equations that are created when equations (30) and (31) of Ref.[25] are combined with appropriate boundary conditions.
(9) Ting |
(10) Ting |
(11) Ting |
(12) Ting |
Below we expanded and simplify eqn. (10) into an eigenvalue problem that we can solve by using IMSL scientific subroutines. Combining equations (11) and (12) from Ref.[25] along with the rotated transformed stiffnesses, yields
(67) |
(68) |
Substituting eqn (68) into eqn (10) of Ref.[25] yields
(69) |
Note that PL exists in every term of the coefficient matrix in eqn. (69) who's determinate yields a cubic characteristic equation in PL2. Although it is possible solve for the roots of this characteristic equation, numerical solutions to this type of problem yield more accurate eigenvalues. Hence eqn.(69) is further developed into an eigenvalue problem where PL are complex eigenvalues and their corresponding complex eigenvectors, νi,L, where the ,L subscript associates the eigenvector with the appropriate eigenvalue (L=1,2,3). The L is not a free indice and the comma does not imply differentiation
The eigenvalue format is
(70) |
It is possible to expand equation (69) into the form of equation (70). First we rewrite equation (69).
(71) |
where we can write equation (71) in abbreviated notation.
(72) |
Multiply equation (72) by the inverse of the coefficient matrix, [A'], yeilds.
(73) |
where [A']-1[A'] becomes the indentity matrix [I], and matrices [B'] and [C'] are premultiplied by [A']-1 and become [B] and [C] respectively.
(74) |
The array [A] in eqn.(70) can be constructed such that
(75) |
where the nonzero components in array [A] of eqn.(70) are summarized below.
(76) |
To show eqn.(70) and eqn.(73) are equivalent requires the following expansion. First substitute eqn.(75) into eqn.(70) where the {Z} eigenvectors are partioned into Xi and Yi components.
(77) |
where
(78) |
Equation (77) can be rewritten interms of Xi and Yi components.
(79) |
Expanding eqn.(79) yields two equations
(80) |
These two equations can be combined into one equation by substituting for Xi, hence eliminating Xi, and rearranging terms.
(81) |
where comparing equations (81) and (73) yield the necessary eigenvalues and eigenvectors
(82) |
Comparing equations (81) and (73) verifies that equation (70) can be used to solve for PL as eigenvalues instead of roots to the characteristic equation derived from the determinate of the coefficient matrix in equation (69). Both eigenvalues, PL, and eigenvectors from equation (81) will be used in the next section to calculate the singular exponent. IMSL routines EVCRG are used to solve equation (70) and EPIRG is used to obtain a performance index for the solutions.
Using Stroh's method Ting et al., Ref.[25], obtained the functional form for displacements and stresses.
(30) Ting |
(31) Ting |
where again we note that ,L does not imply differentiation and L is not a free indice but is used as a subscript to identify eigenvectors associated with the L-th eigenvalue. Here L is assigned 1,2,3 corresponding to the three eigenvalues of eqn.(69).
(26) Ting |
The equations above have been slightly modified. The exponent has been changed from δ to -k and the coefficients, ML and NL have been renamed aL and a'L.
To solve for the exponent singularity, k, we create a systems of 12 linear equations from Ting's equations (30) and (31) and the boundary conditions listed below. Coefficients of these equations are written in terms of k and solved for k as a root to the determinate of the coefficient matrix of these equations. Since the solution is to find only one of the roots between 0.0 and 1.0, this root is first approximated by calculating the determinate as a function ("function-subroutine") of k where a value for k is incremented from 0.0 to 1.0 until the determinate changes sign (+/-). Within this interval a root exists. With an approximation for the root interval other IMSL routines can be used to obtain a more accurate solution for k. Next we construct the twelve equations for the boundary conditions defined in Fig. 38.
The boundary conditions are:
(83) |
From these boundary conditions 12 equations are generated whose coefficients are functions of elastic properties and the exponent on the radius vector whose origin is located at the crack tip where the singularity exists. Below we generate this coefficient matrix.
Boundary condition 1:
(84) |
Boundary condition 2:
(85) |
Boundary condition 3:
(86) |
Boundary condition 4:
(87) |
The blocks defined above are substituted into the schematic shown below.
IMSL routine ZBREN is used to find a solution (root) to the determinate of the coefficient matrix shown in Fig. 30. The ZBREN requires an initial guess where the determinate changes sign in the interval from 0 to 1.0. Hence a function subroutine, F, is setup which calculates the determinate using IMSL routines LFTRG and LFDRG.
With the interactive computer program provided above the student can reproduce the singularities tabulated below. In particular the default values set in the NPIB form of the computer program is configured to calculate the singularity of 2.9388 x 10-2 for θ'=0 and θ=75 with the Pipes-Pagano approximation for elastic properties:
Exponent of r-k |
θ'=0 | θ'=90 | θ'=-θ |
---|---|---|---|
θ=0 | - | 3.3388 x 10-2 | - |
θ=15 | 1.3528 x 10-4 | 3.2814 x 10-2 | 6.4322 x 10-4 |
θ=30 | 2.6286 x 10-3 | 2.8682 x 10-2 | 1.1658 x 10-2 |
θ=45 | 9.6461 x 10-3 | 2.0575 x 10-2 | 2.5575 x 10-2 |
θ=60 | 1.9866 x 10-2 | 1.0519 x 10-2 | 2.3346 x 10-2 |
θ=75 | 2.9388 x 10-2 | 2.6785 x 10-3 | 8.9444 x 10-3 |
θ=90 | 3.3388 x 10-2 | - | - |
T.C.T. Ting and Hoang, Ref.[26], used Stroh's method to solve for the singularity at a ply crack tip acting normal to an interface seperating two dissimilar anisotropic regions, see Figure 40.
It is important to note that singularities coincide with the coordinate system origin. Hence each interface singularity has a unique coordinate system. Since the objective here is to reproduce published results the same coordinate system established by Ting et al. will be used. This choice will require rederivation of an arbitray rotation in the 1-2 plane for orthorhombic symmetry outlined in section 4.1. The derivations below will follow the notation and procedures outlined by Ting et al, Ref.[26]. For this module we have expanded on the solution introduced by Ting et al. and developed the numerical approach to solve for these singularities.
The solution is organized into three parts: 1) transformation of fourth order stiffness tensor, Cijkl, for arbitrary rotation in the 1-3 plane, see eqn 35 Ref.[25], 2) calcuate eigenvalues, pk, and eigen vectors, νi and τij, from equations (9) and (10) Ref.[25], and 3) calcuate exponent, k (singularity) on radius, r, from determinate of equations that are created when equations (21) and (22) of Ref.[26] are combined with the appropriate boundary conditions equation (36) of Ref.[26]. NOTE: The geometry shown in Fig.40 results in a set of boundary conditions different then those used in the previous section 7.1.
(88) |
From these boundary conditions 18 equations are generated. Coefficients are functions of elastic properties and the exponent on the radius vector whose origin is located at the crack tip where the singularity exists. Unlike boundary conditions for the free edge problem of section 7.1, there is an additional notation where + or - appear as superscripts on variables that indicates these variables exist in regions x3 > 0 and x3 < 0 respectively.
Partial results from Ting et al. Ref.[26] are listed below in Table 7 with the Pipes-Pagano approximation for elastic properties:
Exponent of r-k |
θ=0 | θ=15 | θ=30 | θ=45 | θ=60 | θ=75 | θ=90 |
---|---|---|---|---|---|---|---|
θ'=90 | .676635 .500212 .5 |
.676012 .502156 .500222 |
.672441 .505068 .500236 |
.660488 .504813 .500227 |
.630471 .502297 .500314 |
.567537 .500314 .500048 |
.5 |
θ'=75 | .675051 .500189 .433918 |
.665411 .500305 .446234 |
.646234 .500401 .466560 |
.616354 .500401 .485037 |
.571159 .500269 .496300 |
.5 | .499859 .499688 .4322561 |
θ'=60 | .665263 .499839 .379767 |
.642118 .500352 .404849 |
.605628 .500735 .441810 |
.557636 .500673 .477316 |
.5 | .503704 .499085 .430024 |
.498608 .497668 .371402 |
θ'=45 | .643963 .499012 .369508 |
.606474 .500092 .407075 |
.556388 .500639 .456342 |
.5 | .522801 .498293 .443947 |
.515052 .496531 .388842 |
.495691 .495008 .345965 |
θ'=30 | .609394 .498562 .391694 |
.557954 .499986 .442491 |
.5 | .543655 .498192 .445399 |
.558434 .495380 .400299 |
.533668 .493226 .364378 |
.494600 .492252 .339747 |
θ'=15 | .560499 .499279 .438704 |
.5 | .557552 .498554 .444039 |
.592813 .495627 .400073 |
.595139 .492603 .369244 |
.553905 .490548 .349787 |
.497659 .489828 .340358 |
θ'=0 | .5 | .562178 .498904 .441309 |
.609785 .496201 .397243 |
.631434 .493110 .367864 |
.620457 .490626 .350305 |
.566078 .489298 .342401 |
.5 .489038 .341108 |
We provide an interactive computer program module below where you can reproduce these results. The default values in the NPIB form will generate results highlighted in red in the Table above for θ=30 and θ'=60 which uses the Pipes-Pagano approximation for elastic properties. Unfortunately there are unresolved errors in the Fortran computer program (plysing.f) in module15b that predicts singularities listed in Table 7 above.
Singularites predicted for the stress "Free Edge" problem in section 7.1 have been reproduced, listed in Table 6, and can be used in enriched finite elements in Section 8.0 below.