**Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization**

Marco Riva

252 Finite Element Analysis – New Trends and Developments

Nanotechnology, Vol. 2, pp. 114-120.

Method. Sensors 2007, Vol. 7, pp. 697-707.

Saddle River, New Jersey.

pp. 441-451

Lee H. J., Young S. C., Lee Y.P., Jeong K.H. and Kim H.Y. (2007). Deflection of Microcantilever

Li M, Tang H.X. and Roukes M.L.(2007). Ultra-Sensitive NEMS-based Cantilevers for Sensing, Scanned Probe and Very High-Frequency Applications. Nature

Liu Chang. (2006). Foundation of MEMS. Illinois ECE series, Pearson Education, Upper

Loui A., Goericke F.T., Ratto T.V., Lee J., Hart B.R. and King W.P. (2008). The Effect of Piezoresistive Microcantilever Geometry on Cantilever Sensitivity during Surface Stress

Park D. S., Yun D.J., Cho M. W. and Shin B.C. (2007). An Experimental Study on the Fabrication of Glass-based Acceleration Sensor Body Using Micro Powder Blasting

Peiner E., Doering L. and Balke M. (2008). Silicon Cantilever Sensor forMicro/Nanoscale Dimension and Force Metrology. Technical Paper of Microsystems Technology, Vol. 14,

Pramanik C., Saha H. and Gangopadhyay U., (2006). Design Optimization of a High Performance Silicon MEMS Piezoresistive Pressure Sensor for Biomedical Applications,

Rosmazuin A.R., Badariah B., and Burhanuddin Y.M., (2008), Design and Analysis of MEMS Piezoresistive SiO2 Cantilever-based Sensor with Stress Concentration Region for

Saya D., Belaubre P., Mathieu F., Lagrange D., Pourciel J.B., Bergaud C.(2005). Si-Piezoresistive Microcantilever for Highly Integrated Parallel Force Detection

Sone H. , Okano H. and Sumio H. (2004). Picogram Mass Sensor using Piezoresistive Cantilever

Streetman B. G. and Banerjee S. K. (2006). Solid State Electronic Devices 6thEdition. Pearson

Su Y., Evans A.G.R. and Brunnshweiler.(1996). Micromachine Silicon Cantilever Paddles with Piezoresistive Readout for Flow Sensing. Journal Micromechanical

Thundat T., Chen G.Y., Warmack R.J., Allison D.P. and Wachter E.A. (1995). Vapor Detection using Resonating Microcantilever. Anal. Chemical. Vol. 67, pp. 519-521. Vashist K. S. (2007). A Review of Microcantilever for Sensing Applications. Journal of

Yoo K.A., Kim J.H., Nahm B.H., Kang C.J. and Kim Y.S. (2007) Fabrication and Characteristics of Microcantilever-based Biosensor for Detection of the Protein-Ligand

Yu X., Zhang H., Li X., Li T. and Zhang D. (2007). Design of High-Sensitivity Cantilever and Its Monolithic Integration with CMOS Circuits, IEEE Sensor Journal, Vol. 7, pp. 489-494.

Binding, Journal of Physics: Conference Series Vol. 61, pp. 1308–1311.

Biosensing Applications, ICSE 2008 Proceeding, Johor Bahru, Malaysia.

for Biosensor. Japanese Journal of Applied Physics, Vol. 43(7b), pp. 4663–4666.

Applications. Sensors and Actuators A, Vol. 123-124, pp. 23-29.

Prentice Hall, Upper Saddle River, New Jersey

Microengineering, Vol.6, pp. 69-72.

Nanotechnology Online, Vol. 3, June 2007.

by Growing Vapor Bubble. Sensors and Actuators A, Vol. 136, pp. 717–722.

Chemical Sensing. Sensors and Actuators A, Vol. 147, pp. 516-521.

Journal of Micromech. Microeng. Vol. 16, pp. 2060-2066.

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/47325

## **1. Introduction**

The astronomical instrumentation needs high level of image quality and stability. The quality of images processed by an optical instrument can be referred to the size of the spot and/or the point spread function (p.s.f.), while the stability is related to the displacement of the spot centroid during the observations.

The opto-mechanical elements are designed and manufactured in order to have enough stiffness to minimize shape deformations and flexures due to thermo-gravitational loads. Old traditional design philosophy answered to the problem with high thicknesses and related high masses. Heavier glasses means heavier supports and high gravitational dependent misalignment. The technological research is nowadys devoted to the light weighing of opto-mechanical systems either keeping enough stiffness, or actively correcting optical surfaces and/or positions. Complementary to the technological research, the development of powerful numerical tools added to an huge enlargement of CPU computing capacity have been offered a significant improvement into the engineering design enhancing the complexity and efficiency of the design phase.

Optical lens design refers to the calculation of lens construction parameters (variables) that will meet a set of performance requirements and constraints, including cost and schedule limitations. Construction parameters include surface profile types (spherical, aspheric, holographic, diffractive, etc.), and the parameters for each surface type such as radius of curvature, distance to the next surface, glass type and optionally tilt and decenter. The optical design exploits numerical raytracing techniques to maximize the design efficiency. Ray tracing is a method for calculating the path of waves or particles through a system with regions of varying propagation velocity, absorption characteristics, and reflecting surfaces. Under these circumstances, wavefronts may bend, change direction, or reflect off surfaces, complicating analysis. Ray tracing solves the problem by repeatedly advancing idealized narrow beams called rays through the medium by discrete amounts. Simple problems can be analyzed by

©2012 Riva, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. ©2012 Riva, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 2 Will-be-set-by-IN-TECH 254 Finite Element Analysis – New Trends and Developments Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization <sup>3</sup>

propagating a few rays using simple mathematics. More detailed analyses can be performed by using a computer to propagate many rays.

opto-mechanical configuration is prepared starting from the optical design. Then Matlab® implement physical properties into the input file and feeds the MSC-Nastran® solver. The displacements are extrapolated and reorganized in order to be compliant with the Zemax® reference systems. Matlab® runs the Zemax® solver and extract the desired data (Spot radius, EE80, p.s.f., ...). In this way a first order estimation of the Instruments mechanical stability is obtained and can be easily implemented an optimization process based onto the smart modification of the physical and mechanical properties depending onto the Zemax®

Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization

255

Easy modifications can be introduced to improve the accuracy of the algorithm's results. A "qr" based Zernike fitting function has been implemented in order to allow the modeling of optical surface deformations introducing surface errors into the Zemax optical layout. This configuration can be used both for overall system analysis and for detailed object analysis like

The opto-mechanical Finite Element modeling starts from the optical model. Masses, Center of Gravity and Optical Centers of each element are extracted from the optical file. A simplified Matlab® discretization function can be used when the instrument is very simple (i.e. two or three elements); in case of complex geometries the preliminary mesh can be prepared through dedicated software (Femap® , Hypermesh® , ...); this raw discretization usually models a bench or a box that profiles the optical systems through 2D elements. Due to the high interaction level required we decide to omit the automation of the whole bench meshing.

The opto-mechanical elements are then modeled through semi-rgid elements and concentrated masses. The weight of opto-mechanical subsystem is simplified considering double the optical element mass to include the contribution of mechanical mountings. The semi-rigid element (80%) is used instead of the rigid one in order to simulate the finite stiffness of the mountings; the master node is the optical vertex and the slave are the system's CoG and the connection points. Depending onto the geometry and the interface of the instruments, sometimes it is necessary to introduce some reinforcement beams or support trusses that can be easily modeled through 1d beam elements. In some cases, in particular with large optics or

active optical element performances evaluation[11, 13].

**Figure 1.** Numerical Framework

**2.1. FEA**

*2.1.1. Mesh generator*

results.

On the other hand, structural design is nowadays mainly based onto The finite element method (FEM). It is a numerical technique for finding approximate solutions of partial differential equations (PDE) as well as of integral equations. The solution approach is based either on eliminating the differential equation completely (steady state problems), or rendering the PDE into an approximating system of ordinary differential equations, which are then numerically integrated using standard techniques such as Euler's method, Runge-Kutta, etc.

The optimization procedure refers to choosing the best element from some set of available alternatives. In the simplest case, this means solving problems in which one seeks to minimize or maximize a real function by systematically choosing the values of real or integer variables from within an allowed set. This formulation, using a scalar, real-valued objective function, is probably the simplest example; the generalization of optimization theory and techniques to other formulations comprises a large area of applied mathematics. More generally, it means finding "best available" values of some objective function given a defined domain, including a variety of different types of objective functions and different types of domains.

In this chapter we present a possible simplified application of the optimization theory to opto-mechanical design that has been integrated into a multipurpose combined opto-mechanical numerical design process. In particular we will show a single variable optimization routine oriented to minimize mass while keeping the optical displacement below a certain value [12]. In addition considering the general purpose of this book will be briefly shown the modeling techniques used in some application to simulate the performances of functional materials like SHape Memory Alloys and Piezoelectrics. This techniques has been implemented in the optimization algorithm to maximize the efficiency of this devices in the actuation of active Mirrors based onto composite materials.

## **2. Framework**

The integrated design procedure proposed exploits the huge power of numerical computation for the design of instrumentations for astronomy. The framework can be ideally seen as "server to client" communication, where a managing server code feeds input data to computing clients and extracts the desired results. The adopted software are:


The procedure obtains image quality and motion of an opto-mechanical system under thermal and/or gravitational loads. A simplified mesh (1d or 2d elements) of a possible opto-mechanical configuration is prepared starting from the optical design. Then Matlab® implement physical properties into the input file and feeds the MSC-Nastran® solver. The displacements are extrapolated and reorganized in order to be compliant with the Zemax® reference systems. Matlab® runs the Zemax® solver and extract the desired data (Spot radius, EE80, p.s.f., ...). In this way a first order estimation of the Instruments mechanical stability is obtained and can be easily implemented an optimization process based onto the smart modification of the physical and mechanical properties depending onto the Zemax® results.

Easy modifications can be introduced to improve the accuracy of the algorithm's results. A "qr" based Zernike fitting function has been implemented in order to allow the modeling of optical surface deformations introducing surface errors into the Zemax optical layout. This configuration can be used both for overall system analysis and for detailed object analysis like active optical element performances evaluation[11, 13].

**Figure 1.** Numerical Framework

## **2.1. FEA**

2 Will-be-set-by-IN-TECH

propagating a few rays using simple mathematics. More detailed analyses can be performed

On the other hand, structural design is nowadays mainly based onto The finite element method (FEM). It is a numerical technique for finding approximate solutions of partial differential equations (PDE) as well as of integral equations. The solution approach is based either on eliminating the differential equation completely (steady state problems), or rendering the PDE into an approximating system of ordinary differential equations, which are then numerically integrated using standard techniques such as Euler's method, Runge-Kutta,

The optimization procedure refers to choosing the best element from some set of available alternatives. In the simplest case, this means solving problems in which one seeks to minimize or maximize a real function by systematically choosing the values of real or integer variables from within an allowed set. This formulation, using a scalar, real-valued objective function, is probably the simplest example; the generalization of optimization theory and techniques to other formulations comprises a large area of applied mathematics. More generally, it means finding "best available" values of some objective function given a defined domain, including

In this chapter we present a possible simplified application of the optimization theory to opto-mechanical design that has been integrated into a multipurpose combined opto-mechanical numerical design process. In particular we will show a single variable optimization routine oriented to minimize mass while keeping the optical displacement below a certain value [12]. In addition considering the general purpose of this book will be briefly shown the modeling techniques used in some application to simulate the performances of functional materials like SHape Memory Alloys and Piezoelectrics. This techniques has been implemented in the optimization algorithm to maximize the efficiency of this devices in the

The integrated design procedure proposed exploits the huge power of numerical computation for the design of instrumentations for astronomy. The framework can be ideally seen as "server to client" communication, where a managing server code feeds input data to

• Matlab® is the server software: it adapts the inputs for the client codes and evaluates their

• MSC-Nastran® is the FEA client code: it receives from the server code proper models and computes mechanical results (thermo-gravitational displacements, eigen-frequencies, ...) • ABAQUS® is a FEA client code used alternatively to MSC-Nastran® in case of user

• Zemax® is the Raytracing client code: it receives the optical perturbations properly adapted and evaluates the optical performance of the system (image quality, image

The procedure obtains image quality and motion of an opto-mechanical system under thermal and/or gravitational loads. A simplified mesh (1d or 2d elements) of a possible

computing clients and extracts the desired results. The adopted software are:

a variety of different types of objective functions and different types of domains.

actuation of active Mirrors based onto composite materials.

by using a computer to propagate many rays.

etc.

**2. Framework**

outputs.

stability, ...).

defined contitutive laws.

#### *2.1.1. Mesh generator*

The opto-mechanical Finite Element modeling starts from the optical model. Masses, Center of Gravity and Optical Centers of each element are extracted from the optical file. A simplified Matlab® discretization function can be used when the instrument is very simple (i.e. two or three elements); in case of complex geometries the preliminary mesh can be prepared through dedicated software (Femap® , Hypermesh® , ...); this raw discretization usually models a bench or a box that profiles the optical systems through 2D elements. Due to the high interaction level required we decide to omit the automation of the whole bench meshing.

The opto-mechanical elements are then modeled through semi-rgid elements and concentrated masses. The weight of opto-mechanical subsystem is simplified considering double the optical element mass to include the contribution of mechanical mountings. The semi-rigid element (80%) is used instead of the rigid one in order to simulate the finite stiffness of the mountings; the master node is the optical vertex and the slave are the system's CoG and the connection points. Depending onto the geometry and the interface of the instruments, sometimes it is necessary to introduce some reinforcement beams or support trusses that can be easily modeled through 1d beam elements. In some cases, in particular with large optics or

#### 4 Will-be-set-by-IN-TECH 256 Finite Element Analysis – New Trends and Developments Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization <sup>5</sup>

active optical systems, it is important to model the whole optical surface through 2D elements connected to the main bench via 1D beams.

Where *D* is the displacements vector, *K* the stiffness matrix, and *F* is the loads vector. If we combine Equations 1 ans 2, through some algebraical operations, we obtain the general form:

Where *A*<sup>1</sup> is the vector of the maximum amplitudes of the displacements when the *g* vector is in the *xy* plane, *A*<sup>2</sup> is the vector of the maximum amplitudes of the displacements when the *g* vector is orthogonal at the *xy* plane, *α* and *β* are the angles of equation 1 and *φ* is a phase vector of the sinusoids. In order to have the whole displacement we need three known points

Equation 2 can be further detailed introducing the thermal loads simply by adding the *C* · *T* term, where *C* C is a vector constant which ties the temperature *T* at the displacement *D*. If we consider *Td* as the the vector of the displacements due to the thermal loads, Equation 3

Now with only four known points (for four unknown variables: the two amplitudes, the phase and the thermal displacement) we can have the overall displacements and the possibility to

In addition we should note that Equation 4 can be simplified if the known points are wisely

The surface Deformations are processed through a Zernike Fitting algorithm. In precision optical manufacturing, Zernike polynomials are used to characterize higher-order errors observed in interferometric analysis, in order to achieve desired system performance. They are commonly used in active and adaptive optics where they can be used to describe wavefront

to determinate the three unknown variables: the amplitudes and the phase.

1. *α* = 90, *β* = 0, *δT* = 0: gravity vector along *X* axis, no thermal load. 2. *α* = 90, *β* = 90, *δT* = 0: gravity vector along *Y* axis, no thermal load. 3. *α* = 0, *β* = 0, *δT* = 0: gravity vector along *Z* axis, no thermal load.

The most general way to express the Zernike polynomials is in the form:

*im<sup>θ</sup>* =

*n* is chosen; *n* + *m* must be even, and 0 ≤ *m* ≤ *n*. The surface error is defined as[2]:

 *R<sup>m</sup>*

Where the *n* index defines the order of the radial power so an *n* value of 5 would indicate all polynomials whose maximum radial power is *ρ*5. Only certain values for *m* are allowed once

*R<sup>m</sup>*

*R<sup>m</sup> <sup>n</sup>* (*ρ*)*e*

becomes:

decompose the thermal contribution.

4. *g* = 0, *δT* �= 0: only thermal load.

*2.1.3. Zernike fitting*

aberrations.

taken. In fact if the sampling loading condition are:

The algebraic simplification of Equations 1 and 2

*D* = *A*<sup>1</sup> · *sin*(*α*) · *sin*(*β* + *φ*) + *A*<sup>2</sup> · *cos*(*α*) (3)

Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization

257

*D* = *A*<sup>1</sup> · *sin*(*α*) · *sin*(*β* + *φ*) + *A*<sup>2</sup> · *cos*(*α*) + *Td* (4)

*D* = *X* · *sin*(*α*) · *cos*(*β*) + *X* · *sin*(*α*) · *sin*(*β*) + *Z* · *cos*(*α*) + *Td* (5)

*<sup>n</sup>* (*ρ*)*cosmθ*

*<sup>n</sup>* (*ρ*)*sinm<sup>θ</sup>* (6)

Matlab® can read the model mesh and writes the input file for the Finite Element processor modifing both the geometry and the properties (thickness of 2D and section of 1D elements) if needed. This feature is the crucial point that helps the designer to optimize the performances of the system.

Different codes can be used for different application. The interaction with ABAQUS[1] is necessary for smart structures applications thanks to the higher performances of the code in presence of user defined constitutive laws. The routines interacts with NASTRAN[6] for opto-mechanical performances prediction and optimization.

## *2.1.2. Fast extrapolation algorithm*

To evaluate the whole performances of an instrument is necessary to verify the performances within all the loading conditions. In particular it should be possible to find situations where the errors tends to auto-compensate or, on the other hand, are magnified by the interaction with other perturbations.

Here is presented a simple extraction algorithm that can be used to reduce the number of FEA analysis required. The basic idea is to perform less analyses as possible and extracts all the displacement fields in an analytical way. Doing this, one can also predict the displacements in every condition of the gravitational load with or without the thermal expansion.

**Figure 2.** Reference systems

If the *g* vector rotates in the whole 3D space the analytical definition of the displacements is more complex. First of all it is necessary to define how decompose the *g*vector in the three directions *x*, *y* and *z*. We have decided to use that shown in Figure 2. In this way we can simulate the real rotations of the telescope, the *α* angle is the declination and *β* angle represents the rotation of the telescope around the optical axis. Doing this, the components of the g-vector are:

$$\begin{cases} \mathbf{g}\_{\mathbf{x}} = & \mathbf{g} \cdot \sin(\mathfrak{a}) \cdot \cos(\mathfrak{z}) \\ \mathbf{g}\_{\mathbf{y}} = & \mathbf{g} \cdot \sin(\mathfrak{a}) \cdot \sin(\mathfrak{z}) \\ \mathbf{g}\_{\mathbf{z}} = & \mathbf{g} \cdot \cos(\mathfrak{a}) \end{cases} \tag{1}$$

The equation that ties force and displacements can be written as:

$$D = K^{-1} \cdot F \tag{2}$$

Where *D* is the displacements vector, *K* the stiffness matrix, and *F* is the loads vector. If we combine Equations 1 ans 2, through some algebraical operations, we obtain the general form:

$$D = A\_1 \cdot \sin(\mathfrak{a}) \cdot \sin(\mathfrak{\beta} + \mathfrak{\phi}) + A\_2 \cdot \cos(\mathfrak{a}) \tag{3}$$

Where *A*<sup>1</sup> is the vector of the maximum amplitudes of the displacements when the *g* vector is in the *xy* plane, *A*<sup>2</sup> is the vector of the maximum amplitudes of the displacements when the *g* vector is orthogonal at the *xy* plane, *α* and *β* are the angles of equation 1 and *φ* is a phase vector of the sinusoids. In order to have the whole displacement we need three known points to determinate the three unknown variables: the amplitudes and the phase.

Equation 2 can be further detailed introducing the thermal loads simply by adding the *C* · *T* term, where *C* C is a vector constant which ties the temperature *T* at the displacement *D*. If we consider *Td* as the the vector of the displacements due to the thermal loads, Equation 3 becomes:

$$D = A\_1 \cdot \sin(\mathfrak{a}) \cdot \sin(\mathfrak{\beta} + \mathfrak{\phi}) + A\_2 \cdot \cos(\mathfrak{a}) + T\_d \tag{4}$$

Now with only four known points (for four unknown variables: the two amplitudes, the phase and the thermal displacement) we can have the overall displacements and the possibility to decompose the thermal contribution.

In addition we should note that Equation 4 can be simplified if the known points are wisely taken. In fact if the sampling loading condition are:


The algebraic simplification of Equations 1 and 2

$$D = X \cdot \sin(\mathfrak{a}) \cdot \cos(\mathfrak{f}) + X \cdot \sin(\mathfrak{a}) \cdot \sin(\mathfrak{f}) + Z \cdot \cos(\mathfrak{a}) + T\_d \tag{5}$$

#### *2.1.3. Zernike fitting*

4 Will-be-set-by-IN-TECH

active optical systems, it is important to model the whole optical surface through 2D elements

Matlab® can read the model mesh and writes the input file for the Finite Element processor modifing both the geometry and the properties (thickness of 2D and section of 1D elements) if needed. This feature is the crucial point that helps the designer to optimize the performances

Different codes can be used for different application. The interaction with ABAQUS[1] is necessary for smart structures applications thanks to the higher performances of the code in presence of user defined constitutive laws. The routines interacts with NASTRAN[6] for

To evaluate the whole performances of an instrument is necessary to verify the performances within all the loading conditions. In particular it should be possible to find situations where the errors tends to auto-compensate or, on the other hand, are magnified by the interaction

Here is presented a simple extraction algorithm that can be used to reduce the number of FEA analysis required. The basic idea is to perform less analyses as possible and extracts all the displacement fields in an analytical way. Doing this, one can also predict the displacements in

If the *g* vector rotates in the whole 3D space the analytical definition of the displacements is more complex. First of all it is necessary to define how decompose the *g*vector in the three directions *x*, *y* and *z*. We have decided to use that shown in Figure 2. In this way we can simulate the real rotations of the telescope, the *α* angle is the declination and *β* angle represents the rotation of the telescope around the optical axis. Doing this, the components of the g-vector

> *gx* = *g* · *sin*(*α*) · *cos*(*β*) *gy* = *g* · *sin*(*α*) · *sin*(*β*)

> > *<sup>D</sup>* <sup>=</sup> *<sup>K</sup>*−<sup>1</sup> · *<sup>F</sup>* (2)

(1)

*gz* = *g* · *cos*(*α*)

⎧ ⎪⎨

⎪⎩

The equation that ties force and displacements can be written as:

every condition of the gravitational load with or without the thermal expansion.

connected to the main bench via 1D beams.

*2.1.2. Fast extrapolation algorithm*

with other perturbations.

**Figure 2.** Reference systems

are:

opto-mechanical performances prediction and optimization.

of the system.

The surface Deformations are processed through a Zernike Fitting algorithm. In precision optical manufacturing, Zernike polynomials are used to characterize higher-order errors observed in interferometric analysis, in order to achieve desired system performance. They are commonly used in active and adaptive optics where they can be used to describe wavefront aberrations.

The most general way to express the Zernike polynomials is in the form:

$$R\_n^m(\rho)\varepsilon^{im\theta} = \begin{cases} R\_n^m(\rho)\cos m\theta\\ R\_n^m(\rho)\sin m\theta \end{cases} \tag{6}$$

Where the *n* index defines the order of the radial power so an *n* value of 5 would indicate all polynomials whose maximum radial power is *ρ*5. Only certain values for *m* are allowed once *n* is chosen; *n* + *m* must be even, and 0 ≤ *m* ≤ *n*. The surface error is defined as[2]:

257

**Figure 3.** Zernike polynomial basis

$$E = \Sigma (\delta\_i - D\_i)^2 \tag{7}$$

Zemax® supports a number of capabilities under DDE. Each function/item is given a name, that is passed to the Zemax® server using the proper request protocol. Zemax® responds to each item request with requested data. Most of them are passed from Zemax® back to the

Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization

259

In the case of thermal and gravitational performances evaluation, it is necessary an adaptation of the optical model. Coordinate Breaks (CB) are introduced before and after each Optical element paying attention to the transformation sequence. This is crucial to realize a feature that allow the introduction of local optical displacements without modifying the remaining optical layout. Matlab® extracts optical displacements from Nastran (global coordinate system), optical coordinate systems from Zemax® and applies the proper transformation matrices to obtain optical displacements in Zemax® local reference systems. In case of surface deformations, the surfaces of Zemax® file are modified from "Standard" to "Zernike Fringe Sag" in order to allow Matlab® to updates the extra data tables with the Zernike coefficients

After the uploading of mechanical data into the optical file, Zemax® evaluate spot and p.s.f. information through its raytracing engine. If necessary Matlab® request a focusing optimization after having stored the focal plane distance to evaluate the relative focal variation. Finally Matlab® extract from Zemax Spot dimensions (Max and R.M.S. radius) and

The results defines respectively the image quality and stability of the overall optical layout and, if necessary can be passed to an optimization algorithm that manage the mechanical

Whereas optimization methods are nearly as old as calculus, numerical optimization reached prominence in the digital age. Its systematic application to structural design dates to its advocacy by Schmit in 1960. The success of structural optimization in the 1970s motivated the emergence of multidisciplinary design optimization (MDO) in the 1980s. Here will be presented a simplified single variable optimization approach that has been used as starting

The general optimization problem (for example minimization) wants to minimize the

*F*(*X*) (9)

*pj*(*X*) ≤ 0 *j* = 1 : *m* (10)

application (Matlab® ) in a string that must be properly managed.

*2.2.2. Zemax® computations*

obtained from the fitting algorithm.

properties of the system.

**2.3. Optimization**

*2.3.1. Problem statement*

subjected to inequality constraints:

objective function:

point[14].

centroid displacement, eventually of multiple fields.

where *i* is the node number, *δ<sup>i</sup>* is FE displacement of *ith* node, *Di* is the polynomial displacement of the *ith* node. The series of *D* can be symbolically written as *Di* = Σ*cj fij* where *cj* are the coefficients of the polynomial. The best fit coefficients can be found minimizing the error *E* respect to *cj*:

$$\frac{dE}{dc\_{\dot{j}}} = 2\Sigma(\delta\_{\dot{i}} - \Sigma c\_{\dot{j}} f\_{\dot{i}\dot{j}}) f\_{\dot{i}\dot{j}} = 0 \tag{8}$$

This system is then solved finding the best fitting coefficients through an orthogonal-triangular decomposition (Matlab® *qr* function). The first 28 coefficients of the Zernike Polynomial have been taken into account in this application.

#### **2.2. Raytracing**

#### *2.2.1. Zemax® and Matlab® data exchange*

The raytrancing software Zemax® is used to evaluate the optical performances. Zemax® has a very powerful feature which allows another software to establish a communication link to extract lens data. The idea is based onto a program that use Zemax® like a remote application to trace rays through the lenses, and then extracts the data to be sent to other programs for further analysis or computation[17].

The communication between the application and Zemax® is accomplished using Dynamic Data Exchange (DDE). DDE is a protocol defined within the Windows operating system for sharing data between programs. Two programs can establish a DDE link, with one program acting as the "server" and the other the "client". The client generally requests specific data from the server, and the server sends the data back to the client.

Two main function must be used when exchanging data with Zemax that are the link opening and closing. To establish a DDE link with Zemax® , the client program must broadcast a message to all top level windows that includes a reference to the application name, and the topic name. The topic name indicates to Zemax® what data is being requested.

Zemax® supports a number of capabilities under DDE. Each function/item is given a name, that is passed to the Zemax® server using the proper request protocol. Zemax® responds to each item request with requested data. Most of them are passed from Zemax® back to the application (Matlab® ) in a string that must be properly managed.

### *2.2.2. Zemax® computations*

6 Will-be-set-by-IN-TECH

where *i* is the node number, *δ<sup>i</sup>* is FE displacement of *ith* node, *Di* is the polynomial displacement of the *ith* node. The series of *D* can be symbolically written as *Di* = Σ*cj fij* where *cj* are the coefficients of the polynomial. The best fit coefficients can be found minimizing the

This system is then solved finding the best fitting coefficients through an orthogonal-triangular decomposition (Matlab® *qr* function). The first 28 coefficients of

The raytrancing software Zemax® is used to evaluate the optical performances. Zemax® has a very powerful feature which allows another software to establish a communication link to extract lens data. The idea is based onto a program that use Zemax® like a remote application to trace rays through the lenses, and then extracts the data to be sent to other programs for

The communication between the application and Zemax® is accomplished using Dynamic Data Exchange (DDE). DDE is a protocol defined within the Windows operating system for sharing data between programs. Two programs can establish a DDE link, with one program acting as the "server" and the other the "client". The client generally requests specific data

Two main function must be used when exchanging data with Zemax that are the link opening and closing. To establish a DDE link with Zemax® , the client program must broadcast a message to all top level windows that includes a reference to the application name, and the

topic name. The topic name indicates to Zemax® what data is being requested.

*dE dcj*

the Zernike Polynomial have been taken into account in this application.

from the server, and the server sends the data back to the client.

*<sup>E</sup>* <sup>=</sup> <sup>Σ</sup>(*δ<sup>i</sup>* <sup>−</sup> *Di*)<sup>2</sup> (7)

= 2Σ(*δ<sup>i</sup>* − Σ*cj fij*)*fij* = 0 (8)

**Figure 3.** Zernike polynomial basis

error *E* respect to *cj*:

**2.2. Raytracing**

*2.2.1. Zemax® and Matlab® data exchange*

further analysis or computation[17].

In the case of thermal and gravitational performances evaluation, it is necessary an adaptation of the optical model. Coordinate Breaks (CB) are introduced before and after each Optical element paying attention to the transformation sequence. This is crucial to realize a feature that allow the introduction of local optical displacements without modifying the remaining optical layout. Matlab® extracts optical displacements from Nastran (global coordinate system), optical coordinate systems from Zemax® and applies the proper transformation matrices to obtain optical displacements in Zemax® local reference systems. In case of surface deformations, the surfaces of Zemax® file are modified from "Standard" to "Zernike Fringe Sag" in order to allow Matlab® to updates the extra data tables with the Zernike coefficients obtained from the fitting algorithm.

After the uploading of mechanical data into the optical file, Zemax® evaluate spot and p.s.f. information through its raytracing engine. If necessary Matlab® request a focusing optimization after having stored the focal plane distance to evaluate the relative focal variation. Finally Matlab® extract from Zemax Spot dimensions (Max and R.M.S. radius) and centroid displacement, eventually of multiple fields.

The results defines respectively the image quality and stability of the overall optical layout and, if necessary can be passed to an optimization algorithm that manage the mechanical properties of the system.

### **2.3. Optimization**

Whereas optimization methods are nearly as old as calculus, numerical optimization reached prominence in the digital age. Its systematic application to structural design dates to its advocacy by Schmit in 1960. The success of structural optimization in the 1970s motivated the emergence of multidisciplinary design optimization (MDO) in the 1980s. Here will be presented a simplified single variable optimization approach that has been used as starting point[14].

#### *2.3.1. Problem statement*

The general optimization problem (for example minimization) wants to minimize the objective function:

$$F(X) \tag{9}$$

subjected to inequality constraints:

$$p\_j(X) \le 0 \; j = 1:m \tag{10}$$

259

#### 8 Will-be-set-by-IN-TECH 260 Finite Element Analysis – New Trends and Developments Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization <sup>9</sup>

and equality constraints:

$$h\_k(X) \le 0 \; | \; k = 1 : l \; \tag{11}$$

**3.1. Shape memory alloy Finite Element modeling: the Turner micro-mechanical**

Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization

261

To perform Finite Element analysis regarding SMA is necessary to define the proper constitutive law that describe the material behavior. In this paper will be presented the FE implementation of three different constitutive law. A further detailed description of the laws

In physics, a constitutive equation is a relation between two physical quantities (often described by tensors) that is specific to a material or substance, and approximates the response of that material to external forces. It is combined with other equations governing physical laws to solve physical problems, like the flow of a fluid in a pipe, or the response of a crystal to an

During the last ten years researchers have been developed several constitutive laws that can be classified considering the different approaches used in their formulation such as:

The micro-mechanical formulation essentially takes into account the properties of single crystals of the material averaging their behavior over a Representative Volume Element (EVE). Micro-mechanical models have been developed using a thermodynamic approach and evaluating the energy involved during a phase transformation. These models also use

The real benefit of this class of model lies in their ability to predict the real response of the material starting from the lattice parameters for crystalline and crystal and grain level data derived from martensitic transformation. Nevertheless these models are very complex and

The Turner model has been successfully implemented into the commercial code ABAQUS® [10]. This model defines temperature dependent Young modulus and Thermal expansion coefficient. Exploiting some peculiar properties of the software it is possible to define a look up table mapping the variation of material characteristics. As obtained from the calibration of the model (Figures 4 and 5), it has been used the definition of Young modulus and Thermal

Different verification models were performed to select the type of the elements [10]. The comparison underlines how a shell model is precise enough with a high gain in terms of CPU time. Due to this consideration, the final comparison was conducted on a shell-based model

micro-mechanical or macro-mechanical and phenomenological or thermodynamic.

homogenization techniques to derive the overall behavior of the SMA.

require a large number of computational operations.

Expansion coefficient as a function of temperature.

and their application can be found in the references [7].

**approach**

electric field.

**Figure 4.** SMA *E*(*T*)

and side constraints (if applicable):

$$X\_i^l \le X\_i \le X\_i^u \ i = 1 : n \tag{12}$$

where vector *X* is the *n* dimension design variable vector.

In opto-mechanical variable design *X* may consists into the instrument properties (shell thickness, beam diameter, ...) while *F*(*X*) is the overall instrument Mass function dependent to design variables. The inequality constraint *p*(*X*) ≤ 0 is the optical displacement that must be kept within certain specifications.

#### *2.3.2. Iterative procedure*

The optimization algorithm requires an initial set of variables *X*0. In the opto-mechanical optimization we consider as staring point the preliminary mesh defined by the user based onto its own experience. The whole integrated analysis is carried out as shown in this chapter extracting from MSC-Nastran® the overall system mass and from Zemax® the first optical displacement set. A small perturbation is applied onto the design variable (*X*<sup>1</sup> = *X*<sup>0</sup> + *δX*) and the procedure evaluates a second set of Mass and optical displacements. After the initialization the gradient of Mass and Optical displacement functions are evaluate respect to the design variable and the automatic procedure is started.

The optimization routine is launched:


## **3. Modeling techniques**

In the attitude of this book, following will be presented the modeling technique developed for the build up of the Finite Element Models of the Smart Structures. In particular the techniques adopted for the modeling of PiezoComposites and Shape Memory Alloy will be briefly introduced. It will not be described here the phsyical behaviour of the materials, detialed description of the full work and approach can be found in [11] and in [13].

## **3.1. Shape memory alloy Finite Element modeling: the Turner micro-mechanical approach**

To perform Finite Element analysis regarding SMA is necessary to define the proper constitutive law that describe the material behavior. In this paper will be presented the FE implementation of three different constitutive law. A further detailed description of the laws and their application can be found in the references [7].

In physics, a constitutive equation is a relation between two physical quantities (often described by tensors) that is specific to a material or substance, and approximates the response of that material to external forces. It is combined with other equations governing physical laws to solve physical problems, like the flow of a fluid in a pipe, or the response of a crystal to an electric field.

During the last ten years researchers have been developed several constitutive laws that can be classified considering the different approaches used in their formulation such as: micro-mechanical or macro-mechanical and phenomenological or thermodynamic.

The micro-mechanical formulation essentially takes into account the properties of single crystals of the material averaging their behavior over a Representative Volume Element (EVE). Micro-mechanical models have been developed using a thermodynamic approach and evaluating the energy involved during a phase transformation. These models also use homogenization techniques to derive the overall behavior of the SMA.

The real benefit of this class of model lies in their ability to predict the real response of the material starting from the lattice parameters for crystalline and crystal and grain level data derived from martensitic transformation. Nevertheless these models are very complex and require a large number of computational operations.

The Turner model has been successfully implemented into the commercial code ABAQUS® [10]. This model defines temperature dependent Young modulus and Thermal expansion coefficient. Exploiting some peculiar properties of the software it is possible to define a look up table mapping the variation of material characteristics. As obtained from the calibration of the model (Figures 4 and 5), it has been used the definition of Young modulus and Thermal Expansion coefficient as a function of temperature.

#### **Figure 4.** SMA *E*(*T*)

8 Will-be-set-by-IN-TECH

In opto-mechanical variable design *X* may consists into the instrument properties (shell thickness, beam diameter, ...) while *F*(*X*) is the overall instrument Mass function dependent to design variables. The inequality constraint *p*(*X*) ≤ 0 is the optical displacement that must

The optimization algorithm requires an initial set of variables *X*0. In the opto-mechanical optimization we consider as staring point the preliminary mesh defined by the user based onto its own experience. The whole integrated analysis is carried out as shown in this chapter extracting from MSC-Nastran® the overall system mass and from Zemax® the first optical displacement set. A small perturbation is applied onto the design variable (*X*<sup>1</sup> = *X*<sup>0</sup> + *δX*) and the procedure evaluates a second set of Mass and optical displacements. After the initialization the gradient of Mass and Optical displacement functions are evaluate respect

1. Evaluate Mass function gradient as ∇*M*(*Xi*−1)=(*M*(*Xi*−1) − *<sup>M</sup>*(*Xi*−2))/(*Xi*−<sup>1</sup> − *Xi*−2); 2. Update the design variable: *Xi* = *Xi*−<sup>1</sup> − *gM*∇*M*(*Xi*−1) where *gM* is a gain factor properly

3. Evaluate Mass and optical displacement with new variables and update step counter; 4. Is optical displacement below specification? Yes: go to point 5; No: go to point 6

5. Is *<sup>M</sup>*(*Xi*) − *<sup>M</sup>*(*Xi*−1) ≤ *Toll* (where *Toll* is a reference value) i.e. the minimization

6. Evaluate Displacement function gradient as ∇*D*(*Xi*−1)=(*D*(*Xi*−1) − *<sup>D</sup>*(*Xi*−2))/(*Xi*−<sup>1</sup> −

7. Update the design variable: *Xi* = *Xi*−<sup>1</sup> − *gD*∇*D*(*Xi*−1) where *gD* is a gain factor properly

In the attitude of this book, following will be presented the modeling technique developed for the build up of the Finite Element Models of the Smart Structures. In particular the techniques adopted for the modeling of PiezoComposites and Shape Memory Alloy will be briefly introduced. It will not be described here the phsyical behaviour of the materials,

detialed description of the full work and approach can be found in [11] and in [13].

*Xl*

where vector *X* is the *n* dimension design variable vector.

to the design variable and the automatic procedure is started.

converging? Yes: optimization ended; No: go to point 1

set in order to manage iteration and go to point 3

*<sup>i</sup>* <sup>≤</sup> *Xi* <sup>≤</sup> *<sup>X</sup><sup>u</sup>*

*hk*(*X*) ≤ 0 *k* = 1 : *l* (11)

*<sup>i</sup> i* = 1 : *n* (12)

and equality constraints:

and side constraints (if applicable):

be kept within certain specifications.

The optimization routine is launched:

set in order to manage iteration number;

*2.3.2. Iterative procedure*

*Xi*−2);

**3. Modeling techniques**

Different verification models were performed to select the type of the elements [10]. The comparison underlines how a shell model is precise enough with a high gain in terms of CPU time. Due to this consideration, the final comparison was conducted on a shell-based model

and activated via Joule effect. The numerical/experimental comparison shows that, after a short transition, the predicted displacement matches (max. error 20 *μm*) the experimental one

Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization

263

**3.2. Piezoelectric numerical modeling:"Structural-scale" modeling technique for**

Following will be presented the numerical technique adopted to model The Macro Fiber Composites (MFC) that have been developed by NASA Langley Research Center (LARC) in 2000 [16]. The components of the MFC are illustrated in Figure 9. The core is made by aligned piezoceramic fiber included in epoxy resin and joined between two groups of interdigited electrodes (IDE) supported by a Kapton film. The MFC have an overall thickness of 0.3*mm* and the dimensions of the region in which are placed PZT fibers called the active area, can vary from an area of 28 <sup>×</sup> <sup>14</sup>*mm*<sup>2</sup> to 85 <sup>×</sup> <sup>28</sup>*mm*2. The spacing of IDE is 0.5*mm*, while the fibers

The set up of numerical procedure for the performance prediction as the selection of the best technology to couple the structure to the actuators are fundamental to design efficient smart structures. The numerical approach for the design of piezo based smart structures is dealt in this section . The design of smart structure is mainly oriented to the study of authority and then to the stress analysis. With the intent of reducing the design time, a technique able to predict the overall performances of the structure through "light-weight" numerical models

The overall smart structure can be conceptually divided into two sub-system: the structure itself consisting into the composite panel and the MFC actuators. The proposed technique is developed under the ABAQUS® commercial code [4] environment1 and neglect the detailed analysis of stress state adopting essentially two type of elements: S4R four-node shell elements for the host material and C3D20RE twenty-node quadratic bricks (reduced

have a width of 350*μm* and a volume fraction above 85%.

(Figure 7).

**Figure 8.** MFC actuator

**Figure 9.** MFC layers

was developed.

<sup>1</sup> And can be obviously ported to other Finite Element Codes

**MFC**

**Figure 7.** Displacement Turner correlation results

with composite laminate properties. The NiTiNOL actuators were modeled thickening the mesh and considering the material as a "ply" of the lamination sequence. The thermoelastic model receives as input the thermal load on the wires and derives the temperature behavior of all the nodes (Figure 6).

A carbon fiber-reinforced panel with six embedded NiTiNOL wires was modeled and analyzed using the ABAQUS commercial code. The actual panel was manufactured with 12 plies [90/(0)2/90/ + 45/ − 45]*<sup>s</sup>* and its overall dimensions were: 30 x 170 x 1.2 mm. The actuators chosen were OWSME wires (*φ* = 0.38*mm*) trained using standard heat treatments. They were embedded between the 2*nd* and the and 3*rd* plies with a 4% imposed strain. The manufacturing technolgoy is reported in [3, 10]. The panel was constrained on one side and activated via Joule effect. The numerical/experimental comparison shows that, after a short transition, the predicted displacement matches (max. error 20 *μm*) the experimental one (Figure 7).

## **3.2. Piezoelectric numerical modeling:"Structural-scale" modeling technique for MFC**

Following will be presented the numerical technique adopted to model The Macro Fiber Composites (MFC) that have been developed by NASA Langley Research Center (LARC) in 2000 [16]. The components of the MFC are illustrated in Figure 9. The core is made by aligned piezoceramic fiber included in epoxy resin and joined between two groups of interdigited electrodes (IDE) supported by a Kapton film. The MFC have an overall thickness of 0.3*mm* and the dimensions of the region in which are placed PZT fibers called the active area, can vary from an area of 28 <sup>×</sup> <sup>14</sup>*mm*<sup>2</sup> to 85 <sup>×</sup> <sup>28</sup>*mm*2. The spacing of IDE is 0.5*mm*, while the fibers have a width of 350*μm* and a volume fraction above 85%.

**Figure 8.** MFC actuator

10 Will-be-set-by-IN-TECH

with composite laminate properties. The NiTiNOL actuators were modeled thickening the mesh and considering the material as a "ply" of the lamination sequence. The thermoelastic model receives as input the thermal load on the wires and derives the temperature behavior

A carbon fiber-reinforced panel with six embedded NiTiNOL wires was modeled and analyzed using the ABAQUS commercial code. The actual panel was manufactured with 12 plies [90/(0)2/90/ + 45/ − 45]*<sup>s</sup>* and its overall dimensions were: 30 x 170 x 1.2 mm. The actuators chosen were OWSME wires (*φ* = 0.38*mm*) trained using standard heat treatments. They were embedded between the 2*nd* and the and 3*rd* plies with a 4% imposed strain. The manufacturing technolgoy is reported in [3, 10]. The panel was constrained on one side

**Figure 5.** SMA *α*(*T*)

**Figure 6.** Turner Finite Element model

**Figure 7.** Displacement Turner correlation results

of all the nodes (Figure 6).

**Figure 9.** MFC layers

The set up of numerical procedure for the performance prediction as the selection of the best technology to couple the structure to the actuators are fundamental to design efficient smart structures. The numerical approach for the design of piezo based smart structures is dealt in this section . The design of smart structure is mainly oriented to the study of authority and then to the stress analysis. With the intent of reducing the design time, a technique able to predict the overall performances of the structure through "light-weight" numerical models was developed.

The overall smart structure can be conceptually divided into two sub-system: the structure itself consisting into the composite panel and the MFC actuators. The proposed technique is developed under the ABAQUS® commercial code [4] environment1 and neglect the detailed analysis of stress state adopting essentially two type of elements: S4R four-node shell elements for the host material and C3D20RE twenty-node quadratic bricks (reduced

<sup>1</sup> And can be obviously ported to other Finite Element Codes

#### 12 Will-be-set-by-IN-TECH 264 Finite Element Analysis – New Trends and Developments Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization <sup>13</sup>

integration2) for the actuator. The two subsystems are modeled as separated mesh connected together through the \*TIE algorithm which introduce a link between the nodes of the two different surfaces evaluating the distances between the two faces and obtaining an adhesion factor to be applied at each node (Figure 10).

**Figure 13.** Finite Element results of the test specimen

modeled as a voltage equivalent activation4.

**Figure 14.** MFC*EQ* into its MFC counterpart

**4.1. Shape memory actuated deformable mirrors**

<sup>3</sup> The whole substructure was modeled in detail with 0.2*mm* mesh pitch

errors (1*e* − 2%).

**4. Applications**

light-weighted technique are exceptionally close to the detailed ones.

detailed MFC as a sequence of different layers and subcomponents. A previous comparison [4] between detailed models3 confirm in fact, that the displacements resulting from the

Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization

265

The MFC have then been modeled through an homogeneous equivalent piezoelectric plates MFC*EQ* with the same effective coefficients derived from the data-sheets. The dimensions that have been taken into account are the ones of the active area, while the electrodes have been

The MFC*EQ* has been then used for the modeling of simple smart panels with bonded or embedded actuators. The simply bonding of the piezo device onto a panel is simulated through the TIE algorithm and have been compared with the detailed one obtaining limited

This procedure has been used to evaluate optical performances deformable mirrors[7–9] based onto new technologies known as smart structures. Carbon fiber reinforced Mirrors actuated by Shape Memory Alloy wires have been modeled to verify their optical capabilities. Image

<sup>4</sup> the electric field between two digit times the number of digit gives the overall electric field and tension.

#### **Figure 11.** Layout for the structural scale validation.

A simple laminate (Figure 11) 200 × 40*mm* made by three layers, with the piezo in the middle, has been considered to validate the proposed technique. A d.d.p. of 100*V* between the electrodes has been simulated and the transversal motion of the tip has been evaluated through the Classic Lamination Theory (CLT) [4]. Since the laminate is formed by three plates the TIE algorithm is used twice. Due to the fact that Piezo element are quadratic, while shell ones are linear different mesh densities must be adopted as shown in Figure 12

**Figure 12.** Different mesh densities

It is important to pay attention to the constraint conditions, because those nodes are also influenced by the TIE algorithm which connect elements with six degrees of freedom (d.o.f.) per nodes to other with three d.o.f. per nodes. If compared with the CLT, the results coming from those analyses, strongly encourages the use of these techniques in fact the tip deflection error obtained is less then 1*e* − 3%.

The proposed technique can be also adopted for the modeling of smart structures with MFC actuators. For an overall performance evaluation is not necessary to model the whole

<sup>2</sup> the reduced integration is necessary to reduce the so called shear-locking numerical induced effect

**Figure 13.** Finite Element results of the test specimen

detailed MFC as a sequence of different layers and subcomponents. A previous comparison [4] between detailed models3 confirm in fact, that the displacements resulting from the light-weighted technique are exceptionally close to the detailed ones.

The MFC have then been modeled through an homogeneous equivalent piezoelectric plates MFC*EQ* with the same effective coefficients derived from the data-sheets. The dimensions that have been taken into account are the ones of the active area, while the electrodes have been modeled as a voltage equivalent activation4.

**Figure 14.** MFC*EQ* into its MFC counterpart

The MFC*EQ* has been then used for the modeling of simple smart panels with bonded or embedded actuators. The simply bonding of the piezo device onto a panel is simulated through the TIE algorithm and have been compared with the detailed one obtaining limited errors (1*e* − 2%).

## **4. Applications**

12 Will-be-set-by-IN-TECH

integration2) for the actuator. The two subsystems are modeled as separated mesh connected together through the \*TIE algorithm which introduce a link between the nodes of the two different surfaces evaluating the distances between the two faces and obtaining an adhesion

A simple laminate (Figure 11) 200 × 40*mm* made by three layers, with the piezo in the middle, has been considered to validate the proposed technique. A d.d.p. of 100*V* between the electrodes has been simulated and the transversal motion of the tip has been evaluated through the Classic Lamination Theory (CLT) [4]. Since the laminate is formed by three plates the TIE algorithm is used twice. Due to the fact that Piezo element are quadratic, while shell

It is important to pay attention to the constraint conditions, because those nodes are also influenced by the TIE algorithm which connect elements with six degrees of freedom (d.o.f.) per nodes to other with three d.o.f. per nodes. If compared with the CLT, the results coming from those analyses, strongly encourages the use of these techniques in fact the tip deflection

The proposed technique can be also adopted for the modeling of smart structures with MFC actuators. For an overall performance evaluation is not necessary to model the whole

<sup>2</sup> the reduced integration is necessary to reduce the so called shear-locking numerical induced effect

ones are linear different mesh densities must be adopted as shown in Figure 12

factor to be applied at each node (Figure 10).

**Figure 10.** Sketch of the proposed structural scale technique.

**Figure 11.** Layout for the structural scale validation.

**Figure 12.** Different mesh densities

error obtained is less then 1*e* − 3%.

#### **4.1. Shape memory actuated deformable mirrors**

This procedure has been used to evaluate optical performances deformable mirrors[7–9] based onto new technologies known as smart structures. Carbon fiber reinforced Mirrors actuated by Shape Memory Alloy wires have been modeled to verify their optical capabilities. Image

<sup>3</sup> The whole substructure was modeled in detail with 0.2*mm* mesh pitch

<sup>4</sup> the electric field between two digit times the number of digit gives the overall electric field and tension.

#### 14 Will-be-set-by-IN-TECH 266 Finite Element Analysis – New Trends and Developments Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization <sup>15</sup>

quality during refocusing oriented activation has been evaluated. Several runs of this code has been performed, varying the number/position of actuators and the structure lamination sequence in order to have the evaluation of the influence of each ingredient onto the efficiency of the overall system. In this particular application the routine interact with ABAQUS® that is a finite element solver more reliable with user defined Constitutive laws. The aim of the analysis was the evaluation of the authority of the actuators into the variation of the curvature radius of the shell keeping the image quality within acceptable limits.

**Figure 17.** Example of variation of focus position during the activation of a SMA deformable mirror.

Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization

267

**Figure 18.** Starting p.s.f.

**Figure 19.** Starting focused spot

**Figure 20.** Example of Deformed p.s.f.

**Figure 15.** Undeformed SMA actuated Deformable mirror

**Figure 16.** Deformed SMA actuated Deformable mirror

In Figure 17 is shown the overall activation sequence: the first part of the process is dominated by a contraction of the deformable mirror due to the mismatching of the CTE of NiTiNOL respect to CFRP; when the temperature raise up to the *AS* the phase transition starts and the actuators imposes the recovery strain deploying the deformable mirror. The image quality has been evaluated both in terms of focusing spot (Figures 19 and 21) and p.s.f. (Figures 18 and 20)

Several analyses have been performed to evaluare the focusing capabilities respect to the number of actuators and the stiffness of the substrate (i.e. number of composite plies). For the detail and the results of this analyses that are not purpose of this chapter we refer to [8]. As an example here we show how the actuators density modify the efficiency of the smart SMA actuated structure. Thus why thw he comparison between two different configurations that offer similar variations of focal positions is shown. In Figure 22 is displayed the overall focusing variation comparing a 48 plies substrate with 14 actuators (blue) and a 36 plies with 8 actuators (green). This behavior has been crosschecked with the RMS and max spot size variation (respectively Figures 23 and 24).

266 Finite Element Analysis – New Trends and Developments Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization <sup>15</sup> 267 Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization

**Figure 17.** Example of variation of focus position during the activation of a SMA deformable mirror.

14 Will-be-set-by-IN-TECH

quality during refocusing oriented activation has been evaluated. Several runs of this code has been performed, varying the number/position of actuators and the structure lamination sequence in order to have the evaluation of the influence of each ingredient onto the efficiency of the overall system. In this particular application the routine interact with ABAQUS® that is a finite element solver more reliable with user defined Constitutive laws. The aim of the analysis was the evaluation of the authority of the actuators into the variation of the curvature

In Figure 17 is shown the overall activation sequence: the first part of the process is dominated by a contraction of the deformable mirror due to the mismatching of the CTE of NiTiNOL respect to CFRP; when the temperature raise up to the *AS* the phase transition starts and the actuators imposes the recovery strain deploying the deformable mirror. The image quality has been evaluated both in terms of focusing spot (Figures 19 and 21) and p.s.f. (Figures 18

Several analyses have been performed to evaluare the focusing capabilities respect to the number of actuators and the stiffness of the substrate (i.e. number of composite plies). For the detail and the results of this analyses that are not purpose of this chapter we refer to [8]. As an example here we show how the actuators density modify the efficiency of the smart SMA actuated structure. Thus why thw he comparison between two different configurations that offer similar variations of focal positions is shown. In Figure 22 is displayed the overall focusing variation comparing a 48 plies substrate with 14 actuators (blue) and a 36 plies with 8 actuators (green). This behavior has been crosschecked with the RMS and max spot size

radius of the shell keeping the image quality within acceptable limits.

**Figure 15.** Undeformed SMA actuated Deformable mirror

**Figure 16.** Deformed SMA actuated Deformable mirror

variation (respectively Figures 23 and 24).

and 20)

**Figure 19.** Starting focused spot

**Figure 20.** Example of Deformed p.s.f.

#### 16 Will-be-set-by-IN-TECH 268 Finite Element Analysis – New Trends and Developments Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization <sup>17</sup>

**Figure 21.** Example Deformed focused spot

**Figure 24.** MAx spot radius with similar focusing range

**Figure 25.** Undeformed Piezo actuated Deformable mirror

**Figure 26.** Deformed Piezo actuated Deformable mirror

they were already available in the laboratory,

the different configurations and the lamination sequence.

and the curvature radius variation have been verified as functions of the number of actuators,

269

Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization

It has been considered the smallest type of actuator even if the analysis framework set up allow the implementation of several type and dimension of them, simply due to the fact that

Considering this class of actuators the performances are not satisfactory when the MFC are placed onto a single concentric ring. In Figures 27, 28 and 29 the results obtained with a population of 15 actuators evaluated as a technological limit are plotted; the radial

**Figure 22.** Focal position variation with similar focusing range

**Figure 23.** RMS spot radius with similar focusing range

#### **4.2. MFC actuated deformable mirrors**

The same approach (paragraph 4.1) has been used to evaluate Carbon fiber reinforced Mirrors actuated by Piezoelectric MFC actuators.

Analyses have been conducted based on the correlation data obtained by the validation of simplified modeling technique presented in Paragraph 3.2. In particular the surface quality

**Figure 24.** MAx spot radius with similar focusing range

16 Will-be-set-by-IN-TECH

**Figure 21.** Example Deformed focused spot

**Figure 22.** Focal position variation with similar focusing range

**Figure 23.** RMS spot radius with similar focusing range

The same approach (paragraph 4.1) has been used to evaluate Carbon fiber reinforced Mirrors

Analyses have been conducted based on the correlation data obtained by the validation of simplified modeling technique presented in Paragraph 3.2. In particular the surface quality

**4.2. MFC actuated deformable mirrors**

actuated by Piezoelectric MFC actuators.

**Figure 25.** Undeformed Piezo actuated Deformable mirror

**Figure 26.** Deformed Piezo actuated Deformable mirror

and the curvature radius variation have been verified as functions of the number of actuators, the different configurations and the lamination sequence.

It has been considered the smallest type of actuator even if the analysis framework set up allow the implementation of several type and dimension of them, simply due to the fact that they were already available in the laboratory,

Considering this class of actuators the performances are not satisfactory when the MFC are placed onto a single concentric ring. In Figures 27, 28 and 29 the results obtained with a population of 15 actuators evaluated as a technological limit are plotted; the radial

#### 18 Will-be-set-by-IN-TECH 270 Finite Element Analysis – New Trends and Developments Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization <sup>19</sup>

coordinates have been parametrically obtained searching for the best performances that have been obtained placing the actuators at 0.53% of the mirror radius. With this layout the spot size goes out of boundary limit with low radius variations (0.37*mm*) while the max spot radius is always unacceptable.

processed, with the higher number of actuator that are then limited only by technological

Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization

271

**Figure 30.** Focal position variation with 10 and 8 actuators per ring (two ring 60 plies)

**Figure 31.** RMS spot radius with 10 and 8 actuators per ring (two ring 60 plies)

**Figure 32.** Max spot radius with 10 and 8 actuators per ring (two ring 60 plies

analyses that are not purpose of this chapter we refer to [8].

Also in this case have been performed analyses to evaluare the focusing capabilities respect to the number of actuators and the stiffness of the substrate. For the detail and the results of this

aspects, mainly related to the crowding of the surface.

**Figure 27.** Focal position variation with 10 actuators along one ring (60 plies)

**Figure 28.** RMS spot radius with 10 actuators along one ring (60 plies)

**Figure 29.** Max spot radius with 10 actuators along one ring (60 plies)

The performance of the piezo actuated smart deformable mirror are more interesting if we consider two ring of actuators, with the same angular coordinates (Figure 25). In Figure 26 is shown the contour of deformation. The results obtained comparing configurations with respectively 8 and 10 batteries of actuators (Figures 30, 31 and 32) shows the more interesting performances, both in terms of higher radius variation and in terms of better image quality processed, with the higher number of actuator that are then limited only by technological aspects, mainly related to the crowding of the surface.

**Figure 30.** Focal position variation with 10 and 8 actuators per ring (two ring 60 plies)

18 Will-be-set-by-IN-TECH

coordinates have been parametrically obtained searching for the best performances that have been obtained placing the actuators at 0.53% of the mirror radius. With this layout the spot size goes out of boundary limit with low radius variations (0.37*mm*) while the max spot radius

**Figure 27.** Focal position variation with 10 actuators along one ring (60 plies)

**Figure 28.** RMS spot radius with 10 actuators along one ring (60 plies)

**Figure 29.** Max spot radius with 10 actuators along one ring (60 plies)

The performance of the piezo actuated smart deformable mirror are more interesting if we consider two ring of actuators, with the same angular coordinates (Figure 25). In Figure 26 is shown the contour of deformation. The results obtained comparing configurations with respectively 8 and 10 batteries of actuators (Figures 30, 31 and 32) shows the more interesting performances, both in terms of higher radius variation and in terms of better image quality

is always unacceptable.

**Figure 31.** RMS spot radius with 10 and 8 actuators per ring (two ring 60 plies)

**Figure 32.** Max spot radius with 10 and 8 actuators per ring (two ring 60 plies

Also in this case have been performed analyses to evaluare the focusing capabilities respect to the number of actuators and the stiffness of the substrate. For the detail and the results of this analyses that are not purpose of this chapter we refer to [8].

## **4.3. Opto-mechanical design**

Here we show a simple example of the optimization procedure. Within the framework of a feasibility study for a robotic 3m class telescope, equipped with VIS and NIR instruments, able to react to a satellite trigger in less than 50 sec (with a goal of 30 sec): CODEVISIR (Conceptual Design for a VISible and nIR telescope), we exploited the optimization procedure to define the thickness of the main bench of the instrumental suite[15].

(Δ*th* = 0.016*mm*)) within 10 steps. The final thickness is around 25*mm* and the optimized weight of the instrument is quite less than 1.2*Ton* and the final displacement is 2.88*μm*. This procedure produce an overall weight saving of 270*Kg*, starting from a very conservative value.

Integrated FEA and Raytracing in Instrumentation for Astronomy: Smart Structures Evaluation and Structural Optimization

273

The aim of this chapter was to present the design procedure developed to the numerical capabilities of modern softwares and hardwares. The framework is very versatile and can be used simply for performance prediction of optical system or structres and/or for optimization strategies. The Matlab® code works as a "server" that interact with the "client" Finite Element (ABAQUS, Nastran, ...) solver managing the input file and extracting the results; the data are modified in order to be handled by the "client" raytracing software that evaluates the optical performances. The procedure is under development and integration in order to allow a multi-variable optimization. We are planning also to evaluates different optimization algorithms exploiting in particular the sensitivities analyses[5] capabilities of both Finite

**Figure 35.** Image Displacement optimization

**Figure 36.** Thickness and Mass Variation

Element and raytracing softwares.

**5. Conclusions**

The instrument includes seven camera able to cover the wavelength range from the Visible (VIS, 0.4 − 0.9*μm*) to the Near Infrared (NIR, 1 − 2.5*μm*) during the same exposure. This will be allowed by a multichannel imaging configuration that envisages a detector for each photometric band, delivered through a dichroic cascade along the optical path. Two ancillary instruments will complete the instrumentation suite, a visible spectrograph and a photo-polarimeter, fed by rotating the M3 mirror.

**Figure 33.** Codevisir Instrument suite

**Figure 34.** Finite Element simplified Model

The Optimization procedure has been oriented to weight minimizing of the overall system keeping the optical displacement5 below proper requirements. The design variable considered is the thickness of the bench. In Figure 35 can be observed the image displacement behavior calculated as *<sup>D</sup>* <sup>=</sup> <sup>√</sup>Δ*X*<sup>2</sup> <sup>+</sup> <sup>Δ</sup>*Y*2. The requirement was set at one fifth of a pixel (i.e. 3*μm*) and the optimization variable has been modified starting from 35*mm*. The optimization routine converges to a reasonable value (the thickness variation is less than a reasonable value

<sup>5</sup> of one of the seven arms that is also representative of the performances of all the others

(Δ*th* = 0.016*mm*)) within 10 steps. The final thickness is around 25*mm* and the optimized weight of the instrument is quite less than 1.2*Ton* and the final displacement is 2.88*μm*. This procedure produce an overall weight saving of 270*Kg*, starting from a very conservative value.

**Figure 35.** Image Displacement optimization

**Figure 36.** Thickness and Mass Variation

## **5. Conclusions**

20 Will-be-set-by-IN-TECH

Here we show a simple example of the optimization procedure. Within the framework of a feasibility study for a robotic 3m class telescope, equipped with VIS and NIR instruments, able to react to a satellite trigger in less than 50 sec (with a goal of 30 sec): CODEVISIR (Conceptual Design for a VISible and nIR telescope), we exploited the optimization procedure to define the

The instrument includes seven camera able to cover the wavelength range from the Visible (VIS, 0.4 − 0.9*μm*) to the Near Infrared (NIR, 1 − 2.5*μm*) during the same exposure. This will be allowed by a multichannel imaging configuration that envisages a detector for each photometric band, delivered through a dichroic cascade along the optical path. Two ancillary instruments will complete the instrumentation suite, a visible spectrograph and a

The Optimization procedure has been oriented to weight minimizing of the overall system keeping the optical displacement5 below proper requirements. The design variable considered is the thickness of the bench. In Figure 35 can be observed the image displacement behavior calculated as *<sup>D</sup>* <sup>=</sup> <sup>√</sup>Δ*X*<sup>2</sup> <sup>+</sup> <sup>Δ</sup>*Y*2. The requirement was set at one fifth of a pixel (i.e. 3*μm*) and the optimization variable has been modified starting from 35*mm*. The optimization routine converges to a reasonable value (the thickness variation is less than a reasonable value

<sup>5</sup> of one of the seven arms that is also representative of the performances of all the others

**4.3. Opto-mechanical design**

**Figure 33.** Codevisir Instrument suite

**Figure 34.** Finite Element simplified Model

thickness of the main bench of the instrumental suite[15].

photo-polarimeter, fed by rotating the M3 mirror.

The aim of this chapter was to present the design procedure developed to the numerical capabilities of modern softwares and hardwares. The framework is very versatile and can be used simply for performance prediction of optical system or structres and/or for optimization strategies. The Matlab® code works as a "server" that interact with the "client" Finite Element (ABAQUS, Nastran, ...) solver managing the input file and extracting the results; the data are modified in order to be handled by the "client" raytracing software that evaluates the optical performances. The procedure is under development and integration in order to allow a multi-variable optimization. We are planning also to evaluates different optimization algorithms exploiting in particular the sensitivities analyses[5] capabilities of both Finite Element and raytracing softwares.

## **Author details**

Marco Riva

*I.N.A.F. - Osservatorio Astronomico di Brera, Italy*

## **6. References**


**Recent Advances of Finite Element Analysis** 

**in "Civil & Structural Engineering"** 


**Recent Advances of Finite Element Analysis in "Civil & Structural Engineering"** 

22 Will-be-set-by-IN-TECH

[3] Bettini, P., Riva, M., Sala, G., DiLandro, L., Airoldi, A. & Cucco, J. [2009]. Carbon fiber reinforced smart laminates with embedded sma actuators-part i: Embedding techniques and interface analysis, *Journal of Materials Engineering and Performance* 18(5-6): 664–671. [4] DiSanzo, D. [2009]. *Studio di pannelli attivati da mfc: modellazione numerica e prove*

[5] Doyle, K., Genberg, V. & Michels, G. [2002]. *Integrated Optomechanical Analysis*, SPIE

[7] Riva, M. [2010]. *Smart Structures in Instrumentation for Astronomy*, PhD thesis, Aerospace

[8] Riva, M. [2011]. *Smart Structures in Instrumentation for Astronomy*, LAP LAMBERT

[9] Riva, M., Bettini, P., DiLandro, L. & Sala, G. [2009]. Shape memory composite deformable

[10] Riva, M., Bettini, P., DiLandro, L., Sala, G. & Airoldi, A. [2009]. Carbon fiber-reinforced smart laminates with embedded sma actuators-part ii: Numerical models and empirical

[12] Riva, M., DeCaprio, V., Spanó, P. & M.Tintori [2010]. Integrated fi ˛ Anite element analysis and raytracing, oriented to structural optimization, for astronomical instrument design,

[13] Riva, M., DiSanzo, D., Airoldi, A., Sala, G. & Zerbi, F. [2010]. Smart structures for deformable mirrors actuated by piezocomposites, Vol. 7739, SPIE,

[15] Vitali, F., Chincarini, G., M.Zannoni & all. [2010]. A path to the stars: the evolution of the species in the hunting to the grbs, Vol. 7730, SPIE, pp. 77330W–1/ 77330W–14. [16] Wilkie, W., Bryant, G. & et all, J. H. [2000]. Low cost piezocomposite actuator for structural control applications, *7th Annual International Symposium on Smart Structures*

[14] Vanderplaats, G. [2007]. *Multidiscipline design optimization*, VR&D Press.

correlations, *Journal of Materials Engineering and Performance* 18(5-6): 672–678. [11] Riva, M., Bettini, P., DiLandro, L., Sala, G. & Zerbi, F. [2010]. Smart structures for deformable mirrors actuated by shape memory alloy, Vol. 7739, SPIE,

[2] Ahmad, A. [1999]. *Handbook of Optomechanical Engineering*, CRC Press.

**Author details**

**6. References**

Press.

*I.N.A.F. - Osservatorio Astronomico di Brera, Italy*

[1] *ABAQUS keyword reference manual* [2007].

[6] *MSC.Nastran: User guide* [2004].

pp. 77391M–1/77391M–20.

pp. 77391N–1/77391N–16.

Vol. 7738, SPIE, pp. 77380J–2/77380J–9.

*and Materials*, SPIE, Newport Beach.

[17] *ZEMAX: Optical Design Program User's Guide* [2005].

Academic Publishing.

Engineering, Politecnico di Milano.

*sperimentali*, Master's thesis, Politecnico di Milano.

mirrors, Vol. 7288, SPIE, pp. 72881T–1/72881T–12.

Marco Riva

**Chapter 13** 

© 2012 Psarras et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Psarras et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Damage-Tolerant Design of Stiffener Run-Outs:** 

Modern aerostructures are predominantly of semi-monocoque construction characterized by a thin skin and stiffeners. The latest generation of large passenger aircraft also use mostly carbon-fibre composite material in their primary structure and there is a trend towards the utilization of bonding of subcomponents in preference of mechanical fastening. Current design philosophy requires that certain stiffeners are terminated, for example due to an intersecting structural feature or an inspection cut-out. In these circumstances, the loading in the stiffener must be diffused into the skin, leading to complex three-dimensional stressstates. The development and utilization of reliable virtual component testing, in the design of composite aerostructures, can potentially yield significant cost reductions. Such reliability requires a thorough understanding of the damage mechanisms and failure processes in

realistic aerostructures, particularly in critical regions such as stiffener run-outs.

When a stiffener is terminated, the loads which it carries must be transferred to the skin, making the design of the run-out region vital, hence improved design methodologies are required. Several studies of skin-stiffener failure [1-8] have been carried out. Falzon et al. [7, 8] investigated the failure of realistic stiffener run-outs loaded in uniaxial compression. Different stacking sequences and skin thicknesses at the run-out region were tested and a wealth of complexity in the response and subsequent failure was reported. For all tests, failure initiated at the edge of the run-out and propagated across the skin–stiffener interface. It was found that the failure load of each specimen was greatly influenced by changes in the geometric features of these specimens. Falzon and Hitchings [8] used a Virtual Crack Closure Technique (VCCT) [9] to predict the crack growth characteristics of the modelled specimens and reported shortcomings in the quantitative correlation between the predicted

**A Finite Element Approach** 

Additional information is available at the end of the chapter

S. Psarras, S.T. Pinho and B.G. Falzon

and observed failure loads and modes.

http://dx.doi.org/10.5772/50377

**1. Introduction** 

## **Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach**

S. Psarras, S.T. Pinho and B.G. Falzon

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/50377

## **1. Introduction**

Modern aerostructures are predominantly of semi-monocoque construction characterized by a thin skin and stiffeners. The latest generation of large passenger aircraft also use mostly carbon-fibre composite material in their primary structure and there is a trend towards the utilization of bonding of subcomponents in preference of mechanical fastening. Current design philosophy requires that certain stiffeners are terminated, for example due to an intersecting structural feature or an inspection cut-out. In these circumstances, the loading in the stiffener must be diffused into the skin, leading to complex three-dimensional stressstates. The development and utilization of reliable virtual component testing, in the design of composite aerostructures, can potentially yield significant cost reductions. Such reliability requires a thorough understanding of the damage mechanisms and failure processes in realistic aerostructures, particularly in critical regions such as stiffener run-outs.

When a stiffener is terminated, the loads which it carries must be transferred to the skin, making the design of the run-out region vital, hence improved design methodologies are required. Several studies of skin-stiffener failure [1-8] have been carried out. Falzon et al. [7, 8] investigated the failure of realistic stiffener run-outs loaded in uniaxial compression. Different stacking sequences and skin thicknesses at the run-out region were tested and a wealth of complexity in the response and subsequent failure was reported. For all tests, failure initiated at the edge of the run-out and propagated across the skin–stiffener interface. It was found that the failure load of each specimen was greatly influenced by changes in the geometric features of these specimens. Falzon and Hitchings [8] used a Virtual Crack Closure Technique (VCCT) [9] to predict the crack growth characteristics of the modelled specimens and reported shortcomings in the quantitative correlation between the predicted and observed failure loads and modes.

© 2012 Psarras et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Psarras et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Bisagni et al. [10] investigated, experimentally and numerically, the postbuckling response of hat-stiffened composite panels using cohesive zone models and predicted collapse loads which were in agreement with experiments. Krueger [6] used VCCT and a numerically effective shell/3D modelling approach to predict the delamination failure of skin-stiffener run-outs. Camanho et al. [11] implemented a cohesive element to numerically investigate the debond strength of skin-stiffener composite specimens using cohesive zone models and compared with experimental results.

Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach 279

30

*S*  **[MPa]** <sup>10</sup> <sup>20</sup> <sup>20</sup>

*c*

*GIc* **[kJ/m2]** *b*

*GIIc* **[kJ/m2]**

η

*d*

carbon/epoxy pre-preg, with ply thickness 0.25 mm, for the skin and the stiffener, and

50

15 2

(a) (b)

*X* **[MPa]** 

**IM7/8552** 154.1 9.8 4.48 0.34 1572.9 254.6 101.2 0.21 0.61 - **FM300** 2.38 - 0.68 - 61 - 49.8 0.9 2.5 8.0

In previous test results with geometry similar to that in Figure 1a, the specimens failed by unstable debonding of the stiffener from the skin. Therefore, the different configurations in this study were assessed by comparing the energy release rates of the run-outs for a given

All the FE simulations of the parametric study were carried out in ABAQUS [13] and the parameterized models were created in Python [14]. The main model has five different parts, the skin, the adhesive between the skin and the stiffener, the parts of the stiffener and the filler. The material properties for IM7/8552 measured using standard tests and are shown in Table 1. A mesh sensitivity study was carried out to ensure that all results presented are

All parts were modelled with three dimensional hexahedral solid elements C3D8 to accurately capture stresses in the through-thickness direction. Also, solid elements are

*Y* **[MPa]**

**Figure 1.** Designs and dimensions in mm of a) the baseline stiffener and b) the parametric stiffener.

FM300 adhesive film (0.15 mm thick) for the bondline.

1.25

2.5

15 30

10

*Eyy* **[GPa]**

displacement and for several debond lengths.

*Gxy* **[GPa]** *vxy*

**Table 1.** Material properties for IM7/8552, measured in-house, and FM300

100

Material *Exx*

**2.2. The FE model** 

mesh converged.

*2.2.1. Element Idealisation* 

**[GPa]** 

In a recent study by the authors [12], a parametric numerical analysis was conducted to optimize the design of the run-out section to increase the crack growth stability under axial compression. Improved damage tolerance (stable crack propagation) was reported in the modified stiffener run-out design as compared to the baseline configuration. The modified design eventually failed catastrophically by interlaminar delamination, not bondline failure (debonding), which had not been considered in the numerical study. A more detailed analysis of different configurations, which accounts for delamination, was therefore undertaken in the work presented in this chapter. Building on the previous findings, the merits of a C*ompliant* termination scheme are presented.

The research in this paper focuses on stiffener run-outs loaded in compression with a selection of stiffener termination schemes. These schemes are analysed numerically (Section 2) in order to compare the influence of the design on the energy release rate for debonding and delamination. The *Baseline*, *Tapered* and C*ompliant* schemes were manufactured (Section 3) and tested to failure (Section 4). The experimental results are presented in Section 5, and compared to the predictions in Section 6.

The failure mode of the *Tapered* specimen configuration was found to be a combination of interlaminar and intralaminar failure [12]. The *Tapered* stiffeners had an unexpected delamination in the flange between the bottom-most 0o ply and above the 45o ply. Furthermore, the delamination led to an intralaminar failure in the form of a matrix crack across the 0o ply near the filler ends and continued delaminating between the filler and 0o ply, Figure 19a. For this reason, a new finite element (FE) model was created in order to investigate these modes of damage, Figure 19b and c.

## **2. Numerical analysis**

## **2.1. Skin-stiffener configurations**

The main focus of this study was a new configuration to improve the crack growth stability of a stiffener run-out. The structural performance of a baseline skin-stiffener configuration under longitudinal compression, with geometry and dimensions shown in Figure 1a, was compared to that of a modified parametric configuration shown in b. The modified configuration has a widening flange towards the termination end of the stiffener but this added material is offset by the taper of the stiffener web. This results in a stiffener design with a similar overall weight to the baseline design. For the parametric configuration, various values of *b*, *c* and *d* were analysed. The materials used in this study were IM7/8552 carbon/epoxy pre-preg, with ply thickness 0.25 mm, for the skin and the stiffener, and FM300 adhesive film (0.15 mm thick) for the bondline.

**Figure 1.** Designs and dimensions in mm of a) the baseline stiffener and b) the parametric stiffener.


**Table 1.** Material properties for IM7/8552, measured in-house, and FM300

In previous test results with geometry similar to that in Figure 1a, the specimens failed by unstable debonding of the stiffener from the skin. Therefore, the different configurations in this study were assessed by comparing the energy release rates of the run-outs for a given displacement and for several debond lengths.

### **2.2. The FE model**

278 Finite Element Analysis – New Trends and Developments

compared with experimental results.

merits of a C*ompliant* termination scheme are presented.

investigate these modes of damage, Figure 19b and c.

compared to the predictions in Section 6.

**2. Numerical analysis** 

**2.1. Skin-stiffener configurations** 

Bisagni et al. [10] investigated, experimentally and numerically, the postbuckling response of hat-stiffened composite panels using cohesive zone models and predicted collapse loads which were in agreement with experiments. Krueger [6] used VCCT and a numerically effective shell/3D modelling approach to predict the delamination failure of skin-stiffener run-outs. Camanho et al. [11] implemented a cohesive element to numerically investigate the debond strength of skin-stiffener composite specimens using cohesive zone models and

In a recent study by the authors [12], a parametric numerical analysis was conducted to optimize the design of the run-out section to increase the crack growth stability under axial compression. Improved damage tolerance (stable crack propagation) was reported in the modified stiffener run-out design as compared to the baseline configuration. The modified design eventually failed catastrophically by interlaminar delamination, not bondline failure (debonding), which had not been considered in the numerical study. A more detailed analysis of different configurations, which accounts for delamination, was therefore undertaken in the work presented in this chapter. Building on the previous findings, the

The research in this paper focuses on stiffener run-outs loaded in compression with a selection of stiffener termination schemes. These schemes are analysed numerically (Section 2) in order to compare the influence of the design on the energy release rate for debonding and delamination. The *Baseline*, *Tapered* and C*ompliant* schemes were manufactured (Section 3) and tested to failure (Section 4). The experimental results are presented in Section 5, and

The failure mode of the *Tapered* specimen configuration was found to be a combination of interlaminar and intralaminar failure [12]. The *Tapered* stiffeners had an unexpected delamination in the flange between the bottom-most 0o ply and above the 45o ply. Furthermore, the delamination led to an intralaminar failure in the form of a matrix crack across the 0o ply near the filler ends and continued delaminating between the filler and 0o ply, Figure 19a. For this reason, a new finite element (FE) model was created in order to

The main focus of this study was a new configuration to improve the crack growth stability of a stiffener run-out. The structural performance of a baseline skin-stiffener configuration under longitudinal compression, with geometry and dimensions shown in Figure 1a, was compared to that of a modified parametric configuration shown in b. The modified configuration has a widening flange towards the termination end of the stiffener but this added material is offset by the taper of the stiffener web. This results in a stiffener design with a similar overall weight to the baseline design. For the parametric configuration, various values of *b*, *c* and *d* were analysed. The materials used in this study were IM7/8552 All the FE simulations of the parametric study were carried out in ABAQUS [13] and the parameterized models were created in Python [14]. The main model has five different parts, the skin, the adhesive between the skin and the stiffener, the parts of the stiffener and the filler. The material properties for IM7/8552 measured using standard tests and are shown in Table 1. A mesh sensitivity study was carried out to ensure that all results presented are mesh converged.

#### *2.2.1. Element Idealisation*

All parts were modelled with three dimensional hexahedral solid elements C3D8 to accurately capture stresses in the through-thickness direction. Also, solid elements are

capable of modelling several layers of different materials for the analysis of laminated composites which is ideal for our numerical study. In consideration of computational cost, linear order elements were chosen and to ensure higher fidelity results whew required, the models executed using a mesh density based on a ability to define a detailed mesh convergence study. One of the major advantages of solid elements is the feature of defining several ply stack ups within a single element along the stacking direction with an improved element formulation, illustrated in the Figure 2 below.

Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach 281

hourglassing. Even though this particular problem could be partially mitigated by using the hourglass stiffness enhanced element formulation, the use of C3D8R elements in the numerical models was avoided for greater accuracy. Consequently, C3D8 elements were

The skin consists of eight plies, while the stiffener consists of five plies. The stacking sequences used in [1] for the skin and the stiffener, are shown in Table 2. Figure 4 shows the ply orientations for the skin. The thickness of each ply is 0.25mm (double thickness) and the number of integration points is set to three within each ply/element through the thickness.

The number and distribution of elements is a crucial step when using the FE method to solve structural problems. Figure 5 shows the element families that are most commonly used in stress analysis. The main difference among the element families is the geometry type that each of the family represents. As mentioned earlier, an 8-node linear brick element was

used and a course mesh of the modified design is shown in Figure 6.

**Part Stacking Sequence** Skin [45/-45/0/90]s Stiffener (per half section) [0/90/-45/45/0]

**Table 2.** Stacking sequence for the skin and the stiffener

**Figure 4.** Stacking sequence for the skin

*2.2.3. Mesh* 

used in all the models.

*2.2.2. Stacking Sequence* 

**Figure 2.** Solid element with composite layups

All the models were meshed with fully integrated solid elements (C3D8) rather using reduced integration elements (C3D8R), since the reduced integration elements were not able to capture damage variable calculations at all the nodes while delaminating, Figure 3. On further investigating this issue, it was observed that the problems arose because of

**Figure 3.** Contour plot of cohesive damage variable with Hourglass modes

hourglassing. Even though this particular problem could be partially mitigated by using the hourglass stiffness enhanced element formulation, the use of C3D8R elements in the numerical models was avoided for greater accuracy. Consequently, C3D8 elements were used in all the models.

## *2.2.2. Stacking Sequence*

280 Finite Element Analysis – New Trends and Developments

**Figure 2.** Solid element with composite layups

element formulation, illustrated in the Figure 2 below.

capable of modelling several layers of different materials for the analysis of laminated composites which is ideal for our numerical study. In consideration of computational cost, linear order elements were chosen and to ensure higher fidelity results whew required, the models executed using a mesh density based on a ability to define a detailed mesh convergence study. One of the major advantages of solid elements is the feature of defining several ply stack ups within a single element along the stacking direction with an improved

All the models were meshed with fully integrated solid elements (C3D8) rather using reduced integration elements (C3D8R), since the reduced integration elements were not able to capture damage variable calculations at all the nodes while delaminating, Figure 3. On further investigating this issue, it was observed that the problems arose because of

**Figure 3.** Contour plot of cohesive damage variable with Hourglass modes

The skin consists of eight plies, while the stiffener consists of five plies. The stacking sequences used in [1] for the skin and the stiffener, are shown in Table 2. Figure 4 shows the ply orientations for the skin. The thickness of each ply is 0.25mm (double thickness) and the number of integration points is set to three within each ply/element through the thickness.


**Table 2.** Stacking sequence for the skin and the stiffener

## *2.2.3. Mesh*

The number and distribution of elements is a crucial step when using the FE method to solve structural problems. Figure 5 shows the element families that are most commonly used in stress analysis. The main difference among the element families is the geometry type that each of the family represents. As mentioned earlier, an 8-node linear brick element was used and a course mesh of the modified design is shown in Figure 6.

Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach 283

**Figure 7.** Surfaces constraints between the skin (master surface-red) and the adhesive (slave surface-

**Figure 8.** Modified model with BC-1 on the left edge (clamped) and BC-1 on the right edge (U3=1)

Finite element analysis provides an approximate solution and it can only guarantee that equilibrium is satisfied on average over an element. This does not mean that it will satisfy equilibrium over any smaller volume compared to a complete element. As a consequence, equilibrium is enhanced when the size of the element is decreased. Moreover, in regions of stress concentrations, it is necessary to increase the accuracy of the FE solution by either using elements with higher-order shape functions (p-refinement) or by using a finer mesh of

pink)

*2.2.6. Mesh sensitivity study* 

**Figure 5.** Commonly used element families in ABAQUS [13]

**Figure 6.** Course mesh for the modified design with 3D continuum elements

## *2.2.4. Constraints*

Constraints were defined in the *Assembly module* [13] for the initial positions of *instances*. The type of constraint created for this model was a "tie constraint", which allows the user to fuse together two regions even though the meshes created for sectors may be dissimilar [8]. For the current model, sixteen constraints were created between all the parts being in contact. In Figure 7 it can be seen how the master and the slave surfaces are displayed in the model.

## *2.2.5. Boundary conditions*

In the current study, two boundary conditions were defined. The first condition, named BC-1, was a displacement/rotation type boundary condition. All displacements and rotations were set to zero. BC-1 was applied to one edge of the stiffener as shown in Figure 8. The second BC, named BC-2, was a displacement/rotation type BC as well, where all displacements and rotations were set to zero apart from U3, which was set to U3=1. U3 is the displacement in the axial direction and represents the testing machine's displacement.

**Figure 7.** Surfaces constraints between the skin (master surface-red) and the adhesive (slave surfacepink)

**Figure 8.** Modified model with BC-1 on the left edge (clamped) and BC-1 on the right edge (U3=1)

## *2.2.6. Mesh sensitivity study*

282 Finite Element Analysis – New Trends and Developments

**Figure 5.** Commonly used element families in ABAQUS [13]

**Figure 6.** Course mesh for the modified design with 3D continuum elements

Constraints were defined in the *Assembly module* [13] for the initial positions of *instances*. The type of constraint created for this model was a "tie constraint", which allows the user to fuse together two regions even though the meshes created for sectors may be dissimilar [8]. For the current model, sixteen constraints were created between all the parts being in contact. In Figure 7 it can be seen how the master and the slave surfaces are displayed in the model.

In the current study, two boundary conditions were defined. The first condition, named BC-1, was a displacement/rotation type boundary condition. All displacements and rotations were set to zero. BC-1 was applied to one edge of the stiffener as shown in Figure 8. The second BC, named BC-2, was a displacement/rotation type BC as well, where all displacements and rotations were set to zero apart from U3, which was set to U3=1. U3 is the displacement in the axial direction and represents the testing machine's displacement.

*2.2.4. Constraints* 

*2.2.5. Boundary conditions* 

Finite element analysis provides an approximate solution and it can only guarantee that equilibrium is satisfied on average over an element. This does not mean that it will satisfy equilibrium over any smaller volume compared to a complete element. As a consequence, equilibrium is enhanced when the size of the element is decreased. Moreover, in regions of stress concentrations, it is necessary to increase the accuracy of the FE solution by either using elements with higher-order shape functions (p-refinement) or by using a finer mesh of

elements (h-refinement). The goal that a designer needs to achieve is to select the best mesh density which is not prohibitively expensive to run and at the same will provide accurate and acceptable results [15].

Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach 285

a [mm] a [mm]

*b* = 1 mm

*d* = 9

*d* = 4 mm

*b* = 3 mm (Selected)

*b* = 2

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

*a* [mm]

*d* = 6.25 mm (Selected)

*a* [mm]

because of the advantage of using the Abaqus scripting interface to generate automate

The parametric script generates the stiffener run-out models with different crack lengths. The script uses a basic stiffener run-out configuration and changes the design each time according to the parameter values. When a new design is generated, the script propagates

After converging to an optimum mesh density, failure in the stiffeners was captured by consecutive models, propagating a delamination crack or a debonding crack, 1 mm at a time. The results of the parametric study are shown in Fig. 10, where the values of

(a) (b)

0.6

Parameter *d*

0.7

0.8

Normalized

*GT*

0.9

1.0

1.1

Parameter *b*

*a* [mm]

(c) (d)

Normalized

*GT*

0.6 0.7 0.8 0.9 1.0 1.1

*a* [mm]

**Figure 10.** Normalized energy release rates as a function of crack length (a) comparison between baseline stiffener design (Fig. 1a) and selected modified stiffener (Fig 1b with *b* = 3 mm, *c* = 10 mm and *d* =6.25 mm), (b) Influence of parameter *b* on G, (c) Influence of parameter *c* on G and (d) Influence of

*c* = 2 mm *c* = 6 mm

*c* = 10 mm (Selected)

the crack and the energy release rate is calculated for every step.

Baseline Stiffener

repetitive models and execute.

**2.3. Results from modelling** 

0.6

0.6

parameter *d* on G.

0.7

0.8

Normalized

*GT*

0.9

1.0

1.1

Parameter *c*

0.7

0.8

Normalized

*GT*

0.9

1.0

1.1

*2.3.1. Energy Release rate along crack* 

Selected Modified Stiffener

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

In this study, three different meshes were used (a coarse mesh, an intermediate mesh and a fine mesh). 13100 elements were used for the coarse mesh, 18700 elements for the intermediate mesh and 51000 elements for the fine mesh. Table 3 illustrates the three different meshes used, the strain energy and the running time of each model. According to the values of Table 3, the solution is converged. Since the running time for the intermediate mesh was 35 minutes and the results appeared accurate and acceptable, this specific density of elements was selected for the rest of this study.


**Table 3.** Mesh Sensitivity Study results.

Moreover, investigation was undertaken to assess the accuracy of using multiple ply orientations within a solid finite element. Abaqus® provides the option to make a partition through the ply thickness and define the material orientation. As a result, eight partitions were created for the skin and five for each stiffener using material orientations given in Table 2. The results obtained from both approaches were similar, as illustrated in Figure 9.

**Figure 9.** Strain energy with composite lay-up (blue line) and with material orientations (red line)

#### *2.2.7. The Python script*

In order to perform the parametric study, a large number of models were created and a script that generates these models was needed. This script was written in Python [14] because of the advantage of using the Abaqus scripting interface to generate automate repetitive models and execute.

The parametric script generates the stiffener run-out models with different crack lengths. The script uses a basic stiffener run-out configuration and changes the design each time according to the parameter values. When a new design is generated, the script propagates the crack and the energy release rate is calculated for every step.

#### **2.3. Results from modelling**

284 Finite Element Analysis – New Trends and Developments

of elements was selected for the rest of this study.

**(KJ)** 

and acceptable results [15].

**Mesh Strain Energy** 

**Table 3.** Mesh Sensitivity Study results.

*2.2.7. The Python script* 

elements (h-refinement). The goal that a designer needs to achieve is to select the best mesh density which is not prohibitively expensive to run and at the same will provide accurate

In this study, three different meshes were used (a coarse mesh, an intermediate mesh and a fine mesh). 13100 elements were used for the coarse mesh, 18700 elements for the intermediate mesh and 51000 elements for the fine mesh. Table 3 illustrates the three different meshes used, the strain energy and the running time of each model. According to the values of Table 3, the solution is converged. Since the running time for the intermediate mesh was 35 minutes and the results appeared accurate and acceptable, this specific density

> **Number of Elements**

**Coarse** 30150 13100 25 - **Intermediate** 30138 18700 35 0.039 **Fine** 30109 51000 90 0.096

Moreover, investigation was undertaken to assess the accuracy of using multiple ply orientations within a solid finite element. Abaqus® provides the option to make a partition through the ply thickness and define the material orientation. As a result, eight partitions were created for the skin and five for each stiffener using material orientations given in Table 2. The results obtained from both approaches were similar, as illustrated in Figure 9.

**Figure 9.** Strain energy with composite lay-up (blue line) and with material orientations (red line)

In order to perform the parametric study, a large number of models were created and a script that generates these models was needed. This script was written in Python [14]

**Running Time (minutes)** 

**Deviation (%)** 

#### *2.3.1. Energy Release rate along crack*

After converging to an optimum mesh density, failure in the stiffeners was captured by consecutive models, propagating a delamination crack or a debonding crack, 1 mm at a time. The results of the parametric study are shown in Fig. 10, where the values of

**Figure 10.** Normalized energy release rates as a function of crack length (a) comparison between baseline stiffener design (Fig. 1a) and selected modified stiffener (Fig 1b with *b* = 3 mm, *c* = 10 mm and *d* =6.25 mm), (b) Influence of parameter *b* on G, (c) Influence of parameter *c* on G and (d) Influence of parameter *d* on G.

*GT* = *G<sup>Ι</sup>* + *GΙΙ* + *GΙΙΙ*, the total energy release rate, have been normalized to the *GT* of the reference parametric stiffener for 0.5 mm of crack length. Fig. 10a compares the normalised energy release rate for the baseline and selected parametric stiffeners. The negative slope of the *G(a)* curve for the latter indicates stability of crack growth (assuming constant fracture toughness). The influence of parameters *b*, *c*, and *d* is presented in Fig. 10b, c and d. Given the objective of optimising for stability of crack growth, the configuration with *b* = 3 mm, *c* = 10 mm and *d* = 6.25 mm was selected to be carried out for the consequent stages of this study.

Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach 287

Rest of plies

The part descriptions are clearly illustrated in a front view of the skin-stiffener assembly,

Figure 18.

0 plies

Filler

Skin

**Cohesive interaction**

**Figure 13.** Illustration of imposed cohesive properties for debonding mode.

The interface between the skin and the stiffener was connected by using a surface-based cohesive layer, as seen in Figure 12. Cohesive zone modelling is generally used for the numerical simulation of interlaminar failure. Damage initiation is driven by a traction seperation law and the value of the maximum traction, to, Figure 14(a). New crack surfaces are formed when the fracture toughness Gc is equal to the area surface under the tractionseparation curve. Considering the nature of predictions made by the parametric study, the model was analysed with cohesive interaction properties of adhesive FM 300 and listed in

**Figure 12.** The parts of the FE model

Table 4.

**Figure 11.** (a) Front view of failed specimen; (b) Exploded view showing the failed area; (c) Front view of bottom part showing 00 plies ; (d) Bottom view of the upper part showing delaminated 450 plies

Interlaminar and intralaminar failures were observed in the modified run-out stiffener. The stiffener had an unexpected delamination between the 00 ply and the 450 ply. After examination of the specimens it was assumed that the delamination led to an intralaminar failure in the form of a matrix crack across the 00 ply near the filler ends and continued delaminating between the filler and 00 ply. In Figure 11, the failed specimens clearly shows that delaminations and matrix cracks occurred in the 00 ply on both sides of the specimen.

## **3. Modeling the debonding of the stiffener**

The finite element model of the stiffener run-out had the following features:


The part descriptions are clearly illustrated in a front view of the skin-stiffener assembly, Figure 18.

**Figure 12.** The parts of the FE model

286 Finite Element Analysis – New Trends and Developments

study.

*GT* = *G<sup>Ι</sup>* + *GΙΙ* + *GΙΙΙ*, the total energy release rate, have been normalized to the *GT* of the reference parametric stiffener for 0.5 mm of crack length. Fig. 10a compares the normalised energy release rate for the baseline and selected parametric stiffeners. The negative slope of the *G(a)* curve for the latter indicates stability of crack growth (assuming constant fracture toughness). The influence of parameters *b*, *c*, and *d* is presented in Fig. 10b, c and d. Given the objective of optimising for stability of crack growth, the configuration with *b* = 3 mm, *c* = 10 mm and *d* = 6.25 mm was selected to be carried out for the consequent stages of this

**Figure 11.** (a) Front view of failed specimen; (b) Exploded view showing the failed area; (c) Front view of bottom part showing 00 plies ; (d) Bottom view of the upper part showing delaminated 450 plies

Interlaminar and intralaminar failures were observed in the modified run-out stiffener. The stiffener had an unexpected delamination between the 00 ply and the 450 ply. After examination of the specimens it was assumed that the delamination led to an intralaminar failure in the form of a matrix crack across the 00 ply near the filler ends and continued delaminating between the filler and 00 ply. In Figure 11, the failed specimens clearly shows that delaminations and matrix cracks occurred in the 00 ply on both sides of the specimen.

**3. Modeling the debonding of the stiffener** 

1. Skin - single part of (450/-450/00/ 900)S laminate

4. Left - single part of (00/ 900/-450/450) laminate of the stringer

6. Right - single part of (00/ 900/-450/450) laminate of the stringer

2. Filler - made-up of 00 fibres 3. Left - 00 lamina of the stringer

5. Right - 00 lamina of the stringer

The finite element model of the stiffener run-out had the following features:

**Figure 13.** Illustration of imposed cohesive properties for debonding mode.

The interface between the skin and the stiffener was connected by using a surface-based cohesive layer, as seen in Figure 12. Cohesive zone modelling is generally used for the numerical simulation of interlaminar failure. Damage initiation is driven by a traction seperation law and the value of the maximum traction, to, Figure 14(a). New crack surfaces are formed when the fracture toughness Gc is equal to the area surface under the tractionseparation curve. Considering the nature of predictions made by the parametric study, the model was analysed with cohesive interaction properties of adhesive FM 300 and listed in Table 4.

Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach 289

(2)

Cohesive layer Prediction

damage variable, D, and is implemented in the damage evolution model. After the damage initiation, the damage variable monotonically evolves from 0 to 1 on increasing the loading.

For the mix-mode definition of fracture energies, Benzeggagh-Kenane's criterion [16] (BK

( ) *c cc <sup>S</sup> <sup>c</sup>*

*GGG Tns* = + (3)

*G GG Sts* = + (4)

where *Gn, Gs* and *Gt* are the fracture toughness values in the normal and two shear directions respectively which were measured in-house for IM7/8552 and given in Table 4. In this study the value of the BK mode-mixity power parameter, η=*1.6*, was obtained from Maimi [17, 18].

The debonding failure initiation load predicted by the model was compared with the debonding loads predicted by using the energy release rate and the experimental results,

Prediction

Recalling the results from the energy release rate analysis, when the maximum energy release rate across the width of the stiffener was used, the predictions were closer to the experimental results. In addition, when the average energy release rate was used, the difference from the experimental results had an over-prediction of 8.5%. By using the

Also, the debonding growth along the width of the specimen correlates well with the predictions using strain energy release rate across the width of the stiffener as can be seen in Figure 15. Comparing the FE model pattern with the *Tapered* slop on the right, it is clearly observed that there is a correlation between the extent of damaged and the value of

Failure load [kN] 17.7 19.2 17.5 18.7

experiments - 8.5% 1.1% 5.6%

Gmax Prediction

Experiment Gavg

cohesive zone model, the difference with the experimental results was 5.6%.

normalized strain energy release rate across the width of the specimen.

+− =

*<sup>G</sup> G GG G*

*T*

η

*G*

law) was used which defines the energy dissipated as in equation below,

**3.2. Response of the Numerical Model** 

with,

Table 5.

Difference from

**Table 5.** Debonding loads

*n sn*

**Figure 14.** Illustration of (a) Traction-separation law for cohesive zone models (b) Modified law to implement in FEM


**Table 4.** Cohesive interaction properties

### **3.1. Initial linear elastic behaviour**

In order to accurately predict damage initiation, the interaction of traction components was taken into account and the quadratic stress criterion used. This criterion is included in Abaqus and was formulated based on Ye's criterion [13] including the interaction between traction components. According to this criterion, the damage initiates when a quadratic interaction function reaches a value of one. The quadratic interaction function is shown in the following equation.

$$\left\{\frac{t\_n}{t\_n^0}\right\}^2 + \left\{\frac{t\_s}{t\_s^0}\right\}^2 + \left\{\frac{t\_t}{t\_t^0}\right\}^2 = \mathbf{1} \tag{1}$$

where, *t0n , t0s* and *t0t* are the peak values of the stress when the separation is either purely normal to the interface or purely in the first or the second shear direction respectively. In the present numerical study, these values for IM7/8552 are 50 MPa in the normal direction and 100 MPa in the two shear directions[10], see Table 4.

When the initiation criterion is met, the cohesive stiffness degrades at a rate defined by the damage evolution model. The overall damage of the damage zone is represented by a scalar damage variable, D, and is implemented in the damage evolution model. After the damage initiation, the damage variable monotonically evolves from 0 to 1 on increasing the loading.

For the mix-mode definition of fracture energies, Benzeggagh-Kenane's criterion [16] (BK law) was used which defines the energy dissipated as in equation below,

$$\left\{ \mathbf{G}\_n^c + \left( \mathbf{G}\_s^c - \mathbf{G}\_n^c \right) \left\{ \frac{\mathbf{G}\_S}{\mathbf{G}\_T} \right\}^\eta = \mathbf{G}^c \right. \tag{2}$$

with,

288 Finite Element Analysis – New Trends and Developments

implement in FEM

FM300 (adhesive)

Initial linear elastic behaviour

> Damage initiation

Damage evolution

the following equation.

**Table 4.** Cohesive interaction properties

**3.1. Initial linear elastic behaviour** 

100 MPa in the two shear directions[10], see Table 4.

**Figure 14.** Illustration of (a) Traction-separation law for cohesive zone models (b) Modified law to

First shear direction

*Knn in Nmm-3 Kss in Nmm-3 Ktt in Nmm-3*

106 106 106

*N in MPa S1 in MPa S2 in MPa*  50 100 100

0.9 2.5 2.5

In order to accurately predict damage initiation, the interaction of traction components was taken into account and the quadratic stress criterion used. This criterion is included in Abaqus and was formulated based on Ye's criterion [13] including the interaction between traction components. According to this criterion, the damage initiates when a quadratic interaction function reaches a value of one. The quadratic interaction function is shown in

> 222 <sup>000</sup> <sup>1</sup> *nst ns t ttt ttt* ++=

where, *t0n , t0s* and *t0t* are the peak values of the stress when the separation is either purely normal to the interface or purely in the first or the second shear direction respectively. In the present numerical study, these values for IM7/8552 are 50 MPa in the normal direction and

When the initiation criterion is met, the cohesive stiffness degrades at a rate defined by the damage evolution model. The overall damage of the damage zone is represented by a scalar

*Gn in N/mm Gs in N/mm Gt in N/mm* <sup>8</sup>

Second shear direction

(1)

BK law η

Normal direction

$$\mathbf{G}\_T = \mathbf{G}\_n + \mathbf{G}\_s \tag{3}$$

$$\mathbf{G}\_{\rm S} = \mathbf{G}\_{t} + \mathbf{G}\_{s} \tag{4}$$

where *Gn, Gs* and *Gt* are the fracture toughness values in the normal and two shear directions respectively which were measured in-house for IM7/8552 and given in Table 4. In this study the value of the BK mode-mixity power parameter, η=*1.6*, was obtained from Maimi [17, 18].

#### **3.2. Response of the Numerical Model**

The debonding failure initiation load predicted by the model was compared with the debonding loads predicted by using the energy release rate and the experimental results, Table 5.


**Table 5.** Debonding loads

Recalling the results from the energy release rate analysis, when the maximum energy release rate across the width of the stiffener was used, the predictions were closer to the experimental results. In addition, when the average energy release rate was used, the difference from the experimental results had an over-prediction of 8.5%. By using the cohesive zone model, the difference with the experimental results was 5.6%.

Also, the debonding growth along the width of the specimen correlates well with the predictions using strain energy release rate across the width of the stiffener as can be seen in Figure 15. Comparing the FE model pattern with the *Tapered* slop on the right, it is clearly observed that there is a correlation between the extent of damaged and the value of normalized strain energy release rate across the width of the specimen.

Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach 291

**Figure 17.** The refined model that used in the VCCT method

**4.1. Skin stiffener runout configurations** 

The objectives defined in the first part of this study were to create two FE models for the baseline and the modified configuration. The baseline stiffeners failed due to debonding of the skin/stiffener interface, while the modified stiffeners failed by delamination between the 00 and 450 ply interface. Improved damage tolerance (stable crack propagation) was reported in the modified stiffener run-out design as compared to the baseline configuration. The modified design eventually failed catastrophically by interlaminar delamination, not bondline failure, which had not been considered in the numerical study. A more detailed analysis of different configurations, which accounted for delamination, was therefore undertaken. Building on the

The structural performance of three different skin-stiffener configurations – *Baseline* (B), *Tapered* (T) and *Compliant* (C) – under longitudinal compression, with geometry and dimensions shown in Figure 18, was assessed. Compared to the *Baseline* stiffener (Figure 18a), the other two configurations have a widening flange towards the termination end of the stiffener but this added material is offset by the taper of the stiffener web (*Tapered* configuration, Figure 18b). The third configuration includes a curved tape (*Compliant*, Figure

The proposed *Compliant* design was developed by considering potentially beneficial local stiffness variations. This resulted in a stiffener design with a similar overall weight to the *Baseline* design. For the modified configurations, various values of *b*, *c* and *d* (see Figure 18b)

The failure mode of the *Tapered* specimen configuration was found to be a combination of interlaminar and intralaminar failure [12]. The *Tapered* stiffeners had an unexpected delamination in the flange between the bottom-most 0o ply and above the 45o ply.

previous findings, the merits of a compliant termination scheme are presented.

**4. 2nd iteration** 

18c).

were analysed [12].

**Figure 15.** Damage growth pattern compared with energy release rate predictions for tapered specimen

#### **3.3. Modeling the debonding failure using VCCT**

The implemented VCCT method in Abaqus standard, developed by Boeing and Simulia, suggested promising results, especially when the mismatched, tie-constrained, meshes had only 3% error comparing with pairing meshes. As can be seen in Figure 16(a) there was good correlation between the VCCT and the parametric study. On the other hand, the VCCT method could not capture the increase in normalised GT/GC arising from the geometric discontinuity at the edge of the flange. In order to capture the detail at the edge of the flange, the mesh resolution was increased and was biased towards these edges, Figure 17, and the results can be seen in Figure 16(b). Despite the good results, the size of the model and the time needed for execution made the use of VCCT impractical and was not used in further investigations.

**Figure 16.** Comparing the results of the parametric study with the VCCT method (a) along the crack and (b) along the width of the stiffener.

**Figure 17.** The refined model that used in the VCCT method

## **4. 2nd iteration**

290 Finite Element Analysis – New Trends and Developments

**3.3. Modeling the debonding failure using VCCT** 

0 1 2 3 4 5 6 7 8 9 10

and (b) along the width of the stiffener.

0.6 0.7 0.8 0.9 1 1.1 1.2 1.3

0.5

0.7

0.9

Normalised

(a) (b)

*a* [mm]

**Figure 16.** Comparing the results of the parametric study with the VCCT method (a) along the crack

VCCT script results Reference Data

Baseline Stiffener VCCT results

*GT/GC*

1.1

1.3

1.5

Normalised *GT / Gc*

**Figure 15.** Damage growth pattern compared with energy release rate predictions for tapered specimen

The implemented VCCT method in Abaqus standard, developed by Boeing and Simulia, suggested promising results, especially when the mismatched, tie-constrained, meshes had only 3% error comparing with pairing meshes. As can be seen in Figure 16(a) there was good correlation between the VCCT and the parametric study. On the other hand, the VCCT method could not capture the increase in normalised GT/GC arising from the geometric discontinuity at the edge of the flange. In order to capture the detail at the edge of the flange, the mesh resolution was increased and was biased towards these edges, Figure 17, and the results can be seen in Figure 16(b). Despite the good results, the size of the model and the time needed for execution made the use of VCCT impractical and was not used in

02468

Baseline Tapered Compliant

Distance across width [mm]

01234567

Distance along the crack tip [mm]

Baseline Stiffener VCCT script results

VCCT results

Baseline Stiffener

Correlation of debonding damage growth across width of

the stiffener

further investigations.

0.9

1

Normalised

*GT/GC*

1.1

1.2

The objectives defined in the first part of this study were to create two FE models for the baseline and the modified configuration. The baseline stiffeners failed due to debonding of the skin/stiffener interface, while the modified stiffeners failed by delamination between the 00 and 450 ply interface. Improved damage tolerance (stable crack propagation) was reported in the modified stiffener run-out design as compared to the baseline configuration. The modified design eventually failed catastrophically by interlaminar delamination, not bondline failure, which had not been considered in the numerical study. A more detailed analysis of different configurations, which accounted for delamination, was therefore undertaken. Building on the previous findings, the merits of a compliant termination scheme are presented.

## **4.1. Skin stiffener runout configurations**

The structural performance of three different skin-stiffener configurations – *Baseline* (B), *Tapered* (T) and *Compliant* (C) – under longitudinal compression, with geometry and dimensions shown in Figure 18, was assessed. Compared to the *Baseline* stiffener (Figure 18a), the other two configurations have a widening flange towards the termination end of the stiffener but this added material is offset by the taper of the stiffener web (*Tapered* configuration, Figure 18b). The third configuration includes a curved tape (*Compliant*, Figure 18c).

The proposed *Compliant* design was developed by considering potentially beneficial local stiffness variations. This resulted in a stiffener design with a similar overall weight to the *Baseline* design. For the modified configurations, various values of *b*, *c* and *d* (see Figure 18b) were analysed [12].

The failure mode of the *Tapered* specimen configuration was found to be a combination of interlaminar and intralaminar failure [12]. The *Tapered* stiffeners had an unexpected delamination in the flange between the bottom-most 0o ply and above the 45o ply.

Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach 293

Furthermore, the delamination led to an intralaminar failure in the form of a matrix crack across the 0o ply near the filler ends and continued delaminating between the filler and 0o ply, Figure 19a. For this reason, a new FE model was created in order to investigate these

The three different configurations, *Baseline*, *Tapered* and *Compliant* (see Figure 18), were analysed for debonding and delamination growth stability. The results of this analysis are presented in Figure 20, where the values of *GT* = *G<sup>Ι</sup>* + *GΙΙ* + *GΙΙΙ*, the total strain energy release rate, have been normalized to the *GT* of the reference parametric stiffener for 0.5 mm of crack length. Figure 20 compares the normalised strain energy release rate for the *Baseline,* the *Tapered* and the *Compliant* stiffeners. The negative slope of the *G(a)* curve indicates crack growth stability, while a positive slope indicates instability (assuming constant fracture toughness). Recalling the failure modes obtained experimentally [12] the *Baseline* stiffener failed by debonding and the *Tapered* stiffener initially experienced debonding until it finally failed by delamination. This is in agreement with predictions, Figure 20. Consequently, both models were able to correctly describe these experimental results [12]. In addition, the stability analysis for the *Compliant* stiffener predicts that this design will fail stably by

Baseline Tapered Compliant

012345

**Figure 20.** Normalized strain energy release rates as a function of crack length; comparison between *Baseline* stiffener design, *Tapered* stiffener and *Compliant* stiffener with *b*=3 mm, *c*=10 mm and *d*=6.25 mm.

*a* [mm]

modes of damage, Figure 19b and Figure 19c.

Debonding

Delamination

**4.2. Energy release rate along crack** 

debonding, Figure 20.

0

0.2

0.4

0.6

0.8

Normalized

*GT*

1

1.2

1.4

(c) *Compliant* stiffener configuration

**Figure 19.** a) *Tapered* stiffener after testing, b) FE model showing delamination path, and c) FE model of a specimen with boundary conditions.

Furthermore, the delamination led to an intralaminar failure in the form of a matrix crack across the 0o ply near the filler ends and continued delaminating between the filler and 0o ply, Figure 19a. For this reason, a new FE model was created in order to investigate these modes of damage, Figure 19b and Figure 19c.

#### **4.2. Energy release rate along crack**

292 Finite Element Analysis – New Trends and Developments

2.5

15 30 15 2

*d*

(c) *Compliant* stiffener configuration

*<sup>b</sup> <sup>c</sup>* <sup>20</sup>

(a) (b)

(c)

**Figure 19.** a) *Tapered* stiffener after testing, b) FE model showing delamination path, and c) FE model

50

30

<sup>10</sup> <sup>20</sup>

(b) *Tapered* stiffener configuration

(42.5, 11.5)

Delamination area

(30.5, 4.2)

*<sup>b</sup> <sup>c</sup>*

(15, 8)

Displacement

2nd order polynomial

20

(0,0)

*d*

*a*

10

(a) *Baseline* stiffener configuration

30

of a specimen with boundary conditions.

<sup>10</sup> <sup>20</sup>

**Figure 18.** Stiffener design configurations (dimensions in mm).

1.25

100

50

The three different configurations, *Baseline*, *Tapered* and *Compliant* (see Figure 18), were analysed for debonding and delamination growth stability. The results of this analysis are presented in Figure 20, where the values of *GT* = *G<sup>Ι</sup>* + *GΙΙ* + *GΙΙΙ*, the total strain energy release rate, have been normalized to the *GT* of the reference parametric stiffener for 0.5 mm of crack length. Figure 20 compares the normalised strain energy release rate for the *Baseline,* the *Tapered* and the *Compliant* stiffeners. The negative slope of the *G(a)* curve indicates crack growth stability, while a positive slope indicates instability (assuming constant fracture toughness). Recalling the failure modes obtained experimentally [12] the *Baseline* stiffener failed by debonding and the *Tapered* stiffener initially experienced debonding until it finally failed by delamination. This is in agreement with predictions, Figure 20. Consequently, both models were able to correctly describe these experimental results [12]. In addition, the stability analysis for the *Compliant* stiffener predicts that this design will fail stably by debonding, Figure 20.

**Figure 20.** Normalized strain energy release rates as a function of crack length; comparison between *Baseline* stiffener design, *Tapered* stiffener and *Compliant* stiffener with *b*=3 mm, *c*=10 mm and *d*=6.25 mm.

#### **4.3. Energy release rate along the width of the crack tip**

The strain energy release rate along the width of the crack tip was calculated for the *Baseline*, *Tapered* and *Compliant* configurations (Figure 18). Fracture initiation is expected when the *GT* exceeds the fracture toughness *Gc* for a given mixed-mode ratio *GII* / *GT* at each point along the crack tip. In other words, propagation at each point occurs when *GT* / *Gc* >1 [19, 20]. The interlaminar fracture toughness *Gc* can be calculated by using the equation 2 [20]

The value of *Gc* is normalised to the width-average value for the *Tapered* specimen. It can be observed that the trend is similar for the *Baseline* and *Tapered* specimen types but is different in the centre of the *Compliant* stiffener. This is due to the difference in the web of the stiffeners. The curved taper has reduced the normalized strain energy release rate in the centre without affecting the trend in the flange. The maximum value of the energy release rate can be used to predict the load corresponding to the initiation of fracture using

$$\frac{P}{P\_{FE}} = \sqrt{\frac{G\_c}{G\_T}}\tag{5}$$

Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach 295

Experimental failure load [kN]

> 0.34 0.39 16.49<sup>+</sup> −

> 0.16 0.22 17.72<sup>+</sup> −

> 0.16 0.29 18.02<sup>+</sup> −

The tests were carried out in an Instron testing machine, equipped with a 100 KN load cell, at a loading rate of 0.5 mm/min. Load and crosshead displacement were recorded continuously by a PC data logger connected to the load cell and the Instron machine at a sampling rate of 2 Hz. The specimens were aligned by careful measurement in the loading direction to avoid bending. The Imperial Data Acquisition (IDA) program was used to

> Predicted failure load [kN] (% difference with respect to experimental)

Based on *G*avg Based on *G*max

**Table 6.** Failure loads for the different specimen types, as well as the predicted failure loads using Eq. 2.

AE sensors were used to identify and investigate failure, within the specimens, during testing. The AE equipment was manufactured by Physical Acoustic Corporation (PAC) and failure was monitored by AEwin software. Broadband (WD) sensors with an operating frequency range of 100 Hz to 1000 kHz were used and positioned in order to obtain the best

The *Baseline* stiffeners had an average failure load of 16.5 kN while the *Tapered* stiffeners had an average failure load of 17.7 kN and the *Compliant* an average failure load of 18 kN, Table 6. The fracture surfaces for selected specimens are shown in Figure 23, and the load versus displacement curves for selected specimens of the three stiffener designs are shown in Figure 22. The predicted loads (using Eq. 2) match well with the experimental values when

The acoustic emission signals (Figure 22) show that there was an increase in AE activity 0.1 mm before catastrophic failure for the *Baseline* specimen. For the *Tapered* specimen type, the increase in AE emission started about 0.05 mm before catastrophic failure and for the *Compliant* specimen, 0.02 mm. Figure 22 also shows the peak frequency during the tests for all specimen types. It can be observed that there is some very low-energy micro-cracking from the start of the test and this is possibly at the resin pots. The *Tapered* specimen configuration promoted a combination of failure modes including delamination and fibre bridging which preceded catastrophic failure. In addition, the *Compliant* stiffener, according

to AE data and as visually observed (Figure 23c), suffered only from debonding.

16.56 (+0.4%)

17.45 (-1.5%)

18.17 (+0.8%)

(+15.2%)

(+8.2%)

19.93 (+10.6%)

results without affecting the specimens behaviour [21].

the maximum *G* across the width is used, Table 6.

**5. Experiments** 

record load and displacement during the tests.

*Baseline* **Stiffener** 19.00

*Tapered* **Stiffener** 19.17

*Compliant* **Stiffener** 

where *P* is the load at initiation of fracture, *PFE* is the load from the FE model, *Gc* the critical strain energy release rate (Equation 1.26), and *GT* is the strain energy release rate predicted by the FE model as defined previously. Two different predictions for *P* can be made: one using the maximum value of *GT* along the width, and another using the average, Table 6.

**Figure 21.** Normalized *G*T/*G*c across the crack tip for crack *a* = 1 mm for the *Baseline*, *Tapered* and *Compliant* specimens.

## **5. Experiments**

294 Finite Element Analysis – New Trends and Developments

0.6

0.7

0.8

0.9

Normalised

*Compliant* specimens.

*GT / Gc*

1

1.1

1.2

1.3

**4.3. Energy release rate along the width of the crack tip** 

The strain energy release rate along the width of the crack tip was calculated for the *Baseline*, *Tapered* and *Compliant* configurations (Figure 18). Fracture initiation is expected when the *GT* exceeds the fracture toughness *Gc* for a given mixed-mode ratio *GII* / *GT* at each point along the crack tip. In other words, propagation at each point occurs when *GT* / *Gc* >1 [19, 20]. The

The value of *Gc* is normalised to the width-average value for the *Tapered* specimen. It can be observed that the trend is similar for the *Baseline* and *Tapered* specimen types but is different in the centre of the *Compliant* stiffener. This is due to the difference in the web of the stiffeners. The curved taper has reduced the normalized strain energy release rate in the centre without affecting the trend in the flange. The maximum value of the energy release

*c*

Baseline Tapered Compliant

*P G* <sup>=</sup> (5)

*FE T P G*

where *P* is the load at initiation of fracture, *PFE* is the load from the FE model, *Gc* the critical strain energy release rate (Equation 1.26), and *GT* is the strain energy release rate predicted by the FE model as defined previously. Two different predictions for *P* can be made: one using the maximum value of *GT* along the width, and another using the average, Table 6.

02468

Distance across width [mm]

**Figure 21.** Normalized *G*T/*G*c across the crack tip for crack *a* = 1 mm for the *Baseline*, *Tapered* and

interlaminar fracture toughness *Gc* can be calculated by using the equation 2 [20]

rate can be used to predict the load corresponding to the initiation of fracture using

The tests were carried out in an Instron testing machine, equipped with a 100 KN load cell, at a loading rate of 0.5 mm/min. Load and crosshead displacement were recorded continuously by a PC data logger connected to the load cell and the Instron machine at a sampling rate of 2 Hz. The specimens were aligned by careful measurement in the loading direction to avoid bending. The Imperial Data Acquisition (IDA) program was used to record load and displacement during the tests.


**Table 6.** Failure loads for the different specimen types, as well as the predicted failure loads using Eq. 2.

AE sensors were used to identify and investigate failure, within the specimens, during testing. The AE equipment was manufactured by Physical Acoustic Corporation (PAC) and failure was monitored by AEwin software. Broadband (WD) sensors with an operating frequency range of 100 Hz to 1000 kHz were used and positioned in order to obtain the best results without affecting the specimens behaviour [21].

The *Baseline* stiffeners had an average failure load of 16.5 kN while the *Tapered* stiffeners had an average failure load of 17.7 kN and the *Compliant* an average failure load of 18 kN, Table 6. The fracture surfaces for selected specimens are shown in Figure 23, and the load versus displacement curves for selected specimens of the three stiffener designs are shown in Figure 22. The predicted loads (using Eq. 2) match well with the experimental values when the maximum *G* across the width is used, Table 6.

The acoustic emission signals (Figure 22) show that there was an increase in AE activity 0.1 mm before catastrophic failure for the *Baseline* specimen. For the *Tapered* specimen type, the increase in AE emission started about 0.05 mm before catastrophic failure and for the *Compliant* specimen, 0.02 mm. Figure 22 also shows the peak frequency during the tests for all specimen types. It can be observed that there is some very low-energy micro-cracking from the start of the test and this is possibly at the resin pots. The *Tapered* specimen configuration promoted a combination of failure modes including delamination and fibre bridging which preceded catastrophic failure. In addition, the *Compliant* stiffener, according to AE data and as visually observed (Figure 23c), suffered only from debonding.

Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach 297

(a)

(b)

(c)

**Figure 23.** (a) *Baseline* stiffener, (b) *Tapered* Stiffener and (c) *Compliant* Stiffener after failure respectively.

The strain energy release rate analysis yielded good results in the investigation of the runout design influence in debonding/delamination for stiffener terminations. The FE models

**6. Remarks** 

(c)

**Figure 22.** Loads and Peak frequencies versus displacement for a) the *Baseline* b) the *Tapered* and c) the *Compliant* stiffeners. A scale on the right hand side indicates the mode of failure typically associated with these peak frequencies [21].

(a)

(b)

**Figure 23.** (a) *Baseline* stiffener, (b) *Tapered* Stiffener and (c) *Compliant* Stiffener after failure respectively.

#### **6. Remarks**

296 Finite Element Analysis – New Trends and Developments

0

0

0

with these peak frequencies [21].

5

10

Load [kN]

15

20

5

10

Load [kN]

15

20

5

10

Load [kN]

15

20

(a)

Tapered Stiffener

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Dislacement [mm]

(b)

(c)

**Figure 22.** Loads and Peak frequencies versus displacement for a) the *Baseline* b) the *Tapered* and c) the *Compliant* stiffeners. A scale on the right hand side indicates the mode of failure typically associated

Displacement [mm]

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Compliant Stiffener

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Displacement [mm]

Baseline Stiffener

Matrix cracking Delamination

Peak Frequency [KHz]

Peak Frequency [kHz]

Peak Frequecy [kHz]

Matrix cracking Delamination

Fibre / matrix debonding

Fibre failure

Fibre pull-out

Matrix cracking Delamination

Fibre / matrix debonding

Fibre failure

Fibre pull-out

Fibre / matrix debonding

Fibre failure

Fibre pull-out

The strain energy release rate analysis yielded good results in the investigation of the runout design influence in debonding/delamination for stiffener terminations. The FE models accurately predicted the failure loads and failure modes for the specimens tested and the predictions were improved when the distributions of the strain energy release rate across the width was considered. The differences in the predictions using the average and the maximum energy release rates are shown in Table 6 and can be compared to the experimental failure loads. The load-displacement, as well as the peak frequencydisplacement plots (Figure 22), show that the *Tapered* design is slightly more damage tolerant than the *Baseline* one and this improved further with the *Compliant* design. The AE monitoring proved to be valuable in detecting and analysing the failure modes experienced by the specimens.

Damage-Tolerant Design of Stiffener Run-Outs: A Finite Element Approach 299

[4] H. Hosseini-Toudeshky*, et al.*, "Analysis of composite skin/stiffener debounding and failure under uniaxial loading," *Composite Structures,* vol. 75, pp. 428-436, 2006. [5] S. Mahdi*, et al.*, "The mechanical performance of repaired stiffened panels. Part II. Finite

[6] R. Krueger and P. J. Minguet, "Analysis of composite skin-stiffener debond specimens using a shell/3D modeling technique," *Composite Structures,* vol. 81, pp. 41-59, 2007. [7] B. G. Falzon and G. A. O. Davies, "The Behavior of Compressively Loaded Stiffener Runout Specimens – Part I:Experiments," *Journal of composite materials,* vol. 37, 2002. [8] B. G. Falzon and D. Hitchings, "The Behavior of Compressively Loaded Stiffener Runout Specimens – Part II: Finite Element Analysis," *Journal of composite materials,* vol.

[9] R. Krueger, "Virtual crack closure technique: history, approach and applications," *Appl* 

[10] C. Bisagni*, et al.*, "Assessment of the damage tolerance of postbuckled hat-stiffened panels using single-stringer specimens," *AIAA-2010-2696, 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference,* 

[11] P. P. Camanho*, et al.*, ". Fracture analysis of composite co-cured structural joints using decohesion elements," *Fatigue Fract Engng Mater struct 27,* pp. 745- 757, 2004. [12] S. Psarras*, et al.*, "Design of composite stiffener run-outs for damage tolerance," *Finite* 

[13] Simulia, *Rising Sun Mills, 166 Valley Street, Providence, RI 02909-2499, USA, ABAQUS* 

[14] Python, "Python Software Foundation (PSF)," *Wolfeboro Falls, NH 03896-0037, PO Box 37,* 

[15] D. Hitchings, "Finite element modelling of composite materials and structures,"

[16] M. K. Benzeggagh and M. Kenane, "Measurement of mixed-mode delamination fracture toughness of unidirectional class/epoxy with mixed-mode bending apparatus,"

[17] P. Maimi*, et al.*, "A continuum damage model for composite laminates: Part I-

[18] P. Maimi*, et al.*, "A continuum damage model for composite laminates: Part II-Computational implementation and validation" *Mechanics of Materials,* vol. To appear,

[19] R. Krueger*, et al.*, "Panel stiffener debonding analysis using a shell/3D modeling

[20] M. L. Benzeggagh and M. Kenane, "Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending

technique," *Composites Science and Technology,* vol. 69, pp. 2352-2362, 2009.

apparatus," *Composites Science and Technology,* vol. 56, pp. 439-49, 1996.

element modelling," *Composites Part B:Engineering,* vol. 33, pp. 355-366, 2002.

37, 2002.

*6.10,* 2011.

2007.

*USA, Python 2.6.2,* 2010.

*Mech Rev,* vol. 57:109–43, 2004.

12-15 April , Orlando, Florida 2010.

*Elements in Analysis and Design,* vol. 47, pp. 949-954, 2011.

*Cambridge (UK): WOODHEAD PUBLISHING LIMITED,* 2009.

*Composites Science and Technology,* vol. 56, pp. 439-449, 1996.

Constitutive model," *Mechanics of Materials,* vol. To appear, 2007.

## **7. Conclusions**

This study was based on the strain energy release rates for debonding and delamination and successfully predicted the failure loads for the three different specimen types. The predictions were more accurate when the maximum strain energy release rate across the width was used. It can be concluded that the variation of the energy release rate across the width should be considered when stiffener run-outs are designed. AE data recorded during skin-stiffener run-out compression tests proved useful to analyse the failure processes which take place in these specimens. The results show that in the design of skin-stiffener run-outs it is important to consider the possibility of failure modes other than debonding, and that compliant termination schemes offer the possibility of improved damage tolerance.

## **Author details**

S. Psarras\* and S.T. Pinho *Dept. Aeronautics, Imperial College London, SW7 2AZ, United Kingdom* 

B.G. Falzon *Dept. Mechanical and Aerospace Engineering, Monash University, Victoria, Australia* 

## **8. References**


<sup>\*</sup> Corresponding Author

[4] H. Hosseini-Toudeshky*, et al.*, "Analysis of composite skin/stiffener debounding and failure under uniaxial loading," *Composite Structures,* vol. 75, pp. 428-436, 2006.

298 Finite Element Analysis – New Trends and Developments

by the specimens.

**7. Conclusions** 

**Author details** 

and S.T. Pinho

31A, pp. 459-68, 2000.

S. Psarras\*

B.G. Falzon

 \*

**8. References** 

Corresponding Author

accurately predicted the failure loads and failure modes for the specimens tested and the predictions were improved when the distributions of the strain energy release rate across the width was considered. The differences in the predictions using the average and the maximum energy release rates are shown in Table 6 and can be compared to the experimental failure loads. The load-displacement, as well as the peak frequencydisplacement plots (Figure 22), show that the *Tapered* design is slightly more damage tolerant than the *Baseline* one and this improved further with the *Compliant* design. The AE monitoring proved to be valuable in detecting and analysing the failure modes experienced

This study was based on the strain energy release rates for debonding and delamination and successfully predicted the failure loads for the three different specimen types. The predictions were more accurate when the maximum strain energy release rate across the width was used. It can be concluded that the variation of the energy release rate across the width should be considered when stiffener run-outs are designed. AE data recorded during skin-stiffener run-out compression tests proved useful to analyse the failure processes which take place in these specimens. The results show that in the design of skin-stiffener run-outs it is important to consider the possibility of failure modes other than debonding, and that

compliant termination schemes offer the possibility of improved damage tolerance.

*Dept. Mechanical and Aerospace Engineering, Monash University, Victoria, Australia* 

*(Applied Science and Manufacturing),* vol. 35A, pp. 1447-58, 2004.

[1] A. Faggiani and B. G. Falzon, "Numerical analysis of stiffener runout sections," *Applied* 

[2] B. G. Falzon*, et al.*, "Postbuckling behaviour of a blade-stiffened composite panel loaded in uniaxial compression," *Composites Part A (Applied Science and Manufacturing),* vol.

[3] E. Greenhalgh and M. H. Garcia, "Fracture mechanisms and failure processes at stiffener run-outs in polymer matrix composite stiffened elements," *Composites Part A* 

*Dept. Aeronautics, Imperial College London, SW7 2AZ, United Kingdom* 

*Composite Materials,* vol. 14, pp. 145-58, 2007.

	- [21] R. Gutkin*, et al.*, "On acoustic emission for failure investigation in CFRP," *Mechanical Systems and Signal Processing 25,* vol. 1393–1407, 2011.

**Chapter 14** 

© 2012 Roatesi, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Roatesi, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Finite Element Analysis for the Problem of** 

**Tunnel Excavation Successive Phases and** 

Numerical application to geomechanics is an area of research that has grown rapidly since its origins in the late 1960s. This growth led to new methodologies and analysis approaches that are nowadays commonly employed in geotechnical, mining or petroleum engineering

The circular tunnel boundary problem is frequently encountered in rock and soil mechanics, geotechnics and generally in mining and petroleum engineering. Consequently, there is a great interest in solving boundary problems involving the excavation of an underground opening. Moreover, time dependent behavior of geomaterials benefits by a special attention in constitutive and numerical approach. Regarding the constitutive modeling and analytical approach [1] concerns both with vertical and horizontal underground openings excavated in an elasto-viscoplastic rock mass, governed by the constitutive laws which go by his name and for which analytical solution for the displacement and viscoplastic strain are derived; [2] provides analytical results for circular and non-circular openings with thermal effects as well; [3-5] present results in approaching the problem both analytically (convergence-

The viscoplastic approach in FEM has been considered by many authors, in terms either of numerical stability of schemes used in viscoplasticity, see [6-8], or practical applications using FEM solutions in different areas. For instance, tunneling using FEM for viscoplastic materials has been investigated widely, as well, e.g. [9] approaches the problem with outstanding results in computational geomechanics, with special reference to earthquake engineering, in numerical modeling of dynamic soil and pore fluid interaction and earthquake-induced liquefaction and multiphase pollutant transport in partially saturated

**Lining Mounting** 

http://dx.doi.org/10.5772/46168

Additional information is available at the end of the chapter

confinement method) and numerically by FEM.

Simona Roatesi

**1. Introduction** 

practice.
