Random Response of Nonlinear Continuous Systems Citation Whirley, Robert Gary (1988) Random Response of Nonlinear Continuous Systems. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/ezsr-jt33. https://resolver.caltech.edu/CaltechTHESIS:03142013-090423420 Abstract This thesis presents a technique for obtaining the stochastic response of a nonlinear continuous system. First, the general method of nonstationary continuous equivalent linearization is developed. This technique allows replacement of the original nonlinear system with a time-varying linear continuous system. Next, a numerical implementation is described which allows solution of complex problems on a digital computer. In this procedure, the linear replacement system is discretized by the finite element method. Application of this method to systems satisfying the one-dimensional wave equation with two different types of constitutive nonlinearities is described. Results are discussed for nonlinear stress-strain laws of both hardening and softening types. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Mechanical Engineering Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Mechanical Engineering Thesis Availability: Public (worldwide access) Research Advisor(s): Iwan, Wilfred D. Thesis Committee: Iwan, Wilfred D. (chair) Beck, James L. Caughey, Thomas Kirk Culick, Fred E. C. Defense Date: 23 September 1987 Record Number: CaltechTHESIS:03142013-090423420 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:03142013-090423420 DOI: 10.7907/ezsr-jt33 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 7518 Collection: CaltechTHESIS Deposited By: Benjamin Perez Deposited On: 14 Mar 2013 16:16 Last Modified: 16 Apr 2021 22:28 Thesis Files Preview PDF - Final Version See Usage Policy. 2MB Repository Staff Only: item control page

Random Response of Nonlinear Continuous Systems Thesis by Robert Gary Whirley In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 1988 (Submitted September 23, 1987) ii To my parents ©1987 Robert G. Whirley All Rights Reserved iii Acknowledgments First, I would like to thank my parents for their belief in education, and for their encouragement and support for all of my efforts. Next, thanks to the many people at Caltech who have provided friendship and moral support, including Joe and Heather Humphrey, Mike and Gail Chobotov, and many others. A special thanks to Dennis Parsons for his insightful discussions of finite elements and numerical methods in general, and to my office partner, Paul Smith, for introducing me to the fascinating world of personal computers. Thanks also to Professor T. K. Caughey for many interesting discussions on dynamics and computations. Thanks to Karla for understanding. Finally, a special word of appreciation to my thesis advisor, Professor W. D. lwan, who provided continuous guidance and encouragement throughout the preparation of this thesis. iv Abstract This thesis presents a technique for obtaining the stochastic response of a nonlinear continuous system. First, the general method of nonstationary continuous equivalent linearization is developed. This technique allows replacement of the original nonlinear system with a time-varying linear continuous system. Next, a numerical implementation is described which allows solution of complex problems on a digital computer. In this procedure, the linear replacement system is discretized by the finite element method. Application of this method to systems satisfying the one-dimensional wave equation with two different types of constitutive nonlinearities is described. Results are discussed for nonlinear stress-strain laws of both hardening and softening types. v Table of Contents Acknowledgments iii Abstract iv List of Figures viii List of Tables X 1 Introduction and Background 1 1.1 Introduction . . . . 1 1.2 Problem Definition 4 1.3 Background . . . . 7 1.3.1 Introduction 7 1.3.2 Existing Analysis Approaches 8 1.3.3 Brief Review of Research in Nonlinear Stochastic Analysis . 11 1.4 Overview of Current Work . . . . . . . . . . . . . . . . . . . . . . . 15 2 Formulation 17 2.1 Introduction . 17 2.2 Nonstationary Continuous Equivalent Linearization . 18 2.3 Minimization Equations . . . . . . . . . . . . . . . . 20 vi 2.4 Solution of the Minimization Equations 24 2.5 General Discretization 29 2.6 Nonlinear Terms . . . 33 3 Numerical Implementation 36 3.1 Introduction . . . . . . . . 36 3.2 Finite Element Discretization 36 3.3 Equivalent Linear Parameter Basis Functions 40 3.4 Development of the Covariance Equation 43 3.5 Simplifications for Uniform White Noise Excitation . 45 3.6 Solution of the Covariance Equation for White Noise Excitation. 48 . . 3.7 Solution Algorithm for Nonstationary Response . 4 Application 50 52 4.1 Introduction . 52 4.2 Description of the Physical System 53 4.3 Cubic Hardening Constitutive Nonlinearity 54 4.3.1 Introduction . . . . . . . . . . . 54 4.3.2 Relationship to Discrete Systems 55 4.3.3 Minimization Equations . . . . . 57 4.3.4 Finite Element Matrices and Vectors 61 4.3.5 Performance Evaluation . . . . . . . 64 4.3.6 Description of the Example Problem 64 4.3.7 Effect of Type of Rayleigh Damping 66 4.3.8 Effect of Discretization . . . . . . . . 69 4.3.9 Effect of the Cubic Hardening Nonlinearity 71 vii 4.4 5 4o3o10 Engineering Implications o 73 Arctangent Softening Constitutive Nonlinearity 74 4.4o1 Introduction 74 4.4o2 Relationship to Discrete Systems 75 4.4o3 Minimization Equations o o o o o 77 4.4o4 Finite Element Matrices and Vectors 83 4.4oS Description of the Example Problem 83 4.4o6 Results for the Arctangent Softening Material Nonlinearity 84 4.4o 7 Engineering Implications o o o o o o o 0 0 ooooooooo 0 0 0 0 ooooooooo 0 o o 87 0 90 Concluding Remarks 5o1 Su=ary and Conclusions 5o2 Recommendations for Further Work o o o 0 0 ooooooo 0 0 0 ooo 0 0 ooo 0 oooo 0 0 oooo 0 o 90 0 0 0 93 0 References 113 Appendices 117 A Validation Calculations 117 Ao1 Test Problem One: Stationary Response of a One Element System 117 Ao2 Test Problem Two: Stationary Response of a Four Element System 120 Ao3 Test Problem Three: Nonstationary Response of a Four Element System oooooooo 0 0 0 0 ooooooo 0 0 0 0 0 o 0 0 0 oo 0 0 124 0 Ao4 Test Problem Four: Stationary Response of a One Element System with Softening Constitutive Law o o o o o 0 0 0 0 0 0 0 0 0 oo 0 0 128 0 B Evaluation of the Integral for the Arctangent Softening Nonlinearity 131 viii List of Figures 4.1 Nonlinear Shear Beam Under Base Excitation . 95 4.2 Stress-Strain Curve for the Cubic Nonlinearity 96 4.3 Discrete Analog to the Continuous Shear Beam 97 4.4 Nonlinear Single Degree-of-Freedom Oscillator. 98 4.5 Finite Element Model to Represent a SDOF System 99 4.6 Finite Element Model for 2nd Test Problem . .. . . 100 4.7 Discrete Model for 2"d Test Problem with Hardening Nonlinearity 4.8 Discrete Model for 3rd Test Problem with Hardening Nonlinearity. 102 4.9 Comparison of Mean-Squared Displacements vs. Time for 3rd Test 101 Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 4.10 Finite Element Mesh Used for Shear Beam in Example Problems 104 4.11 Comparison of Normalized Response Profiles for Different Types of Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.12 Comparison of Normalized Response Profiles for Different Mesh Sizes106 4.13 Comparison of Normalized Response Profiles for Varying Degrees of Hardening Nonlinearity . . . . . . . . . . . . . . . . . . . . . . . 107 4.14 Comparison of Free End Mean Squared Displacement Time Histories for Varying Degrees of Hardening Nonlinearity . . . . . 108 4.15 Stress-Strain Curves for Arctangent Softening Nonlinearity 109 ix 4.16 Comparison of Normalized Response Profiles for Varying Degrees of Softening Nonlinearity: Early Times . . . . . . . . . . . . . . . . 110 4.17 Comparison of Normalized Response Profiles for Varying Degrees of Softening Nonlinearity: Late Times . . . . . . . . . . . . . . . . 111 4.18 Comparison of Free End Variance Time Histories for Varying Degrees of Softening Nonlinearity . . . . . . . . . . . . . . . . . . . . 112 X List of Tables 4.1 Modal Damping Ratios for Rayleigh Damping Study . . . . . . . . 69 A.1 Su=ary of Results: RMS Relative Displacements y1 • 123 A.2 Ratio of Nonlinear to Linear Restoring Forces . . . . . 124 1 Chapter 1 Introduction and Background 1.1 Introduction Recent developments in computer technology have opened a whole new world to the engineering analyst. Today, problems are routinely solved which only a few years ago were considered prohibitively difficult. Great advances have been made in both the problem size and problem complexity that can be analyzed. These achievements have prompted a reexamination of the assumptions used in the basic formulation of many problems. Whereas previously limitations in the analysis capability were the governing factor, now the limitations rest primarily in the capability to construct accurate mathematical models for engineering systems. One consequence of this progress is a renewed interest in the treatment of uncertainty in problem specification. This thesis presents a technique for the direct incorporation of randomness in applied loads into the analysis of complex nonlinear engineering systems through the theory of random vibration. Most real engineering problems have uncertainty in their definition. One common source of this uncertainty is the structure geometry. For example, often 2 it is difficult to know whether to model an interface within a mechanical system as in full contact, partial contact, or no contact at all. Another frequent source of uncertainty is in the material properties. Such uncertainties can arise from environmental effects, production variations, or a lack of understanding of a material's constitutive behavior. Perhaps the largest source of uncertainty in many analyses is in the specification of the external loads. This uncertainty often arises when the loads result from a physical mechanism so complex that it is best modeled as a random process in time, or a stochastic process. Simplistic assumptions on the time variation of dynamic loads are frequently made to facilitate problem solution, and are often based on inadequate data. Clearly, there is much to be gained from a more realistic treatment of external excitations that are modeled as random processes. The classical deterministic approach to dealing with stochastic excitations is to assume some average or best estimate functions of time for use in the analysis. At best, a parameter study is conducted to examine the sensitivity of the analysis results to variations in the time history of the applied loading. In contrast, a random vibration analysis acknowledges and quantifies these uncertainties, and thereby provides more useful information to the engineer. The analysis can provide not just one solution, but a spectrum of solutions and their relative likelihoods of occurrence. There are two classes of problems which benefit most from a random vibration approach. The first class contains systems where reliability estimates are essential; where a failure would lead to great loss of life or property. Examples of such systems include nuclear power plants, offshore oil platforms, and aerospace structures. The second class contains problems where uncertainty is particularly 3 high, and thus it is especially difficult to construct a meaningful deterministic problem specification. Examples of problems in this category include wind and earthquake response of structures. For engineering analyses falling into one of these two categories, the extra effort required to perform a random vibration analysis is well justified. The remainder of this chapter begins with a brief definition of the problem to be considered in this thesis. Next comes a general discussion of stochastic analysis methods for nonlinear systems. This is followed by a review of previous developments in the field. Several popular techniques used in practice are discussed, each with its strengths and weaknesses. This chapter concludes with an overview of the current work and its relationship to existing techniques. Chapter 2 presents the concepts of the new method developed in this thesis. First, a precise specification of the problem to be solved is given. This is followed by the mathematical development of the nonstationary continuous equivalent linearization technique. A general discretization for the equivalent linear system is then introduced. Finally, incorporation of the nonlinear terms is discussed, and the range of nonlinearities amenable to solution· using the current technique is examined. In Chapter 3, the numerical implementation of the mathematical model developed in Chapter 2 is described. The finite element method is chosen for the discretization, and an algorithm for computation of the equivalent linear parameters on an element-by-element basis is described. Next, the state-space Liapunov equation for the response covariance matrix is derived. Finally, a solution technique is presented for the nonstationary response statistics. Simplifications which arise in the stationary response case are also discussed. 4 Chapter 4 describes the application of this new procedure to a typical class of problems: systems satisfying the one-dimensional wave equation with a nonlinear constitutive law. Specific minimization equations for two example nonlinearities are derived, and their implementation discussed. Validation calculation results are presented and compared with others in the literature. Physical interpretations of additional results are discussed to highlight the influence of each type of nonlinearity. 1.2 Problem Definition This section describes the general nonlinear continuous system to be examined. Consider the continuous nonlinear system in a domain D(x) described by the PDE \7 · TL(a(x), w) + \7 · ,.N(c(x), W, w) + QL(b(x), w)- m(x)w = p(x,t), (1.1) where w = w(x,t) is the dependent variable representing the "displacement" field, a(x), b(x), and c(x) are "stiffness", "damping", and "nonlinearity" parameter fields, respectively, m(x) is the "mass" distribution in the system, p (x, t) is the externally applied load, TL is the linear "stress" operator on the displacements w and "stiffness" function a(x), 5 TN is the nonlinear "stress" operator on displacements w and/ or velocities w and "nonlinearity" function c (x), QL is a linear "damping" operator on the velocities w and "damping" function b(x), and xis a position vector in the domain D(x). The terms in quotes above have been given physical interpretations in the context of structural dynamics. However, the following analysis is in no way restricted to that interpretation. Indeed, many other problems in diverse fields of science and engineering give rise to equations of the form of (1.1) above, and these may benefit as well from the approach presented herein. However, to provide a concrete physical framework, discussion will center on the concepts from structural dynamics indicated in the definitions. Physically speaking, a(x), b(x), and c(x) represent the material properties of the system under consideration. These material properties are permitted to have spatial variations, so the stiffness, damping, and nonlinearity functions depend on x. It is assumed, however, that a, b, and c are known deterministic functions of x. On every point of the boundary oD(x), let the following sets of homogeneous boundary conditions hold: B( 9 l(w) - 0 (1.2) B("l(w) - 0, (1.3) where B(g) and B(n) are sets of geometric and natural boundary conditions which are linear homogeneous operators containing derivatives along and normal to 6 &D(x). Let the initial conditions be given by w(x, 0) = wo(x) (1.4) w(x,O) = wo(x). (1.5) It is useful to digress here and introduce some notational conventions that will be used in the remainder of this thesis. • Boldface terms, such as TL, w, and x in equation (1.1), represent indexed sets. These may be scalars, vectors, or tensors, depending on the numher of indices required to describe the particular physical quantity under consideration. • Capital Latin letters, such as A and C, represent matrices on finite dimensional spaces or linear operators on infinite dimensional function spaces. • Superimposed dots denote differentiation with respect to time, so w= ~~; Any deviation from these conventions will be clearly explained in the context in which it arises. Many systems of engineering interest can be modeled in the form of equation (1.1). For example, consider the axial vibration of a rod with a cubic- hardening constitutive nonlinearity. The equation of motion of such a system is (1.6) This system is cast into the above form by letting the stiffness parameter a represent the elastic modulus E and the nonlinearity parameter c represent the 7 coefficient E', with the operators defined as (1.7) (1.8) Other examples could include torsional vibration of a nonlinear shaft or large amplitude motions of a cable structure. Different nonlinear stress operators could be used to model softening soil behavior or structural response characteristics. Thus, the general form of equation (1.1) allows treatment of a wide class of engineering problems in a unified approach. It is assumed that the nonlinearity is symmetric. This restriction is nee- essary in order that zero mean solutions can be sought for zero mean inputs. Nonsymmetric nonlinearities often induce a constant offset, or drift, into the system response. In order to avoid clouding the present discussion with needless technicalities, this work will treat only the case of symmetric nonlinearities. It is also assumed that the system is linearizable. For the present purposes, it is sufficient to require that the effect of TN becomes small compared to the effect of TL as the response becomes small. Many engineering systems can be formulated such that this is true, therefore this is not a severe restriction in practice. 1.3 Background 1.3.1 Introduction Analysis of linear engineering systems subjected to stochastic excitation has been a subject of interest and study for many years. As a result, this theory is well understood, and is in wide use in engineering practice today. In contrast, the 8 corresponding theory for nonlinear systems is not nearly as well developed. The classes of nonlinearities are diverse, and the behavior of nonlinear systems is much more complicated than that of their linear counterparts. Specific procedures for handling individual nonlinearities have been developed, but little progress has been made in developing widely applicable methods for implementation into general purpose computer programs. As a result, analysis of nonlinear systems subjected to stochastic excitation has been mostly confined to academic and research laboratory settings. It is hoped that the work described in this thesis will help make stochastic analysis of complex continuous nonlinear systems a viable option for the practicing engineer as well. It is useful to briefly consider the classifications of approaches to analysis of systems subjected to stochastic excitation. First, one could classify approaches based on the problem formulation as either linear or nonlinear. Similarly, the formulation could be classed as either spatially discrete or continuous, depending on whether ordinary or partial differential equations are involved. Alternatively, one could consider the solution technique as either exact or approximate. Finally, the desired solution itself could be divided into stationary and nonstationary categories. The technique developed herein is categorized as an approximate approach for the nonstationary solution of a continuous nonlinear system. 1.3.2 Existing Analysis Approaches For each of the above types of problems, specific analysis methods have been developed. These are described in this section, first for discrete systems, then for continuous systems. Discrete system formulations are the most common for stochastic analysis. 9 Extensive results for single degree-of-freedom (SDOF) models are available in the literature. Small multi-degree-of-freedom (MDOF) models are often constructed in an attempt to extend the SDOF results to more complex problems, with varying degrees of success. These models serve well for problems where the lumping of parameters is straightforward, but for many analyses the construction of a meaningful MDOF nonlinear model is a prohibitively difficult task. Solution approaches for discrete systems can be divided into frequency domain solutions, time domain solutions, and Monte-Carlo simulation solutions. Each of these methods is described below. Frequency domain solutions are widely used in the analysis of linear problems. They easily handle complicated excitation spectra and stationary response. Extension to nonstationary response is possible but somewhat complex. Strictly speaking, nonlinear problems are not amenable to frequency domain solution because the technique relies on superposition. Extensions to the theory enable approximate solutions in cases where the nonlinearity is "small." Time domain solutions are used for most nonlinear problems, and for some linear problems with simple excitation spectra. Time domain methods are especially convenient when the nonstationary response is of interest, since ensemble averaging is typically used for the nonstationary statistics. Time domain solutions become quite involved for complicated excitation spectra, however. The most powerful technique available for solution of random vibration problems is direct Monte-Carlo simulation. In this procedure, the digital computer is used to generate a particular realization of the stochastic input, and the deterministic solution is found to this single input using classical analysis techniques. This process is repeated until the ensemble of realizations of the stochastic 10 input represent the statistical distribution to a desired degree of accuracy. Then, statistics are computed across the ensemble of responses obtained from these inputs. Although this method can, in theory, be used to solve almost any random vibration problem, it is very costly in terms of computational resources. In addition, it provides limited insight into the behavior of the system under a different stochastic input. Continuous system formulations have not enjoyed the popularity associated with discrete system formulations for stochastic analysis. This is largely because solutions for continuous systems are much more difficult than for discrete systems. However, many engineering problems are of a continuous nature and do not physically suggest a discrete model. Examples of these problems include aircraft wings and fuselage structures subjected to strong acoustic excitation. Techniques available for stochastic analysis of continuous systems include analytical continuous solution, discretization, and Monte-Carlo simulation. A summary of each of these methods is given below. Analytical continuous solutions are quite difficult to obtain, since partial differential equations replace the more tractable ordinary differential equations of discrete systems. Consequently, few exact solutions are available, and those which do exist are predominantly for linear systems, simple geometries, and white noise excitation. The available exact solutions, while restricted in their applicability, do serve a valuable role as benchmarks for the emerging array of approximate solution methods. An effective approximate approach to solution of the continuous problem is to convert it into a discrete problem. This conversion is classically based on modal decomposition for linear and mildly nonlinear problems. Recent studies 11 have also used finite difference and finite element techniques for the discretization. This approach works very well for linear problems, but the effective discretization of nonlinear continuous systems remains an area of active research. Monte-Carlo simulation can also be used in a varied form for continuous systems. If closed-form deterministic solutions are available, then they can be implemented for calculations on a digital computer. Alternatively, a continuous system model can be implemented on an analog computer and simulation performed with no knowledge of the analytical form of the solution. Monte-Carlo simulations for continuous systems possess the same advantages and disadvantages as described above for discrete systems. 1.3.3 Brief Review of Research in Nonlinear Stochastic Analysis Most methods for the stochastic analysis of nonlinear systems, both discrete and continuous, can be divided into two broad classes: those· which derive a Fokker-Planck-Kolmogorov equation for the transition probability density of the response, and those which deal directly with the stochastic differential equation. The differential equation methods can be further divided into perturbation methods and equivalent linearization methods. In what follows, each of these topics will be examined in more detail. In many cases of practical interest where the excitation is limited to Gaus- sian white noise, it can be shown that the response is a Markov process, and the related Fokker-Planck equation may be derived from the original stochastic differential equation [1,2,3]. As an example, consider the following nonlinear structural 12 response problem in terms of a second order equation of motion: (1.9) where x is a generalized displacement, f(·, ·) is the generalized restoring force which is a nonlinear function of x and/or ~~, and q(t) is stationary Gaussian white noise with constant power spectral density So. Then, assuming a stationary response, the associated Fokker-Planck equation takes the form . ap a . a2p Xax - ax [pf(x, x)]- So ax2 = 0, (1.10) where p(x, :i;) is the joint probability density function of the displacement and velocity response. Successful solution of this equation yields the transitional probability density function of the response process. Caughey [4] has obtained a solution to this equation for certain classes of nonlinearities. Yong and Lin [5] have recently obtained a more general solution applicable to a larger class of nonlinearities. In contrast, general solutions to the nonstationary version of this equation have been found only for linear systems. While many linearization techniques are in use in modern analysis, only two lend themselves well to stochastic analysis: perturbation and equivalent linearization. Each of these approaches will be discussed in some detail in the following paragraphs. Perturbation techniques were applied to random vibration problems by Crandall in references [6,7]. These methods are useful when the nonlinearity in the governing equations is small. This scaling is frequently expressed in the form of a small parameter (the "perturbation parameter") multiplying a nonlinear function. The basic idea in perturbation theory is to expand the solution about the linear solution in a power series in the small parameter. Analysis then 13 proceeds by matching powers of the small parameter, thus yielding a sequence of linear equations to be solved. For example, consider the following equation for a single degree-of-freedom oscillator: X+ 2~w..x + w!x + E17(X, x) = F(t), (1.11) where 17 is a nonlinear function and Eisa small parameter. Assume a solution to this equation as a power series in e (1.12) Substituting this assumed form of solution into the governing equation and matching powers in e yields a set of linear equations for x 0 , Xi> ••• of the form (1.13) (1.14) The solutions to these equations may be combine<;! to compute various statistical properties of the original nonlinear equation. For example, the mean of x(t) is given by E [x(t)] = E [xa(t)J + eE [x1(t)] + · · · (1.15) and the autocorrelation (to first order in e) is E [x(t + r)x(t)] = E [xa(t + r)x0 (t)] + eE [xa(t + r)x 1(t)] +eE [xa(t)x 1 (t + r)]. (1.16) Typically, only the zeroth and first order solutions are calculated. While conceptually simple, this approach can be algebraically cumbersome for continuous systems, and does not easily lend itself to numerical implementation. 14 Equivalent linearization for stationary random vibration was first investigated by Booton [8] and Caughey [9]. !wan and Mason [10] extended the technique to the nonstationary response of discrete systems. Krousgrill and !wan [11] further extended the ideas to the steady state response of continuous systems. The basic idea in equivalent linearization is to replace the given nonlinear system with an "equivalent" linear system. The parameters of the equivalent linear system are chosen to minimize some difference between the original nonlinear system and the equivalent linear system. This procedure leads to expressions for the equivalent linear parameters in terms of the system response statistics. To illustrate the application of discrete equivalent linearization, consider the nonlinear single degree-of-freedom described by the equation x + f(x, :i:) = F(t). (1.17) Assume that an approximate solution can be obtained from the linearized equation x + b.:i: + k.x = F(t). (1.18) The error of linearization is a random process described by e: = f(x, :i:) - b.:i:- k.x, (1.19) which is the difference between equations (1.17) and (1.18). Some measure of this error is then minimized to generate equations for b. and k. in terms of the response statistics. A frequent choice is to minimize the mean square error, thus requiring (1.20) (1.21) 15 Manipulating these equations and assuming that the response process is stationary leads to the following expressions for the equivalent linear damping and stiffness: b _ E[:i;f(x,:i:)] •- E [:i;2] k _ E [xf(x, :i:)] •E [x 2 ] • {1.22) (1.23) Stationary solutions usually require iterative determination of the equivalent linear parameters, while nonstationary solutions avoid iteration if an explicit time integration procedure is used. Equivalent linearization holds two main advantages over perturbation methods. First, equivalent linearization is applicable to more general classes of nonlinearities than are perturbation methods. In particular, hysteretic systems are very difficult to analyze using perturbation techniques. In contrast, many investigators [12,13,14] have successfully applied equivalent linearization techniques to a wide variety of hysteretic systems. Second, the equivalent linearization procedure is better suited to numerical implementation. These observations motivate the use of equivalent linearization in the present effort. 1.4 Overview of Current Work In view of the extreme difficulties in obtaining general solutions to the governing nonlinear equations, a linearization approach is adopted in this work. The original continuous nonlinear system is approximated by a continuous linear systern with time varying coefficients. This replacement is accomplished using a nonstationary equivalent linearization technique developed in this thesis. The equivalent linear partial differential equations are then discretized using the finite 16 element method, and the resulting equations are cast into a first-order Liapunov covariance equation. Finally, a time integration procedure is presented to solve for the evolution of the system covariance response. 17 Chapter 2 Formulation 2.1 Introduction As is evident from the discussion in the previous chapter, there is considerable room for improvement in general methods for obtaining the nonstationary response of nonlinear continuous systems. Exact solutions are rare, particularly for nonlinear problems. Also, analytical closed-form solutions are almost always restricted to simple geometries, making them of limited value for direct application to engineering problems. Monte-Carlo simulation escapes these limitations, but at the expense of considerable computational effort. Thus, an improved approach for obtaining the random response of nonlinear continuous systems is desired. In view of the difficulties in obtaining general solutions for the stochastic response of nonlinear continuous systems, a linearization approach is adopted in this thesis. The original continuous nonlinear system is approximated by a continuous linear time-varying system. This replacement is accomplished using a nonstationary continuous equivalent linearization technique developed in the following sections. 18 2.2 Nonstationary Continuous Equivalent Linearization In this section, the general nonlinear continuous system to be solved is restated. Then, an equivalent linearization procedure applicable to a nonlinear, continuous stochastic system under nonstationary response conditions is derived. Equivalent linearization of a continuous system was first studied by Krousgrill and Iwan [11,15] in the context of steady-state deterministic analysis and stationary random response. In what follows, the technique will be extended to include cases where the excitation and response are assumed to be nonstationary random processes. Consider the continuous nonlinear system in a domain D(x) described by the PDE V · TL(a(x), w) + V · rN(c(x), w;w) + QL(b(x), w)- m(x)w = p(x, t), (2.1) where all terms are as defined in section 1.2. On every point of the boundary oD(x), let the following sets of homogeneous boundary conditions hold: B( 9 l(w) - 0 (2.2) B("l(w) - 0, (2.3) where B(g) and B(n) are sets of geometric and natural boundary conditions which are linear homogeneous operators containing derivatives along and normal to oD(x). Let the initial conditions be given by w(x,O) = wo(x) (2.4) 19 w(x,O) = wa(x). (2.5) The objective of continuous equivalent linearization is to replace the nonlinear PDE, equation(2.1), by some member of a class of linear PDE's for which solutions are easily obtained. The linear replacement equation is chosen such that some difference minimization criteria is satisfied. Since the nonstationary response is of interest, the class of continuous linear systems considered must contain time-varying coefficients. Thus, the equivalent linear system is taken to be of the form '\7 · TL(a, w) + '\7 • TL( a(x, t), w) + '\7 • TL(f3(x, t), w) + QL(b, w) - m(x)w = p(x, t), (2.6) where a(x, t) and f3(x, t) are sets of (unknown) equivalent linear parameter fields, and all other terms are as previously defined. Equation (2.6) is obtained from (2.1) by replacing the nonlinear stress operator rN(c, w, w) with the sum of an equivalent linear stiffness operator rL( a, w) and an equivalent linear damping operator rL(,B,w). All other terms in equation (2.1) are already linear, and therefore are retained in their original form. It is apparent from the form of the equivalent linear equation (2.6) that a characterizes distributed stiffness properties, while f3 represents distributed damping properties of the nonlinear system. Note that the equivalent linear damping is limited to strain-rate related damping. Note also that, even though the equivalent linear parameters are treated as functions of time, their explicit dependence is on the nonstationary response statistics rather than on time. That is, the expressions for the equivalent linear parameters involve the current response statistics (of the linear system), and since the response is considered nonstationary, these response statistics vary with time. 20 It frequently happens in engineering that the region over which loads act is known or easily determined from geometric considerations, while the time history of the applied loading is much more uncertain. Therefore we will herein restrict attention to problems of this class by assuming the applied forcing function p(x, t) to be deterministic in its spatial variation but random, Gaussian, zero mean, and nonstationary in time. Many real excitations result from the combination of a number of different physical processes. According to the central limit theorem of probability theory, under suitable conditions, the sum of a large number of different statistical distributions tends toward a Gaussian distribution. Hence, assuming a Gaussian time dependence for the loading seems like a good first approximation. In addition, many processes of engineering relevance are of sufficiently short duration that they must be considered transient, or nonstationary.. In order to acco=odate these phenomena the statistical parameters of the"Gaussian process are allowed to vary with time. 2.3 Minimization Equations With the form of the equivalent linear system now specified, all that remains is to specify a procedure for determining the equivalent linear parameter fields a(x, t) and ,B(x, t). Since the excitation and response are nonstationary random processes, derivations are cast in terms of a typical sample function, or particular realization, of each random process. Nonstationary statistics are then computed by taking ensemble averages across a collection of sample functions. Let kp(x, t) represent a sample function of the nonstationary stochastic process p(x, t). The 21 sample functions of the response •w and •w then obey V. rL(a, •w) + V. rN(c, •w, •w) + QL(b, •w)- m(x)•w = "p(x,t) (2.7) on 1/(x), and B(g)(kw) - 0 (2.8) B(n)(kw) - 0 (2.9) •w(x,O) - wo(x) (2.10) "w(x, 0) - wo(x). (2.11) on Bll(x), with deterministic initial conditions The approximations to the sample functions obey the equivalent linear equation V. rL(a,"w) + V. rL(a(x, t),"w) + V · rL(,B(x, t), •w) + QL(b,"w) - m(x)•w = •p(x,t) (2.12) on 11 (x), with the same boundary and initial conditions. Strictly speaking, one should use a different notation for the approximate solution obtained from the linearized system and the exact solution to the original nonlinear system. To avoid clouding the discussion to follow with additional notation, this distinction is not made in the notation used in this work. The reader should keep in mind, however, that this distinction does exist. Situations where confusion arises as to which interpretation is intended are clarified in the surrounding text. In order to choose the equivalent linear parameter fields a(x, t) and ,B(x, t), some measure of the error in replacing the original nonlinear system by an equivalent linear system must be developed. Define the "stress difference" function 22 ~(x,t) for the kth sample function as (2.13) This choice of difference definition is motivated by its consistently superior performance for steady-state deterministic problems (Krousgrill [15]) and its reduced continuity requirements when compared to equation differences or strain energy differences. Clearly, the optimum linear replacement system is the one which generates the minimum error. Thus, determination of the equivalent linear parameter fields requires minimization of the difference k ~ (x, t) between the nonlinear system and the equivalent linear system. In general, exact spatial minimization on a point-bypoint basis will not lead to a solution, so an approximate approach is developed below which minimizes an average measure of the difference. Define a spatial averaging operator G.(u) on a continuous field u by G.(u) = { JO(x) u(x, ·) dD(x). (2.14) With this definition, it is easily shown that: • G.(u + v) = G.(u) + G.(v) (Linearity) • G.(u 2 ) > 0 'Vu # 0 and G.(o) = 0 (Positive Definiteness). Using this operator, define a spatially-averaged measure of the stress difference as a function of time by k17(t) = G.(k ~. k~). (2.15) Then, k17 (t) is the spatially averaged value of the norm of the instantaneous stress difference corresponding to the sample function kp(x, t) of p(x, t). Finally, taking 23 the ensemble average of equation (2.15) yields (2.16) It is clear that E (17(t)j provides a measure of the total "stress difference" at any timet resulting from replacing the original nonlinear system (2.1) by the equivalent linear system (2.6). The equivalent linear system will now formally be defined at any time t as that system which causes E [17] to attain its minimum value on the class of solutions to the equivalent linear system. That is, the stress difference ~ is evaluated for solutions of the linear system. Thus, it remains only to find the time and space varying functions a(x, t) and .B(x, t) such that E[17]= minimum Va,j3. (2.17) Using equation (2.16) along with the properties of the expectation operator, (2.17) can be written E [Gz(~ · ~)] =minimum. (2.18) Necessary conditions for the extremization of E[Gz(4 · 4)] are aE [Gz(~ · ~)] aa aE [Gz(~ · ~)] a.a - o, (2.19) - (2.20) 0. Note that these equations are highly symbolic, since a and f3 are really functions, and a complete rigorous argument requires a calculus of variations treatment. The above equations are sufficient for the present purposes, however. Sufficiency of these conditions is established by showing that equations (2.19) and (2.20) yield a minimum, and not a maximum, forE [Gz(~ · ~)]. This result will be established in the next section. 24 Invoking the linearity of the expectation operator and inverting the order of differentiation and expectation, (2.19) and (2.20) become [:a G.(~.~)] _ 0, E[:,BG.(~·~)] - E 0. 2.4 (2.21) (2.22) Solution of the Minimization Equations In view of the fact that a and ,8 are functions of both space and time, general analytical solutions to these minimization equations are not available. As an alternative, a computationally efficient approximate approach for solution of the minimization equations (2.21) and (2.22) is developed below. Approximate the equivalent linear parameter fields a(x, t) and .B(x, t) by the following expansions: M a(x, t) ~ I: 1/>;(x)&;(t) (2.23) i=l M ,B(x, t) ~ 2: 1/>;(x),B;(t) , (2.24) i=l where &;(t) and ,B,(t) are unknown functions of time, and 1/>;(x),i = 1, ... ,M, are linearly independent known functions of x. Using these expansions (2.23) and (2.24), it follows that M TL(a(x,t),w) - TL(:l: 1/>;(x) &;(t), w) (2.25) i=l M - 2: TL(t/>;(x), w)&;(t) (2.26) i=l M - I: T,L(w)a;(t), i=l (2.27) 25 and M TL(p(x,t),w) - TL(I>Mx)fW),w) (2.28) i=l M - .L:>L(cf>;(x),w)P;(t) (2.29) _ i=l M " A L..J' T;L (w)P;(t). (2.30) 1=1 The stress difference A. can now be written in terms of the time-varying coefficients &; and '/:J, as A. = A. (&, ,8, t) (2.31) where (2.32) (2.33) Using the expansions for a and j3 yields the revised form of the Euler equations for the minimization as aE[G.{A.·A.(&,,B,x,t)}) _ a&, aE 0 [a. {A.. A.(&,_B,x,t) }] = O,i=1, ... ,M. a{:J, (2.34) (2.35) As discussed in the previous section, the above equations are necessary, but not sufficient, for minimization of E [G. (A.· A.)]. Following the line of argument of Spanos in [16] and Krousgrill in [15], it is demonstrated below that equations (2.34) and (2.35) do actually produce a minimization of E [G. (A.· A.)]. Let &k and '/:Jk, fork = 1, ... , M, satisfy the minimization relations above. Define another set of parameters ak and Pk by (2.36) (2.37) 26 and let .A(a,,B) (2.38) ..o. - .A( a, ,8). (2.39) ..i - Now expand E [G:(..i · ..i)] in a Taylor series about a and ,8: E(G.(..i·..i)j = E[G.(.A·.A)J M a M a +~a&; E [G,(.A ·.A)] Sa; • + L -. E [G.(.A . .A)] 8(3, i=l a(3, l MM a2 +-I: I: a· -a&E[G.(..o. . ..o.)J sa,s&j 2 i=l i=l a, ' a2 • • E [G.( .A . .A)] 8(3, 8(3i 2 i=lj=l a[3,af3i . MM a2 . • + L: I: . . E [G.(..o. . .A)J s&,of31• i=li=l aa,a{3i 1 M M +- L L A A (2.40) But, since E [G.(.A ·.A)] satisfies equations (2.34) and (2.35), equation (2.40) can be simplified to E(G.(~·~)j - E[G.(.A·.A)J a2 + 2~[; a&,aai E [G.( .A . ..o.)J s&,s&1 1 M M 1 MM a2 +-I: I: . . E [G.( .A· ..o.)J s,B, s,81 2 i=li=l a{3,af3i M M a2 • +I: I: . . E [G,(..o. · ..o.)J s&,of31. i=l i=l aa,a{3i (2.41) Now, recalling that .A is linear in ak and {3k, the higher order derivatives in equation (2.41) can be expressed as a2 E [G,(..o. . .A)J = 2E [G, (a..o. a..o.)] a&,aai a&, . a&i (2.42) 27 (2.43) .a2. E [G.(..o. . ..o.)J = 2E aa;af3; [a. (aa~. a(3; a~)] . a; (2.44) Using equations (2.42) through (2.44) in (2.41) gives (2.45) Next, let two functions be defined as follows: M a..o. M a..o. • - I: a· oa, i=l a; - I:. .5(3,, •=1 a(3, (2.46) (2.47) then equation (2.45) can be written (using the linearity of E [ ·] and G.(·)) E [G.(A · A)j = E [G.(..O. · ..0.)] + E[a. ((ht(oa) + h2(o,8)) · (ht(oa) + h2(o,8))) ]. (2.48) The positive definiteness property of G.(·) requires that (2.49) and therefore the desired result (2.50) 28 follows. This relation assures that the extremum produced by equations (2.34) and (2.35) is indeed a minimum. Up to this point a and {3 have been considered to be sets of unknown parameters. To expedite development of specific minimization equations and simplify notation, assume that a and {3 are sets of one parameter each, i.e., scalars a and /3. The methodology is in no way restricted to this case. Interchanging the order of expectation, differentiation, and spatial averaging in (2.34) and (2.35) and invoking the linearity of Gz gives E[az(~·::)] - 0 (2.51) E[az(~·!A)] = 0. (2.52) Using the definitions (2.27) and (2.30) in (2.51) above gives E [az {(ETf(w)&k(t) + ETf(w)Pk(t)- TN(c,w,w)) · (T,L(w)) }] = 0 (2.53) or rearranging, E [az {ET1L(w) · Tf(w)&k(t)}] + E [az {ET,L(w) · Tf(w)Pk(t)}] = E[az{TN(c,w,w)·T,L(w)}]. (2.54) Equation (2.54) can be written in matrix form as Ail+ c,B = G, (2.55) A;i - E [az hL(w) · Tl(w) }] (2.56) C;i - E [az {T1L(w) · Tl(w) }] (2.57) E [az {TN(c,w,w) · T,L{wJ}] (2.58) where G; - 29 Using definitions (2.27) and (2.30) in (2.52) gives E [a. {(E_rf(w)&k(t) + E_rf(w)tJk(t)- rN(c,w,w)) . (r,L(w))}] =o (2.59) or E [a. {~r/(w) · rf(w)&k(t)}] + [a. {~r/(w) · rf(w)~k(t) }] E E[a.{rN(c,w,w)·r;L(w)}] = (2.60) In matrix form, (2.60) becomes D&+B~=H (2.61) D;,- - E [a. {rl(w) · r/(w) }] (2.62) B;t - E [a. {T;L(w) · T;L(w) }] (2.63) H; E [a. {rN(c, w,w). r;L(w) }]. (2.64) where - Note that C = DT. If all stochastic processes are assumed to be stationary, then the displacement-velocity covariances vanish and it is can be shown that C = D == 0, and thus the two minimization equations uncouple. This result does not carry through for the nonstationary case, however. 2.5 General Discretization In order to develop the minimization equations in a form suitable for numerical implementation, the solution is spatially discretized in terms of deterministic spatial basis functions and stochastic time coefficients. In addition to allowing 30 digital computer implementation, spatial discretization allows solution of prob!ems involving complex geometries which are intractable by other methods. To simplify the notation, assume that the dependent variable, w(x, t), is just a scalar, w(x,t). Equations similar to those which follow are easily developed for cases where w is a vector. Let w(x, t) be approximated using known deterministic spatial basis functions ,Pk(x) and unknown stochastic time coefficients uk(t), k = 1, ... , N. Represent this discretization by N w(x, t) = L ..Pk(x)uk(t). (2.65) k=l This general form accommodates both global discretizations, such as used in modal analysis, and local discretizations, such as used in finite element analysis. Each of these techniques will be discussed further in a subsequent section. Before casting the minimization equations in terms of the discretized solution, it is helpful to note some applicable results. Claim 1 If R and S are linear spatial operators on w, then equation {2.65} implies N E [R(w)] - L R(,Pk)E [uk(t)] (2.66) k=l N N E[R(w) · S(w)] = L:L:R(,Pk) · S(..Pz)E[uk(t)u1(t)]. (2.67) k=ll=l In addition, since w has zero mean in time, N E[w] = L ..Pk(x)E [uk(t)] = 0. (2.68) k=l Furthermore, since the basis functions ..Pk are linearly independent, E[uk(t)] = 0, k = l, ... ,N. (2.69) 31 The following theorem was first stated by the Russian mathematician Kazakov in reference [17], and was popularized in structural dynamics by Atalik and Utku in reference [18). Claim 2 Consider a single-valued function of n variables (2. 70) where • y is a jointly Gaussian random vector process with zero mean, • q(y) is sufficiently smooth so that first partial derivatives with respect to y;, i = 1, ... , n, exist, and " • Jq(y)J < Aexp 'I:,yj, a< 2, for some arbitrary A and any y. ;:1 Then, [YYT] E [Vq(y)], (2.71) a ' ... ' ay,. a] . v T = [ayl ' avz (2.72) E [yq(y)] = E where It is noted that the third requirement in the above claim simply limits how fast q can grow as y grows. As will become clear in the next section, y plays the role of the response state vector, and q plays the role of the nonlinear stress operator. Thus, roughly speaking, the third condition above says that the dependence of the stress on a norm of the state vector should be weaker than exponential growth. This restriction is satisfied by the example nonlinearities discussed in chapter 4, and by many other nonlinearities of engineering interest. 32 Finally, recall that for Gaussian random variables, higher order moments can be expressed in terms of the 2nd order moments [3]; for example: With these results, discrete minimization equations can now be developed by combining the approximation (2.65) with the general minimization equations (2.55) through (2.64). Since it was established that D = CT in the previous section, no explicit expressions for D are given. It is understood that D is easily constructed when needed. Using (2.65) in the definition of T,L(w) yields N T;L(w) = T;L(~= 1/J•u•) k=l N - L Tp('if;.)u.(t). (2.74) k=l Using (2. 74) above in the definition (2.56) gives A;; = - E [a. {(E Tl(.P.)u.(t)) · (~ T,L(.Pl)ul(t))}] E [r JD(x) (t I>l(.P.). T;L(.Pl)) dl'(x) uk(t)ul(t)] k=ll=l N N - L L:E [u.(t)ul(t)] rD(x) T;L(.p.). T/(.Pl) dl'(x). k=ll=l (2.75) } Define (2.76) then equation (2.75) becomes (2. 77) 33 Similarly, using (2.74) in (2.57) gives (2.78) Following the same procedure, using (2.74) in (2.63) gives (2.79) Using the discretization in the expressions for the right-hand-side vectors (2.58) and (2.64) yields G; = E[a.{rN(c,w,ui)·(Erl(tfJk)uk(t))}] N - k=l L lo(x) { E [rN(c,w,ui) · r/(tflk)uk(t)] dD(x) H; - E [a. { rN(c,w,ui) · (E (2.80) r/(tflk)uk(t))}] N - k=)D(x) L r E[rN(c,w,ui)·r/(tfJk)uk(t)]dD(x). (2.81) The above two expressions containing the nonlinear stress operator must now be written in terms of the response covariance matrices. This will be described in the next section. 2.6 Nonlinear Terms In order to calculate the equivalent linear parameters & and ~. it is necessary to express the minimization equations in terms of the response covariance matrix and known quantities. The equations for the A, B, C, and D matrices are already in that form, and need no further manipulation. However, the right-hand-side vectors, given by equations (2.80) and (2.81), contain the nonlinear stress operator 34 TN acting on w and w directly. Following is a procedure to express these vectors in terms of the desired response statistics. First, write equation (2.71) in components as N [7}: aq] . E[y;q(y)) = .L;E[y;y;)E i=l y, (2.82) Recalling the discretization for w and w, it is clear that the dependence of TN on w and w is through the nodal response vectors u and u. This dependence is shown explicitly in the notation (2.83) Now, applying the component form of the Atalik and Utku result (2.82) to equation (2.80) yields G; =f.k=l { l=lf.E [ukuz) E [aa (TN(.,P, u, u) · T,L(~j~k))] j O(x) U! +f:E[ukul)E[a~ (TN(.,P,u,u)·T/('if!k))] dZJ(x) l=l (2.84) U! and similarly equation (2.81) becomes H; = f. { f. E [ukul) E [aa (TN ('if!, u, u) · T;L('if!k) )] k=l j O(x) l=l U! + f.E [ukud E [a~ (TN(.,p, u, u). T/('if!k))] dZJ(x). l=l (2.85) U! Performing the indicated differentiations and noting that T;L (Tjik) is not a function of u or u yields the simplified forms 35 and (2.87) The nonlinear stress operator can frequently be expanded as a power series of u and u, and thus the derivatives of TN can be expressed in terms of higher order moments of u and u. Then, using the property of Gaussian random variables given in (2.73), these higher order moments can be expressed in terms of the 2"d order moments of u and u. Thus, the G and H vectors are expressed in terms of the desired response covariance matrices. 36 Chapter 3 Numerical Implementation 3.1 Introduction In the previous chapter, an algorithm was developed to replace the original nonlin- ear partial differential equation with an equivalent linear time-varying equation. Thus, solution of the original nonlinear system has now been reduced to the solution of an equivalent linear system at each time step. In this chapter, specific procedures are developed for obtaining that solution. 3.2 Finite Element Discretization Thus far, minimization equations have been derived using a general discretization. These equations will now be specialized using the finite element method. The finite element method is widely used to solve solid and structural mechanics problems. The popularity of this technique is due to its ability to handle complicated geometries and inhomogeneous material properties, and its ease of implementation in a general purpose computer program. The ability 37 to solve problems with complicated boundaries stems naturally from the local nature of the interpolations within each element. Also, these local interpolations make it easy to have different material properties in different elements, thus permitting great flexibility in modeling material inhomogeneity. Finally, since complicated problems are treated as an assemblage of simple "elements," the finite element method naturally lends itself to solution by digital computer. These characteristics of the finite element method make it an effective choice for solution of the equivalent linear system, as is described below. Let the domain D(x) be subdivided into N elements J)(•l(x) such that N D(x) = L D(kl(x), (3.1) k=l where J)(kl(x) is the domain of the kth element, and N is the number of finite elements in the domain. Further, let w<kl(x,t), x E J)(kl(x), denote the restriction of w(x,t) to J)(kJ(x), so w<•J (x, t) - w(x, t) Vx E J)(k) (x) Vx rt J)(kl(x). 0 (3.2) Then, from equation (3.1), the solution field is represented as the summation of all of these restrictions, N w(x, t) = L w<kl(x, t). (3.3) k=l Within an element J)(kl(x), let the dependent variable w<•l(x,t) be expressed as the product of spatial basis or shape functions N?l(x) and nodal values u}k)(t) as NEN w<kl(x,t) = L Nj"l(x) u;"l(t), i=l (3.4) 38 where NEN is the number of nodes per element, N?) denotes the finite element shape function for node j of element k, and u)k) refers to the value of wk at the j'" node of the k'" element; this node number is therefore a local one for element k. Thus, the discretization of w(x, t) can be expressed as N N NEN k=l k=l i=l w(x,t) = l::w(kl(x,t) = L L N?)(x) u)kl(t). (3.5) Note that the shape functions N?)(x) are deterministic functions of space and that the nodal values u)k) (t) are nonstationary stochastic processes. The spatial basis functions NY) (x) should satisfy the convergence criteria for finite element shape functions. The sufficient conditions for convergence of the finite element method are described in the book by Hughes [19] and are not repeated here. Special elements have been derived for the solution of 4'" order equations which do not satisfy the convergence criteria but do perform well in practice. Although these criteria are sufficient but not necessary conditions for convergence, special care should be taken when using elements not satisfying the convergence requirements. The finite element methodology described above may now be incorporated into the minimization equations developed in the previous chapter to cast those equations in a more specific form. Using the finite element approximation (3.5) in equation (2.77) and rearranging gives N N A;i = L L NENNEN L L E [ulk)(t)u~"'l(t)] k=l m=l l=l { p=l TL(N/k)) · Tf--(N("')) dD(x). jf)(x) 3 • P (3.6) The above equation can be simplified considerably through use of two 39 observations. First, recall that by definition (3.7) that is, the shape functions for element k vanish everywhere in the domain outside of element k. Second, assume that the chosen reference configuration is stress-free, so {3.8) With these two points in mind, it is now clear that the kernel of the integral in equation {3.6) is nonzero only when N/k) and NJ"'l belong to the same element, i.e., when k = m. Thus, the integral in {3.6) can be replaced by its restriction to element k and the summation over m dropped to get N NENNEN A-. = " ' " ' " ' E(k)(t) L..J L..J •1 L..J k=l 1=1 p=l lp r rL(.M(k)). r.L(N(k))d[)(k)(x) lo<•>( l 1 1 , P • (3.9) X where Ezp is defined as in equation (2.76). Thus, the calculation of the matrices in the minimization equations has been significantly simplified. Applying a similar argument to equation {2.78) and (2.79) and using (3.5) gives {3.10) and N NENNEN B;i = L L L E~k)(t) r k=l 1=1 p=l lp r!-(N/k)). TL(N(k))dD(k)(x). JD<•>(x) 1 ' {3.11) P Incorporating the finite element approximation {3.5) and the observations (3. 7) and {3.8) into the equations for the right-hand-side vectors G and H gives G; = L L L E~;l(t) lD<•l N NEN NEN k=l r=l p=l rhN~k)) · E (x) [---"W arN ] au, dlJ(kl(x) l [JTN] L E;~l(t) <•> rl(N~kl) · E ---:(k) dD{kl(x) k=l r=l p=l D (x) au, N NENNEN +L L {3.12) 40 and (3.13) Comments made in the final section of the previous chapter regarding the nonlinear stress operator apply here also. 3.3 Equivalent Linear Parameter Basis Functions With the minimization equations now specialized using the finite element method, it remains only to make a specific choice for the basis functions used in the approximation of the equivalent linear parameters. Recall the expansions for the equivalent linear parameters: M a(x,t) - I:;</>;(x)&;(t) i=l M ,B(x, t) = I:; </>;(x)_B;(t). i=l Since the replacement linear system is discretized using the finite element method, choose the equivalent linear parameters to be constant within an element. This choice allows the approximation of the equivalent linear parameters to be changed in the same way as the approximation of the equivalent linear solution through refinement of the finite element mesh. Thus, using more elements in a given area increases the accuracy of both the solution and the equivalent linear parameters. This approach is consistent, since sharp gradients in the equivalent linear parameters would be expected to produce sharp gradients in the solution, thus requiring a finer mesh. 41 With this choice, the basis functions ,P;(x) take the form 1.0 \/x E p(il(x) ,P;(x) = 1 (3.14) 0.0 Vx (/_ p{il(x), and M = N = number of elements. The parameters &; and fj, now represent the equivalent linear stiffness and damping coefficients for element i. Note that for this choice of basis function it is easily shown that ,P;(x) and ,P;(x) are linearly independent. It turns out that this choice of basis function leads to a considerable sim- plification of the minimization equations. To see this, recall that, by definition, TL(N(k)) -- TL("'· i l '+',, N(k)) l • (3.15) If the dependence of the operator TL on its parameter is restricted to be such that TL(o, w) = o, (3.16) then since the previous assumption of a stress-free reference configuration gave (from equation 3.8) (3.17) it follows that equation (3.14) implies T/(N?l) - 0 Vi I k _ TL( ,p,, N/')) for i = k. (3.18) (3.19) Thus, T/(Nik)) can be replaced by .,,L(Ni'l) in equations (3.9) through (3.11) and the sum over k dropped without loss of generality. Noting that this argument applies to both stress operators in the integrands, it follows that A,; = B;; = C,; = o for i I i- (3.20) 42 Thus, the particular choice of basis function given above greatly improves the computational efficiency of the minimization procedure. Simplifying equations (3.9) through (3.11) using (3.20) above yields the final form of the minimization equations: NENNEN A··= ~ ~ E(i)(t) ~ " ~ l=l p=l lp NENNEN B-· = ~ ~ E!il(t) 11 ~ l=l ~ p=l lp NENNEN C·· = ~ ~ E!il(t) ~ " l=l ~ p=l lp r TL(N(i)). rL(N(i)) df)(il(x) 1p(il(X ) ' I ' P (3.21) TL(N(i)) df)(il(x) 1rp(il(X ) TL(N(i)). I P (3.22) r (3.23) 1 1 TL(N(i)). TL(N(i)) df)(il(x) 1o<•>(Xl • I • P Similarly, the only nonzero contribution to the right-hand-side vectors given by equations (3.12) and (3.13) occurs when i = k, so these equations become NENNEN G· = ~ ~ E(il(t) { • ~ ~ r=l p=l pr [aN] aUr(•) r.L(N(i)) · E [aN] ~ d!)(il(x) rL(N(;)) · E _r__ df)(il(x) 1o<•>(X) • + NENNEN ~ ~ E(i) (t) { ~ ~ pr 1p(il(x) • r=l p au· r(•) p p=l (3.24) and NENNEN H·1 = ~ ~ E(il(t) { ~ r=l ~ p=l pr 1LD(1 l(x) [aN] r.L(N(i)) · E ~ df)(il(x) 1 p aUr(•) [aN] + NENNEN ~ ~ E\9(t) { rL(N(i)) · E ~ df)(il(x). ~ ~ pr 1o<•>(x) • p a ·(•) r=l p=l (3.25) Ur Thus, this choice of basis functions <;l>;(x) uncouples the minimization equations, even in the case of nonstationary response. The matrix equations for the equivalent linear parameters a and {3 become two scalar algebraic equations for the equivalent linear parameters a; and ~; within element i. The terms in these equations depend only on statistics of response quantities for nodes attached to element i, so the equivalent linear parameters are easily computed during formation of the element stiffness and damping matrices. Thus, the calculation of the equivalent linear parameters can be accomplished in a numerically efficient way. 43 3.4 Development of the Covariance Equation The preceding section concluded development of the linearization procedure. Turn attention now toward the solution of the equivalent linear system. Application of the finite element method, described in an earlier section, to the continuous equivalent linear system (2.6) produces a set of discrete linear equations of the form Mii + C(t)u + K(t)u = f(t) (3.26) where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, and f is the total effective load vector. This discrete system has a dimension equal to the number of degrees of freedom in the finite element mesh for the problem under consideration. Since all external applied loads are assumed to be zero mean Gaussian, the load vector f is also zero mean Gaussian. A result from classical probability theory states that the output of a linear system subjected to zero mean Gaussian input is also zero mean Gaussian. Thus, the displacement vector u and velocity vector u of the equivalent linear system are a zero mean jointly Gaussian process. A Gaussian process is completely described by its mean vector and covariance matrix, so since the mean is zero only the joint displacement-velocity covariance matrix remains to be determined. This covariance matrix contains the variance of the displacement and velocity for all degrees of freedom, plus all of the cross-covariances as well. A procedure to express equation (3.26) in terms of the response covariance matrix is described below. First, the discrete 2nd order structural dynamics differential equation. above is converted to a 1•t order state space equation of higher dimension. Recall the 44 definition of the system state vector s as .~ (: l {3.27) Using this definition, the governing discrete equation (3.26) can be written s =As +F(t), where (3.28) A~ [_:_,K _:_,C is the system matrix, and F ~ M-~r(t) ( l l (3.29) (3.30) is the state space load vector. In the above expressions, I is the identity matrix of appropriate dimension, and it is assumed that the global mass matrix M is nonsingular. Multiplying equation (3.28) by sT, adding the result to its transpose and taking expected values gives the Liapunov equation for the covariance response as e.. (t) = A(t)E>,,(t) + E>,,(t)AT(t) + E [s(t)FT] + E [F(t)sT(t)], (3.31) where (3.32) is the state space covariance matrix. This is a symmetric set of 2n x 2n first order ordinary differential equations for the evolution of the nonstationary 'covariance matrix with time. 45 Thus far, it has only been assumed that the input excitation is zero mean and follows a jointly Gaussian probability distribution. For many engineering problems, it is sufficient to consider the case of white noise excitation, where the power spectrum is uniform across the entire frequency band. Under this restriction, the above equation reduces to a much simpler form. This simplification will be described in the next section. 3.5 Simplifications for Uniform White Noise Excitation Under uniform white noise excitation, the 3rd and 4th terms on the right-handside of equation {3.31) may be replaced with more explicit expressions. First, note that the effective load vector f{t) can now be written f{t) = '"t(t)n(t)g, (3.33) where '"f(t) is a deterministic modulating time function, n(t) is stationary Gaussian white noise with zero mean and constant power spectral density 8 0 , and g is a time-invariant vector which depends on the spatial distribution of the loads. The modulating time function '"Y(t) has been introduced to allow the intensity of the excitation to be scaled as a function of time. 46 Next, use this expression for the load vector in the 3rd term on the righthand-side of the covariance equation to get E [s(t)FT(t)] = "'(t)E [s(t)n(t)J (0 M- 1 g), (3.34) where ( ·) denotes a row vector. Now consider the term E [s(t)n(t)], and recall that s(t) satisfies equation (3.28). Thus, using the principal matrix solution <I>(t), s(t) can be written as the solution to equation (3.28) as s(t) = <I>(t)s(O) + <I>(t) fo' q,- 1 (e)F(e)de, (3.35) where ci> (t) satisfies ib(t) - A(t)<I>(t) (3.36) ci>(O) - I. (3.37) Using the expressions for the load vector (3.33) and (3.30) in this equation yields •l•l ~ O(t)•(O) + • (t) J.' •-•wol<l•l <) 1M~'il ) <e. Multiplying by n(t) and taking expectations gives E [s(t)n(t)J = <I>(t) rt "t(e)<I>- 1 (e) lo l 0 fl M-Ig E [n(e)n(t)] de. 13·38) (3.39) Next define the autocorrelation function for the stationary Gaussian random process n(t) by R,.,.(t- e) = E [n(e)n(t)], so (3.40) 47 For the case where n(t) is stationary random Gaussian white noise with zero mean, the autocorrelation function takes the form Rnn(t- e) = 21rSa8(t- e), (3.42) and thus equation (3.41) reduces to (3.43) Using this expression (3.43) in equation (3.34) gives the desired explicit form Now, note that the 4th term on the right-hand-side of the covariance equation is simply the transpose of the 3rd term: (3.45) Observing that the right-hand-side of equation (3.44) is symmetric, it follows immediately that Finally, using the simplified expressions (3.46) and (3.44) in the general covariance equation (3.31) gives the reduced form for white noise excitation as 48 3.6 Solution of the Covariance Equation for White Noise Excitation In the previous section, the Liapunov equation for the evolution of the system covariance matrix 0 under uniform white noise excitation was derived as (3.48) where 0 is the state-space covariance matrix, A is a 2n X 2n nonsymmetric system matrix containing the n X n mass, damping, and stiffness matrices, So is the constant power spectral density of the input, "f(t) is a deterministic modulating time function, M is the global mass matrix, and g is a constant vector which depends on the spatial distribution of the loads. Note that A in equation (3.48) depends on time since the equivalent linear system contains time-varying coefficients. The covariance equation (3.48) can now be integrated in time to find the nonstationary covariance matrix 0(t). Algorithms for performing this integration numerically on a digital computer are discussed below. Most time integration algorithms fall into one of two classes: implicit or explicit. An implicit time integration scheme requires solution of a matrix equation at every time step. In contrast, explicit algorithms do not involve a matrix 49 decomposition or simultaneous equation solution, but do require smaller time steps to maintain stability and accuracy. Algorithms for the implicit solution of the time-dependent Liapunov equation are complicated and numerically intensive, so to keep computational complexity to a minimum an explicit algorithm is selected here. Also, the size of the covariance matrix grows rapidly with prob!em size, so the time integration algorithm should be chosen such that computer high-speed memory requirements are minimized. Low order explicit algorithms are very memory efficient, thus allowing solution of much larger problems in a given memory size than higher order explicit schemes. The simple modified Euler algorithm is an explicit, second order accurate integration procedure which is well suited to solution of the Liapunov equation in a restricted memory space. Since it is a second order scheme, only a small number of large matrices must be stored in memory at any one time. In addition, since no matrix decompositions are involved, it is simply and efficiently implemented in a vector processing supercomputer environment. The modified Euler algorithm is conditionally stable, with stability assured for a sufficiently small time step D..t. The details of this procedure are described in reference [20]. The stationary response can also be obtained from equation (3.48). Under stationary conditions, the covariance matrix 0 is independent of time, therefore e = 0. The equation to be solved is now the time-independent Liapunov equation (3.49) Solution techniques for equations of this form are discussed in references [21,22]. The remainder of this discussion will focus on techniques for obtaining the nonstationary solution to equation (3.48). 50 3. 7 Solution Algorithm for N onstationary Response The complete procedure for obtaining the random response of a nonlinear continuous system may now be described. As a convenient notation for description of the algorithm, let ek denote 0(tk), i.e., the covariance matrix at the kth time step. Also, let an over bar denote "predicted" quantities in the predictor-corrector modified Euler integration scheme. It is assumed that the finite element method is used to discretize the equivalent linear system, and that the equivalent linear parameters are updated in every time step. The algorithm begins at step k, time tk, when ek is known, and computes the updated covariance ek+l as follows: 1. For each element • compute the equivalent linear parameters for the next time step based on the current response ek • form and assemble the element mass, damping, and stiffness matrices and load vector at time tk 2. Add additional physical damping contribution, if desired 3. Compute predicted covariance 0k+ 1 from (3.50) 4. Repeat steps 1 and 2 with ek and tk replaced with 0H 1 and tk+1 51 5. Compute corrected new covariance 0k+ 1 from (3.51) Note that the above procedure can be rearranged slightly for computer implementation to produce a more memory efficient scheme. Often the physical mechanisms producing damping in a given system are not well understood, and the addition of damping into the mathematical model is based on experience or experimental data rather than rigorous derivation. One form of damping widely used in structural calculations is Rayleigh damping, where the damping matrix is constructed as a linear combination of the stiffness and mass matrices: (3.52) As was indicated in the algorithm, an arbitrary damping contribution is easily introduced into the model with this formulation. This allows greater flexibility in creating realistic models than a procedure which allows only mathematically based damping contributions. This completes the specification of the numerical implementation and solution procedure. The next chapter describes the application of this technique to the solution of several problems of engineering interest. 52 Chapter 4 Application 4.1 Introduction This chapter describes the application of the newly developed solution method to two specific nonlinear problems. Before beginning discussion of the particular nonlinearities to be studied, it is appropriate to consider the objectives of such a numerical investigation. The numerical illustrations described in this chapter are chosen with three goals in mind. The first and obvious motive is to exercise the new technique and gain facility in its application. Second, it is important to validate the accuracy of the new method by comparing it against other available solution techniques. Finally, it is of engineering interest to assess the influence of typical nonlinearities on the nonstationary stochastic response of an example system. The above objectives all point toward selection of a simple nonlinear system for which at least some solutions exist, and for which the physical interpretation of results is clear. The nonlinear system chosen for study is a one-dimensional continuum described by the wave equation with a nonlinear constitutive law. The class of 53 systems modeled by the wave equation is one of significant engineering importance since it contains many common physical phenomena, such as axial vibration of a rod, propagation of plane waves in a continuum, and vibration of a shear beam. The one-dimensional wave equation has been widely studied, both as a continuous system and as a discretized system, using various techniques. This knowledge base provides benchmark solutions against which the current technique can be judged. In addition, the present understanding of phenomena which occur in this system under deterministic excitation permits insightful examination of the results obtained under random excitation. 4.2 Description of the Physical System In order to allow a physical interpretation of the example, it is helpful to choose a particular conceptualization of the one-dimensional wave equation for study. The problem chosen for study is a nonlinear shear beam under random base excitation. As shown in Figure 4.1, the beam lies along the x axis. The transverse displacement in an inertial coordinate frame is denoted by v(x, t), while the displacement relative to the base of any point x at timet is denoted by w(x,t). The beam is assumed to have constant density p and constant cross sectional area Ac. The base acceleration ij(t) is taken to be a random process. For this example, tf(t) is assumed to be Gaussian white noise with zero mean and constant power spectral density S0 • The equation of motion for this system is easily shown to be ar ax= pv, (4.1) ar pw. + pq..()t' ax= (4.2) or in terms of relative displacement 54 where r is the shear stress, and all other terms are as defined above. Using the chain rule and assuming that r depends on x only through the shear strain e:, the governing equation of motion can be written as 2 '()a w _ ..() T c 8x2 = pw + pq t ' (4.3) where a prime denotes differentiation of a function with respect to its argument. The nonlinearity arises from the material constitutive relationship, in which the shear stress is assumed to be nonlinearly related to the shear strain through a constitutive law of the form r=r(e:), (4.4) where the shear strain e: is defined as aw ax e:=-. (4.5) Two different forms of the nonlinear constitutive law are investigated in this chapter. The first example is a hardening stress-strain law, where the rnaterial tangent modulus increases with increasing strain. The second example is a softening stress-strain law, where the material tangent modulus decreases with increasing strain. Each of these cases is discussed in detail in the following sections. 4.3 Cubic Hardening Constitutive Nonlinearity 4.3.1 Introduction The first material nonlinearity considered is the cubic hardening constitutive law. In this material model, shear stress is related to shear strain by the equation (4.6) 55 where e represents shear strain as defined in equation (4.5), G is an elastic shear o is a nonlinearity parameter. A plot of the stress-strain curve for this constitutive relation is given in Figure 4.2 for various values of o. This modulus, and type of constitutive relation, with fi > 0, represents a so-called "hardening" or "stiffening" material because the incremental material tangent modulus increases with the level of strain. This type of material constitutive behavior can occur in porous media such as foams and in some rubbers. For small strain levels, this constitutive relation can be used with o < 0 to approximate the softening behavior found in many ductile metal systems during yielding. However, for larger strains, the material tangent modulus can become negative, thus leading to unstable behavior which may not accurately represent the physics of the problem being modeled. Therefore, except when otherwise explicitly stated, it is assumed that 0 > 0. Substitution of the constitutive law (4.6) and the definition of shear strain (4.5) into equation (4.3} yields the equation of motion in terms of displacement relative to the base as (4.7) This equation will be solved using the newly developed nonstationary continuous equivalent linearization technique. 4.3.2 Relationship to Discrete Systems As discussed in the opening paragraphs of this chapter, it is desirable to study a nonlinear system for which some benchmark solutions already exist. In this section, it is shown that the above cubic constitutive nonlinearity for a one di- 56 mensional continuum is directly analogous to the much studied cubic force nonlinearity of discrete dynamical systems. Consider now a discrete analog to the continuous shear beam described in the preceding section. Let this model contain lumped masses m; and nonlinear restoring forces/;, with the displacement of mass i relative to the base denoted by w;, for i = 1, ... , N. This system is shown in Figure 4.3. The nonlinear restoring force is taken to be cubic, so (4.8) where k; is the spring constant and lJ is the discrete force nonlinearity parameter. The relation between this discrete nonlinear system and the previously defined continuous nonlinear system will now be derived. Let the masses be spaced a distance D.x apart, so the total length of the beam L is given by L = N D.x, where N is the number of discrete masses. Thus, taking the origin of the x-a.xis at the base of the beam, w; = w(iD.x). ( 4.9) Using this expression in equation (4.8) gives f .= k·D. ' ' 3 [w(x) - w(x- D.x) liD. 2 (w(x) - w(x- D.x) ) x D.x + x D.x ] . (4.10) Now, assume that w varies linearly from x;_ 1 to x;. This is consistent with the use of finite elements with linear shape functions for solution of the equivalent linear continuous system. Under this assumption, aw w(x)- w(x- D.x) ax = D.x (4.11) and using the definition of shear strain (4.5), the restoring force f; can now be written (4.12) 57 In this one-dimensional case, the restoring shear force /; can be converted to the restoring shear stress r simply by dividing by the cross-sectional area A., so the effective stress-strain relation for the discrete system is obtained from the previous equation as 1' = 23] • k;t:.x - [e: + ot:.x e: Ac (4.13) Making the equivalences (4.14) and G = kt:.x, A, (4.15) the effective constitutive law for the discrete system now agrees with that of the continuous system (4.16) Thus, many of the results in the literature for the discrete cubic nonlinear system may now be related to the finite element solutions of the continuous nonlinear system with cubic constitutive nonlinearity through use of relations (4.14) and (4.15). 4.3.3 Minimization Equations In order to apply the nonstationary continuous equivalent linearization method to the example shear beam described earlier in this chapter, the general minimization equations developed in chapter 3 must be specialized for this system. First, the operators used in the development of the general theory are identified for the current problem. Next, derivatives of the nonlinear stress operator appearing in the minimization equations are evaluated using these specific definitions. Finally, 58 the remaining terms in the minimization equations are evaluated for the current system. Recall that the general equation of motion for the shear beam was given in equation (4.2) as ar pw- + pq""(t.) ax= (4.17) In this one-dimensional system, the general divergence operator "il· reduces to simply %z. Using the definition of shear strain (4.5) in the constitutive law (4.6) gives the stress in terms of the relative displacement as 3 aw aw r(w)=G-+8G( ax ) ax (4.18) The linear stress operator is now easily identified as (4.19) and therefore the nonlinear stress operator is (4.20) Thus, the linear and nonlinear stress operators are conveniently obtained from the constitutive law and kinematic equation. Next, derivatives of TN appearing in the minimization equations (3.24) and (3.25) must be evaluated using the above operator definitions. Since the finite element method will be used to solve the equivalent linear time-varying system, the operator definitions (4.19) and (4.20) above are first expanded using the finite element approximation. Recall that, in the finite element method, the dependent variable w is approximated within element i by NEN wUl(x,t) = ,L N?)(x) u;n(t), j=l (4.21) 59 where all terms are as defined in section 3.2. Substituting this discretization for w into the definition (4.20) gives the nonlinear stress operator within element i as ?; NEN&N(i) ) ; ulil(t) rN(w(i))=W ( 8 3 (4.22) Differentiating this expression with respect to the nodal value u~i) gives 8rN (NEN 8N(i) . )2 8N(i) -w = 38G L - • ul'l(t) - • . 8 8 8Ur •=1 X (4.23) X Also, note that since the cubic hardening nonlinearity depends only on displacements, 8rN aUr.(i) = 0. (4.24) The squared term in equation (4.23) can be expanded to give 8rN [NEN NEN 8N(i) 8N(i) . . ] 8N(i) - . = 38G ~-·-u<•l(t)u<•l(t) - • 8u~•) q=1 •=1 8x ax q ' ax ' L L (4.25) hence, taking expected values yields E &N(i) NEN NEN &N(i) BN(i) a N] r(i) = 38G-•- L L ,_E [u~il(t)uiil(t)]. [Bur 8x q= 1 •= 1 8x 8x _ q_ _ _ (4.26) Thus, these derivatives of the nonlinear stress operator are now expressed in terms of the element response covariance matrices, as desired. The only remaining term to evaluate in the minimization equations is the linear stress operator acting on a finite element shape function, such as r/(NJil). Recall that, by definition, (4.27) Using the choice of functions ¢; discussed in chapter 3 and the definition of the linear stress operator (4.19), it follows that (i) rL(N(i)) = 8NP • p ax . (4.28) 60 The complete set of general minimization equations may now be evaluated for the one-dimensional shear beam with cubic constitutive nonlinearity. Recall that the minimization equations for element i are (4.29) and (4.30) Using the expressions for the linear stress operator (4.28) and the nonlinear stress operator (4.26) in the general minimization equation terms (3.21) through (3.25) yields A(i) = . L L E [u<•l(t)u<'l(t)] r. NEN NEN l=l p=l P I NEN NEN . j /)!•l(x) . B<•l = I: I: E [uf'l(t)u<'l(t)] l=l c<•l = p=l p NEN NEN . 1 . f)(•l(x) { L L E [uf'l(t)u~'l(t)] lr<l l=l aN<•l aN<•l p=l /) • (x) dD('l(x) (4.31) ---p-dD(il(x) ax ax (4.32) _ l_ _ _P_ ax ax aN(i) aN(i) - 1 aN(i) aN(i) _ l_ _ _P_dD(il(x) ax ax (4.33) NEN NEN NEN NEN Q(i) = 38G L L L L E [u~•lu~•l] E [u~•lu~•l] r=l X 1 p=l q=l s= 1 aN(i) aN<'l aN<•l aN<'l _ _ P_ _ _ r ___ , ___ q_ dD (") ' (x) /)(il(x) ax ax ax ax (4.34) NEN NEN NEN NEN H(i) = 38G L L L L E [u~•lu~•l] E [u~•lu~•l] r=l X p=l q=l s:;:::l aN<•l aN(i) aN<•l aN(') (") (x). __ P_ _ _ r ___ , ___ q_ dD' /)(i)(x) ax ax ax ax 1 (4.35) This completes the specialization of the minimization equations for the shear beam with cubic constitutive nonlinearity. 61 4.3.4 Finite Element Matrices and Vectors Using the minimization equations derived in the previous section, the original nonlinear system is replaced by an equivalent linear continuous system. The replacement system is then solved using the finite element method. This section describes the specific form of the equivalent linear system and the derivation of the appropriate finite element matrices and vectors. The general form of the equivalent linear system is given in equation (2.6). Applying the operator definitions (4.19) and (4.20) from the previous section, this general form reduces to (3 (4.36) w(x, t) = o on aD.,, (4.37) a ( Gaw) a ( a aw) a ( aw) _ ax ax +ax ax +ax ax = pw_+ pq(t). The appropriate boundary conditions are and (4.38) where aD., denotes that portion of the boundary where w is specified, a[)T denotes portions of the boundary where r is specified, and TT is an applied shear stress. The finite element equations will be derived using Galerkin's method. Multiplying the partial differential equation (4.36) by a weighting function W and integrating over the domain gives fo(x) W [:x (G~:) +:x (a~:)+ :X ((3~~) - pw - pq(t)] d[) (x) = o. (4.39) Similarly, the Neumann boundary condition (4.38) becomes fooT W [(a+ a)~:+ (3~~- rT] d(aDT) = o, (4.40) 62 where iJD = iJD., + iJDT· (4.41) Integrating equation (4.39) by parts gives the weak form of the governing equation as aw aw aw aw OX OX aw aw _ _ ] + W pw + W pq(t) dD(x) lr [- G - +-a,-+ -0 {30 P(x) OX uX X X aw + Waaw + W/3aw] d(oD) = 0. - /, [WG0X oL' OX 0X (4.42) The domain of the second integral can be divided into oDT + oD.,. Using the boundary condition (4.40) then gives aw aw aw aw aw aw _ _ ] - G- + - a - + - {3- + w pw + w pq(t) dD(x) 8X 8X 8X 8X 8X 8X lr [ L'(x) - ( lavr r WrTd(oDT) [ - lavw W (G +a) aw ax + f3 aw] ox d(oD.,) = 0 · (4.43) Next, divide the domain D(x) into elements p(kl(x), so D(x) = I::D(kl(x). (4.44) k Using this domain decomposition in equation (4.43) gives (4.45) where it is noted that G, a, and f3 will in general vary from element to element, and either aD~kl or oD~kl or both could be zero for a particular element (such as an interior element). 63 Let N(k) represent a vector of finite element shape functions for the nodes of element k, and let u(k) be the vector of nodal values of w(k) (x, t) for nodes of element k. This gives the finite element approximation, in vector form, as (4.46) To obtain an equation for a typical global node p not on aD,., take the weight function as w(•l = pw(•l = L NJ"l, (4.47) kp i.e., the weight function for node p is taken as the sum, over all elements connected to node p, of the shape function in each element corresponding to node p. This choice yields the Bubnov-Galerkin formulation, and leads to symmetric matrices. Now, using this choice of weighting function W in equation (4.45) gives the following equation for global node p: "'{!r ~ .. [aN(k) aN(k) (aN(k))T _P_(G +a) (aN(k))T - - u (k) + _P_f3 - - U. (k) o<•>(x) ax ax ax ax + pNJ"l(N("))Tii(k) + pNJ"lq(t)] d[)("l(x) - f <•> NJ"lrTd(aoJ."l)} = o, laoT (4.48) where the summation is again over all elements k connected to global node p. Equations for all nodes not on a[),., i.e., not on the boundary where w is specified, can be assembled into Mii + Cli+ Ku = -fq +fr. (4.49) where the corresponding element quantities are given by (4.50) 64 c(k) K(k) h - h (G+a)-- - - aN(k) ( aN(k)) T (3-- - dD(k)(x) o<•l(x) ax ax aN(k) ( aN(k)) T ax o<•l(x) f(k) T f(k) q - ax dD(k)(x) Iaao<•l N(k)rT d(aD(k)) T h (4.51) (4.52) (4.53) T pq(t)N(k) dD(k) (x) Df'l(x) (4.54) and homogeneous geometric boundary conditions have been assumed. Note that for the one-dimensional shear beam under consideration, there are no applied shear tractions, so f, = 0. 4.3.5 (4.55) Performance Evaluation In order to evaluate the performance of the newly developed nonstationary equiv- alent linearization technique on problems involving a hardening constitutive nonlinearity, a series of three examples are solved for which solutions are obtainable by other techniques. The results from the continuous equivalent linearization method are found to agree well with the benchmark solutions. Details of these validation calculations are presented in Appendix A. 4.3.6 Description of the Example Problem With the performance of the new method now demonstrated, it is appropriate to turn attention to the solution of some additional problems. The efforts of this section are directed toward gaining physical insight into the nonstationary response of a continuous shear beam with a cubic hardening constitutive nonlinearity. To this end, several sets of calculations are performed to illustrate the physics of the 65 system and the characteristics of the solution technique. First, since it is well known that an undamped linear second-order system subjected to white noise excitation exhibits unbounded response variance, it is of interest to examine the influence of the form of damping added on the system response. For conciseness, only two classes of Rayleigh damping are explored here: mass proportional damping and stiffness proportional damping. To isolate these effects, this set of calculations is performed on a linear shear beam. It is emphasized that the new method is not restricted to Rayleigh damping in any way, and that many other techniques could equally well be used to incorporate damping into the model. Next, in order to assess the effect of mesh refinement, a set of calculations is performed for the same problem with varying mesh sizes. Again, to prevent the introduction of spurious effects, this analysis is done for a linear system. Finally, a set of calculations is done to investigate the influence of the cubic constitutive nonlinearity on the nonstationary stochastic response. The nonlinear results are compared with the corresponding linear results, and the observed effects are discussed. In order to facilitate comparison of results from the different sets of cal- culations, the same basic problem is used throughout this section. The system is a continuous shear beam, with a cubic hardening constitutive nonlinearity, subjected to Gaussian white noise base excitation. A diagram of the continuous system is shown in Figure 4.1, and the corresponding finite element model is shown in Figure 4.10. The beam is taken to have uniform properties throughout its length, with GL~· = 3.0 and pA,L' = 1.0. The power spectral density of the Gaussian white noise base excitation is 8 0 = 2.0. The type and amount of damping added as well as the value of the nonlinearity parameter vary for each set of 66 calculations. 4.3.7 Effect of Type of Rayleigh Damping It is a well-known result that an undamped linear system subjected to Gaussian white noise excitation exhibits an unbounded response variance. This is physically reasonable since the power spectral density of white noise is constant over the entire frequency range, thus yielding infinite energy content. In spite of this property, white noise remains an extremely useful analytical construct because most dynamic systems respond predominately to energy in narrow frequency bands, and the white noise model can often yield a good approximation to the energy input in those frequency regimes. To preclude problems with infinite response magnitudes, a small amount of damping is added to the mathematical system model for analysis under white noise excitation. All real engineering structures possess some amount of damping, so this addition is not artificial. The objective of this set of analyses is to explore the influence of the form of the damping added on the calculated responses. The mechanism of damping in real engineering structures is complex and not well understood. Therefore, often it is necessary to base the addition of damping on a limited amount of experimental data. This motivates the choice of a simple method for incorporating damping into the analytical model. One very popular procedure for the addition of damping to structural response calculations is Rayleigh damping, wherein the damping matrix C is taken as a linear combination of the stiffness and mass matrices: (4.56) The significance of choosing mass proportional or stiffness proportional damping 67 is best seen by examining the behavior of the modal damping ratios for each case. Consider the discrete equations resulting from application of the finite element method to a linear continuous system and the addition of Rayleigh damping: (4.57) Transforming the above equation into modal coordinates yields a set of modal equations of the form + Ct + CzW;2) Z; + W;2Z; = g; (t ) , 0 Z; 00 ( (4.58) where z. is the ,~h modal coordinate, w; is the natural frequency of the i 1h mode, and g; (t) is the i 1h modal forcing function. Thus, the fraction of critical damping in mode i, for Rayleigh damping, is given as (4.59) As is evident from the above equation, mass proportional damping leads to a fraction of critical damping i; which decreases with mode order (w; increasing), while stiffness proportional damping yields a fraction of critical damping which increases with mode order. Therefore, mass proportional damping leaves the higher modes more lightly damped, while stiffness proportional damping more strongly damps the higher modes. Frequently, mass proportional damping is referred to as external or air damping, and stiffness proportional damping is also known as internal damping. 68 In the following calculations, two limiting cases of Rayleigh damping are examined. In the first case, the damping matrix is taken as proportional to the mass matrix, so c2 = 0.0. The constant of proportionality c1 is chosen such that the first mode fraction of critical damping .;-1 = 20%. This yields c1 = 0.0748. In the second case, the damping matrix is taken as proportional to the stiffness matrix, so c 1 = 0.0. Again, it is desired to have the first mode fraction of critical damping as 20%, and this leads to c2 = 2.14. The resulting fractions of critical damping in other modes are shown in Table 4.1, where it is noted that the natural frequencies used are those of the discretized system produced by the finite element method. The results of the calculations are shown in the top, middle, and bottom plots of Figure 4.11 as comparisons of the normalized displacement variance profile at early time, mid-time, and late-time, respectively. For these calculations, early time corresponds to about one one-hundredth of a signal propagation time across the structure, while mid-time corresponds to one-half of a signal transit time. Late time corresponds to essentially stationary conditions for the two damped cases, at time equal to approximately six signal transit times. In each plot, the solid line shows the undamped results, while the dashed line shows the mass proportionally damped results, and the dot-dashed line gives the results for stiffness proportional damping. In order to allow comparison of the profile shapes of the damped and undamped responses, each curve has been normalized by its peak value. Examination of the figure clearly reveals the influence of the form of damping on the evolution of the response profiles. The response profiles of the undamped and mass proportional cases are very similar in shape, while the signal front spreads much more rapidly in the stiffness proportional case. As Table 4.1 69 mode mass prop . stif. prop. k . 0748 2wk t(2.14)wk 1 20% 20% 2 6.7% 59.7% 5 2.3% 173% 10 1.3% 317% Table 4.1: Modal Damping Ratios for Rayleigh Damping Study shows, modes higher than about 3 are heavily damped for the stiffness proportiona! case, and therefore die out quickly. A superposition of the first two or three modes cannot describe the sharp break in the response profile, so the break appears to spread. In contrast, the higher modes remain lightly damped in the mass proportional case, so the signal profiles are very similar to the undamped case. 4.3.8 Effect of Discretization Whenever a continuous problem is discretized using the finite element method, the question of mesh refinement arises. In order to assure that the mesh used for these example problems adequately represents the physical behavior of the system, a parametric study of mesh size is conducted. The base excited shear beam problem shown in Figure 4.1 is analyzed using uniform meshes of 7, 14, and 28 elements. To allow isolation of the effect of mesh size upon the results, the problem is taken as linear and no damping is added. Since this problem is undamped, linear growth in the response variance is expected, but the spatial 70 response profiles and variance growth rates can still be directly compared between the three mesh sizes. The results of these calculations are presented in Figure 4.12. To aid comparison of the response profiles, all responses at a given time are normalized by the peak value of the large mesh result. The large mesh value was chosen since it is considered the most accurate. In each plot, the solid line represents the large mesh (28 elements) result, the dashed line represents the medium-sized mesh (14 elements) result, and the dot-dashed line represents the small mesh (7 elements) result. The early time response profiles at the top of the figure show some influence of the varying mesh size. The signal has crossed one element, and the variation of element size is thus reflected in the shape of the response profile. This effect is much more pronounced between the coarse and moderate meshes than between the moderate and fine meshes. The fine mesh plot shows that the signal is beginning to reach the second element, and thus both the moderate and fine meshes have captured this portion of the profile. The mid-time results shown in the middle graph of Figure 4.12 also illustrate close agreement between the moderate and fine mesh results. Both of these meshes have successfully captured the shape of the profile near x = 0, while the coarse mesh does not have sufficient resolution to approximate the nearly zero tangent at x = 0. The lower graph in Figure 4.12 shows that the late-time response profiles agree quite well among all of the meshes. By this time, the sharp gradients in the profile have become smooth, and all three meshes yield satisfactory results. Based on these results, the 14 element mesh is judged adequate for the current problem. 71 4.3.9 Effect of the Cubic Hardening Nonlinearity The above results serve to give confidence to the modeling and analysis procedures to be used in the next set of calculations. In this study, the spatial response profiles and variance time histories of the cu hie hardening nonlinear shear beam are directly compared against the corresponding linear results. These analyses clearly illustrate the physical influence of the nonlinearity on the response characteristics. The results of these calculations are shown in Figures 4.13 and 4.14. In order to retain bounded solutions without distorting the response profile, mass proportional damping has been added such that the first mode of the linear part of the system contains 20% critical damping. To clearly illustrate the effect of the nonlinearity, the response profiles in Figure 4.13 have been normalized by the peak value of the linear profile. In these plots, the solid line shows the linear profile, while the dashed line, the dot-dashed line, and the dot-dot-dashed line represent nonlinear results with nonlinearity parameter values of 6 = 0.25, 0.50, and 1.00, respectively. The early time results of the top graph of Figure 4.13 show very little difference between the linear and nonlinear systems. This is reasonable, since the responses are still small, and for very small strains the contribution of the cubic term in the nonlinear constitutive equation is negligible. It should be recognized that the constant value of the response profile in the unaffected portions of the beam just reflects the variance of the base motion, since the upper sections are still at rest in absolute coordinates. The middle graph of Figure 4.13 shows that a marked difference in the response profiles has developed by mid-time (recall that this corresponds to about 72 one-half of the signal propagation time through the linear beam). The linear result (solid line) shows that the profile of mean squared relative displacement gradually bends over to the constant value in the unaffected section of the beam, while the nonlinear results (broken lines) show the development of a sharp "corner" or break at that location. This effect arises from the cubic hardening nonlinearity in that the incremental stiffness, and therefore the wavespeed, increases with strain level. Thus, strain increments at larger strains propagate faster than strain increments at lower strains, creating a sharp front in the response profile during the nonstationary portion of the response. It is also evident from the figure that this sharp front propagates faster than does the corresponding linear signal. As the bottom graph in Figure 4.13 shows, this effect largely disappears by late time. After the signal has crossed the structure several times, the response profile settles down a smooth shape very similar to that of the linear system. A marked difference in the variance growth rates is seen, however. The nonlinear system has scarcely reached a displacement variance 10-20% as large as that of the linear system. This is to be expected, however, since the nonlinear system becomes stiffer as the response level increases, thus attenuating the rate of growth of the variance. Thus, even with damping present, it is clear that the stationary mean squared displacement of the nonlinear system is substantially less than that of the linear system. For comparison, the time history of the variance at the free end of the beam is plotted in Figure 4.14 for the linear case and three differing degrees of nonlinearity. The peak ratio of nonlinear "restoring stress" to linear "restoring stress" under stationary conditions is .45, .65, and .92 for the 6 = 0.25, 6 = 0.50, and 6 = 1.0 cases, respectively. This figure graphically illustrates the attenuation 73 of displacement variance with increasing nonlinearity. Thus, the effect of even these moderate nonlinearities is clearly evident in the results. In summary, this set of analyses has shown that the nonlinearity has a strong effect on the shape of the response profile during the nonstationary portion of the response, and on the magnitude of the displacement variance during all portions of the response. Some engineering interpretations of these results are discussed in the next section. 4.3.10 Engineering Implications It is interesting to briefly consider some implications of the above results for engineering practice. First, the faster propagation of the signal in a hardening nonlinear medium implies that a given point on the structure would experience motion earlier for a nonlinear medium than for a corresponding linear medium. This earlier arrival could then result in an increased duration of shaking at that location in the structure. Next, the development of the sharp break in the nonlinear response profile implies a concentration of deformation and locally larger strains in the nonlinear structure than in a similar linear structure. Thus, a structure made of a hardening material would have to withstand more severe local deformation without failure than would a structure made of a material that is linear up to failure. Finally, the hardening nonlinearity served to strongly attenuate the magnitude of the displacement variance under stationary conditions. Thus, such a material would be a good choice for a shear member with maximum displacement design constraints, if sufficient deformation tolerance is included. As the above discussion illustrates, the hardening constitutive nonlinearity leads to both good and bad characteristics in the structural response, and these must be carefully assessed in the design process. 74 4.4 Arctangent Softening Constitutive Nonlinearity 4.4.1 Introduction This section investigates a nonlinearity that is, in some sense, opposite to the one previously studied. The arctangent nonlinear constitutive law represents materials which soften, rather than harden, with increasing levels of strain. Many real materials fall into this category, such as most metals, composites, and many soils. Thus, the study of a softening nonlinearity has direct application to many areas of engineering practice. The nonlinear softening constitutive law is defined by the equation 2 r(c) = --r tan~ 1 (?r Goc) -- ?rY 2-ry (4.60) with ow ax c=-, (4.61) where • c is the shear strain, • ru is the "yield stress" of the material, and • G 0 is the initial elastic shear modulus. Recall that tan- I ( 0) ""' 0 for small 0, so for small strains the above stress-strain relation reduces to the linear elastic relation r(c) = G 0 c. (4.62) 75 Plots showing the stress-strain behavior produced by equation (4.60) for various values of 01 are given in Figure 4.15. '• Using the constitutive law (4.60) in the general equation of motion (4.3) yields the governing PDE for this nonlinear shear beam as Go G 1 + ( 2!:2 "-ll.!. ., ) 2 2 8 w aX = pw.. + pq"(t ) . 2 (4.63) The solution to this equation will be obtained using the nonstationary continuous equivalent linearization method developed in the previous chapters. 4.4.2 Relationship to Discrete Systems The continuous nonlinear shear beam with an arctangent softening nonlinearity has an analogous discrete representation. In order to allow results from the nonstationary continuous equivalent linearization procedure to be compared with existing results for discrete spring-mass-damper systems, this analogy is developed below. Consider once again the discrete analog to the one-dimensional shear beam. Recall that this model contains lumped masses m; and nonlinear restoring forces /;, and the displacement of mass i relative to the base is denoted by w;, for i = 1, ... , N. This system is shown in Figure 4.3. For convenience, define the relative displacement y; as the displacement of mass i relative to mass (i- 1), so (4.64) With the discrete system now defined, the form of the nonlinear restoring force is derived from the stress-strain relation given in equation (4.60). First, replace the derivative in the strain-displacement relation (4.61) by its backward 76 difference approximation w(i~x)- w((i- l)~x) e~ ~X (4o65) . Let w; = w(i~x), so (4o66) or in terms of the relative displacement y;, Yi (4o67) E:~-0 ~X Next, substitute this approximation for the strain c: into the constitutive law (4o60) to get r(c:) = ~r tan- 1 (~Go 7r y .2_Y•) 0 2 Ty ~X (4o68) Consider now that the spacing between two lumped masses in the discrete model is ~x, and therefore the restoring force is constant over that interval. For the spring between masses i and (i -1), this restoring force is obtained by multiplying equation (4o68) by the cross-sectional area over which the shear stress acts, A. : 2 Go ) J; = -A.ry tan _ 1 (11" -~y; 7r 2 Ty'-l.X o (4o69) Define now a "yield force" as the product of the yield stress and crosssectional area (4o70) then equation (4o69) becomes (4o71) Defining the initial stiffness k 0 by k G0 A. o=~, (4o72) 77 the discrete restoring force expression takes its final form /; = ;fytan- 1 (~~>•). (4.73) Thus, equation (4.73) represents a discrete nonlinear restoring force analogous to the nonlinear continuum constitutive law (4.60). The only approximation made in this derivation was in the replacement of the derivative in the straindisplacement relation with its first-order backward difference approximation. If finite elements with linear displacement interpolation functions are used in the continuous equivalent linearization solution, then strain is constant over an element and this relation becomes exact. 4.4.3 Minimization Equations In order to apply the nonstationary continuous equivalent linearization method to this problem, it is first necessary to evaluate the terms in the minimization equations for the particular constitutive law under consideration. To this end, first the operators TN and TL are defined for the arctangent softening nonlinearity. Next, expressions for the minimization equation right-hand-side terms G; and H; are developed. Finally, the remaining terms in the minimization equations are evaluated. From the preceding comments on the small-strain behavior of this constitutive law, it is clear that the linear stress operator TL is given by (4.74) Now, recall that one of the initial assumptions in the formulation of the continuous equivalent linearization technique was that the system was linearizable, and this was operationally defined as requiring the effect of TN to become negligible 78 compared to rL as the response became small. Since the arctangent nonlinearity reduces to Gas for small responses, it is necessary to subtract this component out in defining TN; thus (7r 2 r tan- 1 - Gas) rN = - - -Gas, 7r y 2 T (4.75) " and the total constitutive law is the sum of the linear and nonlinear components: (4. 76) For the nonlinearity under consideration, it is more convenient to evaluate the nonlinear terms in the minimization equations directly from their definitions in section 2.5 than to use the technique developed in section 2.6. Using equations (2.58) and (2.64), along with the equivalent linear parameter basis functions given by equation (3.14), the right-hand-side terms of the minimization equations for finite element i can be written as (4. 77) and H; = { . E [rN(c,w,w) · r/(w)] dl)(il(x). Jv<•J (x) (4.78) Substituting in the above definition of the linear stress operator (4.74), these equations can be written (4.79) and H; = { . Jl)(•l(x) E [rN (s) t] dl)(il(x). (4.80) As before, it is now necessary to express G; and H; in terms of the element response statistics. This manipulation is described in the following paragraphs. 79 Temporarily consider only the integrands of equations (4.79) and (4.80). If these integrands can be written in terms of the element response covariance matrix, then G, and H; can be easily computed using numerical integration. Substitute the definition of the nonlinear stress operator (4. 75) into the integrand of each of these equations (4.79) and (4.80) and rearrange to get E[rN(e)e] =E[;ryetan- 1 (~~:e)] -GoE[e 2 ] ( 4.81) and ( 4.82) Clearly, the second term in each of (4.81) and (4.82) is easily written in terms of the element response covariance matrix, so concentrate attention on the first term in each of these equations. Recall that e and e are spatial derivatives of the displacement and velocity response processes, and since the replacement system is linear and subjected to Gaussian input, then these approximate response quantities are jointly Gaussian distributed also. Using this fact that e and e are jointly Gaussian distributed, the first term in equation (4.81) can be written as e)] Joo 7r 2 [etan- 1 (7r G 0-r.E - 2 Ty 7r 1 2 -r.etan- = -oo ('l2r G e) p(e)de, 0 --- (4.83} Ty where p( e) represents the Gaussian probability density function of e, ( vz:;ra, -= p e:) = ,o 1 e 2al!' (4.84) and a. is the standard deviation of e:, (4.85) As discussed in reference [24], integrating by parts and using the result -a2z2 ~ " 100 1e +X dx = -ea erfc(a), 2 0 ---:2 2 (4.86) 80 equation {4.83) becomes {4.87) where {4.88) and erfc{a) denotes the complementary error function of a. Note that the right-hand-side of equation (4.87) contains only deterministic functions and u•. It will now be shown that Ue is easily expressed in terms of finite element shape functions and the element response covariance matrix. Recall the definition of shear strain c: as {4.89) and the finite element approximation for w within element i, NEN w(il(x,t)., L N~i)(x)u~l(t), (4.90) k=l where NEN denotes the number of nodes per element. Substituting equation {4.90) into (4.89), squaring the result, and taking expected values gives (4.91) A similar argument shows {4.92) and NEN NEN aN(i) aN(i) E[c:i] = L L k=l 1=1 _k- 1 -E [u1;l(t)ul;'(t)]. ax ax (4.93) 81 The correlation coefficient p between the strain c and the strain rate f: can now be defined as E [cf:] a$ai p=--. (4.94) Thus, the statistics of the of the strain and strain rate are now expressed in terms of the required element response covariance matrix. The final step is to express the integrand of the H; expression in terms of the response statistics of c and f:. Then, equations (4.91) through (4.93) above can be used to get H; in terms of the desired element response covariance matrix. Proceeding as for G;, the first term in equation (4.82) is written in terms of the joint probability density function p(c:, f:) as (4.95) where (4.96) It is shown in Appendix B that pursuing this approach leads to the result (4.97) where a is as defined in equation (4.88). Thus, the integrands for both G; and H; are now expressed in terms of the element response covariance matrix, as required. The final minimization equations for the softening nonlinearity can now be summarized. Note that the left-hand-side terms of the minimization equations 82 depend only on the form of the linear stress operator, and that the form of the linear stress operator .,-L is the same for the arctangent softening nonlinearity as it was for the cubic hardening nonlinearity. Recall that the minimization equations for element i are given by (4.98) and (4.99) and the left-hand-side terms are given by (4.100) B(i) = . fz . L L E [u(•l(t)u(il(t)] NEN NEN l;l p;l l BN(i) aN(i) _ l_ _ _P_ dD(il(x) ax ax l)(oi(x) p . fz aN(i) aN(i) L L E [ul'l(t)u(il(t)] . - 1---P-dD(il(x). l;l p;l p v<•l(x) ax ax (4.101) NEN NEN C(i) = (4.102) Using equations (4.81) and (4.87) in (4.79) yields f'ia.e" erfc(a)- G0 a~] dD(il(x), lv<•l(x) Y; G; = { . [ 2 (4.103) and similarly, using equations (4.82) and (4.97) in (4.80) gives H; = { lv<•l(x) [ f'i.,..pa,e•' erfc(a)- G pa,a,] dD(il(x). y; 0 (4.104) Equations (4.98) through (4.104), along with equations (4.91) through (4.94) for the statistics a., a,, and p in element i, form the complete set of element level minimization equations for the arctangent softening constitutive nonlinearity. 83 4.4.4 Finite Element Matrices and Vectors The problem under consideration is still a one dimensional system satisfying the wave equation, so the form of the equivalent linear replacement system is unchanged from that described in section 4.3.4. As a result, the finite element matrices and vectors are identical to those described in that section, except that the equivalent linear parameters a and f3 are evaluated using the new minimization equations derived in the previous section. To assure that the newly developed minimization equations do yield an accurate solution for problems involving an arctangent softening nonlinearity, a simple test problem is solved for which a solution already exists. Agreement between solutions obtained using two different methods is quite good. Details of this problem are presented in Appendix A. 4.4.5 Description of the Example Problem This section investigates the effect of the arctangent softening constitutive nonlinearity on the stochastic response of a system satisfying the one-dimensional wave equation. Several calculations are conducted with varying degrees of material nonlinearity. The observed response characteristics are discussed, and comparisons are made with the response characteristics of the cubic hardening system discussed earlier in this chapter. Finally, some engineering implications of these results are discussed. As in previous sections, the example system will be physically conceptualized as a nonlinear shear beam subjected to Gaussian white noise base excitation. This structure is shown in Figure 4.1. The beam is taken to have uniform properties throughout its length, with G1~·· = 3.0 and pA,L' = 0.1, where L' denotes 84 the length of a finite element. The finite element model used in these example calculations is shown in Figure 4.10. It contains 14 elements with linear displacement interpolation functions. Mass proportional damping is added to the assembled discrete equations to yield a finite stationary displacement variance without excessively distorting the response profile. The applied base acceleration is given unit power spectral density, and is suddenly applied to an initially quiescent beam at t = 0. The degree of nonlinearity, as measured by the ratio of the initial shear modulus to the yield stress Qo., is varied for each analysis. The way in which the '• solution becomes established is then studied for each case. In addition, the effect of the nonlinearity on the magnitude and shape of the response profile is also examined. The results of these calculations are discussed in the following section. 4.4.6 Results for the Arctangent Softening Material Nonlinearity To study the effect of the nonlinearity on the response, three calculations with varying degrees of nonlinearity are examined. In the first problem, denoted "NL 1," the parameter value Qo. = 1.0, while in the second problem, denoted '• "NL 2," ~ = 2.0, and in the third problem, denoted "NL 3," a value of Qo. = 3.0 'Y 'fy is used. For comparison purposes, a fourth problem with a linear stress-strain law is also analyzed. Clearly, as the ratio Qo. ,, increases, the yield level decreases (since Go is constant) and the problem becomes more nonlinear. The results of the four calculations are presented in Figures 4.16 through 4.18. To graphically illustrate the effect of the nonlinearity, all profiles are normalized by the peak value of the corresponding linear profile. In each figure, 85 the solid line represents the linear result, while the dashed line, the dot-dashed line, and the dot-dot-dashed line represent the NL 1, NL 2, and NL 3 results, respectively. The graphs in Figures 4.16 and 4.17 show the mean squared relative displacement response profiles for the four analysis cases at four fixed values of time. On these graphs, time is measured in multiples of the signal propagation time across a linear beam with the same properties. This transit time is given by k. where c is the shear wave propagation speed in the linear material, given by •• ,...,- v c __ ./Gpo' (4.105) and L is the length of the beam. The top graph in Figure 4.16 shows the response profiles at a time corresponding to one-half of the signal transit time across a linear beam,i.e., t = 0.5~. One effect of the nonlinearity is already clear: the signals are propagating more slowly for the more nonlinear cases. This effect is easily explained, however, since the material tangent modulus decreases with increasing response levels, thus decreasing the effective wavespeed as the response builds up. At this early time, there is no significant difference in the magnitude of the response between the linear and nonlinear cases. The bottom graph in Figure 4.16 shows the response profiles at a time equal to six linear signal transit times across the structure. The propagating front has smoothed out, and the response profiles are nearing their stationary shapes. The nonlinear profiles show a much sharper gradient near the base of the beam than does the linear profile, and this effect increases with increasing degree of nonlinearity. Thus, most of the deformation is concentrated near the base, with the remainder of the beam participating little in the response. At time equal to 12 linear signal transit times across the structure, the 86 top graph in Figure 4.17 shows that the concentration of the deformation at the base of the structure has become even more pronounced. This effect is slightly visible in the NL 1 case, becomes prominent in the NL 2 case, and is extreme in the NL 3 case. For the two most nonlinear cases, only the first 3 or 4 elements participate significantly in the response. At time equal to 40 linear signal transit times across the structure, the bottom graph in Figure 4.17 shows that the profile shapes are very similar to those in the previous graph. The response magnitudes for the two most nonlinear cases have increased slightly, but the shapes remain largely unchanged. Figure 4.18 shows the time history of the free end variance for each of the four analysis cases. As expected, the most nonlinear case attains the highest response level, with the other cases correspondingly lower. The linear case leads to a lower response level than does any of the nonlinear cases. This is expected, since the softening behavior of the constitutive nonlinearity tends to increase response levels over those of the linear system, which has a constant stiffness. This graph also shows that attainment of stationarity is delayed due to the nonlinearity. Again, this is a result of the softening behavior of the system, which slows down the propagation of signals at higher response levels, thus requiring more time for a fixed number of traversals of the beam. For the parameter values used in this set of analyses, ~ "" 2.5, so t = 0.5~, 6.0~, 12.0~, and 40.0~ approximately correspond to actual problem times of t = 1.25, 15, 30, and 100, respectively. As seen in Figure 4.18, the linear response becomes stationary near t = 10, while the mildest nonlinear case NL 1 becomes stationary near t = 35, and the remaining two nonlinear cases NL 2 and NL 3 reach stationarity at t = 60 and t = 80, respectively. Based on these observations, it appears that the concentration of 87 deformation near the base of the beam, as observed in Figure 4.18, continues as stationarity is reached. Engineering Implications 4.4.7 The above results suggest several important observations for engineering applications involving a softening material. First, in contrast to the hardening material, the softening material attenuates the propagation speed of the signal, and thus a given location in the structure will experience motion later than it would in a corresponding linear structure. Also, this slower propagation speed lengthens the time required for the system to reach stationarity, thus underscoring the importance of treating the response of a softening system as a nonstationary random process. Perhaps the most important characteristic observed in the response of a shear beam composed of a softening material is the concentration of the deformation near the base. As previously pointed out, this effect begins early in the nonstationary response, and persists as the response becomes stationary. From a structural engineering point of view, a uniform structure composed of this type of material requires substantially more ductility near the base than in the middle and upper sections if failure is to be prevented. This result has implications for structural design in a seismically active area. For example, if a large shear structure (such as a tall building) is to be located in an earthquake-prone area, it is important to design sufficient deformation tolerance into the lower floors to ensure that the building would safely survive an earthquake event. Since this study considered uniform stiffness and strength throughout the structure, the results for an actual building may vary if these conditions are not satisfied. 88 Alternatively, in a geophysical or earthquake engineering context, one might conceptualize the shear beam as a column of soil excited by stochastic motions of underlying bedrock during an earthquake. Often, the objective of such an analysis is to develop a "transfer function" relating bedrock motions to surface motions (represented by the free end of the shear beam). In this interpretation, the above results suggest that the "effective mass" for a softening material is much greater than for a linear material, since after only about 6 transit times the majority of the column is deforming very little and is simply acting like a large mass supported on a small spring. In this context, "effective mass" is taken as the mass of a single degree-of-freedom oscillator which represents the frequency transfer characteristics of the soil column. This increased effective mass would tend to shift the spectrum of the surface motions lower in frequency than the bedrock motions, and attenuate higher frequency inputs. For earthquake resistant design of large buildings, it is important to know the frequency range of possible inputs. Since tall buildings tend to have lower natural frequencies than short buildings, this shift of the ground motion frequency content toward lower frequencies leads to a more severe design environment than would be generated by a linear ground motion analysis. Allowance for this effect in the structural design therefore results in a more seismically safe structure. Again, it is emphasized that the current example assumes uniform stiffness and strength throughout the column. If actual conditions are substantially different from those considered here, then the magnitude of this frequency shifting effect may be altered. In summary, the presence of a softening material nonlinearity can markedly alter the system response characteristics, and in many cases these nonlinear effects render a linear analysis unconservative for design. Thus, a nonlinear analysis of 89 such systems is well justified. 90 Chapter 5 Concluding Remarks 5.1 Summary and Conclusions In chapter 1, the problem is specified and a general form of the governing nonlinear partial differential equation is given. This is followed by some general background on stochastic analysis, and a brief discussion of existing analysis approaches. The development of the general nonstationary continuous equivalent linearization procedure is presented in chapter 2. The nonlinearity is restricted to be symmetric and linearizable, and this leads to the general form of the equivalent linear continuous system. Although the equivalent linear parameters are functions of space and time, their explicit dependence is on the response statistics, and since these statistics depend on time in nonstationary problems the equivalent linear system is given time-varying coefficients. A system difference minimization criteria is defined as the difference in stress between the nonlinear and the equivalent linear systems. Minimization equations are then derived which give the equivalent linear parameters in terms of the current response statistics. Finally, an efficient numerical solution procedure is developed for the minimization equations. 91 Chapter 3 introduces a numerical implementation of the method derived in chapter 2. First, the finite element method is applied to the equivalent linear system to derive a set of discrete linear ordinary differential equations. Next, basis functions are specified for the equivalent linear parameters which alleviate the need to solve a large coupled system of linear equations, and allow determination of the equivalent linear parameters on an efficient element-by-element basis. The discrete equations are then transformed into a Liapunov equation for the response covariance matrix. Simplifications to this equation are shown to arise under conditions of white noise excitation. Chapter 3 concludes with a summary of the complete solution algorithm for the nonstationary response of a nonlinear continuous system. In chapter 4, the new method is applied to the solution of the one-dimensional wave equation with two types of constitutive nonlinearities: a hardening stressstrain law and a softening stress-strain law. As a physical conceptualization of this equation, the base-excited shear beam is chosen for study. It is demonstrated that results from the nonstationary continuous equivalent linearization method agree well with those of other methods on benchmark problems for which other solutions were available. The concept of the response profile is introduced for studying the way in which the stationary random response becomes established. This profile is shown as a plot of mean-squared relative displacement vs. x at a fixed time. The influence of the type of damping and mesh size are discussed in terms of this response profile. Later in chapter 4, a parameter study is described which considers systems with varying degrees of hardening constitutive nonlinearity. It is found that sharp gradients develop in the response profiles during the nonstationary phase of the 92 solution, and that this becomes more severe with increasing degree of nonlinearity. Engineering implications of these results are discussed. It is noted that this sharp response gradient implies a local concentration of deformation, and thus requires more ductility in the system to preclude failure than would a corresponding linear system. In this context, ductility is used in the sense of ability to withstand deformation without failure. Therefore, a linear analysis for design purposes would be unconservative for systems with this type of nonlinearity. Also, these sharp gradients are seen to disappear as the solution became stationary, thus highlighting the importance of considering the nonstationary portion of the response in an analysis. In the last part of chapter 4, a parameter study is described which considers systems with varying degrees of softening constitutive nonlinearity. It is observed that, as the response progresses, deformation becomes more concentrated near the base of the beam, and the remainder of the structure deforms very little. Engineering interpretations of these results are discussed in both a structural and an earthquake ground motion context. The most important conclusion for a structural system is that this concentration of deformation would require more ductility near the base of a softening nonlinear structure than in a corresponding linear structure, thus again making a linear analysis unconservative for design purposes. In an earthquake ground motion context, one might think of the shear beam as a column of softening soil excited by bedrock motions. Often, the objective of such an analysis is the definition of the frequency transfer characteristics of this soil column. The resulting "transfer function" can then be used to relate bedrock motions from a geological analysis to surface ground motions for an earthquake response calculation. In this application, the present results indicate 93 that the effective mass in these frequency transfer function calculations should be substantially larger than would be indicated by a linear analysis. Thus, in areas where the soil shear response is of a softening character, the surface ground motions will have more low frequency content than did the bedrock motions. In conclusion, the newly developed nonstationary continuous equivalent linearization procedure is seen to be a powerful tool for obtaining the random response of a nonlinear continuous system. It has been shown to perform well on two widely different types of nonlinearities. It is expected that the new method will produce useful engineering results for many nonlinear systems where the actual response process is nearly Gaussian. Although the numerical solution of random vibration problems is quite computationally intensive when compared with deterministic analyses, this method could handle meshes involving several thousand degrees of freedom using present supercomputers. This capability is expected to grow rapidly in the coming years given the current rate of progress in supercomputering power. The numerical implementation described herein couples the new method with the power of finite element analysis to allow nonstationary response solutions for a broad range of complex engineering problems. 5.2 Recommendations for Further Work In the course of an investigation such as this, it is inevitable that many more interesting topics are uncovered than can be explored by the investigator in the time available. The following paragraphs suggest some possible applications and extensions to the present work arising from these observations. One interesting application of the method developed in this thesis would be systems exhibiting hysteretic behavior. In one dimension, a hidden variable 94 model, such as that proposed by Iwan and Asano in reference [12], could be adapted for use with the present technique. In higher dimensions, hysteretic models for stochastic analysis are not well developed, and much remains to be done in this area. Another interesting application is the area of structural finite elements, i.e., elements with rotational degrees of freedom. This class of elements includes beams, plates, and shells, and allows the modeling of many structures of interest to mechanical and civil engineers. When combined with realistic nonlinear material models, these elements could allow construction of a very general analysis tool for the nonstationary stochastic response analysis of a wide variety of systems. A third area for further work involves coupling the techniques developed herein with a damage model. The objective here would be a unified structural analysis for both the stochastic response and the resulting damage level. Clearly, much work remains to be done before a thorough understanding of the response of nonlinear systems to stochastic excitation is achieved. 95 X ---- 71 I I W(X,t) 1 }--v (x,t> I I I I I I I I I I I I I I I I I I I I I 1--QI I t I I I I (·) (·) ....................................... ..... "''"""'''''''' Figure 4.1: Nonlinear Shear Beam Under Base Excitation L 96 STRESS-STRAIN CURVE FOR CUBIC HARDENING NONLINEARITY 5.0 4.0 3.0 ii '~ ..13 2.0 !:: t.O a: Iii 0.0 .... -t.o ~ .J %'/' IID z ?: / ' ' /// -2.0 // / -3.0 I -4.0 I -!1.0 -3.0 ---LINEAR DELTA•0.2 - - - DELTA•0.5 I I - - - - - OELlA•LO I 'i ' ' -s.o 0.0 LO 2.0 STRAIN Figure 4.2: Stress-Strain Curve for the Cubic Nonlinearity 3.0 ~ ~ ;;~ q (t) \ \ \ '\. _li f2 , '"" A~ f.I = NONLINEAR -~ m, :.L \\\ _li I RESTORING . m2 1---'v "( •~A. l" FORCE .") (o) N m, f-'v •i• •.A Figure 4.3: Discrete Analog to the Continuous Shear Beam I 'f. AlA 16xl ~ _l• mN co _, 98 WI r- f ..L : ' "' '. q ( t) m ":"' c -) (· I 8 8 " ,,,,,, Figure 4.4: Nonlinear Single Degree-of-Freedom Oscillator 99 tX X = 1.0 X DENOTES A FREE NODE ® DENOTES A NODE WITH CD SPECIFIED KINEMATIC BOUNDARY CONDITION ® DENOTES ELEMENT N w X= 0.0 q( t ) Figure 4.5: Finite Element Model to Represent a SDOF System 100 r I.LL.LLJ 5 G4 X DENOTES A FREE NODE ® DENOTES A NODE WITH SPECIFIED KINEMATIC 4 G3 DENOTES ELEMENT N ® 3 BOUNDARY CONDITION ® @ G2 ® 2 Gl CD w q( t ) Figure 4.6: Finite Element Model for 2" 0 Test Problem 101 w3 . ?~ ~ ~ q (t) w2 w,----1 c, ..:.J m, /.. ~ .. -,.,f, •) (• c2 ~ C3 • ..t• m2 't2 8 (· ~r- .r"/ . f3 m3 :'\ (. f.I = NONLINEAR RESTORING FORCE Figure 4. 7: Discrete Model for znd Test Problem with Hardening Nonlinearity 102 w3 w2 w,~ .,c 3 .:r ~2 ' ..:1 c, ' q ( t) ":'1 -:.- .. vf m, ""t I o) m2 • ...t. ." f ,. . 2 (• 0 {· ·yr3 •• ~A_ {:) [.) (o) ~ """' '"'"" f.I = NONLINEAR RESTORING m3 . FORCE Figure 4.8: Discrete Model for 3rd Test Problem with Hardening Nonlinearity 103 ••• ~ i 1.0 0.1 i I ~ i 0.1 -· o.• - - - CTS 1!1 liN 0.2 - 0.0 ... ~ a i o.o 10.0 - .... .... .... DIID'!TI! 1!:1 l.IN ..... .... ••• 0.1 I I ••• I -· ••• - - - CT!I El LIN - 0.2 0.0 0.0 10.0 ••• iB i I I - .... .... .... ..... .... ••• ••• 0.1 ~ -· ~ I DISCAI!TE El LIN 0.4 - - - CTS El LIN ••• - 0.0 o.o so.o a.o ,_ a.o - DIICAETE Ell UN .... ••• Figure 4.9: Comparison of Mean-Squared Displacements vs. Time for 3rd Test Problem 104 r x = I 4.0 I5 @) 14 @ x = I 2.0 I3 @ 12 @ x = I 0.0 x I I @ DENOTES A FREE NODE 10 ® DENOTES A NODE WITH SPEC I FlED KINEMATIC X= ® 8.0 9 BOUNDARY CONDITION 8 @ ® G) DENOTES ELEMENT N X= 6.0 7 ® 6 X ® = 4.0 5 4 X = 3 2.0 0 0 ® 2 X = 0 1-w 0.0 q( t ) Figure 4.10: Finite Element Mesh Used for Shear Beam in Example Problems 105 l.O ("" 0.1 ~ i! 0.1 0.7 i 0.1 I 0.1 g ~ I o.• 0.3 0.1 O.l I i I I EARLY TDE tooO.Oi 0-/CJ - - - IHI..... NABS PROP. -·- o.o 0.0 0.2 o.• MIIIIW..IZED DISTANCI ,_,.. BAll UI/U o.o 0.0 l.O 0.1 i 0.1 i 0.7 / ••• I ••• g i •••••• / / / / ---- /// O.l -~ 0.0 0.2 MIIIIW..IZED DISTANCE P'RDM IAil txiU ••• ••• ••• ••• i; 0.1 i! 0.7 I I IHI-., MASS PADP. STIF. PROP. ••• ••• u' / ••• I ••• ~ ,.,.,...------- Mm-TIME t-o.SIL.IC) 0.1 o.o /'' /'' l.O ./ ; I &TIF. PROP. o.• LATE TDe t-e CL1CJ 0.3 .......... 0.1 - - ' - IHI-., ••• ••• o.o MIIIIW..IZED OISTANC! ,.,... MSr - · - STIF. PADP . uvu ••• ••• ••• l.O Figure 4.11: Comparison of Normalized Response Profiles for Different Types of Damping 106 1.0 a ~ 0.1 i 0.7 i 0.1 I 0.5 ~ I T I I I o.o I I !III o.• 0.3 DilLY TINt: t-G.Ot 11./CJ tl i 0.2 0.1 - - l.ulll MESH - - MEDIUM MESH - - - IMAU MESH o.o a.o o.a o.• o.o 0.0 1.0 MOANALJ%!0 DISTANCE I'RDM .,._ 1.0 0.1 i a i 0.1 0.7 i 0.1 I 0.0 " o.• B N I 0.3 - - .......... 0.2 - 0. I - NI!DilM N!BH - - - IMAU. NI8H 0.0 0.0 0.2 o.• 0.1 0.1 1.0 ......... UEU DJITANC! """" 8&Sl!: IXIU 1.0 0.0 ~ i 0.1 0.7 i 0.1 I o.o ~ I 0.4 __ ....,_ .... 0.3 0.1 - 0.1 0.0 - MEDIUM NEIH - - - atAU.. MIEIH o.• 0.0 0.2 IGIIIALJZ!D DISTANCE , _ . .,._ o.o 0.1 1.0 Figure 4.12: Comparison of Normalized Response Profiles for Different Mesh Sizes lOi 1.0 ••• ~ ••• ~ I ~ ~ I 0.7 0.1 ••• ••• EAfL"'' fiNE t.ooO.Di G../Cl o.• -LINEAR ••• - - - NIH-JHEAA 1.501 - - NDtLltEAA 1.251 - - - - · NIJN.,JNENI U. OJ 0.1 o.o o.o o.z ~ZEQ DlST~ o.• ""'* 1!1ASE ••• o.o 1.0 1.0 ••• 0.0 0.7 0.0 o.o o•• 0.3 --LINEAR - o.z - NDtL.ltEAR (.25) - - - HOHL.lHEAA 1.501 - - · - · HOt&.ltEAA It. Ol 0.1 0.0 o.o 0.2 •.• ......UZ!D DlllloHCI , . . . . . . . ••• o.o 1.0 1.0 L.Al! TIME t ... DJCI 0.-1 ~ il 0.1 - - L I..... 0.7 HOtLitEAA 1.251 - - IDLJNEMJ J,!KJ) a I i ~ I - - - - NllfLINEAR IS.OI ••• o.o o.• 0.3 - o.z 0.1 o.o •.• 0.2 ·-·==·;:,:,;;:::=:::=-.:::=:::::::::::::::::••• ••• •.• 1.0 ..-....uza DISTANCE , _ . . . . Figure 4.13: Comparison of Normalized Response Profiles for Varying Degrees of Hardening Nonlinearity 108 CUBIC NONLINEAniTY FnEE END VAn l!HE ltiSTOnY 1.0 0.9 ~ .... c O.B nla: 0.7 atn ~,. ~ . ~ 0.6 LINEAR - NOIILIIIEAR I. 25) - - - NONLitiEAR (.50) ----- llotiLitlEAR II. oJ - 0.5 0.~ :0: a: c z 0.3 0.2 - ~ _.... -----------------......---------------____--~1'/ -_,.~--:::-----------------------/ o .I '"- 0.0 0.0 10.0 20.0 30.0 ~0.0 50.0 l!HE Figure 4.14: Comparison of Free End Mean Squared Displacement Time Histories for Varying Degrees of Hardening Nonlinearity 109 STRESS-STRAIN CURVE FDR AIICT ANGENT SOFTENING NOILINEAIIITV 1.!1 1.0 0.!1 -____________ __ _ ------------- :,::::. ,~ :;" --- ... 0.0 -o.!i -1.0 -1.!1 -3.0 --... --------------·::::::.----~ ------ -- - -2.0 LINEAR ..... -1.0 (G/TAUY)•1. 0 (G/TAUY)•2.0 - · · - · (G/TAUYl•3.0 0.0. 1.0 2.0 STRAIN Figure 4.15: Stress-Strain Curves for Arctangent Softening Nonlinearity 3.0 110 R!SPONSE PFIDFILE AT T•0.5ft./C:I ~ ~ a 1.0 --~ o.s y ./~ .<'0 /1; /I; //; //; !/; !/; //; //; 0.0 ~ 0.7 .;w 0.8 ~ ~ % 0 .• ~ ~ o.< ~ z 0.] 0.2 .... ... ' Lltll&R HL I !1~, /~ 0.1 0.0 0.2 0.' 0.0 NOAMALIIEO DISTANCE FRON BASE (X/\.) 0.8 0.8 '' RESPONSE PFIDFJLE AT T•e lL/C) ~ 0.0 y 0.0 ~ 5 M a g .; w ~ ~#- / 7.0 0.0 -·--- .t' 0.0 ••• !' ' 2.0 _,-· ,/ ~· - - LJN(o\R HL 3 __.-·--· --·-· ./·/' i/ 1.0 .... .... -·-·---------- ,1 I II .// 3.0 o.o --- / % ~ ~ ~ z .. -·-----··---------------. . . . ...................·· - ---- --- - o.o 0.2 o.• NOANALllED DISTANCE FROM BASE 1XILI 0.8 0.8 1.0 Figure 4.16: Comparison of Normalized Response Profiles for Varying Degrees of Softening Nonlinearity: Early Times 111 RESPOHSI PROFIL! AT T•il L/1:1 I " ~ ~ i § i te.o __.. .. t•.o _....... -· 12.D (' fD.D L1).JE4R ~· / o.o ••• II / ___. - ......... !/ ii A.O 2.0 o.o -·---------------------------------- li/ il - - a.o :: -------------------------- -- ----- - -o.2 o.• 0.8 0 • 1.0 HOAMALlZED OlSTAHC! FROM lASE (X/LI ~ 2!.0 i ,--/----- 20 0 ~ 1 .._. ~· '.I _.....--------------------------- to.o 0.0 LlNI!AIII 1 15.0 !.0 --------------------------------- I.I I./ ii " o.o _/- - - -- o.z "'CIAMALlZ!D DISTANCE F'FIDM lAS! IX/U --- -----e.• ••• 0.8 1.0 Figure 4.17: Comparison of Normalized Response Profiles for Varying Degrees of Softening Nonlinearity: Late Times 112 FREE END VARIANCE TIME HISTORIES 25.0 20.0 / __ /_.---· --- --- ------- -------LI"EAR NL I NL 2 _/-15.0 ,/ 10.0 / ,1' ' / I / ,/ I / ,/ ,/ ;/__'.,.-- 5.0 0.0 0.0 ·- / ilL 3 / 20.0 _/_..--- - ~- -- ------------------- ---~0.0 60.0 BO.O 100.0 TIHE Figure 4.18: Comparison of Free End Variance Time Histories for Varying Degrees of Softening Nonlinearity 113 References [1] Caughey, T. K., "Derivation and Application of the Fokker-Planck Equation to Discrete Nonlinear Dynamic Systems Subjected to White Noise Excitation," Journal of the Acoustical Society of America, 35, (1963), p. 1683. [2] Caughey, T. K., "Nonlinear Theory of Random Vibrations," Advances in Applied Mechanics, 11, (1971), p. 209. [3] Lin, Y. K., Probabilistic Theory of Structural Dynamics, McGraw Hill, 1967. [4] Caughey, T. K., "On the Response of a Class of Nonlinear Oscillators to Stochastic Excitation," Proc. Colloq. Intern. du Centre National de Ia Recherche Scientifique, 148, (1964), p. 393. [5] Yong, Y., and Lin, Y. K., "Exact Stationary Response Solutions for Second Order Nonlinear Systems Under Parametric and External White Noise Excitation," Journal of Applied Nfechanics, 54, (June, 1987), p. 415. [6] Crandall, S. H., "Random Vibrations of Systems with Nonlinear Restoring Forces," Proceedings of the International Symposium on Nonlinear Vibrations, No. 1, (1961), p. 306. [7] Crandall, S. H., "Perturbation Techniques for Random Vibrations of Nonlinear Systems," Journal of the Acoustical Society of America, 35, (1963), p. 114 1700. [8] Booton, R. C., "The Analysis of Nonlinear Control Systems with Random Inputs," IRE Transactions in Circuit Theory, 1, (1954), p. 32. [9] Caughey, T. K., "Response of a Nonlinear String to Random Loading," Jour- nal of Applied Mechanics, 26, (1959), p. 341. [10] Iwan, W. D., and Mason, A. B., "Equivalent Linearization for Systems Subjected to Non-Stationary Random Excitation," International Journal of Non- linear Mechanics, 15, (1980), p. 71. [11] Iwan, W. D., and Krousgrill, C. M., "Equivalent Linearization for Continuous Dynamical Systems," Journal of Applied Mechanics, 50, (June, 1983), p. 415. [12] Asano, K., and !wan, W. D., "An Alternative Approach. to the Random Response of Bilinear Hysteretic Systems," Earthquake Engineering and Struc- tural Dynamics, 12, (1984), p. 229. [13] Caughey, T. K., "Random Excitation of a System with Bilinear Hysteresis," Journal of Applied Mechanics, 27, (1960), p. 649. [14] Baber, T. T., and Wen, T. T., "Stochastic Response of Multistorey Yielding Frames," Earthquake Engineering and Structural Dynamics, 10, (1982), p. 403. [15] Krousgrill, C. M., "A Linearization Technique for the Dynamic Response of Nonlinear Continua," Ph. D. Thesis, California Institute of Technology, Pasadena, California, September, 1980. 115 [16] Spanos, P-T. D., "Linearization Techniques for Nonlinear Dynamical Systems," Ph. D. Thesis, California Institute of Technology, Pasadena, California, September, 1976. [17] Kazakov, I. E., "Statistical Analysis of Systems of Multidimensional Nonlinearities," Automation and Remote Control, 26, (1965), p. 458. [18] Atalik, T. S., and Utku, S., "Stochastic Linearization of Multi-degree of Freedom Nonlinear Systems," Earthquake Engineering and Structural Dynamics, 4, (1976), p. 411. [19] Hughes, T. J. R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, 1987. [20] Press, W. H., Flannery, B. P., Teukolsky, S. A., and Vetterling, W. T., Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, 1986. [21] Jameson, A., "Solution of the Equation AX+ XB = C by the Inversion of an M X MoranN x N Matrix," SIAM Journal of Applied Mathematics, 16, (1968), p. 1020. [22] Muller, P. C., "Solution of the Matrix Equations AX+ XB = -Q and STX + XS = -Q, SIAM Journal of Applied Mathematics, 18, (1970), p. 682. [23] !wan, W. D., and Yang, I-M., "Application of Statistical Linearization Techniques to Nonlinear Multidegree-of-Frcedom Systems," Journal of Applied Mechanics, 39, (June, 1972), p. 545. 116 [24] Yang, 1-M., "Stationary Random Response of Multidegree-of-Freedom Systems," Ph. D. Thesis, California Institute of Technology, Pasadena, California, June, 1970. [25] Gradshteyn, I. S., and Ryzhik, I. M., Table of Integrals, Series, and Products, Corrected and Enlarged Edition, Academic Press, 1980, p. 338. 117 Appendix A Validation Calculations A.l Test Problem One: Stationary Response of a One Element System As a first check on the accuracy of the nonstationary continuous equivalent linearization method, a problem with only one degree-of-freedom is solved for the stationary mean-squared displacement. The stationary value is calculated as the long-time asymptotic value of the nonstationary displacement variance time history. Obviously, carrying out the nonstationary solution for a long time is a very inefficient way to calculate the stationary solution, but it serves as a good check on both the solution algorithm and the program coding. This test problem represents the simplest of the nonlinear problems considered. A schematic diagram of the discrete system is shown in Figure 4.4. The equation of motion for this nonlinear oscillator is mw + cw + f(w) = -mij(t), (A.l) 118 where (A.2) and • w is the displacement of the mass relative to the base, • m is the mass, • c is the viscous damping, • k is the linear spring stiffness, • 8 is the nonlinearity parameter, and • ij(t) is the applied base acceleration. For this example, ij is taken as Gaussian white noise with zero mean and constant power spectral density 8 0 • Dividing through the above equation by m yields the standard form (A.3) with the linear natural frequency Wn given as k w2n= - - m (A.4) and the fraction of critical damping 1 given as I= c , 2vkm (A.5) The solution used for comparison is an approximate closed-form analytical result developed by Iwan and Yang for the SDOF Duffing oscillator under stationary random Gaussian white noise excitation. In reference [23], they have 119 given an approximate solution for the mean-squared relative displacement under stationary conditions as (A.6) Yang [24] has shown that the above estimate is within 7.5% of the true RMS stationary solution for arbitrarily large nonlincarities, so this formula provides a good benchmark for evaluation of the new solution method. The finite element model corresponding to this simple system is shown in Figure 4.5. For a one element model of unit length (L' = 1), a unit stiffness k = GL~· and a unit mass m = tpA,L' were used. Damping corresponding to 10% of critical was added to yield a finite stationary displacement variance. A unit power spectral density of base acceleration S 0 was used as input, and the nonlinearity parameter was chosen as 8 = 1.0. The influence of this nonlinearity is discussed below. The degree of nonlinearity present in this model can be assessed by examining the ratio of the restoring force provided by the nonlinear term to that provided by the linear term at the rms displacement response level Uw, where (A.7) A simple calculation shows that this ratio of the nonlinear restoring force to the linear restoring force at the rms displacement response level is given by 8u~. Substitution of the above model parameters into the analytical approximation (A.6) yields an approximate displacement variance E [w 2 ] of 2.127. The nonstationary equivalent linearization calculation yields a stationary variance of 2.128. At this response level, the nonlinear restoring force is 2.1 times as large as the linear restoring force, therefore this system is severely nonlinear. The differ- 120 ence in variance of relative displacement between the two stationary solutions is less than 1.0%, well within the accuracy requirements for almost any engineering analysis. A.2 Test Problem Two: Stationary Response of a Four Element System The second validation problem presented also addresses only the stationary solution, but for a mesh with multiple elements. Thus, this problem represents the next logical step in increasing complexity. The mesh used here consists of 4 elements with linear shape functions. A diagram of the finite element mesh is shown in Figure 4.6. It is important to note a modeling detail at this point. The 4th element in the finite element mesh of Figure 4.6 is used only to contribute mass to node 4. A simple element mass lumping scheme is used which apportions one-half of the mass of an element to each of the 2 nodes of that element. Thus, if a node borders on two elements each having a total mass of m, then that node receives !m from each element for a total nodal mass of m. On the other hand, if a node borders only one element, then it only receives half as much mass, so its total nodal mass is !m. If it is desired that each active node receive equal mass, then an additional element with no stiffness can be added beyond the last active node, and the new end node made inactive. The presence of element 4, therefore, leads to equal amounts of mass being lumped at nodes 2, 3, and 4, allowing the analogous discrete model to also have equal masses. The degree of freedom associated with node 5 is deleted, and since G4 = 0.0, element 4 contributes no stiffness to the 121 problem. This modeling technique allows the total mass lumped at any node of the finite element model to be easily adjusted. The corresponding discrete model for this test problem is shown in Figure 4.7. It consists of three equal masses interconnected by three nonlinear springs and three linear viscous dampers. The configuration of the dampers in this discrete model is not exactly analogous to the mass proportional damping used in the finite element model. The proper configuration is shown in Figure 4.8. The error incurred in using the discrete model of Figure 4. 7 is that a small contribution to the right-hand-side load term from the damping coefficients times the input base velocity is neglected, i.e., is not present in the finite element calculations. However, an analytical estimate of the stationary variance of the system of Figure 4. 7 is available, so it is interesting to compare these results with those of the numerical procedure. The restoring force for the i 1h nonlinear spring is computed from J;(y;) = k;(y; + 8y[), (A.8) where k; is the stiffness coefficient of the i 1h spring, 8 is the nonlinearity parameter (assumed constant for all springs), and y; is the relative displacement between the i1h mass and the (i- 1)'1 mass. Once again, this system will be subjected to zero mean Gaussian white noise base excitation. The benchmark solution for this problem is an approximate analytical result also developed by Iwan and Yang in reference [23]. They gave the following expression for the variance in the relative displacement between mass i and mass (i- 1) : 1 [ 1- \!I 1 + 1 2 tri5So a 2 . = ---= -- ] . y, 615 1 0.1k; (A.9) 122 Once again, it has been shown by Yang in reference [24] that this approximate analytical solution remains within 7.5% of the true RMS displacement response solution for arbitrarily large nonlinearities. In the finite element model, the parameter values were chosen to vary between elements to verify that capability. The parameters used were Gl~' = 3.0, Gl~" = 2.0, and Gl~' = 1.0. Rayleigh damping is added to the model, with the damping matrix taken as proportional to the mass matrix, so C = 0.10M, (A.lO) where C denotes the global damping matrix and M denotes the global mass matrix. From these parameter values for the finite element model, the discrete model parameters are obtained. The spring stiffness values are unequal, with k1 = 3.0, k 2 = 2.0, and k3 = 1.0. The three masses m; are all taken as unity, and the damping coefficients c; are chosen as 0.1 for all three dampers. The power spectral density of the input noise process is S 0 = 2.23. The performance of the new solution method on this benchmark problem is examined for both a mild nonlinearity (5 = 0.1) and a more severe nonlinearity (b = 1.0). The calculated stationary values of the mean-squared relative displacements E [yfl are compared against the approximate analytical results in Table A.l. For comparison, the ±7 .5% values for the analytical solution are also shown, since the true solution lies within these bounds. To gain some insight into the size of the nonlinearity, it is helpful to examine the restoring force/; evaluated at the RMS relative displacement ay,· As a measure of the nonlinearity, consider once again the ratio of the nonlinear restoring force to the linear restoring force at the rms displacement response level. It is easily shown that this ratio is given by bay, 2 • Values of this ratio for each nonlin- 123 Severe Nonlinearity Mass Anal. Mild Nonlinearity No. Est. Anal. Camp. Anal. Camp. Low (-7.5%) 2.92 - 1.75 - 2.72 2.83 1.63 1.69 High ( +7.5%) 2.52 - 1.51 - Low (-7.5%) 3.29 - 1.85 - 3.06 3.05 1.81 1.82 High (+7.5%) 2.83 - 1.67 - Low (-7.5%) 4.00 - 2.33 - 3.72 3.52 2.17 2.02 3.44 - 2.01 - 1 2 3 High ( +7.5%) Table A.1: Summary of Results: RMS Relative Displacements y; ear spring are shown in Table A.2 for both the mildly nonlinear and the severely nonlinear cases. Comparison of the computational and analytical results shows that the nonstationary equivalent linearization technique performs quite well, yielding results that are well within the accuracy requirements for engineering analyses even for a severely nonlinear system. 124 Ratio of Non linear Spring Force Spring to Linear Spring Force Number Mild Nonlinearity Severe Non linearity 1 0.74 2.7 2 0.94 3.3 3 lA 4.7 Table A.2: Ratio of Nonlinear to Linear Restoring Forces A.3 Test Problem Three: N onstationary Response of a Four Element System The final validation calculation presented in this section considers the nonstationary response of a multiple element shear beam subjected to suddenly applied base excitation. The response is calculated from zero initial conditions out to its stationary value. This problem tests all features of the new algorithm and program, including the solution accuracy during the nonstationary portion of the response. The finite element mesh used for this problem is identical to that shown in Figure 4.6 for the second test problem. The analogous discrete model for this system is shown in Figure 4.8. It is similar to the model used in the previous test problem, except that the dampers are now connected between each mass and the base, rather than between each mass and ground. This modification allows the discrete system to exactly represent the mass-proportional Rayleigh damping used in the continuous system solution. The nonlinear spring restoring forces are calculated from equation (A.8), 125 and the input base acceleration is taken as zero mean Gaussian white noise, suddenly applied at t = 0. The reference benchmark for this test is a discrete nonstationary equivalent linearization solution using the method developed by Iwan and Mason in reference [10]. Their technique gives an expression for each element of the equivalent linear stiffness matrix in terms of the current response covariance matrix. The resulting set of discrete equivalent linear equations are next transformed to first order state space equations, and finally to a Liapunov equation for the state space covariance matrix. The application of their technique to the above system of nonlinear ordinary differential equations is briefly outlined below. The equations of motion for this discrete system are m1tii1 + c1w1 + k1 [w1 + Swi]- kz [(wz- w1) + S(wz- w1) 3] = -m 1 q(t) (A.ll) mzWz + CzWz + kz [(w 2 - w1) + 8(w 2 - w!) 3] - k3 [(w3- wz) + S(ws- wz) 3] = -mzq(t) m3tii3 + c3w3 + k3 [(w3- Wz) + 8(w3- wz) 3] = -m3q(t), (A.12) (A.13) where w; denotes the displacement of the i'h mass relative to the base. This set of equations can be combined into a matrix equation of the form Mw + h(w, •v) = -q(t)Mf:, (A.14) where f: is defined by (A.15) 1 126 Now, define the state vector s as (A.16) Using this notation, the second order differential equation (A.14) can be written as a first order state space equation of higher dimension as s = h(s) + ij(t)r, (A.17) r=(o)· i· (A.18) where r is defined as Following the procedure described in section 3.4, equation (A.17) above can be written as a first order matrix Liapunov equation of the form (A.19) where Q is the state space covariance matrix given by (A.20) 8 0 is the power spectral density of the white noise process q(t), and A is the equivalent linear system matrix given, in components, by ah A;=~· uy; (A.21) This equation is integrated numerically using a 4th order Runge-Kutta method to find the time history of the covariance matrix Q. For this analysis, the same stiffness, mass, and damping values were used as for the previous test case. The power spectral density of the input base acceleration for this test was So = 6.8. In order to test the method on a severely 127 nonlinear problem, the nonlinearity parameter 1! was chosen as 1.4. Calculations analogous to those described for the previous test problem show that the nonlinear term provides from 1.5 to 2.2 times as much restoring force as does the linear term at the stationary response level. Thus, the nonlinearity in this model is quite severe. The results of the calculations are presented in Figure 4.9 as plots of the time history of the normalized relative displacement variance at locations corresponding to the three masses of the discrete model. For each of these plots, the actual displacement variance was normalized by the peak value of the discrete variance. In these graphs, the solid line shows results calculated using the new nonstationary continuous equivalent linearization method, while the dashed line represents results calculated using the discrete equivalent linearization theory. The top graph corresponds to mass .l, the middle graph to mass 2, and the bottom graph to mass 3. In the plots, the two lines are coincident in most places. Thus, these graphs illustrate very good agreement between the discrete and continuous linearization methods. Based on the results of the above three validation calculations, the new nonstationary continuous equivalent linearization method is judged to produce reasonably accurate response statistics for nonlinear hardening systems. The long-time asymptotic stationary results calculated with the new method agree well with stationary statistics calculated from an approximate analytical solution to an analogous discrete system. The nonstationary results agree well with those calculated using a discrete equivalent linearization method. 128 A.4 Test Problem Four: Stationary Response of a One Element System with Softening Constitutive Law In view of the fact that most of the solution procedure for the arctangent softening nonlinearity is identical to that used for the cubic hardening nonlinearity, many of the test problems used for the cubic hardening system apply here also. Therefore, only the changes required for the arctangent softening nonlinearity are subjected to validation calculations in this section. In order to verify the formulation of the arctangent softening nonlinearity and the resulting minimization equations described above, a one-element (singledegree-of-freedom) test problem is solved for its stationary solution. This solution obtained using the nonstationary continuous equivalent linearization method is compared against an approximate analytical result derived for a stationary discrete SDOF system. The results of this comparison are discussed below. The equation of motion for the base excited nonlinear SDOF oscillator is mw + cw + f(w) = -mxo(t), (A.22) where m is the mass, cis a damping coefficient, f(w) is the nonlinear restoring force, and x0 (t) is the specified base acceleration. The nonlinear restoring force f (w) is defined by (7rkow) 2 f ytan _1 - - , !( w ) =7r 2/y (A.23) where /y is the yield force, and k 0 is the initial stiffness. The benchmark solution used for comparison is an approximate result developed by Yang and Iwan in reference [23]. They derived the following expres- 129 sion for the mean squared displacement of an arctangent softening SDOF system under Gaussian white noise base excitation: E [ 2] w = 2(7rSom) (cfvF 2 ' (A.24) where S0 is the power spectral density of the input white noise process. In reference [24], Yang shows that the RMS displacement calculated from equation (A.24) is always within 11.4% of the exact solution, so this expression offers a ' useful benchmark for verification of the solution procedure for the softening non- linearity. For this comparison, a one element model was constructed for solution by the nonstationary continuous equivalent linearization method. The parameters used in this one element model of unit length are a£~' = 4.0, pA.L• = 2.0, and !2o. = 2.0. The stationary solution is obtained as the large time value of the •• • nonstationary solution. While not efficient, this method serves as a useful check on the late-time accuracy of the procedure. Using relations derived in a previous section, the analogous SDOF discrete system parameters are found to be k 0 = 4.0, m = 1.0, and ~ = 2.0. Damping corresponding to 5% critical (based on the initial stiffness) was added to both systems. Each model was subjected to zero mean Gaussian white noise base excitation with unit power spectral density. The calculated variance in relative displacement using the numerical procedure is 101.9, while the analytical approximation (A.24) yields a displacement variance of 123.2. This represents a discrepancy of 9% in RMS relative displacement, which is within the accuracy band of the approximate solution. Therefore, the continuous equivalent linearization procedure agrees quite well with this analytical solution. This small problem exercises all 130 of the new minimization equations derived for the arctangent softening system, so this good agreement between the analytical and numerical solutions indicates that the theory and numerical implementation are performing as expected. From the above results for two types of constitutive nonlinearities, it is believed that the new continuous linearization technique produces solutions of sufficient accuracy to be useful in practical engineering analysis. 131 Appendix B Evaluation of the Integral for the Arctangent Softening Nonlinearity In the course of deriving the minimization equations for the arctangent softening nonlinearity, it is necessary to evaluate the integral I= ! !"' oo -oo 1 (11" G0 c:) 2 -r.itan- - p(c:, i)didc:, -oo 7C 2 ry (B.l) where P(e:,i) = 2 7rU~Ui ~l - p2 exp [ (B.2) and E [c:i] p=--. (B.3) Oe(Ji Define some constants: (B.4) 132 (B.5) C"' = '"' C4 O'e(Je -- (B.6) 2Ty (B.7) 2(J =7rC1 1rGo 2ry (B.S) cs = --. Using these constants, the problem can now be written (B.9) or (B.lO) Define the i: integral as 11 by I 1-- J e:e. -..L(''-L<) "' ·J d.s. oo "> (B.ll) -co Using reference [25], entry 3.462.6, it follows that (B.12) Define two additional constants (B.13) (B.l4) so equation (B.12) becomes (B.15) Next, substitute this equation (B.l5) back into the original double integral equation (B.lO) to get the remaining integral (after some manipulation) as (B.l6) 133 Define c8 as 1 - -2 - Cg = (B.17) CG, c2oe so equation (B.16) becomes (B.18) From reference [24], it is known that ! 00 0 0 xtan- 1 (x) e-a"x" dx = - 11' 2a 2 -oo erfc(a). (B.19) Letting x = cse, (B.20) 2 - Cg (B.21) defining a = """"":)'", C"5 and using the above result (B.19) yields c,-C7- 1l' e a' er f'c (a ) . I =c52 ?.... a 2 (B.22) Simplifying out the intermediate constants c1 through c8 , it is easily shown that equation (B.22) becomes (B.23) I= /f;rypa,e"'erfc(a), where a= VZrv 1CG 0 a, (B.24) . Finally, from the definitions in chapter 4, Got:)] [2 I=E -ryitan- 1 (1l' --\2 1C (B.25) , Ty so the desired result is Goc:)] = y(2-:;/vpa,e" erfc(a). c:tan _ (1l'-27 2 . E [-r 11'" 2 1 y This equation appears as equation (4.97) in chapter 4. (B.26)