Aug 20, 2012 - many types of contemporary high-performance military aircraft. .... methods of repairing by patching damaged regions [16] or, for self-...

0 downloads 0 Views 2MB Size

Dynamic Data Driven Application Systems for Monitoring Damage in Composite Materials Under Dynamic Loads by P. T. Bauman, J. T. Oden, E. E. Prudencio, S. Prudhomme, and K. Ravi-Chandar

The Institute for Computational Engineering and Sciences The University of Texas at Austin Austin, Texas 78712

Reference: P. T. Bauman, J. T. Oden, E. E. Prudencio, S. Prudhomme, and K. Ravi-Chandar, Dynamic Data Driven Application Systems for Monitoring Damage in Composite Materials Under Dynamic Loads, ICES REPORT 12-37, The Institute for Computational Engineering and Sciences, The University of Texas at Austin, August 2012.

Dynamic Data Driven Application Systems for Monitoring Damage in Composite Materials Under Dynamic Loads

P. T. Bauman, J. T. Oden, E. E. Prudencio, S. Prudhomme, and K. Ravi-Chandar.

Institute for Computational Engineering and Sciences (ICES) The University of Texas at Austin

August 20, 2012

Contents 1 Introduction

1

2 Physical Example of Damage Identification 2.1 Material Model for Damage . . . . . . . . . . . . . . . . . 2.1.1 Inhomogeneous Specimen . . . . . . . . . . . . . . 2.1.2 Indirect Measurement of Damage . . . . . . . . . . 2.1.3 Three-Phase Homogenized Model for Conductivity 2.1.4 Percolation Model for Conductivity . . . . . . . . .

. . . . .

3 3 5 7 7 8

3 Bayesian Methods 3.1 The General Bayesian Filtering . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 The Extended Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . .

11 12 13 14

4 Experiments

15

5 Integration 5.1 Using the Generated Experimental Data . 5.2 Damage Model Implementation . . . . . . 5.3 Kalman Filter Integration . . . . . . . . . 5.4 A Simple Filtering Problem . . . . . . . . 5.4.1 The Forward Problem Perspective . 5.4.2 The Bayesian Filtering Perspective

16 16 16 16 17 17 18

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . .

. . . . . .

. . . . . .

6 Future Work

24

References

25

i

1

Introduction

This document summarizes the status of ongoing work on the development of a stochastic data-driven system for prediction of material damage in composite materials common to many types of contemporary high-performance military aircraft. The projects is developed around several key ideas and technologies, the principal components being as follows: 1. The physical problem under study is the behavior of thin composite structural components under dynamic (generally oscillatory) loads, such as are common in aircraft fuselages and wings; 2. Continuum damage theories are used to model the onset and evolution of damage, generally in the form of micro-crack densities; 3. A Dynamic Data-Driven Application System (DDDAS) is to be developed to monitor and control the extent of damage accumulated in laboratory experiments performed on composite specimens; and 4. A Bayesian framework is developed for defining, updating, and quantifying uncertainties in the model, the experiment data, and the target quantities of interest. The underlying philosophy of DDDAS is to dynamically connect computational models of the evolution of physical phenomena of interest with experimental systems that deliver relevant data in near-real time so as to allow feedback to control outputs to meet a set of objectives (see [7, 8, 9, 10, 11, 12, 13]). In the present study, the computational models are based on finite element models of highly nonlinear material damage theories of the type used in contemporary fatigue analysis, fracture mechanics, and structural mechanics. These typically involve material parameters that exhibit uncertainties. On the experimental side, an experimental program for monitoring damage is set up which itself is based on a secondary mathematical model with its own uncertainties. This is manifested in the form of a network of Carbon Nano-Tubes (CNT), covering a surface of the test specimen, on which an electric charge is supplied. Local damage is attributed to changes in the induced field and the local current due to micro-crack developed in the substrate material. Thus the system itself must be calibrated and validated, and the inherent uncertainties in data be factored into a statistical analysis of the validation of the full system. The current generation of high-performance military aircraft and the next generation of civil transport aircraft are made almost entirely of composite materials. With repeated loading cycles or accidental overloading such as foreign-object impact events, composite materials accumulate damage that degrades their performance and eventually leads to failure. Aging of this type of aircraft is likely to be quite different from the metal-skin aircraft for which there is almost one century of experience in characterizing fatigue, multisite damage, and other aging phenomena. Therefore, it is essential to monitor and track the evolution of damage in these composite materials. In the early days of damage monitoring, the focus was on nondestructive evaluation (NDE) using inspection techniques such as visual observations, ultrasonic or x-rays, imaging and other modalities to reveal the damage state. These techniques were used during periodic interruptions of service of the vehicle. Recent advances have altered the methodology to 1

embed sensors or sensor arrays into the structure and perform continuous acquisition of data that are to be interpreted in terms of evolving damage. This is called structural health monitoring [22, 24]. The ability to obtain continuous updates in the health of the structure provides a large collection of data that can be coupled within the framework of DDDAS in order to extract appropriate conclusions regarding structural durability and reliability. The type and evolution of damage that appears in composite materials and structures are complex functions of the composition, architecture, geometry and loading history. Furthermore, such damage occurs at multiple scales; visible delaminations of the macro (or structural) scale, transverse matrix cracks, fiber-matrix interfacial failure, fiber breaks, etc. occur at different length scales and result in the degradation of the overall performance. Numerous types of sensors, such as embedded fiber-Bragg sensors [29], distributed CNTs (see [6, 24]) have been developed and used in evaluating damage. Recently, methods of patterning CNTs into addressable arrays with a spacing of a few micros have also been developed [22]. These distributed sensors can be embedded while processing the composite material or pasted as an add-on at critical points in the structure. In this study, we use CNTs embedded inside the polymer matrix; as shown in [26], a small volume fraction - 0.5 wt.% - of multi-wall CNTs creates a percolating electrically conductivity and can easily be detected. These conductivity measurements can be performed at a local level and form the basis of distributed dynamic data that will be derived to drive the simulations. The report is structured as follows. Section 2 presents some physical models for damage and electrical conductivity, followed by Section 3, where we briefly discuss general Bayesian filtering, as well as Kalman and Extended Kalman filters. In Section 4 we explain the experiments to be soon performed in our project. Section 5 then gives an overview of the integration efforts under way, including the explanation of a simple example problem we consider as a beginning exercise. We conclude the report on Section 6, with a brief description of the work to be done in the next 12 months.

2

2

Physical Example of Damage Identification

As a prototypical problem of composite damage evolution let us consider a sheet of a carbonfiber epoxy composite that is commonly used in aerospace applications; the matrix will be infused with a fine dispersion of multiwall CNTs in order to provide an electrically conducting network. A specimen of width W (either with or without a circular hole of radius a) is subjected to monotonically increasing load up to failure. In the conventional method of experiments the quantity of interest (QoI) is the change in compliance of the specimen, which is to be interpreted as the progression of damage or the reduction in the load-carrying capacity of the structure. In the DDDAS approach, two key measurements are added in this process; first, the initial distribution of microstructure is obtained by spatially mapping the electrical conductivity. Second, the specimen is loaded monotonically to create and extend microstructural damage. At different stages in this damage process, the specimen is to be interrogated in order to determine the spatial variation in the conductivity, and mapped to damage. While the primary QoI is still the overall change in compliance with loading, these measurement techniques yield additional significant dynamic data on the evolution of the microstructure in the form of damage. Therefore, embedded in these measurements is the basis for postulating a damage growth law and one must be able to extract this through techniques of DDDAS. So, how does DDDAS enter into the picture? Two strong links between the dynamic data collection and the dynamic simulations are suggested. The first is a one-way coupling: the damaged material is simulated using continuum damage models. This simulation is driven by the dynamic data on the evolution of the microstructure and strain that is obtained as a function of the load obtained from the experiment. The resulting stress field can be used to identify and predict the next increment of damage. We anticipate that interpretation in terms of the local deformation and stress fields determined through dynamic data-driven simulations of damage can significantly reduce the scatter typically observed in such assessments, and allow tighter confidence intervals. The second link is a two-way coupling: methods of repairing by patching damaged regions [16] or, for self-healing of composites, through the embedding of vascular networks for infiltration of “healing agents” (this has been proposed in [28]).

2.1

Material Model for Damage

We begin with a general description of a damage model and then specialize it to the case of uniaxial stress. Consider a body Ω shown in Figure 1. We assume that the body is linearly elastic, with an elasticity tensor C. In addition, the material may accumulate damage in the form of microcracks, and such damage is represented by a scalar variable D = D(x), x ∈ Ω. The stored energy or work function W is the sum of the elastic energy ψe (, D) and the dissipated energy ψD (D), W (, D) = ψe (, D) + ψD (D), (1) where the elastic energy can be written as 1 ψe (, D) = Cijkl (D)ij ij . 2 3

(2)

Figure 1: Example of a body Ω used in our damage modeling. Damaged regions are shown as brighter subdomains, in which increased microcack densities due to dynamic loading are observed. Note that the linear elastic behavior is represented as Cijkl (D), allowing for the elastic response to be damaged. The stress is obtained as σstress,ij (, D) =

∂ψe = Cijkl (D)kl . ∂ij

(3)

In addition, an appropriate model has to be prescribed for damage, obeying specific restrictions: first, damage is irreversible. Second, the change of stress work with damage must be strictly nonnegative. These result in the inequalities D˙ ≥ 0, ∂ψe dψD ∂ψ = + ≥ 0, −G ≡ ∂D ∂D dD

(4) (5)

where −G is the change in energy due to damage accumulation and is equivalent to the energy release rate in fracture mechanics. The two equations above can be combined into the following constraint: ∂ψe dψD ˙ D = 0. (6) + ∂D dD Numerous models for damage have been developed for different processes such as microcracking, creep, etc. See [4, 5, 21] and others for an overview. For purposes of illustration, we will use the Krajcinovic model [20] applied to a uniaxial problem; we anticipate that this will be appropriate for the microcracking induced damage in the specimens to be tested. For 4

this model, the strain-dependent evolution of damage is given as follows: ( s d if = m and d > 0, c dD = 0 if < m or d < 0,

(7)

where 0 and s are material constants and m is the maximum strain experienced by the material up to the present time. Note that damage is reversible; this implies that whenever ≤ m , D remains constant. The damage attained at any strain level is itself written as s+1 m , (8) D= c where m is the maximum strain achieved in the life of the specimen. It can be shown that ψD =

s+3 D. s+1

For implementation into one-dimensional models, it is easier to write the damage evolution equation directly in terms of the strain. Therefore, if we apply a load F , the strain is ( F , ≤ m , K(m )L (9) = F , > m , K()L and the total elongation in the specimen of length L is ∆ = F/K(m ), where K(m ) denotes stiffness. Note that when damage accumulates during > m , Equation (9)-(b) is implicit and then we must set m = max(). If the specimen is homogeneous, and subjected to inhomogeneous strains, equations (7)-(9) are to be used at each material point to track the local evolution of damage, stiffness, and hence the strain. It should be noted that maximum strain, m , the stiffness, K(m ), and the damage parameter, D, are all equivalent measures of the damage; once we know one, the others can be determined, but we need a diagnostic tool for the determination of one of these quantities. As discussed in Section 4, direct measurement of strain fluctuations, and indirect measurement of resistivity changes can be used to identify damage evolution experimentally. An example of the stress-strain response for the Krajcinovic model is shown in Figure 2. 2.1.1

Inhomogeneous Specimen

To fix ideas, we consider a specimen of composite material modeled as a one-dimensional strip of length L. We consider the specimen to have inhomogeneities and therefore spatial fluctuations in the material properties. We take the specimen to be made of N segments of equal length li = L/N, i = 1, 2, . . . , N . Each segment may have different initial damage Di , and may develop additional damage during loading differently. Following the notation above for damage, strain, and stiffness, we define the arrays D = (D1 , D2 , . . . , DN ); ε = (m1 , m2 , . . . , mN ); K = (K1 , K2 , . . . , KN ).

(10)

¯ ε¯ and K, ¯ and covariance These arrays can be taken to be random variables with means D, matrices PD , Pε , and PK . The strain in each segment needs to be determined by application 5

Figure 2: Idealized 1D stress-strain relation showing loss in stiffness due to material damage. P of the above model, but the total strain is N i=1 i . The following illustrates how to proceed (k) (k) (k) in step k + 1: given D , ε and K , find D(k+1) , ε(k+1) and K (k+1) under an increment of load F (k+1) = F (k) + dF . The material model given above can be used to evaluate the strain variation and damage evolution. For each segment i, use Equation (9) to determine the strain: (k+1) (k) (k+1) F ≤ mi , i K(km )L , (k+1) i = (11) (k+1) (k+1) (k) F (k+1) , i > mi . K(i

)L

(k+1)

Note, again, that Equation (11)-(b) is implicit for i . Damage evolution is also obtained as (k+1) (k+1) mi = max(i ). (12) Next, we suppose that we have the measurement of strain from a digital image correlation (k) that represents an independent evaluation: let us call this measurement ei . It is important to recognize that we may not always have a direct measurement of the strain fluctuations (k) and would need other measurements in such cases. We compare the measurement ei with (k+1) (k+1) the prediction i to get an update for damage in each segment and determine mi and hence ε(k+1) . For accomplishing this comparison and updating the damage variables, we use the Bayesian methods described in Section 3.

6

2.1.2

Indirect Measurement of Damage

Now let us assume that a direct assessment of the local strain is not available or possible; (k) i.e., measurement of ei is not available. Another measurement is needed in order to provide an update to the prediction in Equation (12); as described above, distributed sensors can be embedded while processing the composite material or pasted as an add-on at critical points in the structure. In this project, we use CNTs embedded inside the polymer matrix. Changes in the microstructure by the formation of matrix cracks result in a change in electrical conductivity and can be detected with ease. These conductivity measurements can be performed at a local level. The strategy is to consider the specimen as a resistor network, with varying resistances along the length; as the material damages, the resistivity increases. Then we use an inverse analysis (resistor network/percolation model) to map the resistivity to damage. The details of this strategy are discussed in this subsection. As in the mechanical part, we consider the specimen to have inhomogeneities and therefore fluctuations in the resistance. We take the specimen to be made of the same N segments of equal length li = L/N, i = 1, 2, . . . , N . Each segment may have different initial concentrations and distributions of CNTs and therefore the initial resistance can be written as: R = (R1 , R2 , . . . , RN ).

(13)

What is the relationship between R and the damage, D? This can be obtained either through a calibration experiment or through a model of the CNT network. It is expeditious to begin with an experimental calibration and develop the model in parallel. The experimental task is in progress and is discussed in Section 4. Two types of models for this set up are considered: homogenization models based on micromechanics and statistical models based on percolation. These models are discussed briefly in the following two subsections. 2.1.3

Three-Phase Homogenized Model for Conductivity

The following three-phase model is considered in [15]: let the medium be a mixture of two −1 −1 materials - matrix and inclusion with conductivities σcond,m = Rm and σcond,in = Rin respectively. The inhomogeneous material is modeled as a spherical inclusion of radius a, surrounded by a concentric shell of matrix material with outer radius b, and surrounded fur−1 ther by a medium with the effective conductivity σcond,eff = Reff of the mixture. A potential φ = E0 r cos θ is prescribed at the outer boundary of the body such that a homogeneous electric field E = ∇φ = E0 ey (where ey is the unit vector in the y-direction) is established. Since the medium is heterogeneous, and the current is constant, j = σcond E, E varies spatially, but we require that the energy in the three phase composite spheres be equal to the energy in a homogeneous medium with equivalent conductivity. Noting that the potential in each region must satisfy the Laplace equation, we can write the solution as follows: r ≤ a, φin (r, θ) = Ain r cos θ, 2 φ(r, θ) = φm (r, θ) = (Am r + Bm /r ) cos θ, a ≤ r ≤ b, (14) φeff (r, θ) = E0 r cos θ, r ≥ b. Enforcing the interface conditions φin (a, θ) = φm (a, θ)

and

σcond,in · ∂r φin (a, θ) = σcond,m · ∂r φm (a, θ) 7

at the interface between the inclusion and matrix, as well as φm (b, θ) = φeff (b, θ)

and

σcond,m · ∂r φm (b, θ) = σcond,eff · ∂r φeff (b, θ)

at the interface between the matrix and the effective medium, yields the following four equations for the undetermined parameters Ai , Am , Bm , and E0 : Ain a = Am a +

Bm , a2

σcond,i Ain = σcond,m Am − Am b +

2σcond,m Bm , a3

Bm = E0 b, b2

(15) (16) (17)

and

2σcond,m Bm = σcond,eff E0 . (18) b3 These are four homogeneous equations in the four unknowns and, therefore, for nontrivial solutions, the determinant of the coefficient matrix must vanish. This yields the conductivity of the effective medium: c σcond,eff = σcond,m 1 + , (19) (1 − c)/3 + σcond,m /(σcond,i − σcond,m ) σcond,m Am −

where c = (a/b)3 is the volume fraction of the inclusion in the matrix. This homogenized estimate of the conductivity works quite well when there is a uniform dispersion of the inclusions in the matrix. In the case where there is a dilute concentration of carbon nanotubes in the matrix, this may still be an acceptable approach for the determination of the effective conductivity/resistivity. In fact, this analysis can be used to discretize the initial specimen geometry into rectangular grids as illustrated in Figure 3, determine the spatial distribution of resistivity by comparison to experiments and thereby identify the initial spatial distribution of CNTs in the specimens. However, when the CNTs are broken through applied loading, there is effectively no change in the volume fraction of fibers, but only a change in the connectivity, and hence of the conductivity; therefore, this homogenization method cannot be used to identify conductivity changes arising from damage. In order to determine the microcracking induced changes in resistivity, we consider next a percolation model that permits incorporation of connectivity changes. 2.1.4

Percolation Model for Conductivity

The basic idea of the percolation model is as follows: we consider the CNT filled composite as a resistor network, with infinite conductivity for the CNT and zero conductivity for the resin matrix. So, the effective conductivity is simply due to forming connective paths along the CNTs. This model is ideal in dealing with the situation in which the volume fraction of the inclusion does not change, but the connectivity does, due to breaking of the CNT links through microcracking. The governing equation is ∇ · (σcond (x)∇φ) = 0, 8

(20)

Figure 3: Discretization of a specimen geometry into rectangular grids, in order to determine an initial effective conductivity/resistivity profile. where σcond (x) is the spatial variation of resistivity. We develop a two-dimensional percolation model for capturing the effect of the CNTs. Therefore, the domain of interest, shown in Figure 1, is discretized into a rectangular grid with spacing dx and nodes at xi . For interior nodes, Equation (20) can be written in a point-wise form as X σcond,ij (φi − φj ) = 0, (21) i,j

where σcond,ij = dx [σcond (xi ) + σcond (xj )] /2 is the effective conductivity of the i−j th link and the sum is over all links connected to the node i. In order to represent the CNT infiltrated into an insulating polymer matrix, we shall take the conductivity to be 0 or 1 depending on whether the link is in the polymer or the CNT. This is the basic idea behind the bond percolation model. In the numerical implementation, the effective conductivities of the links, σcond,ij are assigned randomly to 0 or 1, with the appropriate volume fraction of CNTs and the effective conductivity of the grid is calculated as the number fraction of percolating paths from the top to the bottom. Then row average conductivity is calculated to determine the resistance in the N segments corresponding to Figure 3. 9

In step 0, we consider the ith segment of the specimen with an effective conductivity (0) (0) σcond,i = 1/Ri , to be composed of a random network of CNTs that provides the appropriate conductivity; in step 1, we can use the measured resistivity and break a certain fraction of the network to alter the effective conductivity in the model; the number of links broken can be taken to be a measure of damage in the material. This process is continued until final failure of the specimen.

10

3

Bayesian Methods

We are interested in statistically assessing how the state of a system of physical volume Ω ⊂ R3 evolves with time t ∈ [0, +∞), so that we can inform ourselves for potential decisions to be taken about the system, and/or for potential control actions to be taken on the system. The state will be indicated by the vector θ ∈ Rnθ , for some fixed positive integer nθ > 0 (N in Subsection 2.1.1). The initial state of the system is specified by the probability density function (PDF) π(θ (0) ). The (eventual) control will be indicated by the vector u ∈ Rnu , for some fixed positive integer nu > 0. In order to assess the system state, we collect measurement data d(0) , d(1) , d(2) , . . . at instants 0 6 t(0) < t(1) < t(2) < . . ., and then apply Bayes’ formula [17] πpost (θ (k) |d(k) ) =

πlike (d(k) |θ (k) ) · πprior (θ (k) ) . πdata (d(k) )

(22)

Measured data will be indicated by the vector d ∈ Rnd , for some fixed positive integer nd > 0. In (22), πpost (θ|d) is the posterior PDF defining the Bayesian update of the prior information embodied in πprior (θ). The likelihood PDF πlike (d|θ), a PDF of d, encapsulates assumptions about the discrepancy between the values d that are measured and the values that can be computed with the computational model we have at our disposal. And the term Z (k) πdata (d ) = πlike (d(k) |θ (k) ) · πprior (θ (k) ) dθ (k) is the normalization value (for a given d(k) ) that makes (22) a PDF. For the DDDAS of interest here: • Ω = some material piece that might be subjected to stresses and consequent damages; • d(k) = electric potential measured at nd positions of the system at instant t(k) ; and • θ (k) = spatial distribution of damage (a vector of size nθ ) inferred by statistical inversion. Once we infer such spatial distribution, the control decision might be (in case the material is part of an airplane in service, for example): • to instantly apply healing (an example of control) to a damaged region, and/or • to instantly update a flight maneuver plan (another example of control) in order to diminish the possibility that any further damage happens to the system, and/or • update the flight computer so that it has up-to-date information before taking a maneuver decision. However, many uncertainties are involved in the process of assessing the system state: • the data is measured only at some points of the system, • the measured data has noise, 11

• the computational model that maps the spatial distribution of electric potential to the spatial distribution of conductivity does not capture reality perfectly, and • the computational model that maps the spatial distribution of conductivity to the spatial distribution of damage is an imperfect characterization of reality.

3.1

The General Bayesian Filtering

If no control is applied, Bayes’ formula (22) is applied every time new data is collected, in order to update one’s knowledge about the system state. Implicit in this understanding is the equality πprior (θ (k+1) ) = πpost (θ (k) |d(k) ). However, if control is applied, then: (a) [Prediction step] One might need to use an evolution equation in order to assess the consequences of such control in the state of the system. A possible discrete form of such evolution equation is θ (k+1) = f (k+1) (θ (k) , u(k+1) , w(k) ),

(23)

where f (k+1) (·, ·, ·) is an evolution function and w denotes the state noise. Once the new state θ (k+1) is predicted, one can also predict the next measurement to be obtained at t(k+1) . A possible discrete form of such prediction is given by the output equation z(k+1) = g(k+1) (θ (k+1) , v(k+1) ),

(24)

where g(k+1) (·, ·) is an output function and v denotes the output noise; (b) [Correction step] Then, finally, one can actually measure data d(k+1) at t(k+1) . The comparison between z(k+1) and d(k+1) in the likelihood PDF will then allow one to statistically update the predicted state θ (k+1) using Bayes’ formula (22). The steps (a)-(b) continue as long as we want to assess and control the state of our system. Such continuing process can be represented by the equation πpost (θ (k+1) |d(k+1) , d(k) , . . . , d(0) ) =

πlike (d(k+1) , d(k) , . . . , d(0) |θ (k+1) ) · πprior (θ (k+1) ) . πdata (d(k+1) , d(k) , . . . , d(0) )

(25)

If one assumes (i) that the system state follows a first-order Markov process, that is, πstate (θ (k+1) |θ (k) , . . . , θ (0) ) = πstate (θ (k+1) |θ (k) ),

(26)

and (ii) that, given a current state, the measurements are independent of previous measurements, that is, πlike (d(k+1) |d(k) , . . . , d(0) , θ (k+1) ) = πlike (d(k+1) |θ (k+1) ), (27) then (25) can be rewritten as [14, 18, 27] πpost (θ (k+1) |d(k+1) , d(k) , . . . , d(0) ) =

πlike (d(k+1) |θ (k+1) ) · πprior (θ (k+1) |d(k) , . . . , d(0) ) . πdata (d(k+1) |d(k) , . . . , d(0) ) 12

(28)

In (28), the term πprior (θ

(k+1)

(k)

(0)

|d , . . . , d ) =

Z

πstate (θ (k+1) |θ (k) , d(k) , . . . , d(0) )·πpost (θ (k) |d(k) , . . . , d(0) ) dθ (k)

involves the evolution Equation (23), while the term πlike (d(k+1) |θ (k+1) ) involves the output Equation (24). The recursive nature of (28) is clear, in the sense that one makes the transition πpost (θ (k) |d(k) , . . . , d(0) )

πpost (θ (k+1) |d(k+1) , d(k) , . . . , d(0) ).

−→

(29)

The PDF transition (29) can be broken into two intermediate transitions, namely the prediction and correction steps mentioned above. The prediction step relates itself to the PDF transition πpost (θ (k) |d(k) , . . . , d(0) ) −→ πprior (θ (k+1) |d(k) , . . . , d(0) ), (30) while the correction step relates itself to the PDF transition πprior (θ (k+1) |d(k) , . . . , d(0) )

πpost (θ (k+1) |d(k+1) , d(k) , . . . , d(0) ).

−→

(31)

The Bayesian filtering procedure (28) is usually intractable computationally. However, it can be simplified with further assumptions, as discussed in the next two subsections. Before proceeding, though, we list two extra important notations. Given a PDF π(θ), we will denote ˆ = arg max π(θ). θ

(32)

θ

Also, we will denote by N (µ, C) a (multivariate) Gaussian distribution of mean µ and (co)variance (matrix) C. It should be noted that hypotheses (26) and (27) are always assumed to be valid throughout this report.

3.2

The Kalman Filter

ˆ (0) , as well Let us suppose the following five assumptions hold ∀ k > 0, for known vector θ as known matrices A(k+1) , B(k+1) , H(k+1) , P(0) , Q(k) , and R(k+1) : f (k+1) (θ (k) , u(k+1) , w(k) ) = A(k+1) · θ (k) + B(k+1) · u(k+1) + w(k) , g(k+1) (θ (k+1) , v(k+1) ) = H(k+1) · θ (k+1) + v(k+1) , ˆ (0) , P(0) ), θ (0) ∼ N (θ w(k) ∼ N (0, Q(k) ), and v(k+1) ∼ N (0, R(k+1) ).

13

(33) (34) (35) (36) (37)

Then the general process (28) can be substituted by the following five steps, assuming that ˆ (k) , P(k) and d(k+1) are known: θ ˜ (k+1) = A(k+1) · θ ˆ (k) + B(k+1) · u(k+1) , θ ˜ (k+1) = A(k+1) · P(k) · A(k+1) T + Q(k) , P −1 (k+1) (k+1) (k+1) T (k+1) ˜ (k+1) (k+1) T (k+1) ˜ K = P ·H · H ·P ·H +R , ˆ (k+1) = θ ˜ (k+1) + K(k+1) · d(k+1) − H(k+1) · θ ˜ (k+1) , and θ (k+1) ˜ ˆ (k+1) = I − K(k+1) · H(k+1) · P . P

(38) (39) (40) (41) (42)

Under the Kalman filter assumptions, the system state will always be Gaussian. So, at each time step it suffices to determine the mean (41) and the covariance (42). Substeps (38)-(39) correspond to the prediction step (30), while substeps (40)-(42) correspond to the correction step (31). The matrix K in (40) is called the Kalman gain.

3.3

The Extended Kalman Filter

In case the system is not linear nor Gaussian, one possible approach is to linearize both the evolution function f (k+1) (·, ·, ·) in (23) and the output function g(k+1) (·, ·) in (24). The general process (28) can then be substituted by the following five steps: ˜ (k+1) = f (k+1) (θ ˆ (k) , u(k+1) , 0), θ (43) T T ˜ (k+1) = A(k+1) · P(k) · A(k+1) + W(k+1) · Q(k) · W(k+1) , P (44) −1 ˜ (k+1) · H(k+1) T · H(k+1) · P ˜ (k+1) · H(k+1) T + V(k+1) · R(k+1) · V(k+1) T K(k+1) = P (45) , ˆ (k+1) = θ ˜ (k+1) + K(k+1) · d(k+1) − g(k+1) (θ ˜ (k+1) , 0) , and θ (46) (k+1) ˆ (k+1) = I − K(k+1) · H(k+1) · P ˜ P , (47) where

(k+1)

(k+1)

∂fi ˆ (k+1) , u(k+1) , 0), (θ ∂θj

Aij

=

(k+1) Wij

∂f ˆ (k+1) , u(k+1) , 0), = i (θ ∂wj

(48)

(k+1)

(49)

(k+1)

(k+1) Hij

∂gi ˜ (k+1) , 0), (θ = ∂θj

and

(50)

(k+1)

(k+1)

Vij

=

∂gi ˜ (k+1) , 0). (θ ∂vj

(51)

ˆ (0) , P(0) , Q(k) , and R(k+1) , ∀ k > 0. Above, obviously, it is assumed that one knows θ We intend to apply the Kalman filtering, as well as the extended Kalman filtering, to the models listed in Section 2, using the parallel MPI/C++ software libraries libMesh [19] and QUESO [23]. 14

4

Experiments

The experimental plan for the DDDAS project is as follows: monotonically increasing uniaxial tensile loading will be applied on epoxy specimens infiltrated with carbon nanotubes. The weight percent of CNTs will be varied in the range from 0.1wt% to 0.9wt% in steps of 0.2wt%. This range of CNTs is the percolation threshold range for increasing conductivity and therefore should provide adequate sensitivity in the experimental measurement of resistivity. The diagnostic methods include two independent measurements. First, the spatial variation of the strain (x) will be measured using Digital Image Correlation (DIC) method [25]. In this method, the specimen is initially decorated with a random speckle pattern. This pattern is photographed at high spatial resolution at different stages of the loading. By correlating the changes in the speckle pattern between the images at different stages, the strain field over the two dimensional region of interest can be determined. This will then be averaged over rectangular regions indicated in Figure 1 and row averaged to obtain the variation (x). The fluctuations in strain that are obtained through this set of measurements will correspond to determining damage ε(k) at every load step k. Second, the specimens will be coated with conducting paint lines at discrete intervals as indicated in Figure 3 to facilitate the measurement of the resistivity R. A potential φ is applied between the ends of the specimen and the voltage drops across eight different lines are measured simultaneously, and at different overall strain levels. The potential across the different x locations is directly related to the resistivity and therefore this measurement yields R. The hardware and software necessary for acquisition of this resistivity data have been assembled, and the system is ready to be used. The main technical challenge remaining is in the preparation of specimens with a uniform dispersion of the CNTs within the epoxy without agglomeration. This will be accomplished by the manufacturer, Designed Nanotubes, Austin, TX, who specializes in functionalizing carbon nanotubes with special properties and dispersing them in the matrix. We expect to have the material for characterization by the end of September. A specialized mold has been designed and built for the purpose of manufacturing uniaxial tensile test samples that are 3 mm thick, 10 mm wide, and 150 mm long. The uniaxial tests will be performed in an Instron test machine with multiple diagnostics to enable determination of the stress-strain response as well as damage development. Therefore, the experiments will provide a direct measurement of the fluctuations in both the strain field and the resistivity; with this data, it is possible to calibrate the damage model first and to establish the connection between the resistivity measurements and the damage. At any load step k, the following experimental data are to be generated: the load Q(k) , the (k) (k) strain variation ei , and the resistivity Ri , i = 1, . . . , N . These are to be passed on to the Bayesian analysis to obtain predictions for the strain and resistivity during the next step.

15

5

Integration

To implement the DDDAS system, the mathematical model for composite damage, the Kalman filter algorithm, and the experimental data generated must be integrated into a cohesive infrastructure. Much of this infrastructure is, and will be, expressed in software. We detail both the software that we generate as well as existing software and tools that we leverage to accomplish building the DDDAS infrastructure.

5.1

Using the Generated Experimental Data

As discussed in Section 4, two primary sets of experimental data are key to the DDDAS under investigation. The first is the spatial strain field at each stage of loading. The strain measurements are gathered infrequently due to the level of involvement in processing the experimental data to obtain strain fields. Once collected, these data are used to construct a model for inferring resistivity in the composite specimen based upon an associated level of damage. The second set of data are the measurements of voltage. As described in Section 4, these data are collected at various lengths along the specimen. These are an integral part of the DDDAS as we use these measurements to compare against predicted measured values of voltage. The error between the predicted and measured values of voltage are then used to perform the Bayesian update (31).

5.2

Damage Model Implementation

We leverage the libMesh library [19] for a parallel, C++ implementation of the damage model described in Section 2. As the project progresses, this will enable potential development of higher dimensional models. Through libMesh, we will also leverage well established linear and nonlinear solver packages, such as PETSc [2, 1, 3], to solve the resulting systems of equations that arise in the damage model.

5.3

Kalman Filter Integration

With the experimental data generated and a standalone application for predicting damage implemented, the final component to be integrated is the implementation of the Kalman Filter algorithm discussed in Section 3.2. There are two primary components that are required by the filtering algorithm: the prediction of the system state θ (k+1) and the prediction of the measurement z (k+1) . In the system under study, θ (k+1) corresponds to the predicted state of damage while z (k+1) is the predicted voltage, based on the predicted state of damage. To elaborate, given the loading at state k + 1 of the composite specimen, the damage model is used to predict the material damage at k + 1, θ (k+1) . This step corresponds to (23) in the filtering problem and the damage model corresponds to the function f in the same equation. The resulting damage prediction for state k + 1 is then used to predict the measurement of the voltage, as described in Section 2, at state k +1. In the case at hand, this is actually a two step procedure. The first step uses the empirically determined relationship between damage and resistivity to determine resistivity in the specimen. Then, using Ohm’s 16

law, the voltage field can be computed using the known applied voltages and the inferred resistivity field. Both the predicted damage and the predicted voltage at step k + 1 constitute the prior information for the Bayesian update, when the posterior for state k + 1 is then computed. This cycle continues until a critical damage state is detected, at which point the process is ordered to stop. Note that, in the current work, the Kalman Filter cycle does not interact directly with the experimental system. Instead, the data is accumulated and the specimen is loaded to a critical state. The data is then used by the Kalman Filter procedure as though it were being gathered in real time to provide proof-of-concept of the ability of this system to prevent failure of the material.

5.4

A Simple Filtering Problem

The example in this Subsection 5.4 is for clarifying the main concepts only, so the team can keep a good communication of ideas, as well as a good level of mutual understanding of the computational steps involved. The example is not aimed to be physically realistic: the models and equations used here fit only the main purpose of exercising the computational flow of a recursive filtering process. The domain is Ω = (0, 1) ⊂ R. The material is assumed to be homogeneous. Points x = 0 and x = 1 are constantly subjected to electric potentials of 0 Volts and 5 Volts, respectively. An external force F : Ω × [0, +∞) → R acts on the material. It is given by πt 1 . (52) F (x, t) = · x · (1 − x) · sin 2 16 The electric potential is measured at the 49 points x = 0.02, 0.04, . . ., 0.98, shown in Figure 4. For computational purposes, the physical domain is divided into 50 elements, each associated with a locally constant resistivity Ri , 1 6 i 6 50. The map between damage Di and conductivity σcond,i at the i-th element is assumed to be σcond,i =

1 = (1 − Di ), Ri

1 6 i 6 50.

(53)

The timely evolution of damage at the i-th is assumed to be governed by the equation (the equivalent to the evolution Equation (23)) |F (xi , t)| D˙ i = , 5

1 6 i 6 50.

(54)

where xi is the middle point of the i-th element. 5.4.1

The Forward Problem Perspective

As the external force is applied, we simulate how the damage, the resistance, and the voltage spatial distributions evolve with time, as shown by Figures 5, 6, 7 and 8, for instants t = 0, (k) t = 9, t = 59, and t = 79, respectively. At each time step t(k) , the damage Di is computed

17

0V

5V Φ

Φ

0.02

Φ

48

R 2

R 1 0

Φ

2

1

49

R 49 0.04

0.96

R 50 0.98

1

x

Figure 4: Schematic representation of the simple problem of Subsection 5.4. (k)

with Equation (54), and the corresponding resistance Ri (53). Defining φ0 = 0, φ50 = 5, and

is then computed with Equation

φ50 − φ0 I (k) = P50 (k) , i=1 Ri one can then compute the corresponding spatial distribution of electric potential through the equation (k) (k) (k) φi = φi+1 − Ri+1 · I (k) , 1 6 i 6 49. Each figure shows the external force (52), the computed spatial distribution of damage, the computed spatial distribution of resistance, and the computed spatial distribution of electric potential. From Figure 5 (t = 0) until Figure 8 (t = 79), it is interesting to observe that the damage gradually increases, up to the point that complete damage (D = 1) happens in the middle of the domain, where the force F (x, t) usually has higher magnitude (and the magnitude is the main driver of damage evolution according to Equation (54)). As intuitively expected, domain regions with higher damage experience more difficulty to support the flow of electric current due to microcracks that accumulate. That is, domain regions with higher damage show higher resistance. Indeed, resistance exhibits higher values in the middle portion of the domain, where the damage is usually bigger. When the damage achieves value 1, the resistance goes to +∞ (according to Equation (53)), electric current ceases to flow, and the potential abruptly decreases from 5V on the right to 0V on the left, as shown in Figure 8. 5.4.2

The Bayesian Filtering Perspective

As the system evolves, we monitor it by first measuring the spatial distribution of voltage from time to time. After each measurement, we infer for the resistance spatial distribution, using Ohm’s law, and then infer for the damage spatial distribution, using Equation (53). The measurements obtained and models used are uncertain, and the inferences are statistical, 18

following Bayes’ formula (22). Once the damage is inferred, one might decide to take control actions: for instance, heal a certain region, or diminish the load.

19

Figure 5: External force and spatial distribution of system state (damage, resistivity, and φ = electric potential) at t = 0.

20

Figure 6: External force and spatial distribution of system state (damage, resistivity, and φ = electric potential) at t = 9.

21

Figure 7: External force and spatial distribution of system state (damage, resistivity, and φ = electric potential) at t = 59.

22

Figure 8: External force and spatial distribution of system state (damage, resistivity, and φ = electric potential) at t = 79.

23

6

Future Work

We will develop and implement the models for predicting the damage state in a onedimensional composite, given a loading state, and for determining resistivity, given the (just predicted) damage state. These models will replace the simple models used in Subsection 5.4, and will be implemented on MPI/C++ computational environments. In parallel to the development and implementation of these models, we will also complete the numerical experiments of Subsection 5.4.2. Regarding experimental data, specimens are currently being manufactured. Once they are obtained, data collection will begin in our laboratory. With the models and actual data at hand, we will apply the Kalman and the extended Kalman filtering algorithms of Section 3 in order to statistically assess the dynamic evolution of damage in the laboratory specimens.

24

References [1] S. Balay, J. Brown, , K. Buschelman, Victor Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang. PETSc users manual. Technical Report ANL-95/11 - Revision 3.3, Argonne National Laboratory, 2012. [2] S. Balay, J. Brown, K. Buschelman, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang. PETSc Web page, 2012. http://www.mcs.anl. gov/petsc. [3] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith. Efficient management of parallelism in object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, Modern Software Tools in Scientific Computing, pages 163–202. Birkh¨auser Press, 1997. [4] J. L. Chaboche. Continuum damage mechanics: Part I - General concepts. Journal of Applied Mechanics, 55:59–64, 1988. [5] J. L. Chaboche. Continuum damage mechanics: Part II - Damage growth, crack initiation and crack growth. Journal of Applied Mechanics, 55:65–72, 1988. [6] T. W. Chou, L. Gao, E. T. Thostenson, Z. Zhang, and J. H. Byun. An assessment of the science and technology of carbon nanotube-based fibers and composites. Composites Science and Technology, 70:1–19, 2010. [7] F. Darema. DDDAS: Dynamic Data Driven Applications Systems. http://www.nsf. gov/cise/cns/dddas. [8] F. Darema. Dynamic Data Driven Applications Systems: A New Paradigm for Application Simulations and Measurements, volume 3038 of Lecture Notes in Computer Science. Springer, 2004. [9] F. Darema. Dynamic data driven applications systems: A new paradigm for application simulations and measurements. In Marian Bubak, Geert Dick van Albada, Peter M. A. Sloot, and Jack J. Dongarra, editors, Computational Science - ICCS 2004, volume 3038 of Lecture Notes in Computer Science, pages 662–669. Springer Berlin / Heidelberg, 2004. [10] F. Darema. Grid computing and beyond: The context of dynamic data driven applications systems. Proceedings of the IEEE, 93(3):692–697, 2005. [11] F. Darema. Characterizing dynamic data driven applications systems (dddas) in terms of a computational model. In Gabrielle Allen, Jaroslaw Nabrzyski, Edward Seidel, Geert van Albada, Jack Dongarra, and Peter Sloot, editors, Computational Science ICCS 2009, volume 5545 of Lecture Notes in Computer Science, pages 447–448. Springer Berlin / Heidelberg, 2009.

25

[12] F. Darema and M. Rotea. Dynamic data-driven applications systems. In Proceedings of the 2006 ACM/IEEE conference on Supercomputing, SC ’06, New York, NY, USA, 2006. ACM. [13] F. Darema and H. E. Seidel. Report of the August 2010 Multi-Agency Workshop on InfoSymbiotic/DDDAS. The Power of Dynamic Data Driven Applications Systems, 2011. [14] A. M. Fraser. Hidden Markov Models and Dynamical Systems. SIAM, 2008. [15] Z. Hashin and S. Shtrikman. A variational approach to the theory of the effective magnetic permeability of multiphase materials. Journal of Applied Physics, 33:3125– 3131, 1962. [16] Hexcel Corporation. Composite Repair. technology-manuals, 2011.

http://www.hexcel.com/Resources/

[17] E. T. Jaynes. Probability Theory: The Logic of Science. Cambridge University Press, 2003. [18] A. H. Jazwinski. Stochastic Processes and Filtering Theory. Dover, 1998. [19] B. S. Kirk, J. W. Peterson, R. H. Stogner, and G. F. Carey. libMesh: a C++ library for parallel adaptive mesh refinement/coarsening simulations. Engineering with Computers, 22(3):237–254, 2006. [20] D. Krajcinovic and G. U. Fonseka. The continuous damage theory of brittle materials, Parts 1 and 2. Journal of Applied Mechanics, 48:809–824, 1981. [21] J. Lemaitre and R. Desmorat. Engineering damage mechanics. Springer, 2005. [22] NASA. Carbon Nanotube-Based Structural Health Monitoring Sensors, 2011. NASA Tech Briefs, LAR-16475-1. [23] E. E. Prudencio and K. W. Schulz. The parallel C++ statistical library QUESO: Quantification of Uncertainty for Estimation, Simulation and Optimization. In M. Alexander et al., editor, Euro-Par 2011 Workshops, Part I, volume 7155 of Lecture Notes in Computer Science, pages 398–407. Springer-Verlag, Berlin Heidelberg, 2012. [24] J. Rausch and E. M¨ader. Health monitoring in continuous glass fibre reinforced thermoplastics: Tailored sensitivity and cyclic loading of CNT-based interphase sensors. Composites Science and Technology, 70:2023–2030, 2010. [25] M. A. Sutton, W. J. Wolters, W. H. Peters, W. F. Ranson, and S. R. McNeill. Determination of displacements using an improved digital correlation method. Image and Vision Computing, 1:133–139, 1983. [26] E. T. Thostenson and T. W. Chou. Carbon nanotube networks: sensing of distributed strain and damage for life prediction and self healing. Adv. Mater., 18:2837–2841, 2006. 26

[27] G. Welch and G. Bishop. An Introduction to the Kalman Filter. http://www.cs.unc. edu/~welch/kalman/kalmanIntro.html, 2006. UNC-Chapel Hill TR 95-041. [28] S. R. White, N. R. Sottos, P.H. Geubelle, J. S. Moore, M. R. Kessler, S. R. Sriram, E. N. Brown, and S. Viswanathan. Autonomic healing of polymer composites. Nature, 409:794–797, 2001. [29] G. Zhou and L. M. Sim. Damage detection and assessment in fibre-reinforced composite structures with embedded fibre optic sensors review. Smart Materials and Structrures, 11:925–939, 2002.

27