An inverse model and mathematical solution for inferring viscoelastic properties and dynamic deformations of heterogeneous structures

Accepted Manuscript An inverse model and mathematical solution for inferring viscoelastic properties and dynamic deformations of heterogeneous structures Qinwu Xu, Hongyu Zhu PII: DOI: Reference:

S0045-7825(15)00365-5 http://dx.doi.org/10.1016/j.cma.2015.11.012 CMA 10756

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date: 22 January 2015 Revised date: 7 September 2015 Accepted date: 6 November 2015 Please cite this article as: Q. Xu, H. Zhu, An inverse model and mathematical solution for inferring viscoelastic properties and dynamic deformations of heterogeneous structures, Comput. Methods Appl. Mech. Engrg. (2015), http://dx.doi.org/10.1016/j.cma.2015.11.012 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An inverse model and mathematical solution for inferring viscoelastic properties and dynamic deformations of heterogeneous structures Qinwu Xua,b,*, Hongyu Zhub a

Civil, Architectural and Environmental Engineering; b Institute for computational Engineering and Sciences, the University of Texas at Austin, Austin, TX 78712 USA Email: [email protected], [email protected]; *Phone: +1512-232-1701; Fax: +1 (512) 475-7314

ABSTRACT We proposed an adjoint-based inverse method and mathematical solutions using Lagrange multiplier theorem for inferring dynamic viscoelastic properties of heterogeneous structures to improve computation reliability and efficiency. Existing methods such as the one using finite-difference approximated gradient may not be efficient and accurate enough for considering complex situations of the coupled effects of dynamic loading, material deformation memory, and structural heterogeneity. The proposed method derives both gradient and Hessian mathematically to satisfy the first-order necessary and second-order sufficient optimal conditions, which has great advantages compared to other approaches. We also proposed robust numerical algorithms for accurate and fast computations, including a regularization to control reasonable parameter ranges, a Newton’s method with Hessian function to determine search direction for fast convenience, and a modified Armijo rule to find a stable step length. We developed a Galerkin time-domain finite-element method for numerical solutions of resulting partial differential equations, and validated the model for a layered structure. Results indicate that the proposed method has greatly improved computation accuracy and speed as compared to existing approaches. It may be applied to a broad range of heterogeneous materials and structures at different length and time scales. Key words: Viscoelastic material, optimization, finite elements, dynamics, structures Nomenclature u: displacement : displacement discretized at nodes : variation of : incremental variation of : test function or Lagrangian multiplier : joint discretized at nodes : differential of : Incremental variation of : material model parameter : material parameter vector discretized at nodes : pre-defined material model parameter : variation of : incremental variation of : relaxation modulus : dynamic modulus : elastic modulus of springs of Maxwell model : viscosity of dashpots of Maxwell model

: Young’s modulus : relaxation modulus matrix : lagrangian g: gradient, H: hessian, H: hessian discretized at nodes : loading b: body force : stress tensor M: mass matrix; C: damping matrix : shape function matrix : strain-displacement matrix, : gradient, : divergence, : for all

1

1. Introduction Dynamic modulus

equals to the ratio of stress over strain under vibratory conditions (i.e.,

).

is a characteristic of material viscoelastic (VE) property, which can be decomposed into the storage modulus (the real or elastic part) plus the loss modulus (the imaginary or viscous part) as . For viscoelastic materials, stress and displacement are dependent on time and temperature [1][2]. Heterogeneous structures and material systems (multiple components or layers with different material properties) displaying the VE property exist widely in the contexts of engineering and human life at different length scales, including human tissues of skin [3] and arterial walls [4], memory foam, and polymer structures [5] as shown in Figure 1. Material viscoelastic property plays a crucial role in a subject’s responses and performance under loading effects. For example, viscoelastic plaques cause blockages of arterial lumen [6] due to the viscoelastic flow. Proper viscoelastic properties of memory foam mattress offer optimal comfort for sleeps. Asphalt concrete (AC) is one of the mostly used construction materials and display typical VE property. AC pavement-a multilayered heterogeneous structure built on soil foundation, typically exhibit rutting (permanent displacement) [7] and fatigue cracking performance under repeated loading, which is dependent on the AC material’s VE property [1]. In order to safely perform material and structural design of the heterogeneous structures and material systems, a fundamental understanding of the VE properties and dynamic deformation of materials is important. a )

b ) Skin Dermis

Subcutaneous fat

Normal artery

c )

d )

Plaque

Blocked artery

Figure 1. Heterogeneous viscoelastic media at different length scale: a) skins [3]; b) heart’s artery wall [9]; c) memory foam mattress (www.suppemattre.net); d) layered polymer molecular structure [5] all with permissions. Inverse computations using the response data collected from nondestructive tests can predict or estimate material properties—without directly measuring them. More often these inverse computations generally involve two stages for the inverse computation of material properties: 1) simulating structural responses (e.g., deflections) under a load source to emulate nondestructive tests; and then 2) minimizing the error function of the differences between response simulations and measurements—or observations—using a mathematical optimization algorithm to estimate appropriate material properties. The second stage uses the response outputs from 1) for numerical approximation such as calculating gradients with finite difference method (FDM) without interactions on the partial differential equation (PDE) of the first stage of response modeling. The neural network approach has also been proposed, in which it establishes a large database at stage 1), and then inverts model parameters in stage 2) based on the database using the system identification or pattern search method [8]. Thus, for stage 1), a separate product such as software may be used to perform the forward modeling and provide response results. This is referred as a two-stage approach in this study to differentiate the proposed approach discussed later. Inverse computations of material viscoelasticity and deformations have been studied in engineering and science including biology and the polymer sciences. Catheline et al. [10] used the wave-propagation-based experimental

2

data with an inverse method to estimate VE properties of a membrane. Brigham [11] used the SMARS approach (a combination of a classical random search algorithm plus the neural networks-based method) to invert the VE material properties of solids immersed in fluids. Lei and Szeri [12] and Zhao et al. [13] developed a two-stage approach to invert the material hyperelastic and VE properties of biomaterials, where the forward responses were computed using ABAQUS and the inverse computation was based on the Levenberg-Marquardt algorithm embedded in MATLAB. Araújo et al. [14] used a gradient-based optimization technique with the Gauss-Newton algorithm (presented by [15]) and a minimization error function (based on the difference between experimental vibration data and finite element (FE) modeling results) to invert the linear VE properties of a sandwich structure. Sims et al. [16] used an inverse FE method to determine the VE properties of subcutaneous fat by combining the use of ANSYS software for forward computations and MATLAB program for the inverse computation with the nonlinear least square algorithm. Giavazzi et al. [17] inverted the VE parameters of skins with the FSQP (Feasible Sequential Quadratic Programming) algorithm coupling with a FE modeling of response (method presented by [18]). Araújo et al. [19] studied a constrained minimization problem, in which a gradient-based optimization technique was employed to invert five parameters of a VE fractional derivative model. Magnuson et al. [20] backcalculated creep compliance of AC material using a mathematical power function with three model parameters. Scarpas and Blaauwendraad [21] inverted the complex modulus of AC materials using the fourparameter Burgers model. Levenberg [22] used the min-max optimum approach for inverse computation of a twoparameter VE model—elastic modulus combined with viscosity for the pavement structure. Varma et al. [23] presented a genetic algorithm for inverse analysis of VE properties of AC material. Wang and Lytton [8] accounted for the dynamic loading effects in inverting materials’ elastic moduli. Liang and Zhu [24] used the dynamic analysis to invert the creep compliance and fatigue parameters of AC materials. Base on this study’s thorough literature review, the research problems for inverting viscoelastic material properties are summarized as follows: i) Typically, a two-stage approach—response forward modeling and then inverse computation—is used for inverse computation of material properties, including those discussed above (e.g. the FE forward modeling + FSQP inverse computation combined method [17]). For the stage of inverse computation, different approaches have been used including the neural network method and Newton method. For the gradient based method, the gradient of an objective function with respect to the model parameters (e.g. modulus) is often approximated using FDM [25] as follows: (1) where

is the simulated deflection response and

is the variation of model parameter value.

In this method, is approximated based on the modeled response results of stage 1. For this method each gradient requires one forward response modeling, which is much more computationally expensive than the method to be developed in this research. A modified newton-Raphson method was developed to compute gradients of all model parameters in one forward response modeling [25], but it is primarily used for relatively simple inverse problems with a small number of model parameters. For Neural network method, it builds a large database to include response results at variable model parameters for inverse search purpose. However, for the complex problem it is very difficult to build such a large database to consider numerous parameter variations and combinations of a number of model parameters. ii) The inverse computation of dynamic VE properties of heterogeneous structures involves complexities due to the coupled effects of dynamic feature, time dependency, memory effect of material displacement, and structural heterogeneity. Therefore, Existing studies have often used simplified methods such including the simplified material model and numerical inverse procedures. These models include at least the power law, sigmoidal

3

function and the four-parameter Burgers model. The mathematical power law model and sigmoidal function may offer a good fit to the experimental data, but they lack details to describe the physical mechanism. The Burgers model—in a relatively simple format with only a few model parameters—is unable to represent the full creep compliance and dynamic modulus master curve at a wide range of frequency [1]. In addition, the complex model conditions may be simplified in the inverse computation. For example, although modeling of the heterogeneous structure under applied loading is a dynamic VE problem, either the material nonlinear behavior or the dynamic loading effect is approximated using the linear elastic approach or in less accurate ways including examples discussed early. These simplifications may not fully capture the accurate material behaviors and properties. iii) With an increased number of model parameters, the process of computation becomes more numerically expensive when using two-stage approaches. The motivation of this research is to develop an inverse model and mathematical solution for imaging VE properties and dynamic deformation of heterogeneous materials and structures, thereby improving the reliability and efficiency of existing methods. This study proposed an adjoint-based method and applied the Lagrange theorem-based mathematical framework to a new problem. This study has formulated the inverse problem to infer the dynamic viscoelastic material properties of heterogeneous structures from the displacement observations. The main mathematical contribution is to solve the dynamic viscoelastic material inverse problem of heterogeneous structures within the Lagrange theorem-based framework. Moreover, the resulting parameters obtained from the inversion have been compared with some experimental data, which, to some extent has validated this work.

2. Proposed Inverse Method and Procedures The proposed method is adjoint based method. The unknown material model parameters are connected to measured responses of the heterogeneous structure through the partial differential equations (PDEs). The solution of inverse problem through the solution of PDEs is becoming increasingly feasible with more powerful computers available. We have extended the Lagrange theorem-based mathematical framework to solve a dynamic viscoelastic material problem within heterogeneous structures. The challenging of the solution lies in the time dependency, dynamics, material deformation memory effect, and structural heterogeneity.

2.1. Inverse problem and model domain The VE material was modeled by the generalized Maxwell (GM) model, one of the most popular VE models, made up of parallel spring and dashpot series as plotted in Figure 2a. The model domain includes both the VE and elastic materials in which each material is projected to a layer. The stress of the dashpot ( ) is dependent on the strain rate as ( is strain). The relaxation modulus , relaxation shear and bulk moduli and for the VE materials in Prony series are expressed as follows [26]: (2a)

(2b) (2c) where , , spring element and the viscosity of the

is the modulus at the infinite time of the first spring; is the modulus of the is the number of terms of the spring-dashpot series (e.g. N is 7 or greater) and is dashpot element.

4

a) VE lay er

Input: material model

b) Forward modeling Output: responses

E1

Elastic layer Output: material properties

Input: responses

Inverse computation Figure 2. Forward and inverse models: a) material model and structure under loading and b) model domain and its resposnes. As shown in Figure 2 the material model parameters are input/output and the measurement of responses is output/input accordingly for the forward/inverse problem. The material properties are considered either VE using the GM model, or elastic with modulus value of . No damages such as plastic deformation and cracking are considered in this study. It is reasonable to assume that no damages occur within a short time period under a normal loading in the nondestructive testing for healthy monitoring purpose. The inverse problem can, therefore, be formed as follows: 1) the given information includes the input loading source, measurements of response pulses, and the geometry and boundary conditions of the model domain (see Figure 2b) and 2) the unknown information are the model parameters of each material ( and , and ) within the heterogeneous structure to be inverted (see Figure 2a). According to the stress equilibrium, the strong form of the governing state equation of the model domain can be formed as follows: (3) where is the body force; is the damping of structural system and materials; stress tensor; is displacement response at time ; is a space domain and According to the material viscoelasticity, the stress tensor

is material density; is a time period.

is

can be expressed as follows [26]: (4)

where is the deviatoric strain tensor; identity matrix.

is the hydrostatic strain scalar at time vairable

and is the

It can be rewritten in a reduced formula as follows: (5)

5

where

is the relaxation modulus that can be express as follows [30]: (6)

The governing state equation is subjected to the natural boundary condition or stress equilibrium at the loading surface as follows: (7) where

is the loading area;

is normal direction and

is the loading (see Figure 2).

The essential boundary condition satisfies the following: (8) where

is the constrained displacement on the boundary area

for

with M boundaries.

The response measurements can be defined as follows: (9) where is the displacement response measurements (or called observations) on the loading surface area (see Figure 2b). Following the Galerkin method, by applying a test function on both sides of the governing state equation and Green’s equation, and then integrating with time and space domains, the weak form of the governing state equation can be formed as follows:

(10) where

is the test function.

2.2. Langrangian Formulation with Optimal Conditions Inverse computation is an iteration procedure to find the optimal material properties that minimize the differences between response observations and the modeling results. The objective function (a scalar) is defined as follows: (11) where is the simulated response (e.g., displacement) and (see Figure 2).

is the observation at the surface area of

Multiple solutions may exist for the inverse computation [25], and thus, here a regularization term is proposed to penalize or limit the values and range of inverted material properties as follows:

6

(12) where is a regularization parameter; is the material model parameter; is the preferred material model parameter range for penalization that and and are the lower and upper bounds of material property values that may be determined from either engineering knowledge or experimentally based database. With space independent material properties, the formula can be rewritten as follows (e.g. for the validation example discussed later in this study): (13) Different material properties may give almost the same ―observations‖ (the problem is ill-posed), the regularization can help pick preferable material properties, e.g., would prefer results close to m*. The cost function is formed as a sum of the objective function and the regularization term as follows: (14) The cost function is subjected to the governing state equation as a PDE constraint. Based on the Lagrangian multiplier theorem, the Lagrangian can be formed as a sum of the cost function and the governing state equation weak form, which is used to find the stationary point for inverse computing material model parameter , as follows:

(15) where (displacement), the Lagrangian.

(test function or Lagrangian multiplier), and

are three ―independent‖ variables of

To find material model parameter using the Lagrangian framework, the optimality conditions are satisfied for the inverse procedure. Here, both the first order necessary and second order sufficient conditions were imposed for optimality. The first order necessary condition is also called KKT (Karush–Kuhn–Tucker) condition for nonlinear programming, which states: suppose that is a local minimizer and is continuously differentiable in an open neighborhood of , then the gradient is zero ( ) [27]. The second order sufficient condition states: suppose that the Hessian matrix ( ) is continuous in an open neighborhood of ; if and is positive definite then is a strict local minimizer of [27]. The first order necessary condition finds the stationary point, including those local stationary points where the tangents are ideally zero. The second order sufficient condition is to assure that the stationary point is the strict

7

local minimum. In contrast, most existing methods for inverting material properties primarily only calculates the gradient to satisfy the first order necessary condition [25].

2.3. Proposed Inverse Procedure The inverse computation is an evolving procedure in which the material model parameters are updated at each iteration step until reaching an optimal point at which the objective function is small enough. Figure 3 shows the proposed flow chart of the numerical inverse computation procedures. The inputs include the model geometry, loading sources, boundary conditions, and initial guessed material model seeds. With the inputs the first step is to calculate responses of displacements. Applying the first order variation of Lagrangian with respect to , the simulated response of displacement can be determined from the following equilibrium (the same as the governing state equation):

(16) where

is the variation of .

Different iteration stop criteria have been applied for inverse computation. RMS (root mean square) is a statistical measure that has often been used as an iteration stop criterion for inverse problems in engineering [20][25]. RMS can be expressed in a norm form as follows:

(17)

where is the specific simulated response values at the is the observation or measured response value.

time step (total

) and

location (total ), and

Other techniques included the gradient norm based method has been proposed as follows [28]: (18) where

is the initial gradient at the beginning and is the gradient.

In this study a

norm based iteration stop criterion similar to RMS rule is proposed as follows: (19)

where

is a small percentage number.

For numerical examples studied in this paper, an error of =5% for deflection response would produce minimal changes of material property to be accurate and acceptable for the engineering analysis and design purpose. However, one shall note that α=5% may not warrantee the exact solution, and thus this value may be changed for different engineering problems. Generally a lower α value would achieve a higher accuracy, but require more computation cost.

8

If the iteration stop criterion satisfies during the iteration the computation stops and outputs material model parameters as the inverse results; otherwise, the test function, gradient vector and Hessian matrix are calculated in order to invert the new material model parameters (see Figure 3).

Geometry, load, model seeds Simulating responses

Measured responses

<5%

Yes

No Solve test function

Inverted parameters

Calculate gradient

Search direction

Step length Update model parameters Figure 3. Flow chart for the proposed inverse computation procedure. To determine the test function, applying the first order variation of Lagrangian with respect to , the first-order adjoint equation is formulated as follows in order to determine :

(20) where

is variation of deflection variable

that is calculated from Equation (16).

Given the calculated and , applying the first order variation to the Lagrangian with respect to material parameter , the gradient g can be determined from the following first order decision equation:

(21)

9

where

is a variation of

, and determined as follows:

. can be dismissed on both sides. Thus, the strong form of the gradient can be

(22) In summary, the calculation of gradient involves two linear system solutions (one for the first order state governing equation (16) the other from the first order adjoint Equation (20)). Finding the Hessian matrix also involves three computation steps with three equilibriums. These three equilibriums are based on the second order variations of lagrangian with respect to , , and , respectively, as follows:: (23)

(24)

(25) From Equation (23) the incremental deflection can be calculated, and then from Equation (24) the incremental test function is computed given value. Given the calculated values of , , and , the hessian matrix is determined from Equation (25). These computations also involve two linear system solutions (Equations (24) and (25)). The vast majority of the computation costs are spent on the computation of gradient vector and Hessian matrix since it involves four linear system solutions in total for ach iteration step (two for the gradient and another two for the hessian computation). With the gradient and hessian matrix calculated, the inverse search direction at the iteration step can be determined according to Newton’s method that is derived from the Taylor’s theorem as follows: (26) where

and

are the Hessian matrix and gradient computed at the

iteration step, respectively.

The hessian matrix is used to improve convergence speed as compared to other methods which only calculate gradient to satisfy the first order optimal condition [25]. The Wolfe conditions were proposed for performing the inexact line search in an efficient way to compute the search direction and step length as follows: (i)

(27a)

(ii)

(27b)

10

where is displacement vector; is the material model parameter computed at the iteration step; is the transpose of ; are two constants such that and is the step length as used to constrain the ―diversion‖ (over- or under-estimated search direction) toward a more stable inverse computation. Inequality (i) in Equation (27a) is known as Armijo rule, and the inequality (ii) in Equation (27b) is known as curvature condition. The Armijo rule is to ensure that decreases the parameter value sufficiently, and the inequality (ii) is to ensure that the gradient has been reduced sufficiently [27]. In particular, if the two inequalities hold when the cosine of the angle between from zero, then the gradient . research is set at and inequality of satisfies.

and the gradient

(

) is bounded away

value is usually set quite small while value is much bigger [27]. In this for numerical implementation. To maintain a descent search direction, the

The inequality (ii) of Equation (27) is modified in this study as a strong Wolfe condition on the curvature as follows: (28) Consequently, the Armijo rule is also modified in this study to determine a stable step length following inequality:

from the

(29) The initial step length is set as a value proportional to for numerical practice purpose. Then at each iteration step, is reduced to the half of the value at the previous iteration step until the modified Armijo rule is satisfied. With the calculated and values, the material model parameter values at the iteration step are inverted or updated as follows: (30) Given and the modeled response values, the based value can be recomputed for evaluation, and this iteration continues until is lower than the target value (e.g. 5%) when the updated material model parameters are considered acceptable inversion results. The inverted material model parameters include the GM model parameters ( and ) for the viscoelastic material(s) and elastic moduli ( ) of the elastic material(s) within the heterogeneous structure. For the VE properties, the master curve of dynamic modulus in the frequency ( ) domain can be transformed from the inverted relaxation modulus through the Laplace transformation as follows [29]: (31) where variable and

is relaxation modulus; .

is angular frequency;

is Laplace transformation;

is the transform

The storage and loss modulus (real and imaginary part) at the frequency domain can be expanded as follows: (32)

11

(33) where

is the frequency in Hz.

By assigning a wide range of frequency to these formulas, the master curve of dynamic modulus ranging from the minimum modulus at the rubbery stage to the plateau modulus at the glassy stage can be determined.

2.4. Numerical solution method The numerical solution involves computations of gradient vector, Hessian matrix, search direction and step length, and lastly the inversion of material model properties. The computation involves complexity as considering the time dependency, displacement memory effect of material viscoelasticity, and the dynamics feature within a heterogeneous structure. The vast majority of computation costs and complexity rely on the calculations of gradient and hessian by solving the first and second order PDEs with triple integrals of space domain, sub-time domain of material viscoelasticity and the major time domain. Accordingly in this study a time-domain finite element (FE) method is developed for solving those resulting PDEs. The material and structural responses are solved at the forward time steps from the first order governing state equation, and details were already presented in another article [30]. The test function values are solved at the backward time steps (from the time end of to the time beginning of zero) due to the time dependency of dynamics and material viscoelasticity. The solution methods of the second order differentials are similar to that of the first order differentials. Mathematical derivations and numerical solutions of gradient and Hessian for this proposed dynamic VE material inverse problem are developed and presented in the Appendix A with details. Different linear system solvers have been proposed including the direct and iterative ones. Iterative linear solvers have used different algorithms including the preconditioner, parallel computation, and multigrid scalable method, etc. The main focus of this work is establishing a systematic framework for the treatment of the inverse medium problem of dynamic viscoelastic material properties. In this study a linear system solver based on the Cholesky factorization method with iterative refinement is used to solve the positive definite system of linear equations in band storage mode. This research primarily deals with the 2-D (two-dimensional) problems where using the factorization with iterative refinement is affordable, although it is necessary to use scalable algorithms for treating large-scale 3-D (three dimensional) problems. The chief advantage of the Cholesky factorization with iterative refinement method lies in its robustness [31] and accuracy [32]. However, this method may have poor memory scalability for large system problem, and thus it would be very useful to use other methods such as the scalable multigrid Krylov method for solving the 3-D large-system problems as the future research. The scalable Krylov iterative solver preconditioned by multigrid are intrinsically suitable for parallel computing of the 3-D large scale problem that requires massive amount of computing resources [33]. Multigrid used as preconditioner within a Newton-Krylov solver is algorithmically scalable as its chief advantage, and it yields faster computation speed than the diagonal preconditioned conjugate method [34]. Parallel multigrid based solver with Krylov subspace methods have also been developed for solving complex physical problems such as for simulating the 3-D groundwater flow using either finite element [35] or finite difference discretization [36]. It can scale to thousands of processor cores for solving billions of unknown parameters [37]. Multi-level domain decompositions may also be included [33]. However, its convergence could be worse for ill-conditioned problems [38]. A computer program written in Intel Visual Fortran® (IVF) language was developed for the entire numerical computation following procedures presented in Figure 3. The code is stable and throughout, we used double precision for the numerical simulations. The factorization method with iterative refinement produces more accurate and stable solutions if the double working precision is applied [39]. The inputs include the dynamic

12

loading time histories, the model geometry information, and measured structural responses. Two core modules were designed in the computer coding: 1) R-G-H (response-gradient-hessian) Module for the computations of response, gradient, and Hessian and 2) Inverse Computation Module, which computes the search direction and step length based on the R-G-H results with a dynamic linker to Module 1, and then inverts material model parameter values at each iteration step. The program outputs the inverted material model parameter values at the ending.

3. Validation The inverse model was implemented to a multilayered flexible pavement-soil structure under dynamic loading pulses for validation. Flexible pavement is one of the mostly used infrastructures, accounting for 90% of the paving areas in the US including highways, airports, and parking lots. This heterogeneous structure includes the AC as the surface layer, and the aggregate or stabilized base and/or subbase layers sitting on the soil foundation or earth as shown in Figure 4. The property of each layer is different than each other within the heterogeneous structure. AC is a typical VE material and its property is temperature dependent [1][40][41], and here it is modeled by the GM model. Other materials can be considered elastic [42].

Figure 4. FWD dynamic loading test on a multilayred pavement structrue. A dynamic falling weight deflectometer (FWD) equipment was used to produce a loading pulse on the surface by dropping a mass on the plate sitting on the top of the structure, which is used to emulate the vehicle loading. The loading pulse time history was recorded by a loading sensor. Several geophones were set at variable distances from the loading center to 1.5m to measure the vertical deflection responses. The Poisson’s ratios of AC/base/subbase/soil are given as 0.35, 0.35, 0.4, and 0.4 as common values and they have relatively minimal effects on structural responses [30][43][44].

3.1. Model domain Because the model domain has a much larger size (e.g., a 400 m2 parking lot) as compared to the loading area (a circle with a radius of 0.15 m and area of 0.07 m2), the model domain can be regarded axisymmetric as shown in Figure 5. Wave reflection at boundaries may affect numerical simulation results and thus a relatively large model size (i.e. 168.64 m in radius and 147.71 m in depth) is adopted so that the stress wave is unable to reach the boundaries at the ending of loading time. Alternatively, perfectly-matched-layers [45] may be used to accurately account for the truncation boundaries, without causing any reflections, as is done in [46][47][48] for subsurface

13

imaging of soils. The 8-node two-dimensional ring element was used to model the space domain (see Figure 5). The essential boundaries conditions are expressed as follows: (34a) (34b) (34c) where

is the symmetric line,

and

are the boundary area.

The measured vertical displacements or observations at the geophone zone of

are expressed as: (35)

where is the axisymmetric line of the model; is the boundary area of the far field; is the boundary area at the model bottom; is the geophone position area and are geophone measured vertical displacements.

Geophones

B

SB

147.71 m

AC

S

168.64 m Figure 5. Axisymmetric finite element model. Very often the material property of the same layer is considered space independent for the pavement structure, and thus the material parameter vector can be expressed as follows: (36) where, course; unbolded

, and are model parameters of the GM model of the AC material; is elastic modulus of base is elastic modulus of subbase course and is elastic modulus of soil foundation. Note that the is discretized values of , the same for and as corresponding to and .

The FE method is developed for the numerical solution. The isoparametric eight-node ring element is used for the model body of the space domain, and the two-node surface element is used for the surface traction area under

14

the FWD loading. Variable mesh is used for the model, in which the loading and displacement measurement area is meshed with fine size and coarser elements are used at the far field (see Figure 5). The forward response modeling is based on the dynamic VE model [30]. Sensitivity analysis is conducted to determine the minimum element and node sizes that achieve accurate simulations (i.e. responses have minimum changes when increasing element size). As a result, a minimum element number of 492 (41x12 in radius and depth) and node number of 1583 are used for the model. The inverse computation is performed following the proposed numerical procedure discussed earlier.

3.2. Forward modeling of response Extensive validations were performed for the forward modeling of responses by comparing its simulation results with that of the analytical solution and commercial software with details presented somewhere else[30][43]. The elastic modelling results were validated with the analytical solution at ideal boundary conditions. The dynamic elastic, static VE, and dynamic VE modeling results were validated by ANSYS – one of the most popular commercial software for FE analysis. Using the same mesh but different numerical algorithms, the proposed forward model attains almost identical results with that of ANSYS simulations (e.g. 0.1% or less) as shown in Figure 6 for the four-layered structure [30].

Relaxation modulus of AC - MPa E(t) - MPa

1400

a)

1200

AC HMA Base Subbase Soil

11.18 cm 17.78 cm 30.48 cm Infinite

3 2568 Mpa 2400 kg/m 213 MPa 1600 kg/m3

264 MPa 1800 kg/m3 212 MPa 2000 kg/m3

1000

800

600

400 0.001

0.01 Time (second) - logarithmic scale

0.1

15

Coding d1 Coding d3 Coding d5 Coding d7 ANSYS d2 ANSYS d4 ANSYS d6 Loading (kN)

b)

5.E-04

Deflection (m)

4.E-04

Coding d2 Coding d4 Coding d6 ANSYS d1 ANSYS d3 ANSYS d5 ANSYS d7

45 40 35 30 25

3.E-04

20 2.E-04

Loading (kN)

6.E-04

15 10

1.E-04

5 0.E+00 0

0.02

-1.E-04

0.06 0

0.04

-5

Time (second)

Figure 6. Dynamic viscoelastic forward model validation of a four-layered strucrture: a) structural and material properties and b) simulation results of deflections vs. that of ANSYS software [30]. (Note: “Coding di” and “ANSYS di” are the modeled vertical deflections on the surface of the distance of 0, 20.32, 30.48, 45.72, 60.96, 91.44, and 152.40 cm by the developed coding and ANSYS, respectively).

3.3. Validation of gradient and hessian It is crucial to attain accurate computations of gradients for accurate and robust inverse computation of material properties. For validation purpose, the gradients of the objective function with respect to the material model parameters are computed using the discrete method. The cost function can be written as follows based on Taylor’s theorem: (37) where is the cost function; and is the remainder.

is material model parameter vector;

is the material property variation vector

The gradient can be approximated as follows: g=

(38) (39)

where

and

A very small

is the ith material model parameter and gradient, respectively. value shall be used for higher accuracy (i.e. 0. 1%

).

The hessian matrix can be estimated following the Taylor’s theorem as well:

16

(40) where H is the hessian matrix. The individual element value of the Hessian matrix can be further expanded as follows:

(41) The validation of gradient and hessian values will be included in the next section.

3.4. Validation of inverse results Both theoretical and experimental validations are proposed. The theoretical validation is proposed in the following steps: 1) give the material properties as true values for the multilayered structure; 2) simulate the structural response of deflections under a given loading pulse using the powerful commercial FE software (i.e. ANSYS) and the simulation results are treated as observations; 3) invert material model parameters given the observations and loading pulse and 4) compare and as the validation results. In this evaluation, the inverse modelling has used the same mesh size as that of the forward modelling. This is fair since mesh size affect the forward modelled responses and then the inverted material properties. Experimental validation approach is proposed in the following steps: 1) conduct the dynamic loading tests on the field structure and measure deflection pulses with several geophones placed on different distances (variable distances can capture the effect of wave propagation); 2) invert material properties using the measured deflection pulses; 3) the same materials are extracted from field and their properties are measured from laboratory tests and 4) the inverted material properties are compared with that of laboratory measurements for evaluation. To minimize the potential chance of inverse crime, three strategies have been included in the model development and validation: 1) the regularization term is used to pick up the preferable material properties since different model parameter values may attain the same observations; 2) for the theatrical validation observations are produced by another independent commercial software, and realistic material and structural properties are used as from the experimentally based database of engineering projects [49] and 3) experimental validation was performed by comparing the inverted material properties with laboratory measurement values.

3.4.1. Theoretical validation A dynamic loading pulse was applied on a three-layer pavement structure with given geometries and material properties with information presented in Table 1. All material and structural properties fall within reasonable ranges of the engineering design according to the experimentally based database [49]. Given the loading pulse input and material properties, deflection responses at variable distances from 0 to 152.4 cm were simulated using the ANSYS software, which are called observations or ―measurements‖. With the same loading pulse input and observations, material properties are inverted and then compared to the given values. The peak loading and deflections at each distance have been used for inverting material elastic modulus where all materials are assumed elastic [8]. Table 1 listed the elastic modulus values of each material including the AC at

17

the loading frequency of 33 Hz. Two cases with and without regularization term are studied to invert the elastic modulus of each material using the measured peak loading and deflection. The regularization has applied a predefined material modulus range as listed in Table 1 (the last row) according to the experimentally based database of Federal Highway Administration’s Long Term Pavement Performance (LTPP) program [49]. The simulated deflection basin at those seven distances and inverted material modulus values are plotted in Figure 7. Results indicate that both cases yield almost identical observations as that of true values (see Figure 7a). However, the case without regularization has significantly different modulus values than the true values, e.g. the modulus of AC and base is 16.01% lower and 23.79% higher than that of true values, respectively. The case with regularization has improved the inversion accuracy since it tries to pick up the preferable material properties. Table 1. Structural and material parameters Layer

AC

Base

Soil

Thickness (cm)

19.84

17.78

Infinite

3

2,400

1,441

1,500

1500 1000-3500

125 35-350

115 35-350

Density (kg/m )

Deflection (mm)

True modulus (MPa) Regularization range (MPa)

0.50 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00

a) TRUE Without regularization With regularization

0

50

1600

150

200

b)

1400

Elastic modulus (MPa)

100 Distance (cm)

1200 TRUE

1000

Inverted: with regularization

800

Inverted: without regularization 600 400 200 0

AC 1

Base 2

Soil 3

18

.06

Ds-0 Ds-8 3.5E-04 Ds-12 Ds-18 3.0E-04 Ds-24 gap Figure term: a) deflection basins and b) inverted 2.5E-047. Elastic inversion with and without regularization Ds-36 Ds-60 modulus versus true values. Dm-0'' 2.0E-04 Dm-8'' Dm-12'' 1.5E-04 Dm-18'' Dm-24'' and loading pulses instead of the peak 1.0E-04 When the AC material is considered viscoelastic, the full deflection Dm-36'' values5.0E-05 are used for the inverse computation. The AC relaxation modulus has a master curve as time Dm-60''

gap

dependent. Figure 8a shows the simulated vertical deflections or deformations (with the final inverted material 0.0E+00and observations at variable distances, indicating a close match though some small gaps exist. Figure properties) 0.07 b) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 8b presents -5.0E-05 simulated max deflections at the time of 0.021 s (only portion of the full model is plotted). Figure 9 Time (second) plots the simulated deformation-loading loop, in which the max deflection has a time lag to the peak loading with o a phase angle of about 25 due to the wave propagation effects. Transient Viscoelastic

a) 4.5E-04

Ds-0 Ds-8 Ds-12 Ds-18 Ds-24 Ds-36 Ds-60 Dm-0'' Dm-8'' Dm-12'' Dm-18'' Dm-24'' Dm-36'' Dm-60''

06

Transient Elastic

4.0E-04

Deflection (m)

Ds-0 Ds-8 Ds-12 Ds-18 Ds-24 Ds-36 Ds-60 Dm-0'' Dm-8'' Dm-12'' Dm-18'' Dm-24'' Dm-36'' Dm-60''

gap

4.0E-04

Deflection (m)

3.5E-04 3.0E-04 2.5E-04 2.0E-04 1.5E-04 1.0E-04 5.0E-05

b)

Ds-0 Ds-8 Ds-12 Ds-18 Ds-24 Ds-36 gap Ds-60 Dm-0'' Dm-8'' Dm-12'' Dashed lines: simulations Dm-18'' solid lines: observations Dm-24'' Dm-36'' Dm-60''

3.7 3.2 2.8

d) -5.0E-05

m m 2.4

m 2.0

m

1.7

1.4 1.1 8.5

0.0E+00

0.07

Max deflections

m

m

15 .7 m

m m

m 0.021 s

0.00

0.01

0.02

0.03

Time (second)

0.04

0.05

0.06

0.07 7.1 m

Figure 8. Dynamic viscoelastic inversion: a) Inverted deflections vs. observations ( is deflection at distance of 0, 20.3, 30.5, 45.8, 61.0, 91.4, and 152.4 cm or 0, 8, 12, 18, 24, 36, and 60 inch, and b) simulated max deflection (only part of them model area is presented).

19

0.05

32.44

35.79

37.77

38.62

38.97

39.22

38.94

37.01

16.17

0.55 2.42

38.94

Simu-0 cm Simu-20.3 cm

32.44

Simu-30.5 cm

28.13

-1.01

-0.90

-0.69

-0.37

-0.02

0.41

0.99

1.57

Simu-152.4 cm

-0.94

Simu-91.4 cm

-0.95

Simu-61.0 cm

1.78

7.86

Simu-45.8 cm

1.47

4.96 0.55 2.42

-0.05

37.77

35.79

23.43 18.87 14.77 11.11

0.00

0.05

38.62

0.10

38.97

0.15

39.22 Deformation (mm)

-0.46-0.42 0.40 -0.44 -0.41 -0.32 0.35 -0.18 0.30 0.02

0.00

5

37.01

0.20

0

0.25

5

4.96

32.30

0.00 0.00

0

7.86

24.84

0.02

5

1

32.30

8.45

0

23.43 18.87 14.77 11.11

5

-1.02 -1.01 -1.06 -1.06 -1.02 -0.95-0.76 -0.41

0

28.13

0.11 0.14 0.87 3.32

0.11 0.14 0.87 3.32

24.84

0.05

5

16.17

0.02

8.45

0.00 0.00

0.88

0.00

0.35

0

Figure 9. Simulated deformation-loading radar chart (Simu-x cm: simulation at the distance of x cm). Figure 10 presents the inverted dynamic moduli of AC at frequency domain versus true or given values. Results indicate that the inverse results of (real modulus) are close to the given or true values. Inversed values (imaginary part) oscillate somehow, which can be explained by the limitation of the GM model in Prony series as it may produce non-smooth master curves for the AC material especially for the loss modulus due to the discrete spectrum of the formula [1]. In contrast, the existing method using the FDM for gradient approximation without Hessian calculation [25] produces less accurate inversion values which are lower than true values at the high frequency range. The proposed method also reduces iteration times (five vs. seven for this case study). Inverted elastic moduli of base/soil materials are close to given values, too (see Table 1). The inverted modulus of base and soil is 113 and 127 MPa, respectively, having reasonable agreements with the true values. Regarding the computation times, the proposed method solves total 1,200 linear systems ( , iteration times linear systems of each iteration total time steps) which is significantly lower than that of the existing method with total 7,140 linear systems ( ).

20

6000

a)

True values

Inverted with existing method 4000

Inverse values

3000

2000 1000 0 1.E-03 -1000

True values Inverted with new method

b)

1000

Inverted with new method

Dynamic moduli (MPa)

Dynamic moduli (MPa)

5000

1200

True values

Inverted with existing method Inverse values

800

True values

600

400 200

1.E-01

1.E+01

1.E+03

1.E+05

1.E+07

0 1.E-03

1.E-01

Reduced frequency (Hz)

1.E+01 1.E+03 Reduced frequency (Hz)

1.E+05

1.E+07

Figure 10. Inversion results vs. true values: a) storage modulus and b) loss modulus. For the exact Newton method of smooth problems the iteration number is mesh independent [50]. Meshindependent iteration algorithms with inexact Newton method have also been developed [51]. The developed numerical procedure with the Newton iteration for the dynamic viscoelastic problem of solids (2nd order derivative exists) is independent on mesh as illustrated in Table 2. Results indicate that variable meshes (e.g. node number ranging from 818 to 6531) almost attain the same iteration number of 5 (including the initial step using seeds as the 1st iteration). The case with an element mesh of 53 24 (radius by height) has iterated 6 times as its value (for iteration stop) at the 5th iteration is around the borderline of target value 5%. Table 2. Iteration number at variable mesh. Mesh 31x9 41x12 53x24 64x33

Element No. 279 492 1272 2112

Node No. 818 1583 3971 6531

Iteration No. 5 5 6 5

value1 1.13% 2.98% 1.61% 2.95%

1

Figure 11 plots the gradient values of AC, base, and soil. Note that is the summation of gradient values of each parameter of the GM model for AC material as (e.g. total 15 values for a GE model with 7 parallel series of elastic springs and viscous dashpots). Results indicate that gradient values drop rapidly with the increase of iteration number. The soil material has the highest gradient value mainly due to its lowest modulus value while significant impact on the cost function (

). Table 3 lists the gradient values

th

attained at the 5 iteration as compare to that of the discrete estimation. Results show that the computed gradients are fairly close to that of discrete approximations where a very small (i.e. 0.001E) is used.

21

-2.7E-02 Seeds

Gradient (mm2/MPa)

-2.2E-02

1st iteration 2nd iteration

-1.7E-02

3rd iteration 4th iteration

-1.2E-02 -7.0E-03 -2.0E-03

AC

Base

Soil

Figure 11. Gradients values of different materials at different iteration number. Table 3. Gradient values of adjoint method versus discrete approximation (10-4 mm2/MPa) Material AC Base Soil

Adjoint based -1.3724

Discrete approximation -1.3755

-1.2291 -5.6726

-1.2233 -5.7344

Table 4 lists the Hessian matrix element values of the diagonal as compared to that of the discrete estimation. Results have shown reasonable match between simulations and estimations. The errors of Hessian are slightly higher than that of gradients due to its 2nd order discretization with additional numerical errors added up. Table 4. Hessian element values in the diagonal (10-4 mm2/MPa) Parameters Adjoint based Discrete approximation

E0

E1

E2

E3

E4

E5

E6

E7

Ebase

Esoil

-1.4400 -1.4414

-7.9193 -7.1339

-1.7263 -1.7281

-4.7952 -4.7969

-1.4544 -1.4543

-1.7280 -1.7272

-1.1520 -1.1521

-2.4480 -2.4161

-2.5515 -2.6034

-3.1406 -3.1386

(2) Continued Parameters

2

3

4

5

6

7

Adjoint based

-1.5921

-3.3433

-1.7831

-5.2893

-1.2247

-1.2577

-7.0446

Discrete approximation

-1.5922

-3.34437

-1.7809

-5.2825

-1.2261

-1.2735

-7.0554

Note:

and

are for the GM model’s elastic springs and dashpots of the AC material.

3.4.2. Experimental validation FWD test was conducted at Florida State for a four-layer flexible pavement structure with a thick AC surface course (20.8 cm). Figure 12 presents the measured loading and deflection pulses using the DynaTest equipment. The vertical deflections ( to ) on the top of the structure were measured at the distances of 0, 20.3, 30.5, 30.5, 45.7, 70.0, 91.4, 121.9, and 152.4 cm. A sudden ―drop‖ of the loading pulse at around 0.02 s was observed, which was induced by the buffer effect of the dropping plate of the DynaTest equipment. Laboratory frequency sweep tests were conducted to determine dynamic moduli of the AC material at six frequencies (0.1, 0.5, 1, 5, 10,

22

and 25 Hz) and five temperatures (-10, 4.4, 21.1 and 54.4oC). Using the temperature-time superposition rule, the measured moduli were shifted to that at the reduced frequency for the reference temperature of 21.1oC. The master curve of laboratory measured dynamic moduli values ( ) at the reduced frequency was fitted using the sigmoidal function model. Figure 13 presents the inverted dynamic moduli master curve of versus laboratory measurements and the constructed sigmoidal function curve. Results indicate that for this case the vast majority of the inverted dynamic moduli values of the AC material are close to experimental results.

90

Input loading [-0.67, 28.8kN]

D1

80

D2

m)-6m) Deflection (x10

70

D3

60

D4

50

D5

40

D6

30

D7

20

D8

10

D9

-10 -20

0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0.022 0.024 0.026 0.028 0.03 0.032 0.034 0.036 0.038 0.04 0.042 0.044 0.046 0.048 0.05 0.052 0.054 0.056 0.058

0 2

6

10 14 18 22 26 30 34 38 42 46 50 54 58 Time ( 10-3 second)

Time (second)

Figure 12. Input loading and deflection pulses (

is the deflectyion at the

distance).

23

Dynamic Modulus (MPa) (MPa)

40050 35050 30050 25050 20050 15050

10050 5050 50 1.E-06

Unit:

- MPa,

E0 E1 ŋ1 E2 ŋ2 E3 ŋ3 E4 ŋ4 E5 ŋ5 E6 ŋ6 E7 ŋ7

MPa.s

63 4421.15506 0.30324 4829.8416 0.01716 4188.46613 13.29398 1.30E-06 3820.31484 6720.428 1.00E-06 4917.77846 385.29108 5575.51507 7855.42130

1.E-04

Sigmoidal model Inverse results C -10.0 ooC oC 4.4 oC oC 21.1 oC oC 37.8 oC oC 54.4 oC

1.E-02

1.E+00 1.E+02 Frequecny f (Hz)

1.E+04

1.E+06

Lab test

1.E+08

Figure 13. Inverse results versus measurements of asphalt concrete. Table 5 presents the laboratory-measured moduli and the inverted ones of unbound materials (base, subbase and soil). Laboratory tests were conducted on 30 samples of each material. Results indicate that the inverted values are closer to that of the max values of the laboratory measurements, while greater than the median values. The difference may be explained by the following reasons: 1) the laboratory test conditions of axial loadings applied on a small specimen was very different than the in-situ conditions of a large pavement-soil structure under the vehicle loading; 2) in laboratory tests variable axial loading levels were applied, resulting in variable moduli values (see Table 5), but it is unclear which one could more accurately emulate the in-situ situation and 3) FWD tests were performed in a different date than laboratory testing, during which the moduli of unbound materials might change due to repeated loading effects and the changing of environmental conditions (e.g. moisture). Table 5. Measured and inverted unbound material modulus in MPa Lab test Base Subbase Soil

Min 94 63 48

Max 413 409 114

Mean 233 192 82

Inverted 428 359 201

3.5. Computation time issues The computation costs primarily rely on the solutions of the linear systems. To calculate the gradient and Hessian, for each iteration step the Lagrangian optimization method involves only four linear system solutions (two for gradient: response and test function, and another two for Hessian: incremental response and test function) as independent on model parameter number. In comparison, the two-stage approach will require the same number

24

of linear system solutions as the number of model parameters for each iteration of the inverse computation using response outputs in order to attain reasonable accuracy. For example, the gradient-based finite difference algorithm requires one forward linear system solution of response to estimate the gradient value of each model parameter (see Equation (1)). The proposed method may also reduce iteration times (e.g. 5 versus 7 for example 1) as it calculates Hessian for accelerating the convergence speed. Table 6 presents the computation time or number of linear system solutions for only one iteration step of a four-parameter elastic media inverse problem and an 18-parameter dynamic VE inverse problem (e.g. Example 2). Results show that the developed method has dramatically improved the computation speed for the dynamic VE problem as compared to the finite difference algorithm (total time of 400 versus 1800), regardless that it also improves computation accuracy (see Figure 10). Table 6. Computation time comparison for one iteration Model parameters Approach Linear systems Time steps Total time

4 Two-stage 4 1 4

18 Proposed 4 1 4

Two-stage 18 100 1800

Proposed 4 100 400

However, the Lagrangian optimization method involves much more complex mathematical derivations. When dealing with relatively simple inverse problems with a small model parameter number such as for the elastic inverse problem with only a few model parameters, the two-stage method is sufficient and accurate enough [25], in which the analytical solution instead of finite element solution can also be used for faster forward computation of responses.

4. Conclusion We developed an adjoint-baseded inverse model and mathematical solution based on the Lagrange multiplier theorem to invert VE properties and dynamic deformation of heterogeneous structures to improve computation reliability and efficiency. We proposed robust numerical algorithms for fast and accurate calculations, and a Galerkin finite-element time-domain method to solve the resulting PDEs. Validation results on the layered structure under dynamic loading tests indicate that the proposed model can successfully invert the master curves of dynamic moduli at variable temperatures for accurate modeling of material properties and dynamic deformations. Compared to the existing approaches as often being used such as the finite-difference-based gradient method, the newly developed method improves computation accuracy and speed. The developed method may be applied to a broad range of heterogeneous materials and structures for the dynamic viscoelastic material property inverse problem by changing the length and time scales. The numerical implementation in this study is primarily for a 2-D model domain, and it can be extended to the 3-D domain where a scalable multigrid Krylov solver is recommended for linear system solutions as the future research.

Acknowledgement The experimental data were supported by the US Department of Transportation (DOT) Federal Highway Administration’s LTPP program. We greatly appreciate Professor Jorge Prozzi at Civil Engineering for his important advice on using LTPP data, Professor Georg Stadler’s (currently at New York University the Department of Mathematics of the Courant Institute) patient guidance and Professor Omar Ghattas’ critical comments on inverse computation both at the Institute for Computational Engineering and Sciences (ICES), Dr. Arash Fathi at the Structural Mechanics group (currently postdoctoral researcher at ICES) and Dr. Tobby Isaac at the Department of Computational and Applied Mathematics for their very useful inputs, all at the University of Texas at Austin. The first author studied at Civil Engineering and worked at the ICES afterwards at the University

25

of Texas at Austin where the main work and additional analysis were done, respectively, before he moved to the state of Minnesota as a Computational Mechanics Engineer in May 2015.

26

Appendix A: numerical solution method for PDEs of dynamic viscoelastic inversion problem of heterogeneous structures We have developed a time-domain finite element (FE) method for numerical solutions of resulting PDEs of the dynamic viscoelastic inverse problem of heterogeneous structures. The core of numerical computation is to solve the gradient and hessian. With numerical values of gradient and Hessian, the search direction and step length can be calculated for inverting material properties as discussed in in the main text. The three key variables in the Lagriangian include displacement , test function , and material model parameter m. In the FE solution, they are discretized to that at FE nodes as , , ( is the shape function matrix, , , and m are displacement, test function, and material model parameters discretized at FE nodes of elements, respectively). Substitute the discretized form of these three variables into the Lagrangian equation (15), the Lagrangian weak form can be re-expressed as follows:

(A1) where

is the shape function matrix to discretize the loading at the surface area

damping matrix; matrix.

is the mass matrix and

;

is the

is a spatially discretized relaxation modulus

(A2) where

is the strain-displacement matrix and

is the relaxation modulus matrix.

With this Lagrangian weak form, the gradient vector ( ) and Hessian matrix ( ) at FE nodes, can be computed with details presented in the following sections. Note that all PDEs are all with respect to the discretized forms of and in the numerical solution. A1. Computation of Gradient Vector The computation of gradient includes three steps of consequence (see Equation (20) to (22)), and the numerical solutions for the specific dynamic VE problem of heterogeneous structures are developed as follows. A1.1 State equation to determine displacement Applying the first order variation of the Lagrangian with respect to the test function , the governing state equation weak form can be formed as follows:

27

(A3) is arbitrary and can be dismissed on both sides

, thus the following satisfies: (A4)

where The displacement tangent is discretized as follows based on the explicit Euler rule: (A5) (A6) where

is the time step, and is the sub-time step.

Substitute (A5) and (A6) into (A4) and then discretize it in time to attain the final weak form:

(A7) From this weak form the dynamic VE displacement at time steps of can be solved. A timedomain finite element method with fast and robust numerical algorithms was developed for this solution and already presented elsewhere [30][43]. Alternative methods for solving dynamic viscoelastic response of the structural system include the model reduction method [52]. A1.2 Adjoint equation to determine test function Space discretization The numerical solution of the test function for the dynamic VE problem of a heterogeneous structure was developed step-by-step in the following manner. Using the Galerkin method, apply the first order variation of the Larangian equation (A1) with respect to displacement variable , the weak form of the adjoint equation can be achieved as follows:

28

(A8) In this weak form the 2nd term is defined as: (A9) In this

expression the first integration is derived as follows:

(A10) Substitute (A10) into the

expression and re-derive it to a new formula as follows: (A11)

This double integration of with respect to to an equivalent integration with respect to

and then and then

(see Figure A1a) can be converted (see Figure A1b) as follows: (A12)

Figure A1. Integration of 2-D time field of and notations are renamed to each other to arrive at the new form of

. with notation change only: (A13)

Define

for substitution, and according to the integration by parts,

can be re-derived as

follows:

29

(A14)

Define

as a dynamic inertia integral.

the integration by parts (given

can be re-derived as follows according to

):

(A15) Substitute Equations (A13) to (A15) into (A8) to re-derive the weak form of the first order adjoint equation as follows:

(A16) According to the chain rule and integration by parts, the time integration of the third term in the adjoint equation above can be derived as follows:

(A17) Substitute Equation (A17) into (A16), and dismiss the arbitrary

term on both sides

: (A18)

As a tangent,

can be approximated using the explicit Euler method as follows: (A19)

Given a very small time step length, this approximation would be accurate enough. Define the VE stiffness matrix at the sub-time step where as follows:

for each of

time steps

30

(A20) Define the displacement misfit vector as follows: (A21) where is measured displacements or observations discretized to FE nodes, and displacements at FE nodes as determined from the first order state equation.

is simulated

Substitute Equations (A19) to (A21) into (A18) to arrive at a briefed equilibrium: (A22) where

is the current time step for time

The numerical solution of following.

.

is a key for the successful solution of this equilibrium as described in the

Firstly when the Poisson’s ratio is considered constant often as the case when material temperature has no changes during the loading period, the shear and bulk relaxation moduli can be determined from the relaxation modulus as follows: (A23a) (A23b) Substitute Equation (A23) into Equation (6), the relaxation modulus matrix can be re-expressed as follows for numerical solution purpose: (A24) where

is a viscoelasticity shape function matrix.

is expressed as follows [30]: (A25) Substitute Equation (A24) into Equation (A20), the VE stiffness matrix can be derived as follows:

(A26)

31

Thus, given

and relaxation modulus

and VE stiffness matrix

can be numerically calculated.

Time Discretization and Linear System The Houbolt finite-difference method (FDM) is adopted for the time discretization of test function, velocity and acceleration since it is less dependent on time step length and unconditionally stable [53], as follows: (A27a) (A27b) where is a negative or backward time step length that equals to forward calculation of response .

) as opposite to that in the

As results, a relatively large time step length (e.g. 0.001 s) can be used to achieve satisfied accuracy based on sensitivity analysis. Substitute Equation (A27) into Equation (A22) and rearrange the equilibrium of weak as follows:

(A28) where

is

solution at time step of as

.

Equation (A28) can be re-arranged as a linear system in the following: (A29) (A30)

(A31) where equation.

is the dynamic VE stiffness matrix and

is the reaction vector for the first order adjoint

By assembling all elements of each material including both the VE and elastic materials, the global linear system can be formed as follows: (A32)

32

(A33) (A34) where is global stiffness matrix of the model domain, a 2 matrix with as the total FE node number; is global elastic stiffness matrix, a 2 matrix, and its element values at the VE materials are assigned with zeros while it is only functional for the elastic materials and is the 4th order elasticity tensor of elastic materials. To solve the global linear system, the Cholesky factorization with iterative refinement method is employed. As shown in Figure A2 the test function is computed for at the reversed order of time s steps. The total time is broken into time steps, and the calculation starts from the loading period end of which is assigned as the 0th time step for the calculation of , and then the time step of is solved in a backward order, until the last time step at the time of zero ( ) is reached. The time order for solution is opposite to that for the forward solution of displacement discussed early.

Figure A2. Time discretization for the solution of test function. A1.3 Decision equation to determine gradient Given the displacement and test function calculated at FE nodes, the gradient can be determined in the following. Weak Form of Gradient The first order decision or control equation is formed by applying the first order variation of the Lagrangian with respect to material model parameter to achieve the weak form as follows:

(A35)

33

where

is an identity matrix if m is considered constant or space independent.

The damping matrix C is determined following the Rayleigh damping model as follows: (A36) where

are coefficients; M is the mass matrixand

is the elastic stiffness matrix.

Thus, the derivative of C can be calculated as follows: (A37) is arbitrary and thus it can be dismissed on both sides of Equation (A35). Therefore, the gradient vector determined as follows:

is

(A38) Define the relaxation modulus derivative as follows: (A39) is a second rank matrix with a size of material model parameter number). For each

with regard to the

(

is the total FE node number and

is the total

material model parameter: (A40)

where

is the

material model parameter,

Define the VE stiffness matrix derivative at the

. time step as follows: (A41)

Discretize the displacement tangent of the sub-time step following the Euler method as follows: (A42)

Break the time period

into

steps, applying explicit Euler rule to discretize

the time integration, and then apply the trapezoidal rule for the time integration of formula as follows:

which is taken out of to achieve the discretized

(A43)

34

Substitute the time discretization Equations (A41) to (A43) into Equation (A38) to arrive at the final weak form as follows for the solution of gradient:

(A44) where

,

are the displacement and test function at FE nodes as solved in the earlier sections.

Solution of Gradient of Individual Material Parameter within Heterogeneous Structure One of the key elements in the gradient Equation (A44) is computed in the following means.

(VE stiffness matrix derivative), which is

For the GM model of the VE materials within the heterogeneous structure, substitute (A24) into (A41), the numerical solution of is derived as follows: (A45) where

is the viscoelasticity shape function matrix (see (A25).

i) for the elastic modulus of the first spring at infinite time

: (A46)

Substitute (A46) into (A44), the gradient of

is determined as follows: (A47)

ii) for other elastic moduli of those springs distributed in the GM model:

(A48)

Therefore, the time integration of

can be calculated as follows:

35

(A49)

Substitute Equations (A49) into (A45) and then into (A44) to arrive at the final form of gradient for the elastic modulus of springs ( ):

(A50)

where

is the elastic modulus of the

spring element

.

iii) for the viscosity of those dashpots: For the

dashpot viscosity of the GM model

), the gradient of Equation (A44) can be reduced to:

(A51) where: (A52)

(A53)

Substitute Equation (A53) into Equation (A51), and then rearrange it to arrive at the final discretized form for the gradient of viscosity as follows:

(A54) iv) for the elastic materials within the heterogeneous structure:

36

The gradient of the elastic modulus of the elastic material of the heterogeneous structural system ( , is the total number of elastic materials). The gradient is calculated as follows: (A55) where is the elastic modulus of the elastic material of the structural system; range for penalization and is the elastic stiffness matrix derivative.

is predefined moduli

is derived as follows: (A56)

where

is the 4th order elasticity tensor.

Applying time discretization (the same as those discussed above), the gradient of the modulus of the material in the final work form can be determined as follows:

elastic

(A57) Thus, the gradient for each material model parameter within a heterogeneous structure is solved by far. A2 Second Order Variational Method Computing Hessian Matrix The PDE solutions of Hessian matrix of the dynamic VE problem of heterogeneous structures are developed and detailed in the following. Computation of Hessian matrix involves three steps in sequence (second order state, adjoint, and decision equations) and two linear system solutions following equations (23) to (25). A2.1 Incremental state equation to determine incremental deflection The second order variation of lagrangian with respect to Equation (A3)) with respect to

is equal to the first order variation of

(see

as follows: (A58)

The second order variation with respect to follows:

is equal to the first order variation of

with respect to

as

(A59)

37

The second order variation of lagrangian with respect to respect to as follows:

is equal to the first order variation of

with

(A60) where

:

(A61) Substitute Equation (A61) into (A60):

(A62) Sum Equations (A59), (A60) and (A62), and then zero it as follows:

(A63) is arbitrary and can be dismissed on both sides

. Hence, the following equation satisfies:

(A64)

where

is a matrix with size of

The time integration of

,the ith column of the matrix is

for

.

is derived as follows:

(A65)

38

Substitute Equation (A65) into (A64):

(A66) Discretize the velocity and acceleration following the Houbolt method for the incremental displacement: (A67a) (A67b) Discretize the tangent of incremental displacement following the explicit Euler rule: (A68) Substitute equation (A67) and (A68) into equation (A66) and then discretize the integrations:

(A69)

where

is the VE stiffness matrix and

is the first

order derivative of VE stiffness matrix. Equation (A69) can be rearranged and reduced to a global linear system as follows: (A70) )

(A71)

(A72) where is a stiffness matrix; matrix and

is a reaction vector; is the first order derivative of VE stiffness are computed incremental displacements at previous time steps.

39

Now one can solve the global linear system to determine the incremental displacement method discussed in section A1.1.

following the same

A2.2 Incremental adjoint equation to determine incremental test function The second order variation of Lagrangian with respect to Equation (A8)) with respect to as follows:

is equal to the first order variation of

(see

(A73) The second order variation of Lagrangian with respect to respect to as follows:

is equal to the first order variation of

with

(A74) (A75) where

is a mass matrix unit.

The second order variation of Lagrangian with respect to respect to and can be derived as follows:

equals to the first order variation of

with

(A76) Sum these three terms (Equations (A73), (A74) and (A76)) and zero it as follows:

(A77) is arbitrary and can be dismissed on both sides

, and thus the following equation satisfies:

(A78)

40

Discretize the velocity and acceleration of the incremental test function with the Houbolt method: (A79a) (A79b) Discretize the tangent of incremental test function following the explicit Euler rule: (A80) Substitute equations (A79) and (A80) into (A78) to arrive at the equilibrium:

(A81) where

is the VE stiffness matrix;

is the first order derivative of VE stiffness

matrix and

is calculated value at the current time or the

time step.

Equation (A81) can be rearranged as follows for building a global linear system: (A82) (A83)

(A84) where system.

is stiffness matrix for test function linear system and

is a reaction vector for test function linear

The global linear system is solved to determine the incremental test function discussed in section A1.2. One shall note that the solution of , the same as for as opposite to that of .

following the same method , is in a reversed time order

A2.3 Incremental decision equation to determine Hessian

41

Given the incremental displacement and test function calculated above, the Hessian matrix can be determined as follows in sequence. The second order variation of Lagrangian with respect to Equation (A35)) with respect to as follows:

is equal to the first order variation of

(see

(A85) The second order variation of Lagrangian with respect to respect to as follows:

is equal to the first order variation of

with

(A86) The second order variation of Lagrangian with respect to respect to as follows:

is equal to the first order variation of

with

(A87) Sum these three terms to determine the hessian matrix as follows:

(A88) is arbitrary and thus can be dismissed on both sides. is an identity matrix if material properties are space independent. Therefore, the Hessian matrix can be determined as follows in this case: (A89) where

and

is the first and second order derivative of VE stiffness matrix, respectively.

is defined as follows: (A90) Using the same time discretization methods discussed early, equation (A89) can be rewritten as follows as the final work form:

42

(A91)

+ where H is the hessian matrix with a size of inversion computation.

by

and

is the total number of model parameters for the

With the gradient vector and hessian matrix solved, the step length and search direction can be determined as discussed in section 2, and then the material model parameters are inverted or updated at the current iteration step. This procedure repeats until the simulated responses are very close to observations based on the iteration stop criterion when the updated material properties at the current iteration step are turned out as the inverse results.

Reference [1]. Xu Q, Solaimanian M. Modeling linear VE property of asphalt concrete by Huet-Sayegh model. International Journal of Pavement Engineering 2009; 10: 401─422. [2]. Motamed A, Bhasin A., and Liechti KM. Constitutive modeling of the nonlinearlyvVE response of asphalt binders; incorporating three-dimensional effects. Mechanics of Time-Dependent Materials 2013; 17: 83–109. [3]. Cobb III CH. Summer Skin Care: http://www.hughston.com/hha/a_17_3_4.htm

Get

your

vitamin

D

in

moderation.

2004;

[4]. Xu F, Lu TJ, Seffen KA. Biothermomechanics of skin tissues. Journal of the Mechanics and Physics of Solids 2008; 56: 1852–1884. [5]. Tronto J, Bordonal AC, Naal Z, and Valim JB. Conducting polymers/layered double hydroxides intercalated nanocomposites. Materials Science -Advanced Topics 2013; Chapter 1, 1st edition. [6]. Hossain, SS, Hossainy SFA., Bazilevs Y, Calo VM, and·Hughes TJR. Mathematical modeling of coupled drug and drug-encapsulated nanoparticle transport in patient-specific coronary artery walls. Computational Mechanics 2012; 49: 213–242. [7]. Xu Q., LN Mohammad. Modeling asphalt pavement rutting under accelerated testing. Road Materials and Pavement Design 2008; 9: 665-687. [8]. Wang F, Lytton R. System identification method for backcalculating pavement layer properties. Transportation Research Record 1993; 1384: 1-7. [9]. Xiong GC, Yang Y, and Chuan W. Surgical Treatment for Diffuse Coronary Artery Diseases. Chapter 16 in "Artery Bypass"2013; book edited by Wilbert SA. [10]. Catheline S, Gennisson, JL, Delon G, Fink M, Sinkus R, Abouelkaram S, and Culiolic J. Measurement of VE properties of homogeneous soft solid using transient elastography: An inverse problem approach. Journal of Acoustic Society of America 2014; 116: 3734–3741. [11]. Brigham JC, Aquino W, Mitri FG, Greenleaf JF, and Fatemi M. Inverse estimation of VE material properties for solids immersed in fluids using vibroacoustic techniques. Journal of Applied Physics 2007; 101: 023509. [12]. Lei, F., Szeri, A.Z. Inverse analysis of constitutive models: biological soft tissues. Journal of Biomechanics 2007; 40: 936–940.

43

[13]. Zhao R, Wyss K, and Simmons CA. Comparison of analytical and inverse finite element approaches to estimate cell VE properties by micropipette aspiration. Journal of Biomechanics 2009; 42: 2768-73. [14]. Araújo AL, Soares CMM, Herskovits J, Pedersen P. Estimation of piezoelastic and VE properties in laminated structures. Composite Structures 2009; 87: 168–174. [15]. Herskovits J, Dubeux V, Soares CMM., Araujo AL. Interior point algorithms for nonlinear least squares problems. Inverse Problem Science and Engineering 2004; 12: 211–23. [16]. Sims AM, Stait-Gardner T, Fong L, Morley JW, Price WS, Hoffman M, Simmons A, Chindhelm K. Elastic and VE properties of porcine subdermal fat using MRI and inverse FEA. Biomech Model Mechanobiol 2010; 9: 703–711. [17]. Giavazzi S, Ganatea MF, Trkov M, Šuštarič P, Rodič T. Inverse determination of VE properties of human fingertip skin. Materials and Geoenvironment 2010; 57: 1–16. [18]. Yuung HL, Fung-Huei Y, Ching LL. Application of FSQP to Inverse Estimation of the Constitutive Constants and Friction Coefficient in the Nosing Process. Materials Science Forum 2005; 505–507: 685-690. [19]. Araújo AL, Cristovao, Soares MMS, Soares CAM, Herskovits J. Characterisation by inverse techniques of elastic, VE and piezoelectric properties of anisotropic sandwich adaptive structures. Applied Composite Materials 2010; 17: 543–556. [20]. Magnuson AH, Lytton RL, and Briggs R. Comparison of computer predictions and field data for dynamic analysis of falling-weight deflectometer data. Transportation Research Record 1991; 1293: 6171. [21]. Scarpas A, Blaauwendraad J. Spectral element technique for efficient parameter identification of layered media, part III: VE aspects. International Journal of Solids and Structures 2002; 39: 2189–2201. [22]. Levenberg E. Inverse analysis of VE pavement properties using data from embedded instrumentation. International Journal for Numerical and Analytical Methods in Geomechanics 2012; 37: 1016–1033. [23]. Varma S, Kutay E, and Levenberg E. A VE genetic algorithm for inverse analysis of asphalt layer properties from falling weight deflections. Transportation Research Record 2013; 2369: 38-46. [24]. Liang RY, Zhu JX. Efficient computation algorithms for forward and backward analysis of a dynamic pavement system. Computers & Structures 1998; 69: 254-263. [25]. Xu Q, Prozzi JA. A finite-element and Newton-Raphson method for inverse computing multilayer moduli. Finite Elements in Analysis and Design 2014; 81: 57-68. [26]. Petra, N. and Stadler, G. Model Variational Inverse Problems Governed by Partial Differential Equations. ICES Report 11-056, the University of Texas at Austin.Christensen RM. Theory of Viscoelasticity 1982, 2nd edition, Acadenuc Press, Inc. [27].

Nocedal, J, Wright SJ. Numerical Optimization 2006, 2nd edition, Springer-Verlag.

[28]. Petra N., Zhu H., Stadler G., Hughes T.J.R., Ghattas O. An inexact Gauss–Newton method for inversion of basal sliding and rheology parameters in a nonlinear Stokes ice sheet model. Journal of Glaciology 2012; 58: 889-903.

44

[29]. Park SW, Schapery RA. Methods of interconversion between linear VE material functions. Part I—a numerical method based on Prony series. International Journal of Solids and Structures 1999; 36: 1653-1675. [30]. Xu Q, Prozzi JA. A time domain finite-element method for dynamic viscoelastic solution of layered half space. Computers & Structures 2015; 160: 20-39. [31]. Chadwick J.N., Bindel D.S. An efficient solver for sparse linear systems based on rank-structured Cholesky Factorization. arXiv:1507.05593. 2015. [32]. Li S., Gu M., Wu C.J., and Xia J. New efficient and robust HSS Cholesky factorization of SPD matrices. SIAM Journal on Matrix Analysis and Application 2012; 33: 886–904. [33]. Badia S., Martin A.F., Principe J. A highly scalable parallel implementation of balancing domain decomposition by constraints. SIAM Journal on Scientific Computing 2014; 36: C190-C218. [34]. Jones J.E., Woodward C.S. Newton–Krylov-multigrid solvers for large-scale, highly heterogeneous, variably saturated flow problems. Advances in Water Resources 2001; 24: 763-774. [35]. Mahinthakumar K., and Saied F. Multigrid and krylov solvers for large scale finite element groundwater flow simulations on distributed memory parallel platforms. Technical report, ORNL/TM134441, 1197, Oak Ridge National Laboratory, 1998. [36]. Ashby S.F., and Falgout R.D. A parallel multigrid preconditioned conju-gate gradient algorithm for groundwater ow simulations. Technical Report UCRL-JC-122359, Lawrence Livermore National Laboratory, CA, 1995. [37]. Burstedde C., Ghattas O., Stadler G., Tu T., Wilcox L.C. Parallel scalable adjoint-based adaptive solution of variable-viscosity Stokes flow problems. Computer Methods in Applied Mechanics and Engineering 2009; 198: 1691–1700. [38].

Trottenberg U., Oosterlee C., and Schuller A. Multigrid. 1 edition, Academic Press, 2001.

[39]. Nicholas H. Accuracy and Stability of Numerical Algorithms. SIAM: Society for Industrial and Applied Mathematics. University of Manchester, Manchester, England, 2nd edition, 2002. [40]. Xu, Q. and Solaimanian, M. Modeling temperature distribution and thermal property of asphalt concrete for laboratory testing applications. Construction and Building Materials 2010; 24: 487–497. [41]. Xu Q and Mansour Solaimanian. Measurement and evaluation of asphalt concrete thermal expansion and contraction. Journal of Testing and Evaluation 2008; 36: 507. [42]. Xu Q, Prozzi JA. Static versus VE wave propagation approach for simulating loading effects on pavement structures. Construction and Building Materials 2014; 53: 584-595. [43]. Xu Q. Development of a computational method for inverting dynamic moduli of multilayer systems. Dissertation, the University of Texas at Austin 2014. [44]. Xu Q, Chang GK. Experimental and numerical study of asphalt material geospatial heterogeneity with intelligent compaction technology on roads. Construction and Building Materials 2014; 72C: 189198. [45]. Fathi A, Poursartip B, and Kallivokas LF. Time-domain hybrid formulations for wave simulations in three-dimensional PML-truncated heterogeneous media. International Journal for Numerical Methods in Engineering 2015; 101: 165-198.

45

[46]. Kallivokas LF, Fathi A, Kucukcoban S, Stokoe KH, and Ghattas O. Site characterization using full wave form inversion. Soil Dynamics and Earthquake Engineering 2013; 47: 62-82. [47]. Fathi A, Kallivokas LF, Poursartip B (2015) Full-waveform inversion in three-dimensional PMLtruncated elastic media. Computer Methods in Applied Mechanic and Engineering 2015; 296: 39–72. [48]. Fathi A. Full-waveform inversion in three-dimensional PML-truncated elastic media: theory, computations, and field experiments. Dissertation, the University of Texas at Austin, 2015. [49]. US Department of Transportation, Federal Highway Administration. Long-term pavement performance program, https://www.fhwa.dot.gov/research/tfhrc/programs/infrastructure/pavements/ltpp/ [50]. Deuflhard P., and Bornemann F. Scientific Computing with Ordinary Differential Equations. 1st edition, Springer-Verlag, New York, 2002. [51]. Brown P.N., Vassilevski P.S., and Woodward C.S. On mesh-independent convergence of an inexact Newton multigrid algorithm. SIAM Journal on Scientific Computing 2003; 25: 570–590. [52]. Zghal S, Bouazizi ML, Bouhaddi N, Nasri R. Model reduction methods for VE sandwich structures in frequency and time domains. Finite Elements in Analysis and Design 2015; 93: 12–29. [53]. Chung J, Hulbert GM. A family of single-step Houbolt time integration algorithms for structural dynamics. Computer Methods in Applied Mechanic and Engineering 1994; 118: 1-11.

46

Highlights We developed an inverse model and mathematical solutions based on Lagrange-multipliertheorem It inverts dynamic viscoelastic properties of heterogeneous structures with displacement observations The developed method improves computation speed and accuracy of existing approaches