Nonlinear Response of Structures in Regular and Random Waves by Arthur William Lipsett B.Sc, The University of Alberta, 1974 M.Sc, The University of Alberta, 1976 A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of D O C T O R O F P H I L O S O P H Y in The Faculty of Graduate Studies Department of Civ i l Engineering We accept this thesis as conforming to the required standard The University of British Columbia August 1985 © Arthur William Lipsett 1985 In presenting this thesis in partial fulfilment of the requirements for an ad-vanced degree at The University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Civ i l Engineering The University of British Columbia 2324 Main Mal l Vancouver, British Columbia Canada V6T 1W5 August 1985 ii Abstract The problem of the dynamics of a flexible offshore structure in either a reg-ular or random sea is considered in this thesis. A simple single degree of freedom model of the structure is assumed and the relative velocity formulation of the Morison equation is used to describe the fluid force. The resulting equation of motion is a nonlinear ordinary differential equation with either harmonic or stochastic forcing depending on the wave de-scription. Solutions are obtained for regular deterministic waves by numerical integration, various linearization methods and a new perturbation method developed in this thesis. The numerical solution is used to assess the accuracy of each of the approximate solution methods. Of these, the perturbation method is found to give the best approximation to the numerical solution over the complete frequency range of interest. For random seas the response spectrum and the mean square response are obtained by various linearization methods, the method of equivalent linearization, and by the new perturbation method. The perturbation method and the method of equivalent linearization are very similar in that they both yield the same values of effective damping. Comparison of the results obtained by a numerical simulation method with the results of the perturbation method and the widely used method of equivalent linearization shows that the perturbation method gives a better estimate of the response mean square value than does the method of equivalent linearization. For all of the approximate solution methods that are discussed it was found that the use of Hermite polynomials to represent the solution is very effective in obtaining various expected values required in the computational procedure. In addition to the average response statistics, such as the response mean square value, the probability density of the response is also considered. It is well known that the response of a linear system to Gaussian forcing is itself Gaussian. The wave force given by the Morison equation is non-Gaussian and therefore the response is also non-Gaussian but of unknown form. The hypothesis that for a linear equation, the probability density of the response is of the same form as the probability density of forcing, even for the case of non-Gaussian forcing, is investigated and verified using the results of numerical simulations. Design considerations of interest which follow from the response probability density are also discussed. iii Contents Abstract i i Figures v Tables vi i Symbols vi i i Acknowledgement x i 1 Introduction 1 2 Deterministic Wave Forcing . . . . . . . 11 2.1 Equation of Motion 13 2.2 Parameter Values 17 2.2.1 Damping 17 2.2.2 Inertia and Drag Coefficients 18 2.2.3 Keulegan-Carpenter Number 20 2.2.4 Reynolds Number . . 21 2.2.5 Mass Ratio . . 21 2.2.6 Frequency Ratio 21 2.3 Solution Methods 22 2.3.1 Numerical Integration 22 2.3.2 Linearization Methods . 23 2.3.3 Perturbation Method 28 2.4 Results 32 3 Random Wave Forcing 49 3.1 Time Dependent Damping Linearization 51 3.2 Equivalent Linearization 58 3.3 Perturbation Method 62 3.4 Numerical Simulation 66 3.5 Response Spectra and Mean Square Results 71 3.6 Probability Density of Response 89 3.6.1 Constant Drag Damping 93 4 Summary and Conclusions 102 4.1 Regular Waves 103 4.2 Random Waves . 105 4.2.1 Response Spectrum and Mean Square Value . . . . . 105 4.2.2 Response Probability Density 106 4.3 Recommendations for Further Study 107 References 108 Contents jy Appendices I l l A l Review of the Solution of Linear ODE's with Harmonic Forcing . I l l A2 Review of Linear Wave Theory 114 A2.1 Wave Spectra 118 A3 Review of Linear Spectral Theory 121 A4 Hermite Polynomials 128 A5 Convolution Spectral Calculations 135 V Figures Fig. 2.1 Single degree of freedom model 12 Fig. 2.2 Inertia coefficient as a function of Re and K (Sarpkaya and Isaacson 1981) 19 Fig. 2.3 Drag coefficient as a function of Re and K (Sarpkaya and Isaacson 1981) 19 Fig. 2.4 Frequency response for a = 0 and e = 0.05 34 Fig. 2.5 Frequency response for a = 0 and e = 0.15 . 3 4 Fig. 2.6 Frequency response for a = 0 and e = 0.25 . . . . . . 35 Fig. 2.7 Frequency response for a = 0 and s = 0.50 35 Fig. 2.8 Frequency response for perturbation method with e = 0.05 . . 3 9 Fig. 2.9 Frequency response for perturbation method with e = 0.15 . . 3 9 Fig. 2.10 Frequency response for perturbation method with s = 0.25 . . 40 Fig. 2.11 Frequency response for perturbation method with e = 0.50 . . 40 Fig. 2.12 Effective damping for deterministic forcing 41 Fig. 2.13 Time history of the response for frequencies O = 1/3 and fl = 1 . 43 Fig. 2.14 Frequency response for e = 0.05 and a = 0.2 44 Fig. 2.15 Frequency response for e = 0.05 and a = 0.5 44 Fig. 2.16 Frequency response for s = 0.05 and a = 0.8 . . . . . 4 5 Fig. 2.17 Frequency response for e = 0.05 and a — 1.2 . . . . . 4 5 Fig. 2.18 Effective damping for various current ratios 47 Fig. 2.19 Static offset for various current ratios 47 Fig. 3.1 Pierson-Moskowitz spectra for Vw = 11.6 . 7 4 Fig. 3.2 Pierson-Moskowitz spectra for Vw = 21.8 75 Fig. 3.3 Pierson-Moskowitz spectra for Vw = 31.6 76 Fig. 3.4 System transfer function for various effective damping values . . 77 Fig. 3.5 Response spectra for a = 0, Vw = 11.6, e — 0.25 80 Fig. 3.6 Response spectra for a = 0, Vw = 21.8, e = 0.05 . . . . . 81 Fig. 3.7 Response spectra for a = 0, Vw = 31.6, e = 0.05 82 Fig. 3.8 Response spectra for a = 0, Vw = 31.6, e — 0.25 82 Fig. 3.9 Relative velocity and wave spectra for Vw = 31.6 and e = 0.25 . . 85 Fig. 3.10 Response spectra for Vw = 21.8 and a = 1.0 88 Figures yj Fig. 3.11 Response spectra for Vw = 31.6 and a = 1.0 . . . . . . 88 Fig. 3.12 Response probability density e = 0.05 96 Fig. 3.13 Response probability density e = 0.15 96 Fig. 3.14 Response probability density s = 0.25 97 Fig. 3.15 Response probability density e = 0.50 97 Fig. 3.16 Response peak probability density for a = 0,Vw = 31.6, and e = 0.15 100 Fig. A2.1 Definition sketch for linear wave theory 115 Fig. A2.2 Pierson-Moskowitz wave height spectrum Vw — 31.6 . . . .. 120 Fig. A5.1 Convolution spectra for harmonic wave 137 vii Tables Table 2.1 Solution methods for deterministic forcing 33 Table 2.2 Static response for various nonlinearity and a = 0 . . . . 36 Table 2.3 Static response for various current ratios and e = 0.05 . . . 46 Table 3.1 Solution methods for random forcing . 72 Table 3.2 Standard deviation of response for a = 0 79 Table 3.3 Effective damping for a = 0 84 Table 3.4 Standard deviation of response for e = 0.05 86 Table 3.5 Effective damping for e = 0.05 89 Table 3.6 Simulation results for second and fourth response moments . . 98 viii Symbols A constant in P - M spectrum = 0.0081, constant Ak = amplitude of kth sinusoidal component B = constant in P - M spectrum = 0.74, constant B{nun2) bispectrum c = linearization constant c(nltn2M = trispectrum drag coefficients = inertia coefficients D — cylinder diameter H wave height HP) = frequency response function Hek(z) = Hermite polynomial of order k with argument z J Jacobian K Keulegan-Carpenter number L = wavelength M = number of sinusoidal components R(T) correlation function R(TI,T2) bicorrelation function tricorrelation function Re = Reynolds number sift) spectral density function T = period U water particle displacement un discrete value of U at t = nAt vw = wind speed for P - M spectrum a, b — constants, coefficients c = structural viscous damping per unit length, coefficient filter coefficients d water depth, constant e = mean square error, constant f force, function fm = force maxima 9 = function, constant, acceleration of gravity h function M O = impulse response function k = structural stiffness per unit length, wavenumber m = cylinder mass plus added mass per unit length m 0 = cylinder mass rep unit length m 3 third moment = fourth moment = probability density function of x Symbols ix p(x, x') = joint probability density of (z, x') p(x, x',x") = joint probability density of (x, x',x") q — function r = relative displacement, zero mean relative displacement f = constant part of relative displacement f = relative displacement s = vertical coordinate measured upwards from seabed t = time, thickness u = zero mean water particle displacement w = zero mean unit variance variable = '"o/°'rJ w(t) = white noise function x = absolute cylinder displacement x m = response maxima z = vertical coordinate, zero mean unit variance variable = u'/ = phase angle, velocity potential $ = constant spectral density of white noise rpi, tp2 = parameters u = frequency fl = dimensionless frequency = U/OJ^ < > expected value Symbols x superscripts * *r c subscripts D eff f k, v I m N u x 0,1,2 a dimensional quantity differentiation with respect to t convolution spectrum of order r complex conjugate average over one period drag effective value forcing index inertia maximum value natural value in water refers to a variable based on water particle kinematics refers to a response variable terms in a perturbation series xi Acknowledgement The author wishes to thank Dr. M . Isaacson for his advice and assistance in the preparation of this thesis. The financial support of the Natural Sciences and Engi-neering Research Council of Canada is gratefully acknowledged. CHAPTER 1 Introduction The complete understanding of the interaction between waves and offshore structures is essential for the safe and economical design of such structures. Due to the expansion of the search for oil and natural gas into harsher environments, offshore struc-tures are now being designed and constructed for increasing depths of water. For example, a guyed jacket platform was recently built in 300 metres of water in the Gulf of Mexico. With structures operating in deep water the effects of dynamic wave-structure interaction can be important and may even govern the design. Structural dynamics is most important when the frequency range of the exciting forces, in this case the wave frequencies, are in the range of one or more of the natural frequencies of the structure. In such cases large, potentially damaging, resonant responses can occur. At the water "depths now being con-sidered for offshore structures there is a strong possibility of these resonance conditions occurring due to the matching of the wave frequencies and the natural frequencies of the offshore structures. Obviously the effects of dynamic wave-structure interaction must be one of the major considerations in the design process. The recent comprehensive book of Sarpkaya and Isaacson (1981) reviews all of the pertinent aspects of wave forces on offshore structures. They indicate a number of areas which are still in the research stage, one of which is the dynamic response of flexible offshore structures in random waves. The analysis of this problem would not be difficult if it could be assumed that a linear rela-tionship exists between the excitation and the response. Then the powerful linear theory l Introduction / 1.0 2 of spectral analysis could be successfully used for a complete solution. However, as is usu-ally the case, linear theory is only an approximation valid in a limited range, usually for the case of small excitation and response. In the design situation the conditions of large excitation and large response are often of paramount importance and the effects of various nonlinearities cannot be ignored. Some of the more important nonlinear aspects of wave forces on offshore structures are: 1. The wave force on a slender offshore stuctural element is usually obtained using a modified form of the well known Morison equation. The resulting wave force is nonlinear in both the structure's velocity and the wave velocity. 2. Wave particle kinematics are often computed using a nonlinear wave theory. It is especially important to use a nonlinear wave theory to find the wave kinematics for extreme waves in shallow water. 3. For bottom founded structures there is a nonlinear interaction between the structure and the foundation. Often the foundation can be modeled by nonlinear springs and dampers. 4. Under extreme conditions the structure might experience plastic behavior at some locations which implies not only nonlinear behavior but history dependent behavior as well. In this thesis attention is focused on the first nonlinearity mentioned above, the nonlinearity present in the equation used to compute the fluid force on offshore struc-tures and structural elements. Airy or linear wave theory will be used to compute the wave particle kinematics. No consideration is given to the effects of the soil structure interaction and it will be assumed that physical properties of the structure are known constants and a linear structural model is assumed to apply. For fixed slender offshore structures wave forces are usually computed using the Morison equation. In their landmark paper Morison, O'Brien, Johnson and Schaaf (1950) assumed that the wave force on a structure is given as the sum of an inertia force Introduction / 1.0 3 and a drag force. For the usual case of a cylindrical section of diameter D the force per unit length is where p is the fluid density, dU/dt and cPU/dt2 are the undisturbed velocity and acceler-ation of the flow respectively and Q and Cm are the drag and inertia coefficients. The inertia force is analogous to that on a body in a uniformly accelerating flow of an ideal fluid, while the drag force is analogous to that on a body in a steady flow of a real fluid with a fully developed wake formation behind the body. Measurements of the forces in these reference flow situations leads to the specification of the drag and inertia coefficients. However, due to the differences between the fluid flow in waves and the two reference flows mentioned above, the drag and inertia coefficients for wavy flow are not the same as those given in the reference flows. Numerous experiments, summarized by Sarpkaya and Isaacson (1981), have shown that the two coefficients depend mainly on the flow Reynolds number, Keulegan-Carpenter number and the relative roughness. The Morison equation in the form given above has been extensively used in the design of offshore stuctures. These structures have been built in relatively shal-low water, that is water depths less than approximately 30 metres, where the effects of wave-structure interaction are insignificant or are easily taken care of with a linearization procedure. In shallow water an offshore structure is quite stiff with natural frequencies usually much greater than the wave frequencies. In deeper water however, offshore struc-tures become more flexible with natural frequencies approaching the wave frequencies and consequently consideration must be given to the possibility of dynamic effects increasing the response. The phrase wave-structure interaction is used to describe the coupling that occurs when the large response of the structure influences the force on the structure which in turn determines the structural response. In order to include the effects of wave-structure interaction, a modified form of the Morison equation which includes the structural veloc-ity is often assumed to apply. One popular method is to assume that the drag force is dU dt (1.1) Introduction / 1.0 4 proportional to the square of the relative velocity between the fluid and the structure. This formulation, called the relative velocity formulation, has been used by Malhotra and Penzien (1970a,b) and Grecco and Hudspeth (1983) among others. Usually the nonlinear-ity is handled by a linearization technique resulting in a greatly simplified analysis. Even though linearization is used, acceptable solutions often result. The second method that can be used to modify the Morison equation to include the effects of wave-structure interaction is the independent flow fields formulation recently described by Laya et al. (1984). This method considers two independent flow situations, the first being a far field flow in which the fluid force is given by the original form of the Morison equation and arises from the wave force acting on a stationary structure. The second is a near field flow caused by the structure oscillating in an otherwise quiescent fluid. In this case the fluid force is again given by the Morison equation written in terms of the structure's velocity and acceleration. Laya et al. (1984) discuss the relative merits of each of these two formulations for wave-structure interaction and qualitatively suggest their respective range of applicability. Irrespective of which formulation is used to modify the Morison equation, the resulting equation of motion will be nonlinear in both the response velocity and the fluid velocity. The incident wave kinematics are assumed to be known and the type of solution method depends on the description of the wave field. Two wave field descriptions are often used for design purposes. The first method considers a deterministic design wave of permanent form. Using linear wave theory to obtain the wave kinematics, this design wave will have a harmonic form with a single fixed period and wavelength. There are no known analytic solutions to the equation of motion even for the case of simple harmonic forcing and recourse must be made to approximate or numerical methods. In this thesis the equations of motion with harmonic forcing are solved by three methods; numerical integration, linearization, and perturbation. As is usually the case the results from the numerical integration will be used to assess the accuracy of the other approximate methods. Introduction / 1.0 5 The second wave description represents the wave as a Gaussian random pro-cess. Due to the assumption of linear wave theory, the wave velocity and acceleration will also be Gaussian random processes. However the nonlinearity in the Morison equation yields a wave force and response which are non-Gaussian. Much research has been con-ducted in the past fifteen years to determine the probabilistic nature of the wave force on a rigid structure using the Morison equation, linear wave theory and Gaussian waves (see for example Borgman 1967). The corresponding probability of the response, which is ultimately the response statistic of interest, has received much less attention. It has been hypothesized that the response probability is of the same form as the probability of the forcing and this hypothesis is considered here. Many solution techniques focus on the mean square response statistics, specifically the response variance and the response spectrum. This is because these statistics are much easier to find than the probability density of the response. For example, for the case where wave-structure interaction is important, mean square solutions have been found using numerical simulation (Shinozuka and Wen, 1971), linearization (Gudmestad and Connor, 1983), and by perturbation (Eatock-Taylor and Ra-jagopalan, 1982). The similarity in methods of solution for harmonic forcing and random forcing is not surprising because a method that is successful in obtaining an approximate solution of a nonlinear differential equation with harmonic forcing can often be adapted to obtain a solution to the same equation with random forcing. A new perturbation method for the random forcing case is developed here using this approach of adapting a solution technique from the harmonic forcing case. The subject of randomly forced nonlinear differential equations has had a very active research period beginning in the early 1960's. Reviews on the subject have periodically appeared in the literature, for example; Crandall (1966), Caughey (1971), Iwan (1974), Smith (1978), Vanmarcke (1979), Spanos (1981), Roberts (1981a,b), Crandall and Zhu (1983), and most recently To (1984) and Roberts(1984). In addition, the books by Lin (1976), Soong (1973) and Nigam (1983) also treat the problem of random forcing Introduction / 1.0 6 of nonlinear differential equations. Due to the extensive nature of the reviews cited above, another comprehensive review is not warranted at this time. However, the topics that are of interest in this thesis will be briefly discussed in what follows. Basically there are two totally different approaches that can be taken in order to try and find solutions to randomly forced nonlinear differential equations. The first approach starts with the differential equation and then modifies a solution technique suitable for deterministic forcing. The most widely used of these methods are the method of equivalent linearization, the perturbation method, and numerical simulation. Each of these methods will be briefly described. By linearizing the nonlinear equation the very powerful spectral analysis methods of linear random vibrations can be used to obtain the solution. If the forcing is known to have a Gaussian probability distribution then due to the linearity of the system, the system response will also be Gaussian. Therefore the response spectral density completely characterizes the probability distribution and the problem is reduced to finding the spectral density, a simple task for a linear system. The drawback to any linearization method is that certain features of the nonlinear solution are not reproduced by the linear solution. Specifically, the probability distribution of the response of a nonlinear system will in general be non-Gaussian even though the forcing is Gaussian. The departure from a Gaussian distribution is most evident in the tails of the distribution, making the accurate prediction of extreme events of the nonlinear process very difficult. Also, the response spectral density of the nonlinear system may have multiple peaks whereas the spectral density of the corresponding linear system may only have a single peak. In spite of these shortcomings the method of equivalent linearization has proved to be a very powerful method to solve randomly forced nonlinear differential equations. Specifically, the mean square value predicted by the method of equivalent linearization has been found to be very accurate, even for cases of very large nonlinearity. The application of the method of equivalent linearization to the problem of wave forces on offshore structures has been investigated by Spanos and Chen (1981). Introduction / 1.0 7 In contrast, the perturbation method for randomly forced nonlinear equa-tions first described by Crandall (1963) only applies when the nonlinear term is very small compared to the other terms in the equation. Basically this method is the extension to the randomly forced case of the well known regular perturbation method. As with most perturbation methods, the nonlinear equation is transformed into a set of recursive linear equations which are then solved by a standard method such as spectral analysis. Practi-cally, the computations become very involved for all but the first or second terms in the perturbation series. Also, usually only the moments of the response can be found directly using this method, with only the lower order moments being easily found. One problem with the regular perturbation method described by Crandall is its relationship to the case of harmonic forcing. It is well known that a regular perturbation method often does not give a uniformly bounded solution for the case of harmonic forcing. Such is the case for the equations being considered here. The idea used to overcome this problem is to introduce additional degrees of freedom into the set of equations whose values are chosen to suppress the terms responsible for the unboundedness. These methods, generally called multiple scale perturbation methods, have been very successful in finding solutions to harmoni-cally forced nonlinear equations. One method, Lindstedt's method, is useful for finding the steady state solution to such problems. It seems appropriate then to extend Linstedt's method to the randomly forced case in order to find the steady state or stationary solution. A major contribution of this thesis is to describe how Linstedt's method can be applied to randomly forced nonlinear equations. It turns out that the first term in Linstedt's pertur-bation series is identical to the method of equivalent linearization while the higher order terms are improvements to the method of equivalent linearization. It also turns out that the perturbation method described here is valid over a much wider range of nonlinearity than is the regular perturbation method of Crandall. Lindstedt's method is applied to both harmonic and random forcing and the experience gained in applying the method for harmonic forcing gives guidance to the implementation of the method for random forcing. Introduction / 1.0 8 It should be noted that the perturbation method described here is applied to the govern-ing equation of motion without any simplification. Eatock-Taylor and Rajagopalan (1982) describe a regular perturbation method for wave forces on offshore structures that starts with a linearized form of the equation of motion. Comparisons between their method and the perturbation method developed here will be made to show where each method is applicable. The final solution technique applicable directly to the differential equation is numerical simulation and will be described here in the context of wave forces on offshore structures. This method, often called Monte Carlo simulation, starts by simulating a possible realization of a random sea surface. The criterion that is used to generate this realization is that the wave height spectral density of the realization approximates one of the commonly used wave height spectra often used in design. Borgman (1969) and Spanos (1983) review the techniques that can be used for this simulation process. The fluid velocity and acceleration are then calculated using linear wave theory and the equation of motion is numerically integrated using the simulated velocity and acceleration records in the Morison equation to obtain the time history of the response. Two scenarios about this simulation procedure can occur. The first is that only one long simulated record is required to completely represent the random character-istics of the sea state. The response statistics of interest are then found by taking averages over the resulting single response history. This implies that the sea state is assumed to be both stationary and ergodic. When the sea state is not ergodic a number of sea states must be simulated and each one applied to the equation of motion and integrated to yield an ensemble of response histories. The response statistics of interest are then computed using averages over the ensemble of response histories. In either case the computational effort is quite significant which serves as the motivation for finding efficient approximate methods. As was the case for harmonic forcing, the numerical simulation results are used to assess the accuracy of the other approximate methods. Introduction / 1.0 9 A totally different approach to the problem can be taken when the random forcing can be assumed to be equivalent to white noise. Then it can be shown that the response is a Markov vector process for which the response probability can be completely described by a transition probability density and an initial condition. An equation for the transition probability density called the Fokker-Planck-Kolmogorov equation, abbreviated here FPKE, can be found from the original equation (see for example Caughey, 1963). Once a solution of the F P K E has been obtained the response statistics of interest can be found in a straightforward fashion. The F P K E is a second order linear partial differential equation for the type of equations of motion considered here. Note that the problem has changed from trying to solve a nonlinear ordinary differential equation to trying to solve a linear partial differential equation, both being quite difficult. Unfortunately there are only a few known solutions to the FPKE in a few special cases. Even the stationary FPKE, whose solution yields the stationary probability density, has only a few known solutions and approximate methods are often required in order to find a solution. Nigam (1983) reviews some of the methods that have been proposed to find solutions to the FPKE. As mentioned above the F P K E approach can only be applied when the forc-ing can be considered to be white noise. White noise is a physically unrealistic model for natural processes because it implies infinite energy for the process. However, in the case of a wide banded process applied to a lightly damped system, the assumption of white noise as a model for the forcing results in very realistic solutions. The usefulness of white noise as a model for real processes in this case is similar to the usefulness of the concept of ideal point loads in the theory of elasticity. If a natural process is not wide banded, as is the case of a random sea state, then the F P K E approach can still be used if an enlarged system of equations is used. The white noise process is first passed through an auxiliary system which shapes the input to the real system. The disadvantage here is that the resulting FPKE is even more difficult to solve. The organization of this thesis follows the above discussion. Chapter 2 de-rives the equations of motion of a single degree of freedom model of an offshore structure Introduction / 1.0 10 which is then solved for the case of harmonic forcing. The equation of motion is solved by numerical integration, various linearization schemes, and a perturbation method. Results are computed for a range of typical parameter values and comparisons are made between the approximate methods and the numerical integration solution. The experience gained in solving the deterministic forcing case is used to guide the solution methods for the randomly forced situation, especially with respect to the perturbation method. Chapter 3 starts with the equations of motion derived in Chapter 2 and obtains solutions for the case of random forcing by various linearization methods, the perturbation method and numerical simulation. The solutions from each method are com-pared using the mean square value and the response spectrum as the main parameters for comparison. The probability density of the response is also discussed, specifically with respect to the hypothesis that the response probability density is of the same form as the probability density of the forcing. Finally, Chapter 4 concludes this thesis by summarizing the results obtained and by putting the newly developed perturbation method in its proper perspective. CHAPTER 2 Deterministic Wave Forcing In this chapter the problem of a single degree of freedom model of an offshore structure excited by a deterministic wave is considered. The flow velocity is assumed to be the sum of a constant current and a sinusoidal time variation. First, the equation of motion of the model will be formulated using the relative velocity and the independent flow fields formulations to modify the Morison equation to account for wave-structure interaction. In either case a nonlinear ordinary differential equation with harmonic forcing results. The equations for the relative velocity formulation are then solved by numerical integration, by linearization methods and by a perturbation method. The effects of the nonlinearity on the response is noted and an assessment of how well each approximate method captures the nonlinear behavior is discussed. The numerical integration solution is assumed to approximate the exact solution for all cases to be considered. Due to the similarity between the independent flow fields equation and the relative velocity equation, the solution methods applicable to the relative velocity formulation could also be applied to the independent flow fields equation. Figure 2.1 shows the single degree of freedom model along with the parame-ters and variables of interest. The equation of motion for this system on a per unit length basis is m ° ^ + C ^ + * * ' = / ; + / i ( 2 > 1 ) where m 0 is the mass per unit length of the cylinder, D is the cylinder diameter, c is the 11 Deterministic Wave Forcing / 2.0 12 F I G U R E 2.1 Single degree of freedom model. structural viscous damping per unit length, k is the structural stiffness per unit length, x* is the absolute cylinder displacement, U* is the wave fluid particle displacement, p is the density of water and V is time. The superscript * denotes a dimensional variable. Morison et al. (1950) have shown that the wave force on an offshore structure can be expressed as the sum of an inertia force, j\ and a drag force, which presently forms the basis for computing wave forces on slender offshore structures. The original form of the Morison equation applies to a fixed structure and for the typical circular section is, on a per length basis dU* ~dF where Cm and Cj are the inertia and drag coefficients respectively. The values of these two coefficients have been the object of much study and they have been found to de-pend mainly on the Reynolds number, the Keulegan-Carpenter number and the relative roughness. The Reynolds number is the ratio of the inertia force to the viscous force and is Re = (dU*/dt*)mDjv where v is the viscosity of water and (dU*fdt*)m is the maximum flow velocity. The Keulegan-Carpenter number K is proportional to the ratio of the maximum wave particle displacement Um, to the diameter of the cylinder and is given by K = (Um/D)2ir. A comprehensive summary of the dependence of C m and Cd on „ „ nD2 „ d2U* 1 n „ dU* ri+rD = P-Tcm-1^+-pDcd— (2.2) Deterministic Wave Forcing / 2.1 IS these parameters is given in Sarpkaya and Isaacson (1981). It is beyond the scope of the present study to consider the effect of variable force coefficients and, as is usually the case, representative values of Cm and Cj will be used. 2.1 EQUATIONS OF MOTION The Morison equation in its original form applies to a fixed structure. To account for the effect of wave-structure interaction, the Morison equation is modified to include the effect of the structure's velocity on the wave force. Presently there are two formulations that have been proposed to consider the structural motion, the relative ve-locity formulation and the independent flow fields formulation. The equations of motion are derived for both formulations but solutions are only obtained for the relative velocity formulation. As will be seen later the similarity of the equations for the two formulations means that solution techniques applicable to the relative velocity formulation can also be applied to the independent flow fields formulation. The wave force postulated by the relative velocity formulation is i * 1 Tin (dU* d x * \ \ d U * d x * + 2 p D C d \ l F - w ) \ l F - - d T (2.3) The inertia force, the first two terms on the right hand side (r.h.s.), is the sum of the inertia force on a rigid body in an accelerating inviscid fluid plus the inertia force due to an accelerating body in an otherwise quiescent inviscid fluid. The term pvD2(Cm — l)/4 is commonly called the added mass coefficient. The drag force is taken to be proportional to the square of the relative velocity between the fluid and the structure, hence the name relative velocity formulation. Note that the fluid force is now a nonlinear function in the response velocity, which is the source of all the computational difficulty. The use of the relative velocity formulation is widespread and usually appears in a linearized form which Deterministic Wave Forciny / 2.1 14 greatly simplifies the analysis. Typical examples of its use are given by Malhotra and Penzien (1970a,b) and Grecco and Hudspeth (1983). The independent flow fields formulation to account for wave-structure inter-action has recently been described by Laya et al. (1984) and they assume a fluid force of the form dU* »£2~ 1, the values of the inertia and drag coefficients tend to values independent of K. For a < 1 the coefficients tend to the harmonic values of Figures 2.2 and 2.3. Often, in the presence of a current, Re and K are computed based on a velocity which is the algebraic sum of the current and the maximum harmonic velocity. Figures 2.2 and 2.3 are then used to find the force coefficients. Details of other methods that have been proposed to account for the effects of a current on the force coefficients are given by Sarpkaya and Isaacson (1981). Deterministic Wave Forcing / 2.2 19 F I G U R E 2.3 Drag coefficient as a /unction of Re and K (Sarpkaya and Isaacson 1981). Deterministic Wave Forcing: / 2.2 20 2.2.3 Kenlegan-Carp enter number The following discussion is for the case of zero current. The Keulegan-Carpenter number K, is proportional to the ratio of the wave particle displacement to the diameter of the cylinder and is given by K — (Um/D)2ir. If K < 2, then a fluid particle will not separate from the cylinder during a one half cycle of the wave motion and no wake will form behind the cylinder. Typically, this condition denotes a "large body" where the effects of diffraction coincidentally become important. Using drag forces as defined above is not appropriate in this case because they are assumed to arise from the pressure drop associated with a well developed wake. For the "large body" problem the fluid is usually assumed to be inviscid and the forces on the body are found on the basis of the Laplace and Bernoulli equations and appropriate boundary conditions for the fluid motion. Details are again given by Sarpkaya and Isaacson (1981). For larger values of K the fluid particles separate from the cylinder eventually resulting in a wake structure of alternately shedding vortices during each half cycle of wave motion. It is for this case that the Morison equation applies. For intermediate values of K the flow patterns are very complex and further study is needed to correctly model the wave forces acting on the structure. The parameter UmQ/D which appears in the governing equations is related to the Keulegan-Carpenter number by the relationship K = (UmQ/D)(2ir/il). Laya et al. (1984) suggest that the relative velocity formulation applies to the situation of a well developed wake, quasi-steady flow and therefore large Keulegan-Carpenter number, K > 10 — 15. By choosing UmCl/D = 2.55 this condition is satisfied for the frequency range of interest here. This value was also used by Blevins(1977). Experimental studies are needed to determine exactly the range of K for which the relative velocity formulation applies and whether or not the independent flow fields formulation is more appropriate for smaller values of K. Determinbtic Wave Forcing / 2.2 21 2.2.4 Reynolds number The Reynolds number based on the maximum wave velocity for the case of zero current is where v is the viscosity of water and is about 1 0 - 6 m 2 / s . Therefore if the diameter is of the order of 1.0 m and the velocity is about 1 m/s, then the Reynolds number for a prototype offshore structure is about 10 6. To design small scale experiments with such large Reynolds numbers in oscillatory flow is very difficult, so that the appropriate values of the inertia and drag coefficients can only be obtained by extrapolation of the small scale data shown in Figures 2.2 and 2.3. 2.2.5 Mass ratio is seen to represent the importance of the nonlinearity in the governing equations. As some solution methods are only valid for small values of the nonlinearity, for example the perturbation method, the range of this parameter is seen to be very significant. Consider, for example, a hollow steel pipe of thickness t and external diameter D. The nonlinearity parameter e can be written for this case as where p„ is the density of steel. Taking t/D = 0.1, Cm = 1.25 and Cd — 1.45 results in a value of e = 0.27. Certainly other combinations of variables can be used to determine a specific value of e, but in what follows the range 0 < e < 0.5 will be considered. 2.2.6 Frequency ratio The final parameter of interest is the ratio of the forcing frequency to the natural frequency of the structure. As structures are being built in deeper waters their natural frequencies approach the predominant wave frequencies and potentially damaging The mass ratio pD2/m or equivalently the parameter e = (pD2/m)Cd/2 Deterministic Wave Forcing / 2.3 22 resonant responses can result. Therefore the dimensionless frequency range of interest is 0 < Q < 2. Of special interest are the frequencies Cl = 1/p, where p is an integer. The case corresponding to p = 1 is the fundamental resonance condition, which occurs when the natural frequency of the structure matches the wave frequency. The other resonant frequencies for p > 1 arise due to the nonlinear nature of the Morison equation. As will be shown later, for a = 0 only the odd values of p are important, with the first two being the most significant. 2.8 will be described. A Runge-Kutta numerical integration scheme is used to assess the accuracy of all other approximate methods. Several linearization methods can be used for the approximate solution of these equations and the relative merits of each will be discussed. Finally, a perturbation method similar to Lindstedt's method will be described. The solution of an inhomogeneous ordinary differential equation consists of two parts: the homogeneous solution which for a damped equation always decays with time, and a particular solution which represents the steady state response. Realistically only the steady state response is of interest and this will be the case for all methods of solution. 2.3.1 Numerical Integration In the absence of an analytical solution, numerical integration usually pro-vides the most accurate representation of the solution to the equation of motion. However, a major drawback of any numerical solution is that it only applies for the specific parame-ters used in the numerical procedure. Therefore the computational effort required to do a comprehensive parametric study is usually prohibitive. The solution procedure used here starts by rewriting equation 2.7 as a pair of first order equations by introducing two new variables defined by 2.3 SOLUTION METHODS In this section a number of techniques for the solution of equations 2.7 and Deterministic Wave Forcing / 2.3 23 resulting in dxl = -2$Nxx -x2 + iU" + e(U' - x') lU'-x1 dt di2 dt where U is given by equation 2.10. Solution of this pair of equations is efficiently achieved using a fourth order Runge-Kutta numerical integration scheme. Specifically, the algorithm given by White (1974) which performs the integration using Gill 's method, was used with a time step of 0.2. Since the steady state solution is of interest, zero initial conditions are prescribed and the integration carried out until an oscillation of constant amplitude is reached. 2.3.2 Linearization Methods As the solution of a linear equation is usually much simpler than that of a nonlinear one, often the first approach to solving a nonlinear equation is linearization. Although certain characteristics unique to the nonlinear solution cannot be reproduced by the corresponding linear solution, often valuable results can be obtained over a limited range of system parameters. It is therefore very important to determine the range of appli-cability of a linearization technique and this will be discussed for each method described. There are many different criteria that can be used to linearize a nonlinear equation. Those that have found widespread use for the problem of wave forces on offshore structures will be described here. For the relative velocity formulation either equation 2.7 or 2.8 can be linearized by various methods. Each of these methods will now be discussed in turn. The drag force term in equation 2.7 contains the nonlinearity to be linearized. Rewriting the drag force as Id = where Deterministic Wave Forcing / 2.3 24 immediately suggests the linearization possibilities. If the structure velocity is very much less than the wave velocity then both the quadratic term and the linear term in the ratio of these velocities can be neglected in comparison to 1. The result is a drag force which is independent of the structure's velocity, in which case equation 2.7 simplifies to x" + 2$Nx' + z = iU" + eU'\U'\ (2.11) The only effect of the motion of the structure using this linearization method is the added mass. Note that the wave force is given by an equation very similar to the original form of the Morison equation and that there is no hydrodynamic damping in this case. Because there is no hydrodynamic damping a very large response will result, especially near the fundamental resonance fi = 1. As has been pointed out by Blevins (1977), this linearization is not very practical but it is instructive to compute the solution to equation 2.11 for comparative purposes. The solution technique is described by Blevins (1977) and starts by expanding the right side of equation 2.11 in a Fourier series. The solution is then also written as a Fourier series with unknown coefficients which are found by equating like harmonic terms. A typical example of this procedure is given in Appendix A l . For completeness, the Fourier series expansion for U' \U'\ where U is given by equation 2.10, is also contained in Appendix A l . It should be noted that this Fourier series expansion contains only odd harmonic terms when a = 0. The even harmonic terms are due to the steady current for non-zero a. Additionally if or > 1 then only the first two Fourier terms remain. Due to the solution procedure, x[t) also contains the same Fourier composition as does the forcing for the various values of a. Gudmestad and Connor (1983) also present the Fourier series expansion for U'\U'\, which is only slightly different from that given in Appendix A l as they define the wave displacement slightly different from the definition in equation 2.10. As the above method is only the first of the relative velocity linearization schemes, the notation R V L l will be used to describe it. Alternatively, this linearization is described as a zero drag damping linearization. Deterministic Wave Forcing / 2.3 25 The second linearization scheme neglects only the quadratic term in the drag force. Using this linearization, equation 2.7 simplifies to The only difference between this linearization and the previous one is the hydrodynamic damping in equation 2.12. This represents a more realistic estimate of the damping and the solution should be accurate subject to the condition that the structure's velocity is small compared to the wave velocity. The effective damping is now time dependent and equal to The solution of linear ordinary differential equations with time dependent coefficients can be found by a number of methods. A perturbation method described by Eatock-Taylor and Rajagopalan (1982) for this equation will be discussed in the next chapter on random forcing. For harmonic forcing the solution is obtained using the Runge-Kutta numerical integration procedure described above. The time dependent drag damping case will be noted as RVL2 in what follows. equations is to use an average value of the time dependent coefficient, reducing the equation to a constant coefficient one. The average value of \U'\ is found by integrating over one period of the wave, or where the overbar represents the average value over one period. Substituting for U from equation 2.10 and integrating gives (2.12) An approximate method that is often used to solve time dependent coefficient a < 1 a > 1 where sin/? = a; 0 < / ? < TT /2. When a = 0 the familiar result \U'\ a-0 — JT D Deterministic Wave Forcing* / 2.3 26 is obtained. Therefore for zero current and constant drag damping the effective damping is The constant drag damping scheme will be denoted RVL3 . Gudmestad and Connor (1983) also present the linearization of the Morison equation using a constant drag damping. However, they propose the linearization of the drag force as [U1 - x')\U' - x'\ = U'\U'- x'\ -x'\U'- x'\ « U'\U'\ - \U'\ x' which underestimates the hydrodynamic damping in equation 2.12 by a factor of two. This linearization does not appear to have a rational basis whereas the linearization which leads to equation 2.12 is mathematically consistent and is believed to be the correct one. This concludes the linearization schemes that can be applied to equation 2.7. The final linearization procedure is a linearization of equation 2.8 for the relative displacement r. The nonlinear drag force term is linearized in the general case as fD =er'\r'\ = a + 6r' where a and 6 are linearization constants. These constants are chosen such that the result-ing linear equation, in some sense, comes as close as possible to the nonlinear equation. One way to do this is to minimize the mean square error of the linearization. This method is commonly called the method of equivalent linearization. The mean square error e, is e= {er'\r'\-{a + br')}2 where the overbar again denotes the average over a period. The linearization constants a and 6 are obtained by minimizing e with respect to a and b or Deterministic Wave Forcing / 2.3 27 which results in a set of equations whose solution gives a and b. The two equations for a and 6 are _ a + br1 = er1 {^l _ _ (2-13) ar1 + br12 = er12 [r1] In order to evaluate the averages, the expected functional form for the relative velocity over one period is required. This can be found by considering the equation resulting from the linearization procedure which is r" + 2(Cjv + b/2)r' + r = (1 - i)U" + 2)} where the phase $ is taken to be zero for convenience. Note that 8 = 0 corresponds to a = 0, the case of zero current. Two cases must be considered, 5 < 1, and 6 > 1. If 6 < 1 then the integrations over one period yield a = — [(1 - 252)/3 + {\b - ASZ) cos/? + (25 2 - \ ) sin 2/? ir i 3 2 2 + -6 cos/? cos 2/3 e2B 7T 5 2<5/? + (282 + - ) cos /? - 6 sin/? o — ^ cos /? cos 2/?j where sin/? = 8; 0 < / ? < TT /2. If 8 > 1 then a = eB2(--82) v 2 y 6 = eB28 Note that in either case the values of the linearization coefficients depend on the response velocity through the parameters B and 8. Therefore an iterative procedure is usually required in order to obtain the correct values of the linearization constants. This can be done by initially choosing a and 6 to be zero and then finding the response velocity. New Deterministic Wave Forcing / 2.3 28 values of a and b are found from the above expressions and the iterations continued until the values of o and b converge to their correct values. If a = 0 then 5 = 0 and the linearization constants simplify to o = 0 _ eBB (2-14) 6 = = 1 7 Thus for a = 0 the effective damping is the familiar equation AB f e " = co a n 0 ^ /o given by Co = % ( l - n 2 ( l - 7 ) ) <*o = - 2 f j v % ^ Deterministic Wave Forcing / 2.3 29 The perturbation method proceeds by expanding both the response r(t) and the damping ratio £/v m a perturbation series in powers of e r = r 0 + erx + e 2r 2 H (2.16) $N = fo + s?l + * fc + • • • The nonlinear damping is expanded in a Taylor series about e = 0 r'|r'| = r4K| + e2ri|rS|+---and substituting equations 2.16 and the above expansion into equation 2.15 yield a set of recursive linear ordinary differential equations: *o + 2 f t r o + r o = co c o s ft* + do s m fit + eo* + 9o = /o(0 (2.17o) r'/ + 2$0r[ + r x = -r{, |^ 0| - 2 ^ = f^t) (2.176) r'i + 2?0^ 2 + r 2 = -2ri |r{,| - 2ftri - 2fcr& = /2(<) (2.17c) These equations are solved in turn, the forcing for the second one depending on the solution to the first. One important feature of this set of equations is the realization that ft. represents the effective damping for the perturbation solution. Initially, ft is unknown and must be found as part of the solution. Note that because the nonlinearity results in increased damping, the expansion of fty in a perturbation series is appropriate in this case due to the way the effective damping appears in the resulting perturbation equations. As is usually the case for a multiple scale perturbation method, the parameters ft and ft are determined by suppressing the so called secular type terms on the right hand sides of equations 2.176 and 2.17c respectively. Therefore in order to obtain the first estimate of effective damping both equations 2.17a and 2.176 must be considered. Equation 2.17a can be solved by the Fourier series method described in Appendix A l and the response is given by ro(0 = ao c o s fit + fy)sm fit + AQt + BQ Deterministic Wave Forcing / 2.3 30 In order to solve equation 2.176 the nonlinear term on the right hand side is expanded in a Fourier series. That is —r'o \ r'Q\ = c 1 0 + ^ ^ C i n cos nQt + dln sin nQt where the Fourier coefficients C j „ and dln are found as usual. The forcing function for equation 2.176 can now be written as /i(0 = c io + (cn - 2^6 0Q) cos Qt + (dn + 2fta 0O) sin Qt + c l n cos nQt + dln sin nQt It is easy to realize that cos Qt and sin Qt represent the nearly secular terms here. If we considered the regular perturbation expansion of a similar undamped equation then as is well known (see for example Kevorkian and Cole 1981) the presence of these terms in fi(t) gives rise to a solution which is unbounded for large time at the fundamental frequency 0 = 1. This is a result of sinf and cost being homogeneous solutions of the undamped equation as well as being forcing terms. Obviously then a regular perturbation solution does not work and the idea used to overcome this problem is to introduce extra degrees of freedom which are used to suppress the terms (i.e. secular terms) responsible for the unboundedness. A perturbation method based on this idea is generally called a multiple scale perturbation method. With the presence of damping, the solution is not unbounded but the terms cos Of and sin Qt represent nearly secular or small divisor terms. Therefore they must again be suppressed using the parameters ft and ft. This is equivalent to ensuring that only r0(t) represents the response at the frequency O. This interpretation is useful in discussing the perturbation method for the random forcing situation. Suppressing the secular terms results in ( c n - 2ft60O) = ( d u + 2fta 0O) = 0 which are two equations for the parameter ft. It can be shown that both equations give the same value for ft ?i = ~ ( « o + 6o)1/a Deterministic Wave Forcing / 2.3 31 The first estimate of the effective damping is now found from the root of the equation tor = ?o + *fc(fo) where the explicit dependence of ft on ft is due to the dependence of a 0 and 6 0 on ft through the solution rQ(t). Equation 2.176 can now be solved using the Fourier series method with r i (0 = °io + X * ° l n c o s s m n ^ n=2 and the coefficients found as usual. Note that the summation starts from n = 2 due to the suppression of the fundamental frequency terms. If there is zero current then all even terms in the series are zero and there will be no static offset. Finally to this order of approximation, the absolute response is x(t) = u-r0{t) -er^t) where r 0(f) and ri(t) are given above. The solution to the next order of approximation proceeds in a similar fashion. First, the nonlinear term on the right hand side of equation 2.17c is expanded in a Fourier series with f2{t) written as / 2 W = c2o + (C21 ~ 2fefy)fl) c o s M + ( d 2 l + 2 & a 0 n ) sin tit 5^( C 2n — 2ft&innil) COS Tltlt n-2 + [d2n + 2ftOi„nft) sin ntlt where c 2 n and d2n are the Fourier coefficients of the nonlinear term in f2(t) and the other symbols have been defined previously. Again f2 1 3 *° D e found by suppressing the fundamental frequency terms in f2(t) resulting in two equations for the single unknown £ 2 . Unfortunately in this case the two equations for £2 do not reduce to one single equation. Therefore it is not possible to suppress completely the terms cosfM and sinflt in f2(t) but we can find the $2 which minimizes these secular type terms. This is achieved by Deterministic Wave Forcing / 2.4 32 minimizing the amplitude of the fundamental frequency terms by finding the ft which satisfies j - [(c 2i - 2ft& 0 n) 2 + ( 0.6 and s < 0.15. It overestimates the solution near the secondary resonance peak at ft = 1/3 but as will be shown later an improved solution results if the perturbation series is extended to 0(e 2 ) . For frequencies less than about 0.6, the perturbation solution is an improvement to the solution from the method of equivalent linearization. Characteristic of an asymptotic series, the perturbation solution to O(e) is a very good approximation of the numerical solution. Figures 2.8 through 2.11 show the improvement in the perturbation method by including terms to 0(e 2 ) . Comparison is again made with the numerical integration so-lution with improvement in the perturbation solution evident for all values of nonlinearity. The solution only changes near the secondary resonance peak where the largest departure from the numerical solution occurs. Even for the largest value of nonlinearity, e = 0.50, the perturbation solution gives a good estimate of the solution over the entire frequency range. As is usually the case, the effort required to obtain the perturbation solution to 0(s2) is significantly greater than the effort required for the 0{s) solution. Figure 2.12 shows the values of effective damping predicted by the method of equivalent linearization and by the perturbation method. Both the first estimate of effective Deterministic Wave Forcing / 2.4 39 0.0 0.2 0.4 0.8 0.8 1.0 12 1.4 1.6 1.8 DIMENSIONLESS FREQUENCY F I G U R E 2.8 Frequency response for perturbation method with s = 0.05. F I G U R E 2.0 0.0 02 0.4 0.6 0.8 1.0 12 1.4 1.8 1.8 2.0 DIMENSIONLESS FREQUENCY Frequency response for perturbation method with e -0.15. Deterministic Wave Forcing / 2.4 40 o 00 W CO o Pu CO w C J CO in o u S - l 5 Q o d EXACT PERT1 PERT2 c=0.25 - i — 0.0 02 0.4 0.6 0.8 1.0 12 1.4 1.8 1.8 DIMENSIONLESS FREQUENCY" Frequency response for perturbation method with e = 0.25. - I— 1— 1.4 2.0 F I G U R E 2.10 o m. CO CM o CO w o E-1 i s CO o o m d O i—i s s 8 © EXACT PERT1 PERT2 \ _ c=0.50 F I G U R E 2.11 0.0 02 0.4 0.8 0.8 1.0 12 1.4 1.8 1.8 2.0 DIMENSIONLESS FREQUENCY Frequency response for perturbation method with e = 0.50. Deterministic Wave Forcing / 2.4 41 o • RVL3 RVL4 P E R U PERT 2 0.0 02 0.4 0.6 0.8 1.0 12 1.4 1.6 1.8 £ 0 DIMENSIONLESS FREQUENCY F I G U R E 2.12 Effective damping for deterministic forcing. damping using terms of O(e) indicated by P E R T l and the second estimate of effective damping using terms of 0(e 2 ) indicated by P E R T 2 are shown in the figure. Identical values of damping are found by the method of equivalent linearization and the perturbation method to O(e). This equivalence in effective damping allows us to put the perturbation method in its proper perspective; that is, it represents the consistent correction to the method of equivalent linearization. Inclusion of terms to 0(e 2 ) only change the estimate of the effective damping hear the second resonance peak 0 = 1/3 where the solution has the largest departure from the numerical result. The decrease in effective damping near 0 = 1.0 can be explained by the way it is computed by the method of equivalent linearization. As can be seen from equation 2.16 the effective damping is proportional to the relative velocity between the structure and the the fluid. The relative velocity is the greatest when the structure is at rest and it is least for the largest structural motion, near the fundamental resonance O = 1.0. Therefore Deterministic Wave Forcing / 2.4 42 the effective damping is largest near ft = 0.0 and smallest near ft = 1.0. Also shown in Figure 2.12 are the effective damping values from the average drag damping method. For convenience the relevant equation for the effective damping by this method is repeated here 2 UmCl This effective damping is a constant, independent of the frequency and is only shown at zero frequency. The average damping method predicts an effective damping that is greater than the method of equivalent linearization and the perturbation method. However in comparing the magnitudes of effective damping for the various methods it must be remembered that they apply to different differential equations and only comparison of the response predicted by each method gives the complete picture. To obtain a feel for the time history of the response Figure 2.13 presents a few cycles of the response for ft = 1.0 and ft = 1/3 for a nonlinearity of e = 0.15 as computed by both the perturbation method to O(e) and by numerical integration. The dashed line with symbols represents the perturbation solution while the solid line is the numerical integration solution for each frequency. The differences in these responses are obvious with the response at ft = 1/3 having two distinct frequency components as might be expected considering the way in which the solution has been found. This figure also shows that the perturbation method accurately predicts the steady state time history of the response. The next set of results is concerned with the effect of a constant current on the response. The response curves plotted in Figures 2.14 through 2.17 represent the solution for current to harmonic velocity ratio a = 0.2,0.5,0.8 and 1.2 respectively. The nonlinearity chosen here is e = 0.05 with the other parameters the same as in the previous case. The solutions are given only up to ft = 1.1 as they show no differences with the case of zero current beyond this value. Quite obvious now is the appearance of another response peak at ft = 1/2. Again this peak can be predicted from the Fourier series expansion of Deterministic Wave Forcing / 2.4 43 o F I G U R E 2.13 Time history of the response for frequencies fi = 1/3 arid Q = 1. the nonlinear term for the case of nonzero current. In addition, a static offset appears in the solution and the responses shown in these figures include the offset. As can be seen, increased current increases the hydrodynamic damping with a resulting decrease in response. In order to obtain the maximum absolute dynamic response, the static responses need to be known and are given in Table 2.3. The values of static offset predicted by the perturbation method will be shown later. Four methods have been used to compute the response with a constant current; numerical integration, perturbation to 0(e), the method of equivalent linearization and the linearization with the average drag damping. The average drag damping method is essentially that given by Gudmestad and Connor (1983) with the proper form for the effective damping used here as previously discussed. The influence of the contribution to the solution of the Q = 1/2 term increases as the current increases, effectively dominating Deterministic Wave Forcing / 2.4 44 0.0 0J2 0.4 0.6 0.8 1.0 DIMENSIONLESS FREQUENCY F I G U R E 2.14 Frequency response for s = 0.05 and a = 0.2. 0.0 02 0.4 0.6 0.8 1.0 12 DIMENSIONLESS FREQUENCY F I G U R E 2.15 Frequency response for e = 0.05 and a = 0.5. Deterministic Wave Forcing / 2.4 45 0.0 0J2 0.4 0.6 0.8 1.0 DIMENSIONLESS FREQUENCY F I G U R E 2.16 Frequency response for e = 0.05 and ct = 0.8. EXACT PERT1 RVL3" RVL4 a=1.2 T — 0J2 -I— 0.6 -I— 0.0 0.4 . 0.8 1.0 DIMENSIONLESS FREQUENCY F I G U R E 2.17 Frequency response for e = 0.05 and a = 1.2. 12 Deterministic Wave Forcing / 2.4 46 Table 2.3 Static response for various current ratios and e = 0.05 a Static response 0.0 0.327 0.2 0.470 0.5 0.734 0.8 1.056 1.2 1.577 the 0 = 1/3 term for a = 0.8. When a > 1 it is easily shown that only the 0 = 1 / 2 and 0 = 1 terms remain. The perturbation method again shows good approximation to the numeri-cal solution for all values of a and again the method of equivalent linearization fails to reproduce the secondary peaks at O = 1/3 and 1/2 as would be expected. However, near O = 1.0 the method of equivalent linearization again is very accurate. The average drag damping method shows erratic behavior for increasing current, underestimating the response at O = 1.0 for low values of a and overestimating it for large values of a. The values of effective damping and static offset computed by the pertur-bation method and by the method of eqivalent linearization are shown in Figures 2.18 and 2.19 respectively. Again the effective damping computed by both methods are iden-tical. Increasing the steady current increases the effective damping as might be expected. When the current ratio a > 1 the effective damping is a constant and does not vary with frequency. Also shown in Figure 2.18 are the effective damping values predicted by the average drag damping method. The average drag damping estimates of the effective damp-ing are greater than those from the perturbation method. As is shown in Figure 2.18, the differences between these values decreases as the current ratio increases. The static offset also increases with increasing current ratio as expected. As shown in Figure 2.19, the offset predicted by the method of equivalent linearization and Deterministic Wave Forcing / 2.4 47 8 d S 3 -OH Q W > i—i O w w o_ 0.0 0J2 0.4 0.6 0.8 DIMENSIONLESS FREQUENCY F I G U R E 2.18 o Effective damping for various current ratios. 0.0 02 0.4 0.6 0.8 1.0 DIMENSIONLESS FREQUENCY F I G U R E 2.19 Static offset for various current ratios. Deterministic Wave Forcing- / 2.4 48 the perturbation method are almost identical. The static offset predicted by the average drag damping method is also shown in Figure 2.19. These are the same as the values from the perturbation method for small frequency. CHAPTER 3 Random Wave Forcing In this chapter attention is focused on the solution of equations 2.7 and 2.8 when the forcing is described as a random process. In general the wave velocity will be assumed to be a Gaussian random process with a mean value of V. That is, the wave velocity can be written as U' = V + v! = (3-4) The evaluation of the expected value in this equation is accomplished by using the expan-sion of \U'\ in an Hermite polynomial series. That is \U'\= | aerf ( a /v^ ) + \[^e~a*/2 j (3-6) 6X = 2cverf(o7\/2) If a = 0, the case of zero current, these expressions reduce to those given by Eatock-Taylor and Rajagopalan (1982). Random Wave Forcing / 3.1 53 The expected value (\U'\) is shown in Appendix A4 to be simply equal to the coefficient 60 of the Hermite polynomial series expansion. That is (|C/'|) = cr1 l (|aerf(a/V2) + y | e - V 2 | which is the same expression given by Gudmestad and Connor (1983) but obtained by a weighted residual technique. Substitution of this expression into equation 3.4 gives the effective damping as teff = to + « V | a e r f ( a / v ^ ) + yj^e""2/2 J (3.7) which reduces to for the case of zero current, a = 0. Similarly an expansion of the drag term U' \ U'\ can be obtained as U' \U'\ = = fc\, {z'(t)z'(t + r)> + w e [(z>{t) £ akHek(z(t + r))} + akHek(z(t))z'(t + r)}] Jt=o it=o + * 2 ( £ «i*«y(*(0) E*kBek{z{t + r))) y=o it=o Random Wave Forcing / 3.1 56 As is shown in Appendix A4, the sum of the second and third expected values in the above equation are identically zero while the other two expected values are easily computed using the orthogonality property of the Hermite polynomials. The resulting correlation function is r=0 Using the following relationships for a stationary random process and taking the Fourier transform gives the desired spectrum S / o ( 0 ) = £ V 0 5 ( O ) + tftf + £^}SAtl) + * 2 £ ^ { 5 . , ( 0 ) } " 7^2 an' The superscript * r denotes a convolution spectrum of order r and 5(0) is the Dirac delta function which has the useful property that for a function g(0) < K O ) 5 ( o - o 0 ) d o = <7(o0) A computationally efficient algorithm to compute convolution spectra based on the fast Fourier transform is given in Appendix A5 . Also shown in Appendix A5 are the convolution spectra of order 2 and 3 applied to a simple harmonic wave. The sum on r in the equation for equals 6 0. The correlation function is = (/i(0/i(' + *)> = ^(^bjHei{z{t))x'S)Y.hHek{z[t + r))x'0(f + r)) j=i k=i Random Wave Forcing / 3.1 57 Again the cumulant discard hypothesis is used to obtain an approximate form for this expected value as * 4 e 2 ( E M W ) ) 5>ffe t (*(* + r))) j=l k-l = 4e2J2b2rrl{Rz(r)YRxlo{r) r = l where the expected value of the Hermite polynomial is found using the relationships in Appendix A4. Now 5^(0) is found by taking the Fourier transform of RJ^T). Finally, using the above results, the response spectrum is given as - 2 „ 2 5,(0) = * 2 a 2 5(0) + | t f , ( 0 ) | 2 { ( 7 2 0 2 + ^ ± ) 0 2 5 „ ( 0 ) + ^ E g [ M 0 ) ] " ( 3 . 1 2 ) +<«'*[£ § ? ( * « ) ' * • . « ] } r = l «' where 7[ ] represents the Fourier transform. Note here that 6(0.) = 0 if O ^ 0 and *, (0) = 1. If the last term in equation 3.12 is neglected, then only the first term in the perturbation series is retained. That is, if the time dependent drag term is neglected, then the constant drag damping response spectrum is given by 5,(0) - e2a2S(U) + | t f , ( 0 ) | 2 { ( 7 2 0 2 + ^ ) 0 2 5 „ ( 0 ) 2 i (3-13) ••"E^fc-H } r=2 «' which corresponds to the case considered by Gudmestad and Connor (1983). As previously stated usually only terms up to and including r = 3 are evaluated in the sum in equation 3.13. Therefore the notation cubic constant drag damping will be used to describe this method. An even simpler estimate of the response spectrum results by neglecting the convolution spectrum terms in equation 3.13. Effectively this means that the drag term U'\U'\ has been linearized resulting in the response spectrum to be given by 5,(0) = s 2 a 2 5(0) + | t f , ( 0 ) | 2 { 7 2 0 2 + ^ ± } o 2 5 „ ( 0 ) (3.14) Random Wave Forcing / 3.2 53 The notation linear constant drag damping is used to describe this method. In equations 3.12 through 3.14 the known input wave height spectrum is simply related to the displace-ment and velocity spectra by equation 3.2. The response mean square value predicted by any of these methods can then be found by integrating the spectrum. That is as described in Appendix A3 . 3.2 EQUIVALENT LINEARIZATION The method of equivalent linearization has been extensively used to obtain mean square solutions to non-linear random vibration problems. A review of the method and its many applications has been given by Spanos (1981). In the context of random wave forces on offshore structures Spanos and Chen (1981) have applied the method to equation 2.8 which is now briefly described. Starting with equation 2.7 and assuming that the velocity is as indicated by equation 3.1 and introducing a relative displacement defined by f = u — x the equation of motion becomes f" + 2$Nr' + e{V + f')\v + ?\ + r=(l- i)u" + 2cNu' + u (3.15) analogous to equation 2.8. Due to the asymmetric non-linearity the solution to this equa-tion is assumed to be of the form f = F + r where r is the constant offset and r is a zero mean random process. Substitution into the equation of motion results in r" + 2$Nr' + s{V + r')\V + r' | + r + r = ( l - i)u" + 2foru' + u (3.16) Random Wave Forcing- / 3.2 59 To linearize this equation, the nonlinear term on the l.h.s. is replaced by C V where C is the linearization constant to be determined by mimimizing the mean square error of the linearization as was done for deterministic forcing. Performing the minimization gives the linearization constant as Cmt{= 9(r'W)dr' J-co where p(r') in this case is the zero mean Gaussian probability density function. Therefore the value of C can be shown to be given by C = 2ecxr, j ^ e r f ^ / v ^ ) + \ / f « " ^ / 2 | where aT< is the standard deviation of r1 and f3 = V/a^ is the ratio of the constant velocity to the standard deviation of the relative velocity. Not surprisingly a similar expression was found for b0 in equation 3.6 using the Hermite polynomial expansion described previously. The only difference here is that 0 ^Hcl{w{t))+2<0Hel{w(t)) + j Hei(w(t))dt^ (3.25) which indicates that f0(t) is simply related to the first Hermite polynomial of to. This fact is important in evaluating expected values involving fo(t). The spectrum of the relative response 5 r (fl) is given by 5r(n) = \Hr(n)\2{sfo(n) +«(s f o f l {Q) + sflfo{a)) + ^(Sfoh(fi) + Shfo{Q) + Sfl(Ci))} analogous to equation 3.11. The transfer function in this case is , i 7 r ( n ) l 2 = (1 - 0 2 ) 2 + ( 2 r o f t ) 2 which is known once the effective damping has been determined. As was previously the case, the effective damping is found by suppressing the secular type terms on the r.h.s. of equations 3.246 and 3.24c. That is, by using ft and & to suppress the secular type terms, the effective damping is then found from the second of equations 2.16. To find ft, the nonlinear term in f^t) is expanded in an Hermite polynomial series in w. Thus /i(0 = ~°\{P + w)\0 + w\- 2err» fcto - *-= -^akHek(w(t)) - 2\) £ 0 Therefore f\(t) simplifies to hit) = " ( " I - 2(vof l)fTe1(ti;(0) " 5>*ffe*M0) (3-27) Now since f0(t) is simply related to Hei(w{t)) by equation 3.25 and since equation 3.24a is linear, r0(t) is linearly related to 5'e1(iu(f)). Analogous to the deterministic case, all of the contribution to the solution due to fTe1(tw(i)) is given entirely by r 0 ( i ) . Therefore the secular type term in equation 3.246 is the term containing jffe1(io(<)) and ft is determined by suppressing this term. That is, ft is using equation 3.9 for aj. That this is the proper procedure in this case is seen when the effective damping is found using the second of equation 2.16 with terms to O(e). That is Co = fjv _ £$i = & + ear,Q | /?erf(/?/V^) + \ / f « " ^ / 2 | which is seen to be identical with the effective damping found using the method of equiva-lent linearization, equation 3.19. The perturbation method and the equivalent linearization method are closely tied together through the effective damping as was also observed for the deterministic forcing situation. Now with an estimate of the effective damping, the terms in the relative response spectrum can be evaluated. The first term Sf0(Vl) is simply related to the wave displacement spectrum due to the linear equation relating / 0 and u. sn(n) = \Hn(n)\2sf0(Q) Random Wave Forcing / 3.3 65 where |#„(n) | is given by equation 3.20. Again equation 3.2 relates the wave displacement spectrum to the known input wave height spectrum. Now consider the 0(e) terms in equation 3.26, 5 ^ ( 0 ) + 5y 1y 0(fi). Because fo(t) is only a function of Hei(w(t)) and fi(t) only contains Hermite polynomials of order 2 and greater, the correlation functions R^f^) and R^f^r) are identically zero due to the orthogonality of the Hermite polynomials (see Appendix A4). Therefore there will be no contribution to the relative response spectrum due to the sum •S'/0/1(O) + 5^y 0 (f l ) . A similar argument can be made to show that the terms Sj0f2(Q) and Sf2f0{tt) also do not contribute to the relative response spectrum. The parameter ft ap-pearing in f2(t) is determined by suppressing the secular type term containing Hei(w(t)). Therefore after suppressing ffe1(to(f)), f2(t) will only contain Hermite polynomials of order 2 and greater and the correlation functions Rf0f2{r) and #/ 2 / 0 (T ) will be identically zero due to the orthogonality of the Hermite polynomials. The value of ft is not determined explicitly here and the effective damping is based on the determination of ft as described above. The final term Sj^Cl) can now be easily evaluated using the Hermite poly-nomial expansion of equation 3.27 for fi(t). The correlation function Rfx{r) is y=2 k=2 r=2 r> where the expected value is found using the expressions in Appendix A4. Taking the Fourier transform gives the desired spectrum »/.<«) = E$K(")}*r r=2 r'0 and again the superscript * r denotes a convolution spectrum of order r. The details of the algorithm used to compute the convolution spectra is given in Appendix A5. It is possible to relate 5^(0) to the wave height spectrum due to the linear relationship between these two terms. Random Wave Forcing / 3.4 66 Actually the absolute response x is of more interest and the absolute response spectrum is analogous to the method of equivalent linearization. The only new spectral terms that need to be found here are Snr(Q) + Sru(Cl) which are found by the same procedure followed by the method of equivalent linearization for similar spectral terms. Using equations 3.19, 3.20 and 3.23 and with some simplification finally results in the absolute response spectrum being given as This expression is almost the same as the expression in equation 3.13 for 5 z(ft) as predicted by the cubic constant drag damping method except that here o~rtQ replaces k) k=l Random Wave Forcings / 3.4 68 The random phase angles oo by the central limit theorem. Borgman (1969), for example, gives more details of this method. A similar method which considers random amplitudes is discussed by Tuah and Hudspeth (1982). One disadvantage of the summing of sinusoids to simulate a wave record is that a periodic record results. This is not a realistic representation of a random sea state. Obviously, the longer the period of the record the better the simulation approaches the stationary ideal. This implies that the frequency increment should be as small as possible thus requiring a large number of harmonic components in order to accurately represent the target spectrum. If a small number of harmonic components are used in the simulation then the statistics of interest must be computed based on ensemble averages. For example Spanos and Chen (1981) use 20 harmonic components and 300 simulated wave velocity records to obtain mean square response values in order to assess the accuracy of the method of equivalent linearization. If extreme statistics were of interest then many more records would be required. Another method that can be used to obtain a wave record compatible with a target spectrum is to use the output of a filter with white noise input. For a linear filter the input-output relationship in the frequency domain is S0(Q) = IHfitltfSiitt) Random Wave Forcing / 3.4 69 where 5,(0) is the input spectrum, S0(ft) is the output spectrum and |/fy(ft)| 2 is the filter transfer function. The spectral density of the input white noise is constant and setting Si(Q) = 1 results in the output spectral density being equal to the transfer function of the linear filter. Therefore the task is to design the filter such that it's transfer function closely approximates the desired target spectrum. Spanos (1983) describes this design process for three different types of filters; autoregressive (AR) , moving average (MA) and autoregressive-moving average ( A R M A ) filters. For example, the time series that results from an A R M A filter is m p where and are the filter coefficients, wn are random numbers which represent the white noise process and At is the sampling time step. The random numbers are obtained from a standard random number generation subroutine and then scaled to give the spectrum of the white noise equal to one over the range 0 < ft < TT /AJ and zero for higher frequencies. To obtain a Gaussian white noise representation, the random numbers are sampled from a Gaussian distribution. The details of this procedure will be discussed later. Similar relationships can be written for the A R and M A filters which are special cases of the more general A R M A filter. Specifically A R : = 0 f > 1 and d^O M A : C£ ^ 0 = 0 A distinct advantage of simulation by filtering is that very long simulated records can be produced before repeating itself. This is a result of random number generators being able to produce on the order of 2 3 0 random numbers in sequence. Therefore it is possible to assume that the simulated record is both stationary and ergodic and that response statistics of interest can be found by averaging over one long response record. For example, Roberts (1977) integrates a nonlinear single degree of freedom system subjected to white noise forcing for 106 time steps and then computes response statistics such as the response Random Wave Forcing / 3.4 70 probability density, in order to assess the accuracy of an approximate method. Another advantage of digital filtering is that the spectral density of the filter is known explicitly in analytical form so that it can be used directly in an analytical procedure if desired. The details of determining the coefficients of an M A filter are given below. This method is relatively straightforward and easy to implement. The simulated time series from the M A filter of order q is Un = U(nAt) = ]jT c„wn_v u=-q with the transfer function of the filter given in terms of the filter coefficients as 15/(0) | 2 = (c 0 + 2 ^ c l / c o s f i i / A * ) 2 i/=i where At = n/Qc is the sampling time step and Cic is the upper frequency cutoff beyond which the spectrum is assumed to be zero. The filter coefficients c„ are found by minimizing the mean square error of the approximation of the filter's transfer function to the target spectrum. This results in the coefficients being given by C " = f t - " / 0 VSitt)cos{—)dtt i / = 0 , l , . . - f g where S(Q) is the target spectrum. Obviously c„ = c_v and the filter is completely described by q+1 coefficients. Spanos (1983) has used an M A filter of order 29 to accurately represent the P - M wave height spectrum. A n example of how well this filter approaches the P - M spectrum is shown in Appendix A2. Once the wave records have been simulated then the equation of motion can be numerically integrated to produce a sample response history. If the input simulation is long enough to be assumed ergodic, then the response statistics are found by averaging over the one response history. For a non-ergodic record, a number of wave records must be simulated and the statistics of interest found by averaging over the ensemble of response records. The method described in Clough and Penzien (1975) is used here to integrate equation 2.7, given simulated time histories of both the wave velocity and acceleration. Random Wave Forcing' / 3.5 71 This method assumes a linear variation of acceleration over each time step and also all system coefficients such as the stiffness and the damping are assumed to remain constant in each time step. The fundamental period of equation 2.7 is Tj = 2ir and therefore a time step of At = O.lTf is used in the integration routine to ensure an accurate result. This time step then fixes the frequency cutoff to Ctc = 5.0. As mentioned above, the white noise is assumed to have a spectrum equal to one over the range 0 < fl < Clc which in this case implies a variance of cr^ = 5.0 for the white noise. However, as given in Clough and Penzien (1975), the spectrum 5^(0) of a white noise time series with time increment At is c ^ 2 A i s i n 2 ( ^ ) , x 5 „ ( 0 ) = a l ~ - ^ L ± (3.29) where a\ is the variance of the set of random numbers used to model the white noise. Integration of equation 3.29 gives the variance of the white noise to be cr2 = 0.77a 2. Therefore to obtain the desired white noise properties, the variance of the random numbers must be a 2 = 5/0.77. Since most Gaussian random number generators produce random numbers with zero mean and unit variance, these numbers must be scaled by >/5/0.77 to achieve zero mean white noise with CT2 = 5.0. Finally, in the simulation process the velocity and acceleration are uncorrected and therefore two M A filters can be used, one to simulate the velocity time history from the velocity spectrum and the other to simulate the acceleration time history from the acceleration spectrum. 3.5 RESPONSE SPECTRA AND M E A N SQUARE RESULTS The results to be given in this section are concerned with the response spec-trum and the mean square value of the response computed by each of the six methods described above. The pertinent equation for the response spectrum and the corresponding notation used as an abbreviation for each approximate method are shown in Table 3.1. The Monte Carlo numerical simulation scheme is used in this section to compute the mean Random Wave Forcing / 3.5 72 Table 3.1 Solution methods for random forcing. M E T H O D Equivalent Linearization Perturbation Linear Constant Drag Damping Cubic Constant Drag Damping Time Dependent Drag Damping Numerical Simulation N O T A T I O N E Q L I N P E R T L I N E A R C U B I C T I M E D E P M O N T E E Q U A T I O N 3.23 3.28 3.14 3.13 3.12 square value of the response in order to assess the accuracy of each of the approximate methods. Numerical simulations are only done for the case of zero current. Spanos and Chen (1981) have compared the results of simulation studies to the method of equivalent linearization for the nonzero current case. For all of the computations the parameters used are Cd — 1.45 , Cm = 1.25 and fty = 0.02 with the wave particle kinematics obtained at the water surface in deep water, s = d = 200m. Four values of the nonlinearity parameter are considered, namely e = 0.05, 0.15, 0.25, and 0.50, along with four velocity ratios a = 0.0, 1.0, 2.0 and 3.0. The P - M spectrum is used throughout and is given in dimensionless form as (see Appendix A2) where B = 0.74, A = 0.0081 and g = 9.81m/s2 is the acceleration of gravity. The diameter of the cylinder is assumed to be D = l m and the natural frequency is = 1 rad/s. These values, in addition to being quite reasonable also mean that all dimensionless lengths and Random Wave Forcing / 3.5 73 times have the same numerical values as do their dimensional counterparts due to the nondimensionalization of equation 2.6. The P - M spectrum is characteristic of a well developed sea state and only depends on a single parameter, the wind speed Vw. Three values of wind speed are consid-ered; Vw = 11.6 m/s, 21.8 m/s and 31.6 m/s. These values were chosen so that the peaks in the velocity spectrum and in the convolution velocity spectrum coincide with the peak in the system transfer function. This relationship is illustrated below for the cubic constant drag damping method defined by equation 3.10a and repeated here for convenience. x" + 2$effx' + x = ^U" + sU'\U'\ = fQ{t) (3.10a) The corresponding frequency domain equation is sx(n) = \Hx(n)\2sf0{u) where S/ 0 (n) is given by 5 / o(0) = e2al6(n) + {72<12 + ^ } 5 . , ( 0 ) + £ 2 £ ^ { M O ) } " a « ' r=2 a « ' Because of the simplicity of the frequency domain equation, the nature of the response spectrum is easily assessed by comparing the nature of the transfer function to the spectrum of the forcing. For example, if a peak in the forcing spectrum coincides with a peak in the transfer function then a corresponding peak in the response spectrum will also occur. If on the other hand a peak in the forcing spectrum occurs where the transfer function is essentially zero then the response spectrum will also be essentially zero. As is seen from the above equation the forcing spectrum in this case depends on the wave velocity spectrum and the velocity convolution spectra. Figures 3.1 through 3.3 show the P - M wave height spectrum, the velocity spectrum and the velocity convolution spectra of order *2 and *3 for each of the wind speeds chosen. Only these convolution spectra are shown as they are the only ones used for computational purposes. Figure 3.4 shows the transfer function for four different values of Random Wave Forcing- / 3.5 74 0.0 0.5 DIMENSIONLESS FREQUENCY 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4 0 DIMENSIONLESS FREQUENCY 2.0 2.5 3.0 3.5 4.0 > 0.0 0.5 1.0 1.5 2.0 2.5 30 35 DIMENSIONLESS FREQUENCY DIMENSIONLESS FREQUENCY 4.0 F I G U R E 3.1 Pierson-Moskowitz spectra for Vw = 11.6. Random Wave Forcing / 3.5 75 DIMENSIONLESS FREQUENCY DIMENSIONLESS FREQUENCY F I G U R E 3.2 Pierson-Moskowitz spectra for Vw = 21.8. Random Wave Forcing / 3.5 76 DIMENSIONLESS FREQUENCY DIMENSIONLESS FREQUENCY F I G U R E 3.3 Pierson-Moskowitz spectra, for Vw = 31.6. Random Wave Forcing / 3.5 77 F I G U R E 3.4 0.2 0.4 0.6 0.8 1.0 \2 1.4 1.6 1.8 2.0 DIMENSIONLESS FREQUENCY System transfer function for various effective damping values. effective damping. The choice of these particular P - M spectra with wind speeds of 11.6m/s, 21.8m/s and 31.6m/s is now evident . That is each of these has a peak near ft = 1.0, the location of the peak of the transfer function, for either the velocity spectrum and/or one of the velocity convolution spectra. As can be seen from Figure 3.4 the amplitude of the transfer function depends greatly on the value of effective damping so that the above effects should only be noticeable for relatively small effective damping values. Note that for low frequencies ft < 0.7, the transfer function has a magnitude of approximately 1.0 and therefore the response spectrum should closely match the forcing spectrum in this frequency range. The peak of the transfer function occurs at or near 0 = 1.0 and it drops off sharply for frequencies H > 1.6. Therefore the response spectrum should also drop off sharply for frequencies Ct > 1.6. A l l of these observations will be seen in the results that follow. The above discussion also applies to the method of equivalent linearization Random Wave Forcing / 3.5 78 and to the perturbation method, but now the forcing function depends on the relative velocity spectrum and it's convolution spectra. As wil l be seen later the peaks of the relative velocity spectra and those of the wave velocity spectra coincide so that the P - M spectra chosen here are the appropriate ones for all the methods under consideration. First consider the standard deviation of the response for the case of zero cur-rent. Table 3.2 shows the values computed by each of the approximate methods described above for each of the P - M spectra and range of nonlinearity parameter. For this case it is observed that the response computed for the time dependent drag damping case is essen-tially the same as that for the cubic constant drag damping case. As such, only the results for the cubic constant drag damping method are shown in Table 3.2. Numerical simula-tions were only performed for one wind speed, Vw = 31.6 m /s. For increasing nonlinearity the response increases because the magnitude of the forcing increases with increasing e, although the effective damping also increases which tends to decrease the response. Of the approximate methods, the cubic drag damping method predicts the largest response, followed in order by the perturbation method, the linear constant drag damping method and finally the method of equivalent linearization. Figure 3.5 shows the response spectra predicted by each method for a = 0, Vw = 11.6 and e = 0.25. Similar spectra are observed for the other values of e. As can be seen from Figure 3.5 there is no apparent difference in the spectra from each method. This conclusion could also have been made based on the response standard deviations predicted by each method as shown in Table 3.2. The shape of the response spectrum could be predicted because the peak in the velocity spectrum for V w = 11.6 coincides with the peak of the transfer function at ft = 1.0. Therefore the response spectrum should only have one peak occuring at ft = 1.0. Since the current is zero, only the *3 convolution spectrum contributes to the response, but because the second peak of the *3 convolution spectrum occurs at ft = 3.0 where the transfer function is essentially zero, there will be no response at this frequency. Random Wave Forcing- / 3.5 79 Table 3.2 Standard Deviation of Response for a = 0. Vw = 11.6 e E Q L I N P E R T L I N E A R C U B I C 0.05 0.248 0.251 0.250 0.253 0.15 0.482 0.489 0.479 0.486 0.25 0.618 0.632 0.610 0.619 0.50 0.808 0.831 0.790 0.804 Vw = 21.8 £ E Q L I N P E R T L I N E A R C U B I C 0.05 0.395 0.414 0.399 0.419 0.15 0.801 0.839 0.817 0.856 0.25 1.086 1.137 1.118 1.170 0.50 1.565 1.635 1.632 1.707 Vw = 31.6 e E Q L I N P E R T L I N E A R C U B I C M O N T E 0.05 0.573 0.629 0.579 0.638 0.715 0.15 1.292 1.386 1.332 1.430 1.512 0.25 1.851 1.966 1.948 2.071 2.082 0.50 2.853 2.979 3.094 3.268 3.006 Response spectra obtained by the different methods are shown in Figure 3.6 for a = 0, Vw = 21.8 and e = 0.05. The spectra predicted by the method of equivalent linearization and the linear constant drag damping method are almost identical, as are the spectra of the perturbation method and the cubic damping method. This is not surprising due to the similarity of the equations for these spectra as discussed above. Notice the contribution to the spectrum near f) = 0.5 as would be expected. The standard deviations given in Table 3.2 confirm the above observations for Vw = 21.8. Random Wave Forcing / 3.5 o-80 PERT 0.0 02 OA 0.6 0.8 1.0 12 1.4 1.6 1.8 2.0 DIMENSIONLESS FREQUENCY F I G U R E 3.5 Response spectra for a = 0, Vw = 11.6, e = 0.25. Figures 3.7 and 3.8 show two response spectra for U = 31.6 for zero current with e = 0.05 and s = 0.25 respectively. The differences in these spectra result from the different values of effective damping for each case. As will be seen later the effective damping for e = 0.05 is about 0.10 while the effective damping for e = 0.25 is about 0.40. The differences in the transfer function for increasing effective damping are evident in Figure 3.4, where there is a significant magnitude at Q = 1.0 for the lower damping value but the transfer function has a magnitude of approximately 1.5 at 0 = 1.0 for the larger damping value. Therefore for the lower damping both the response at the peak of the velocity spectrum and at the second peak of the *3 convolution spectrum give rise to a response, yielding a bimodal response spectrum. The response at Q = 1.0 results from the matching of the peak in the transfer function and the second peak of the *3 convolution spectrum. For the larger damping the transfer function does not significantly contribute at Q = 1.0 so that the response only occurs at fl = 1/3 and a unimodal spectrum results. Random Wave Forcing / 3.5 81 1 r 0.0 0.2 0.4 0.8 0.8 1.0 12 1.4 1.6 1.8 2.0 DIMENSIONLESS FREQUENCY F I G U R E 3.6 Response spectra for a = 0, V w = 21.8, e = 0.05. In order to assess the accuracy of these methods, the Monte Carlo numerical simulation method was used for Vw = 31.6 and each of the nonlinearity values. Two MA filters, each of order 29, were used to simulate velocity and acceleration time histories from their corresponding P-M velocity and acceleration spectra. A total of 500,000 time steps were used in the integration procedure. All statistics of interest were continually updated in the integration scheme in order to save on storage requirements. Table 3.2 shows the standard deviation results obtained from the simulation procedure. As can be seen, the simulation results are slightly higher than all of the approximate methods, with the perturbation method and the cubic drag damping method providing results which are the closest to the simulation results. Further results of the simulation procedure will be discussed in Section 3.6 dealing with the response probability. The reliability of the simulation results were assessed in two ways. First, for each set of parameters at least three different simulations were performed, each using a Random Wave Forcing / 3.5 O PERT EQLIN LINEAR CUBIC" 0.0 02 6.4 0.8 08 1.0 12 1.4 1.6 1.8 Z0 DIMENSIONLESS FREQUENCY F I G U R E 3.7 Response spectra for a = 0,Vw = 31.6, e = 0.05. o cvi-PERT EQLIN LYNEAR CUBIC 0.0 02 0.4 0.6 0.8 1.0 12 1.4 1.8 1.8 Z0 DIMENSIONLESS FREQUENCY F I G U R E 3.8 Response spectra for a = 0, Vw = 31.6, e = 0.25. r Random Wave Forcing / 3.5 83 different seed value in the random number generator. The response variance after 500,000 time steps for the various seeds were not significantly different, as the variances were within 2% of each other. Second, the simulation procedure was performed for 106 time steps for one case as will be described later in the section on the probability density of response maxima. Again the response variance was not significantly different from the variance obtained after 500,000 time steps. With this in mind, some comments about the accuracy of the approximate methods are in order. Taking the Monte Carlo simulation results shown in Table 3.2 as representing the true solution, the perturbation method gives the best estimate of the mean square response for the complete range of the nonlinearity. As can be seen from Figures 3.5 through 3.8, the response spectrum predicted by each method are of similar shape. The only difference is that the perturbation method and the cubic drag damping method usually predict slightly larger spectral peaks. It is felt that the perturbation method gives the best estimate of the response spectrum due to its performance with respect to the mean square values mentioned above. It is interesting to note that the method of equivalent linearization has proven to be a very robust method for the solution of nonlinear random vibration problems. It gives accurate results of the mean square response regardless of the size of the nonlinearity for many problems of interest (see for example Spanos 1981). Also the method is usually found to underpredict the response variance in the range 0 — 20 percent depending on the problem. The results given in Table 3.2 are consistent with this observation. Further, the perturbation method is a consistent improvement to the method of equivalent lineariza-tion, which is another reason for concluding that the perturbation method most closely represents the true solution among the approximate methods. The effective damping predicted by the constant drag damping methods and by the perturbation and equivalent linearization methods are given in Table 3.3 for the zero current case. As expected the effective damping increases with increasing nonlinearity Random Wave Forcing / 3.5 Table 3.3 Effective Damping for a = 0. 84 Vw = 11.6 e E Q L I N + P E R T L I N E A R + C U B I C 0.05 0.048 0.050 0.15 0.101 0.109 0.25 0.154 0.169 0.50 0.292 0.318 Vw = 21.8 e E Q L I N + P E R T L I N E A R + C U B I C 0.05 0.075 0.076 0.15 0.181 0.189 0.25 0.282 0.301 0.50 0.517 0.582 Vw = 31.6 e E Q L I N + P E R T L I N E A R + C U B I C 0.05 0.102 0.104 0.15 0.260 0.271 0.25 0.407 0.438 0.50 0.737 0.856 and increasing wind velocity. The constant drag damping methods predict larger effective damping than the other two methods. Usually, a larger damping corresponds to a smaller response but here, for example, the linear constant drag damping method predicts a larger response than does the method of equivalent linearization. This can be explained by writing the equations for these two cases. For the linear constant drag damping case the governing equation can be shown to be i " + 2$effx' + x = 7 t f " + e\fl is the second order correlation function. Obviously the third moment is found by setting both r x and r 2 equal to zero in these expressions with the result Rx(o,o) =< x3(0 >= m 3 = f f Bjn1,n2)dnldn2 which is analogous to the relationship between the second moment and the response spec-trum. The fourth moment can be found in an analogous manner by defining a third order correlation function as Rz(Tl,T2, TS) =< X(t)x{t + Tx)x{t + T2)X(t + T3) > and a trispectrum Cx(Cli, n 2 ,n 3) as the three dimensional Fourier transform of the third order correlation function. Setting all r's to zero gives the fourth moment as £ . ( 0 , 0 , 0 ) = < x \ t ) >= m 4 = / f f c x(n 1,n 2,n 3)(fn 1dn 2dn3 JQi Jn2 Jn3 Random Wave Forcing / 3.6 92 By defining higher order correlation functions and spectra the expressions for the higher order moments can be found. To show the application of higher order spectra in finding higher order re-sponse moments, consider a linear equation with general forcing f(t) x" + 2$x' + x = f{t) (3.31) This equation is linear and the input-output relationship in the frequency domain is sx(n) = \Hx(n)\2sf(n) as discussed in Appendix A3. The second response moment is then found by integration. A similar procedure can be used to obtain the input-output relations for the higher order spectra which are Bx(nu n2) = + Q ^ H ^ H ^ B ^ , ft2) cx(nun2M = Hi{nx + n2 + C I ^ H ^ H ^ H ^ C ^ Q , , ^ ) where the superscript c denotes the complex conjugate. The response moments are then found by integration. Now it remains to evaluate the bispectrum and the trispectrum of the forcing function in terms of the known input forcing spectrum. This can be done in principle but the details soon become very complicated. For the case of Gaussian forcing the response will also be Gaussian. That is if P(/) = J - e x p ( - ^ ff) \/2ir Borgman's expression is Pif) • V ( o , « + 6 X P l 1 6 * ? ^ ' U l a - tyrft ) where U(Q, t) is the parabolic cylinder function which is related to the modified Bessel functions Ii/i{t2/4) and I_1/±(t2/4). Using the relations found in Abramowitz and Stegun (1965): U(0,t)=1-(nt)1/2 t2 t2 1 / - l / 4 ( 7 ) - / l / 4 ( T ) a 1 t > 0 t < 0 ' - 1 / 4 ( 4 ) + fy(T) Pierson and Holmes (1965) have found expressions for the second and fourth moments of the force probability density, which depend only on two parameters; V i and That is a) = 3V? + $ mif = 10501 + 18V»2V2 + 302 which can be inverted to give (3.33) Random Wave Forcing / 3.6 95 The odd moments are zero here because of symmetry for zero current. By substituting equation 3.33 into the equation for p(f), the desired form of the probability density of the force in terms of / and the moments of / is obtained. Therefore by the above argument the probability density of the response is obtained by replacing / by x and the moments of / by moments of x in the above equation. Thus the fourth response moment is required here and the method outlined above using the higher order spectra can be used. In order to illustrate how well the approximate probability density obtained from the probability density of the force approaches the actual probability density the results of the simulation procedure can be used. Figures 3.12 through 3.15 show the probability densities of the response for e = 0.05, 0.15, 0.25, and 0.50 respectively for the P - M spectrum with Vw = 31.6. The approximate probability density is obtained as described above using the values of the second and fourth response moments computed by the numerical simulation. Also shown are the Edgeworth series form of the probability density and the Gaussian density. The M A simulation procedure was used with a total of 500,000 time steps to integrate the equation of motion. A l l response statistics were computed over the length of a single record. However, the simulations were run for a few different seed values for the random number generator to check that consistent results from the simulation were indeed obtained. Table 3.6 shows the second and fourth response moments, the kurtosis and the maximum and minimum response values computed in the simulation procedure. The approximate probability density is very close to that of the simulations for all values of nonlinearity. As mentioned above, Tickell (1977) has shown the same result using experimental results. The Edgeworth series form of the probability density is accurate for small values of i but shows oscillations for larger response values. If the Wave Forcing / 3.6 £=0.05 F I G U R E 3.12 Response probability density e = 0.05. 55 = o : cu -CO Kb £=0.15 EDGEWORTH GAUSSIAN +""siMULATi6"N - i — -3.0 -6.0 -3.0 0.0 3.0 6.0 X F I G U R E 3.13 Response probability density e = 0.15. Random Wave Forcing- / 3.6 97 E-< i—( CO £ J w-Q'oJ T — E -O -5t <: PQ o w'o. CO o (X co Kb £=0.25 APPROXIMATE EDGEWORTH GAUSSIAN +"" ~S I M~uL ATI ON"~ -6.0 -3.0 T — 0.0 X 3.0 6.0 F I G U R E 3.14 Response probability density e = 0.25. >-• E E-» W'o_ CO -o CO KTo_ £=0.50 APPROXIMATE EDGEWORTH '^"'GAUSSIAN +"" STMULATI ON"" -8.0 -4.0 o.o X " i — 4.0 8.0 F I G U R E 3.15 Response probability density s = 0.50. Random Wave Forcing- / 3.6 Table 3.6 98 Simulation results for second and fourth response moments. £ mix A 4 xmax xmin 0.05 0.511 1.219 1.668 5.34 -5.41 0.15 2.287 25.47 1.871 11.0 -11.5 0.25 4.334 87.90 1.697 14.4 -14.8 0.50 9.035 323.54 0.963 18.1 -17.7 kurtosis were larger than indicated in Table 3.6, the probability density would become negative and have a multimodal characteristic for the larger response values. The Gaussian density underpredicts the response for large response values, which is the reason for trying to find the non-Gaussian response probability density. In the design process, the maximum values of the response are often of paramount importance. Therefore in order to make probability statements about the re-sponse maxima, the corresponding probability density of the maxima is required. Maxima are defined by the conditions z'(i) = 0 and x"{t) < 0 and it is assumed that the statistics of the minima of the response are the same as those of the maxima. Therefore in order to find the probability density of the response maxima, the joint probability density p(x,x',x") is required. This could be obtained by the above hypothesis from the corresponding joint probability density of the forcing p(f, / ' , / " ) . In principle it is possible to find p(f, /', /") using the technique of transformation of variables, but it is often very difficult to directly integrate out the auxiliary transformation variables. This results in an integral expression for the joint density which is of limited use for the purposes of further analysis. With the assumption of a narrow banded process it can be shown (see for example L in 1976) that only the joint probability density p(x, x') is required to make probabilistic statements about the the response maxima. Therefore if p(f, / ' ) were known, then using the above hypothesis, the corresponding joint probability density of the response would be known. Transformation of variables can be used to find p(/, / ' ) in integral form Random Wave Forcing / 3.6 99 with the integration of the auxiliary variable again being quite difficult. However Moe and Crandall (1978) have performed the integration asymptotically for large values of / . The resulting joint probability density of the force can then be used to determine the asymptotic probability density of the force maxima fm as where ^ i = 7 c v and rp2 = e = w(0 Assuming both f(t) and w(t) to be Gaussian processes, implies that the operator § of the second equation is linear. Therefore this equation can be given in the frequency domain as sf(Q) = \H9(n)?$ where $ is the spectral density of the white noise. The form of the equation for Q can be determined by matching the known spectrum of /(f) to the transfer function associated with Q. Finally, the expanded system of equations with white noise forcing is 9{£{x(t)}} = w(t) for which a F P K E can be derived. This approach has often been suggested, most recently by Crandall and Zhu (1984), but to date no results have appeared in the literature. CHAPTER 4 Summary and Conclusions The aim of this thesis is to investigate various methods that can be used to obtain a steady state solution to the problem of regular or random wave loading on a slender flexible member of an offshore structure. The Morison equation, more specifically the Morison equation modified to include the relative velocity between the fluid and the structure, is used to model the fluid force on the structure. This modification accounts for the effects of the motion of the structure on the fluid force. A simple single degree of freedom model of the structure is used to model the structure for ease of comparison of various solution methods. The resulting equation of motion is nonlinear in both the response velocity and the wave velocity and at present there are no analytical solutions available, even for the simplest wave representations. Therefore approximate or numerical methods are required to obtain a solution. As is usually the case, the numerical solution is used to assess the accuracy of the approximate solutions. An alternative modified form of the Morison equation, the independent flow fields formulation, is also briefly discussed. The governing equation in this case is similar to the relative velocity formulation in that it is a nonlinear differential equation in the response velocity with the forcing given by the wave particle kinematics. The methods that can be used to solve the independent flow fields formulation equation are therefore similar to those described here for the relative velocity formulation. Only solutions of the relative velocity formulation are discussed here. 102 Summary and Conclusions / 4.1 10S It should be mentioned at this point that at present there is no theoretical justification for either the relative velocity or the independent flow fields formulation being the correct one to account for the effects of fluid structure interaction in the problem of wave forces on slender offshore structures. Experimental verification is also not presently available. Therefore one important aspect of wave forces on flexible offshore structures that requires further attention is to determine, most likely by small scale experiments, the appropriate equations required to include the effects of fluid structure interaction and the range of parameter values for which they are applicable. Solutions to the relative velocity equation of motion are obtained by various existing methods, such as the method of equivalent linearization, various other linearization methods, and numerical integration for both regular and random waves. A new perturba-tion method, based on Lindstedt's method, is also used to obtain a solution. The equation of motion is solved for typical values of the parameters that characterize the problem of wave forces on offshore structures. Further, the effect of a constant current is also treated for both regular and random waves. As mentioned above, the numerical integration so-lutions are used to assess the accuracy of the various approximate methods. The above solutions refer to the frequency response curves for the case of regular wave forcing and to the response spectrum and the response mean square value for the case of random wave forcing. In addition, for random wave forcing, the probability density of the response is also discussed. The conclusions regarding this thesis are naturally grouped into three dif-ferent categories. These are regular wave forcing, the mean square response and response spectrum for the random wave forcing and finally, the response probability. Each of these will be discussed in turn. 4.1 REGULAR WAVES The governing equation of motion for the case of regular waves is a nonlin-ear ordinary differential equation with harmonic forcing. Solutions of this equation are Summary and Conclusions / 4.1 104 obtained by numerical integration, various linearization techniques which include (a) zero drag damping, (b) constant drag damping and (c) time dependent drag damping, the method of equivalent linearization and the perturbation method. As previously stated the numerical integration solution is used to assess the accuracy of each of the approximate methods. Solutions are obtained for both a simple harmonic wave and a harmonic wave with a constant superposed current. The solutions are presented as frequency response curves giving the maximum response as a function of the wave frequency. The importance of the nonlinear response term is denoted by the parameter e, which also represents the magnitude of the forcing term. Therefore there are two opposing effects that occur for increasing e. First, increasing e implies increasing hydrodynamic damping hich tends to decrease the response, while on the other hand, increasing e increases the forcing which increases the response. For smaller values of s and zero current the frequency response curves show peaks at dimensionless wave frequencies of ft = 1, which is the primary resonance, and at ft = 1/3 which is a secondary resonance. Larger values of e tend to flatten out these resonant peaks. If a constant current is included, another secondary resonant peak occurs at a frequency of ft = 1/2. Of the various linearization techniques, the zero drag damping case is not appropriate in this case because it grossly underpredicts the damping and results in an overprediction of the response, especially at the resonant frequencies. Very little difference was observed between the constant drag damping method and the time dependent drag damping method. However, in comparison with the numerical solution, these methods underpredict the response near the primary resonant frequency. This is due to the fact that the structure's velocity is the same order as the wave velocity, which contradicts the basic assumption of the validity of these methods. The linearization condition that the ratio of the structure's velocity to the wave velocity is small must always be checked when using either of these methods. Summary and Conclusions / 4.2 105 In contrast to the above linearization methods which retain the nonlinear wave velocity term, the method of equivalent linearization linearizes both the nonlinear response term and the nonlinear wave velocity term. Due to this linearization scheme, this method is unable to predict the secondary resonant peaks at all, but approximates the response near the primary resonance peak at ft = 1 very well. A perturbation method based on Lindstedt's method is also described and is shown to be closely related to the method of equivalent linearization. That is, the effective damping, which is the sum of the stuctural and the hydrodynamic damping, estimated by each method is identical. Therefore the perturbation method is seen to represent a con-sistent improvement to the method of equivalent linearization. The perturbation method is able to predict the secondary resonant peaks with improvement in the solution being realized by increasing the number of terms in the perturbation series. The perturbation method also estimates the response near ft = 1 very well. Of all the approximate methods, the perturbation method most closely ap-proximates the numerical solution for the complete wave frequency range. 4.2 RANDOM WAVES 4.2.1 Response spectrum and mean square value In this case the waves are considered to be represented by a Gaussian stochas-tic process. Constant current is included by considering a process with a nonzero mean value. Two of the response statistics of interest are the response spectrum and the mean square value of the response. These are average response statistics, and are considered in this section. Conclusions regarding the probability density of the response and extreme value statistics are discussed in the next section. A numerical integration solution or Monte Carlo simulation is used to com-pute the response mean square value, in order to assess the accuracy of each of the ap-proximate methods. The basic approach taken for each of the approximate methods is to Summary and Conclusions / 4.2 replace the original nonlinear equation of motion with one or more linear equations which are then solved using standard techniques. In this regard the use of Hermite polynomials are shown to be appropriate due to their nice orthogonality properties. The approximate methods are divided into two groups: the method of equiv-alent linearization and the perturbation method; and the linear constant drag damping, the cubic constant drag damping and the time dependent drag damping methods. For the case of zero current there is no appreciable difference between the cubic constant drag damping and the time dependent drag damping methods. Only a small difference between these two methods is noted for the case of a constant current. As is the case for the regular wave forcing, the effective damping estimated by the method of equivalent linearization and by the perturbation method are identical. The perturbation method predicts a better estimate of the response mean square value than does the method of equivalent linearization in comparison with the simulation results. 4.2.2 Response p robab i l i t y density The major drawback of each of the approximate methods is that they are able to find only the average response statistics such as the response spectrum and the mean square value of the response. They are not able to give useful information about the re-sponse probability density and the extreme response values which are of major importance in the design situation. The numerical simulation solution does yield information about the probability density of the response, however at the cost of excessive computational effort. The cubic drag damping method results in a linear equation of motion with the forcing given by the Morison equation. Previous investigations have derived, either exactly or approximately, the probability density of the forcing. This probability density function is non-Gaussian and therefore the response probability density function is also non-Gaussian, but of unknown form. It has been hypothesized that for a linear equation with non-Gaussian forcing, the response probability is of the same form as the forcing Summary and Conclusions / 4.3 107 probability density. This hypothesis is shown to be approximately true using the results of the numerical simulation for both the response probability density and for the probability density of the response maxima. Probablity statements that are useful in the design situa-tion follow from the latter probability density. In order to obtain the response probability, both the second and fourth response moments are required. The second response moment is the mean square value of the response and is easily obtained from the response spectrum while the fourth moment is obtained from the higher order trispectrum. Unfortunately, the computations required to find the fourth moment are quite involved and it remains to find simple effective estimates of the fourth moment of the response. An alternative approach to find the probability density of the response is by using the Markov vector approach. An enlarged system of equations is required in order to model the non-white noise nature of the random waves. This approach has often been suggested for nonlinear random vibration problems, but as yet no solutions using this method have appeared in the literature. 4.3 RECOMMENDATIONS FOR FURTHER STUDY Some specific recommendations for further studies based on the work in this thesis are: 1. Experimentally verify the form of the fluid force equation required to ac-count for the wave-structure interaction for wave forces on slender offshore structures. 2. Extend the analysis presented here to the case of a more realistic multi-degree of freedom model, especially with respect to the perturbation method. 3. Determine a simple effective method to determine higher order moments of the response of a linear system subjected to non-Gaussian forcing. 4. Carry out the Markov vector approach for a simple nonlinear system with non-white noise forcing. JOS References A B R A M O W I T Z , M . and S T E G U N , I .A. (1965) Handbook of Mathematical Functions Dover, 1046p. A T A L I K , T.S. and U T K U , S. (1976) Stochastic Linearization of Mul t i Degree of Freedom Nonlinear Systems Earthquake Engineering and Structural Dynamics Vol . 4, pp. 411-420. B L E V I N S , R .D. (1977) Fiow Induced Vibrations Van Nostrand Reinhold, 363p. B O R G M A N , L . E . (1967) Random Hydrodynamic Forces on Objects Annais of Math-ematical Statistics Vol . 38, No. 1, pp. 37-51. B O R G M A N , L . E . (1969) Ocean Wave Simulation for Engineering Design Journal of the Waterways and Harbours Division (ASCE) Vol . 95, pp. 557-583. B O R G M A N , L . E . (1972) Statistical Models for Ocean Waves and Wave Forces A d -vances in Hydrosciences Vol. 8, pp. 139-181. C A U G H E Y , T . K . (1963) 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 Vol . 35, No. 11, pp. 1687-1692. C A U G H E Y , T . K . (1971) Nonlinear Theory of Random Vibrations Advances in Ap-plied Mechanics Vol 11, pp. 209-253. C L O U G H , R .W. and P E N Z I E N , J . (1975) Dynamics of Structures Mc-Graw Hi l l , 634p. C R A N D A L L , S.H. (1963) Perturbation Techniques for Random Vibrations of Nonlinear Systems Journal of the Acoustical Society of America Vol . 35, No. 11, pp. 1700-1705. C R A N D A L L , S.H. (1966) Random Vibrations Applied Mechanics Surveys pp. 681-689. C R A N D A L L , S.H. and Z H U , W.Q. (1984) Random Vibrations: A Survey of Recent Developments Journal of Applied Mechanics (ASME) Vol . 50, pp. 953-962. D R A P E R , N . R . and T I E R N E Y , D .E . (1972) Regions of positive and unimodal series ex-pansion of the Edgeworth and Gram-Charlier approximations Biometrika Vol. 59, pp.463-472. E A T O C K - T A Y L O R , R. and R A J A G O P A L A N , A . (1982) Dynamics of Offshore Struc-tures, Part 1: Perturbation Analysis Journal of Sound and Vibration Vol . 83, No. 3, pp. 401-416. G R E C C O , M . G . and H U D S P E T H , R .T . (1983) Stochastic Response of Prototype Off-shore Structures Journal of Structural Engineering (ASCE) Vol. 109, No. 5 pp.1119-1138. 109 G U D M E S T A D , O.T. and C O N N O R , J .J . (1983) Linearization Methods and the In-fluence of Current on the Nonlinear Hydrodynamic Drag Force Applied Ocean Research Vol. 5, No. 4, pp. 184-194. H O O F T , J.P. (1982) Advanced Dynamics of Marine Structures J . Wiley and Sons, 345p. I W A N , W.D. (1974) Application of Nonlinear Analysis Techniques Applied Mechanics in Earthquake Engineering A S M E A M D Vol. 8, pp. 135-161. K E N D A L L , M . and S T U A R T , A . (1977) The Advanced Theory of Statistics, Vol. 1 MacMillan, 4th edition 472p. K E V O R K I A N , J . and C O L E , J .D. (1981) Perturbation Methods in Applied Mathemat-ics Springer Verlag 558p. L A Y A , E . J . , C O N N O R , J .J . and S U N D E R , S.S. (1984) Hydrodynamic Forces on Flex-ible Offshore Structures Journal of Engineering Mechanics (ASCE) Vol . 110, No. 3, pp. 433-448. L I N , Y . K . (1976) Probabilistic Theory of Structural Dynamics R. E . Krieger, 368p. M A L H O T R A , A . K . and P E N Z I E N , J . (1970a) Response of Offshore Structures to Ran-dom Wave Forces Journai of the Structural Division (ASCE) Vol . 96, No. 10, pp. 2155-2173. M A L H O T R A , A . K . and P E N Z I E N , J . (1970b) Nondeterministic Analysis of Offshore Structures Journai of Engineering Mechanics Division (ASCE) Vol. 96, No. 6, pp. 985-1003. M O E , G . and C R A N D A L L , S.H. (1978) Extremes of Morison-Type Loading on a Single Pile Journal of Mechanical Design (ASME) Vol. 100, pp. 100-104. M O R I S O N , J.R., O ' B R I E N , M.P. , J O H N S O N , J .W. and S C H A F F , S.A. (1950) The Force Exerted by Surface Waves on Piles Petroleum Transactions of the American Institute of Mining and Metallurgical Engineers Vol . 189, pp. 149-154. N A Y F E H , A . H . (1973) Perturbation Methods J.Wiley and Sons, 600p. N E W L A N D S , D .E . (1975) Random Vibrations and Spectral Analysis Longman, 285p. NIG A M , N . C . (1983) Introduction to Random Vibrations M I T Press, 341p. OCHI , M . K . (1982) Stochastic Analysis and Probabilistic Prediction of Random Seas Advances in Hydroscience Vol . 13, pp. 217-375. P I E R S O N , W . J . and H O L M E S , P. (1965) Irregular Wave Forces on a Pile Journal of the Waterways and Harbours Division (ASCE) Vol . 91, pp. 1-10. R O B E R T S , J .B . (1981a) Response of Nonlinear Mechanical Systems to Random Exci-tation, Part 1: Markov Methods Shock and Vibration Digest Vol . 13, No. 4, pp. 17-28. R O B E R T S , J .B . (1981b) Response of Nonlinear Mechanical Systems to Random Ex-citation, Part 2: Equivalent Linearization and Other Methods Shock and Vibration Digest Vol. 13, No. 5, pp. 15-29. 110 R O B E R T S , J .B . (1984) Techniques for Nonlinear Random Vibration Problems Shock and Vibration Digest Vol . 16, No. 9, pp. 9-14. S A R P K A Y A , T . and I S A A C S O N , M . (1981) Mechanics of Wave Forces on Offshore Structures Van Nostrand Rheinhold, 65lp. S H I N O Z U K A , M . and W E N , Y . K . (1971) Nonlinear Dynamic Analysis of Offshore Structures Proceedings of the First International Symposium on Stochastic Hydraulics C . L . Chiu ed. pp. 507-521. S M I T H , E . (1978) On Nonlinear Random Vibrations Division of Structural Mechanics, Norwegian Institute of Technology, Report No. 78-3, 167p. SOONG, T .T . (1973) Random Differential Equations in Science and Engineering Aca-demic Press, 327p. SPANOS, P.T.D. (1981) Stochastic Linearization in Structural Dynamics Applied Mechanics Reviews Vol . 34, No. 1, pp. 1-8. SPANOS, P.T.D. (1983) A R M A Algorithms for Ocean Wave Modelling Journal of Energy Resources Technology (ASME) Vol . 105, pp. 300-309. SPANOS, P.T.D. and C H E N , T . W . (1981) Random Response to Flow-Induced Forces Journal of the Engineering Mechanics Division (ASCE) Vol . 107, pp. 1173-1190. T I C K E L L , R . G . (1977) Continuous Random Loading on Structural Members The Structural Engineer Vol . 55, No. 5, pp. 209-222. T O , C. W.S. (1984) The Response of Nonlinear Structures to Random Excitation Shock and Vibration Digest Vol . 16, No.4, pp. 9-14. T U A H , H . and H U D S P E T H , R .T . (1982) Comparison of Numerical Random Sea Simula-tions Journa of the Waterways, Port, Coastal and Ocean Division (ASCE) Vol. 108, No. WW4, pp. 569-584. T U N G , C .C . and H U A N G , N . E . (1976) Interactions Between Waves and Currents and their Influence on Fluid Forces Proceedings of BOSS' 76 pp. 129-143. V A N M A R C K E , E . H . (1979) Some Recent Developments in Random Vibrations Ap-plied Mechanics Reviews Vol . 32, No. 10, pp. 1197-1202. W H I T E , F . M . (1974) Viscous Fluid Flows Mc-Graw Hi l l , 725p. APPENDIX A l Review of the Solution of Linear ODE's with Harmonic Forcing In this Appendix the particular solution of the linear equation x" + 2$x' + x = f(t) is given where f(t) can be represented as a Fourier series f(t) = cQ + ^ c n cos nQt + dn sin nfif n=l The Fourier coefficients are found from the usual equations c 0 = r - / f(t)dt cn = - / f(t) cos nCltdt {Al.l) * J-it/n dn = - f{t) sin nClt dt * J-x/n The solution is also written as a Fourier series x(t) = aQ + ^ an cos nQt + bn sin nQt n=l where the coefficients o n and bn are found by substitution and then equating like harmonic terms. This procedure gives a 0 = c 0 _ c w ( l - (nfl)2) - dn2(nQ n ( l - ( n n ) 2 ) 2 + (2cnn)2 d n ( l - ( n n ) 2 ) + c„2^nn n (1 - (nft)2)2 + (2cnft)2 111 Review of the Solution of Linear ODE's with Harmonic Forcing / A1.0 112 This solution is frequently used for the linearization and perturbation methods described in Chapter 2. It should be repeated here that the above solution represents the steady state solution of the equation of motion which is the solution of interest in the context of deterministic wave forces on offshore structures. Due to system damping the homogeneous solution will quickly die out and is not of interest here. The remainder of this Appendix gives the Fourier series expansion of U'\U'\ which is required in the development of the methods of Chapter 2. The wave displacement is given by equation 2.10 and is repeated here for convenience U= E^.at + ^coscit (2.10) so that differentiation with respect to t gives the wave velocity as U' = ^-(a-smQt) Now the Fourier series expansion of t7'|U'j is U'\U'\ = c 0 + c„ cosnfli + dn sin nQt n=l with the Fourier coefficients given by equations A l . l . Two cases need to be considered, a < 1 and a > 1. When a < 1, by appropriately splitting up the integrals, the Fourier series coefficients can be shown to be c°= \ (^ir)2 [(°2 + \ ) w + 4acos/3~ rin2/3] / s in (2-n) /? sin(2 + n)/?\ ^ " \ (2 - n) + (2 + n) ) ) Review of the Solution of Linear ODE's with Harmonic Forcing / A1.0 113 1 fUmQ\2 f, 2 1.4 0 A / s i n ( l - n ) / 3 sin(l + n ) £ \ n—1,3,— / cos(2 - n)/3 cos(2 + n)p\ \ ~~ V ( 2 - n ) + (2 + n) / J where sin/? = a; 0 < / ? < T T /2 . When there is no constant current, a = 0, /? = 0 and the Fourier series expansion reduces to ir\ D J ^-f n(4-n2) as given by Blevins (1977). The second case considers a > 1 where U'\U'\ = U'2 for all t and the Fourier series expansion is easily shown to be given as Note that there are only two harmonic components in this expression along with a constant term. These Fourier series expansions are similar to those reported by Gudmestad and Connor (1983), the only difference being a few sign changes from + to — because they define their velocity as U' = ^ ( a + s inm) which is slightly different from the definition used here. APPENDIX A2 Review of Linear Wave Theory Linear wave theory forms the basis for most random representations of a sea state because of the simple relationships that exist for Gaussian random processes as they pass through linear systems. The governing equation for linear wave theory is the Laplace equation in two dimensions d26 d2A where is the velocity potential. The boundary conditions that apply in linear wave theory are *L = o z * = _ + gn* = 0 z* = 0 (A2.2c) where the two dimensional region to be considered is shown in Figure A2.1 with n* being the dimensional water surface elevation above z* — 0. The linearization assumptions made in the boundary conditions are that equations A2.26 and A2.2c are applied at the still water surface instead of the actual water surface and that both these equations neglect nonlinear terms. Instead of the two free surface conditions above, the following two can be used. 114 Review of Linear Wave Theory / A2.0 115 F I G U R E A2.1 Definition sketch for linear wave theory. Another condition that must be satisfied is {x*,z*,t*) = (x* - ct\z*) which ensures permanence of form or the periodic nature of the wave. The complex solution of this boundary value problem is ,/ * * *v iHgcoshk*(z* + d*) . . where the wave crest f?* = H/2 passes the location x* = 0 at time t* = 0, k* is the wave number related to the wavelength L by k* = 2x/L and u is the wave frequency related to the wave period T by u = 2ir/T. A l l boundary conditions except equation A2.2d have been used to this point. This last boundary condition results in the dispersion relationship between the wave fre-quency and the wave number. Substituting equation A2.3 into equation A2.2d gives — = k*d* tanh k*d* 9 Review of Linear Wave Theory / A2.0 116 which is a transcendental equation for the parameter k*d*. The relationship between the wave kinematics and the surface elevation follow directly from equation Al.lt 1 d, Hiru cosh k*d* and the definition of velocity potential W = a ? = Y Mhk-d- e x p { ' ( * x ~ u t ) ) The desired relationship between the surface elevation and the velocity is obtained by combining these equations to give du* _ coshfc'fV +d*) , dt* ~ U sinh A M * V Applying the non-dimensionalization of equation 2.6 yields , _ ^cosh ks sinh kd ^ where s = z + d is the dimensionless distance upwards from the seabed and the ' denotes differentiation with respect to t. In order to find the spectral relationships of interest, first the correlation function RUI{T) is determined and then taking the Fourier transform gives the spectrum. i ^ r ) = <«'(*)«'(*+ r)> ^9 cosh 2 ks _ , . sinh 2 fed , V ' Taking the Fourier transform gives The general input-output relationship for linear systems finally results in sU.(Q) = |H..,(n)| a5,(n) where |iY t t# v(n)| 2 is the transfer function between the velocity and the wave height spectrum and is obviously given by UWn)|> - « 2 5 ? ^ sinh kd Review of Linear Wave Theory / A2.0 117 Some other relationships of interest are those between the wave acceleration spectrum, the wave displacement spectrum and the wave height spectrum. For stationary random processes it is easy to show (see Appendix A3) that Mft) = n2s„(ft) 5tt»(n) = n25t,,(n) so that the desired relationships are For the method of equivalent linearization the forcing is given in equation 3.16 as /(<) = ( l - 7 ) u " + 2 ^ u ' + u The input-output relationship of this equation in the frequency domain is Su((l) = \HM(Q)\2Sf{Q) where it is easily shown that ' ^ ( n | 2 = (i - [i - 7]n2)2 + (2!Nw Review of Linear Wave Theory / A2.1 118 A2.1 WAVE HEIGHT SPECTRA A popular wave height spectrum often used in offshore engineering calcu-lations is the Pierson-Moskowitz (P-M) spectrum. The analytical form of this spectrum along with some of it's important properties are given below. Also shown are the results of simulating surface elevation records consistent with a given wave height spectrum by digital filtering. The results shown here are applicable to the simulation of any wave record based on the associated target wave spectrum. First the dimensional form of the P - M spectrum is M « ) = ^ « P ( - ^ ) where g is the acceleration of gravity, Vw is the wind speed and the usual values of the constants are A = 0.0081 and B = 0.74. This spectrum is a one parameter spectrum depending only on the wind speed Vm and it applies to a fully developed sea state. Alter-native forms of the P - M spectrum can be written in terms of other parameters such as the significant wave height for example, instead of the wind speed. The dimensionless wave height spectrum is given by w-&«p(-sO where the nondimensionalization of equation 2.6 has been used with fl=_V_ b = Bg* Of importance is the peak frequency flp, at which the maximum spectral amplitude occurs. This is given by O p = (0.86)1/* The spectral moments are Review of Linear Wave Theory / A2.1 119 which can be shown to exist for n < 4 only. In order to evaluate higher order moments the spectrum must be truncated at some upper frequency limit. The moving average (MA) filter can be used to simulate a time history of the surface elevation consistent with a Pierson-Moskowitz wave height spectrum as described by Spanos (1983) and summarized in Chapter 3. The filter is completely described by the filter coefficients given by c, = ^ n c v ^ ( n ) c o s ( ^ ) d n where S(U) is the target P - M spectrum and ftc is the cutoff frequency of the band limited white noise input to the filter. The transfer function of the filter is given by ^ / ( ^ ^ ( c o + E c . c o s ^ ) 2 l/=l c Spanos uses a filter of order q = 29 to simulate a time series with the P - M wave height spectrum used as the target spectrum. This is also done here for Vw = 31.6 and Figure A2.2 shows the target spectrum and the spectrum of the time series. The spectrum of the time series is su(Q) = |^(n)| 25 n(n) where Sn(ft) — 1 is the assumed spectrum of the white noise input to the filter. The agreement between these is excellent. In the numerical simulations described in Chapter 3, the P - M velocity and acceleration spectra are used as target spectra to generate simulated velocity and acceler-ation time histories which are then used to specify the fluid force in the Morison equation. Numerical integration of the equation of motion then yields the response. F I G U R E A2.2 Pierson-Moskowitz wave height spectrum with Vt Review of Linear Spectral Theory / A3.0 121 APPENDIX A3 Review of Linear Spectral Theory Consider a single degree of freedom system with random forcing governed by the equation x" + 2fx* + x = f(t) (A3.1) which often occurs in the linearization methods and the perturbation method described in Chapter 3. The solution of this equation can be achieved in both the time domain via the convolution integral and in the frequency domain. What is to be described here are the details of the solution in both the time and frequency domains and the relationships between them. In the time domain the solution is given by the convolution integral *(*)= f h{t-6)f(6)d9 (A3.2) J — CO where h( ) is the impulse response function of equation A3.1 which is t / , e -' 1 sin v/1 - c21 MO = /r-^ f < 1 (A3.3) v i - r The impulse response function only depends on the system characteristics, specifically the damping and the natural frequency (= 1 in this case). The impulse response function is causal, so that h{t) = 0 if t < 0. That is the effect of the impulse is only felt after the impulse has occurred. Therefore h(t — 0) = Q if 6 > t and the upper limit of integration in equation A3.2 can be replaced by oo without loss of generality. Review of Linear Spectral Theory / A3.0 122 An alternative form of the convolution integral is obtained by letting £ = t—8, which transforms equation A3.2 to x(t)= f ° - f l d * U3.4) Jo Again since the impulse response function is causal, the lower limit of integration can be replaced by — oo without loss of generality. Either of the forms of the convolution integral can be used to obtain the solution given the time history of the forcing. For the case of random forcing, f(t) is not usually known explicitly and direct integration only applies to simulated forcing histories. The computational effort in this case would generally be excessive. Fortunately the solution in the frequency domain is much simpler and more widely used for random forcing. The frequency response function of equation A3.1 is found by setting /(<) = etQt and letting the response be x(t) = H(Cl)etnt where H(Cl) is the frequency response function. These substitutions give 1 tf(n) = I - n 2 + »'2f n 1 - ft2 2cQ (i - n 2 ) 2 + (2c:n)2 (i - n 2 ) 2 + (2cn)2 The magnitude squared of the frequency response function is called the transfer function and is \H[Q)\* = ,2 1 (i - n 2 ) 2 + (2f-n)2 A very simple relationship exists between the impulse response function and the frequency response function which is h(t) = — / H{Cl)eintdCl 1 U3-5) H{ti) = / h(t)e-int dt J—CO which are recognized as a Fourier transform pair. The importance of the frequency domain approach is the simplicity of the resulting frequency domain equations. For example taking Review of Linear Spectral Theory / A3.0 123 the Fourier transform of equation A3.4 results in X(tl) = H(Q)F(Q) where X(f l ) and F(Q) are the Fourier transforms of the x(t) and f(t) respectively and /oo X{Q)eiQt dQ /"/.oo (A3.6) X(Cl) = - J j(t)e~iUldt and similarly for the pair /( i)and F(Cl). The definition of the Fourier transform pair in equation A3.5 applies only to the relationship between the impulse response function and the frequency response function, whereas in general all other Fourier transform pairs are defined by equation A3.6. Two important functions that are extremely important in frequency domain analysis are the correlation function and the spectral density. These functions will be described here for a single variable but the concepts are easily extended to more than one variable (see for example Clough and Penzien 1975). The correlation function of a random process is given by J J x ( * i , « 2 ) = < * ( * i ) x ( t 2 ) > where < > denotes the expected value. If the random process is stationary, as is assumed throughout this thesis, then the correlation function only depends on the time difference T = \ti — <2|- Therefore the correlation function can be written as RX{T) =< X{t)x{t + T) > (A3.7) ==RX{-T) i.e., the correlation function is an even function. This property will be frequently used in what follows. The convolution integral solution for x(t) is now used to find a relation-ship between the correlation function of the forcing and the response correlation function. Substitution of equation A3.4 into equation A3.7 gives /CO f- oo -oo J oo Review of Lin ear Spectra] Theory / A3.0 124 where integration and expectation can be interchanged as they are independent of each other. Using the stationarity property of equation A3.7 gives /OO f oo / M*i)M6)*/(* + * 1-6)^ 1(*6 U3.8) -oo J oo which is the desired relationship between the input forcing correlation function and the output response correlation function. In general the integrations required to obtain Rx(r) from this equation are quite tedious and alternative more efficient methods in the frequency domain are available. The spectral density is defined as the Fourier transform of the correlation function and taking the Fourier transform of equation A3.8 results in the simple relation-ship 5,(0) = \H(Cl)?Sf(Vl) That is, the output response spectral density is the product of the system transfer function and the input spectral density. This relationship is the reason that spectral analysis of linear random processes is so efficient, especially if the forcing is Gaussian. The above equation is only valid for linear systems because the convolution integral, which is based on the principle of superposition, is only applicable to linear systems. If the forcing function can be assumed to be Gaussian, then the response spectral density completely describes the probability density of the response. This rela-tionship is as follows. First, if the forcing is Gaussian and the system is linear then the response will also be Gaussian. A Gaussian probability is completely described by it's first two moments, the mean and the variance. The mean value can be found by taking the expected value of equation A3.4 ax =< *(*) >= ^ MO < /(« - 0 > d£ J — o o The integral can be simply evaluated if the forcing is stationary, because then < f[t—f) >= y,j is a constant independent of £. Then the mean response is HX = vfH{0) Review of Linear Spectral Theory / A3.0 125 where the relationship between h(t) and H(il) in equation A3.5 and the specific form of H(Q) for equation A3.1 have been used. The response variance, or the response mean square value, is found by considering the Fourier transform relationship between the cor-relation function and the spectral density. Rx(T)= r Sx(n)ein* dfl J —OO Setting T = 0 gives the mean square value Rx{0) =< x2(t) >=^= f sx(n)dn so that the mean square value is simply the integral of the spectrum over the entire frequency range. Therefore by knowing the response spectrum, the probability density of the response is then known completely. Some other important relationships for derivatives of stationary processes can be derived using the stationarity property in equation A3.7. Specifically by differentiating equation A3.7 it can be shown that (Clough and Penzien 1975) D2Rx[r) „ , , A d*Rx(r) n , N where Rxi{r) is the correlation function of the derivative process x' for example. Also since Rx(T)= r Sx(Cl)eiQr dCl J —OO /oo sx,{nynT dn -oo RAT)= f sx,,(n)einT dn J —CO it is easy to show that sx,(Q) = tfsx(n) 5 ± « ( 0 ) = n2sx,(n) = n45x(fi) which are very useful relationships between the spectral density of a stationary random process and the spectral density of the derivative processes. Review of Linear Spectral Theory / A3.0 126 Often the spectrum of the sum of two or more stationary processes is re-quired. For example consider q — x + y with the generalization to the sum of more terms being straightforward. Then the correlation function of q is Rq(r) =< q(t)q(t + r) > =<[x(t) + y(t)}[x(t + T) + y(t + T)}> =< x{t)x{t + T)> + < x(t)y(t + r) > + < y{t)x(t + r) > + < y{t)y{t + r) > = RX(T) + Ry(T) + RXY(T) + Ryx(T) Taking the Fourier transform gives the spectrum of q Sq(tl) = SX(Q) + Sy(Q) + Sxy(Q) + syx(n) For both the method of equivalent linearization and the perturbation method described in Chapter 3, the absolute response is given by x = u — r and knowing the equations for u and r and how they are related enables the response spectrum to be found. The details of this calculation will be given here. Based on the above discussion for the spectrum of a sum, the response spectrum is sx(Q) = sa(Q) + sr(n) - (5or(n) + sru(ft)) (AS.Q) The relative response spectrum is related to the wave displacement spectrum by Now it is required to find the relationship between the cross spectral terms and the wave displacement spectrum. These spectra are found in the usual way by taking the Fourier transform of the correlation function. The relationship in the time domain between r and u for the method of equivalent linearization is r" + 2 /oo poo / M £ i ) M & W ' + £ i -- 0 0 •» — 0 0 A similar expression can be obtained for RUT{T). The Fourier transform of Rru(r) + RUR{.T) gives, after some manipulation 5ru(n) + sur(n) = 2[(i - n2)(i - (1 - 7 ) n 2 ) + 4 f e / / ir„fi 2] |jy r(n)|2s t t(n) (^3.11) where the transfer function is given by equation 3.20. Substitution of equations A3.10 and A3.11 into A3.9 gives the response spectrum as 1 1 5,(0) + \Hr(n)\* |fr.(n)p - 2{(i - n2)(i - (I - 7)n 2) + 4 f e / / ^ > ] |frr(n)|25.(n) Now with the definitions of the transfer functions from equation 3.20, the term in the [ ] simplifies to (£,yy — &v) 2- Finally, using the definition of the effective damping in Chapter 3 for the method of equivalent linearization, the response spectrum is 5 x (n) = { 7 2 n 2 + ^ i } n 2 | ^ r ( n ) | 2 5 t t ( n ) where aj is given by equation 3.9 with au> being replaced by oy. By including the possibility of an offset the response spectrum is given by equation 3.22. APPENDIX A4 Hermite Polynomials Hermite polynomials Hek(z), possess the important property that they are orthogonal with respect to a zero mean unit variance Gaussian probability density weight-ing function. That is /OO j ^ -j=e-> '2He,iz)Hek{z) dz = *!6> where 8jk = 1 only if j = k and is zero otherwise. Because of this property it is often useful to use Hermite polynomial expansions for Gaussian random variables or for functions of Gaussian random variables. Specifically it is very easy to evaluate expected values of Hermite polynomial expressions of Gaussian variables due to the orthogonality relationship as will be shown in this Appendix. In what follows some useful properties of Hermite polynomials are discribed. Abramowitz and Stegun (1964) define the Hermite polynomials Hek(z) by the Rodrigues formula Hek{z) = ( - l ) V 2 / 2 ^ e - * 2 / 2 fc=l,2,... The lowest few Hermite polynomials are He0(z) = 1 He^z) = z He2(z) = z2 - 1 Hez{z) = z3 -Zz He^z) = z* - 6z2 + 3 Hes(z) = zs - 10z 3 + 15z 128 Hermite Polynomials / A4.0 129 The expected value of a function of a zero mean unit variance Gaussian random variable z is given by (/(*)) = r -}=e->2/2f(z)dz J-oo V2TC which can be easily evaluated using the orthogonality relationship if f(z) is expanded in an Hermite polynomial series. That is, if / ( * ) = = ! > * * « * ( * ) (A4.1) it=o then (/(*)) = J > f ° ~!i=e-z2/2Hek(z)He0(z) dz = c 0 That is the expected value of f(z) is simply given by the first coefficient c 0 of the Hermite polynomial expansion of that function. The coefficients ck can be found in a number of equivalent ways. Three different methods are described here. 1. Orthogonality Multiplying both sides of equation A4.1 by e~z2/2Hej(z)/y/2n and integrat-ing over the range of z gives /(*)—=-HeJ{z)dz = J2c>< / ff«*(*)-7=-ffey(*)«fe -co V 2 T T ^TQ • ' - C O V2TT Using the orthogonality relationship gives 1 f 0 0 2 / ck = -j=j_J{z)e-> l2Hek{z)dz 2. Least square minimization The coefficients ck can also be found by minimizing the expected value of the squared error of the approximation. That is minimizing with respect to the coefficients gives 9 ((f(z)-^2ckHek(z))2)=0 k=Q dck Hermite Polynomials / A 4.0 130 which can easily be shown to yield the same expression for ck as given above. 3. Weighted residuals The residual R of the expansion process is R = f(z)-Y^ckHek(z) k=Q The coefficients ck can be found by setting the expected value of the weighted residual to zero. That is where Wj = Wj(z) are the weighting functions. If the weighting functions are themselves Hermite polynomials, Wj{z) = Htj{z), then the same expression for the coefficients ck will again result. Often the functions of interest are polynomial functions of z and it is neces-sary to evaluate integrals of the form in order to find the coefficients ck. Some of the first few integrals of this type are indicated here. First for n = 0 which follows directly by substituting Rodrigues formula into the integral. When k = 0 then HeQ(z) — 1 and the integral to be evaluated is This integral cannot be evaluated as an indefinite integral but by setting the limits of integration between 0 and a yields the error function Hermite Polynomials / A 4.0 131 For n = 1 J ze-z^2Hek{z) dz = - (zHek^{z) + Hek_2(z)) e~z^2 k^0,l If k = 0, then j zt-^l2 dz = -t->*'2 and if A; = 1 j * z2e~z*l2 dz = - c*e - a 2 /2 + y | e r f (a/2) 'o Finally, for n = 2 y z2e>2/2Hek(z) dz = - (z2Hek_1(z) + 2zHek_2(z) + 2Hek_z{z)) t^l2 k # 0,1,2 The case of k = 0 corresponds to n = 1, k — 1 given above, and if k = 1 then If * = 2 / z V * I2 dz = -a(a2 + 2)e~a 2/2 + v'Sirerf (a/2) Jo '0Generalizations for arbitrary n could obviously be made if desired. Often a limit of integration is z = 0 so that an expression for Hek(0) is required. The useful identity [ 0 k Hek(0) = I < - I ) * " M K A: even odd can be deduced from Abramowitz and Stegun (1964). Following Eatock-Taylor and Rajagopalan (1982) the expected values of the product of two functions of z can be found using Hermite polynomial expansions. For example if / ( * ( « i ) ) = E « y ^ ( « 3 y(*(*2)) = E 6 * i r e *W2)) Hermite Polynomials / A4.0 132 then the correlation function i2/ f f(£i,* 2) 1 3 In order to evaluate this expected value the joint probability density of z{t{) and z(t2) is required. For z Gaussian, L in (1976) has shown that the joint probability density is P ( z i > * 2 ) = ^ e x p z\ + z\ 2 where zx = z[tx), z 2 = z[h) an^ P = < >• Using the orthogonality relationship for Hermite polynomials the desired expected value is /C O i»CO / P&U 2 2 ) / ( Z l M * 2 ) dzxdZ2 - c o J — C O = £ a r W r The extension to finding correlation functions with more terms is accomplished in a straightforward fashion. For example, in addition to the functions f{z{) and 17(22)> consider the function h{zz) defined by l where 23 = z(t$). It is desired to find the correlation function Rf9h(h,t2,h) =< f{z{h))g{z{t2))h{z{t3)) > In order to evaluate this expected value the trivariate probability density p (2 1 ,2 2 ,2 3 ) is required and is p{zi, z2, zz) = (2^3/2 e x P ( ~ \ & + z% + E E E ( ^ ) r ( ^ « n ^ ) ' ( i 4 4- 2 ) r » t ^ e r + , ( 2 1 ) i Y e r + t ( 2 2 ) i 7 e , + < ( 2 3 ) where =< 2(^)2(^2) > and similarly for the other correlation terms. The desired correlation function is then given by /CO />O0 /• OO / / V{zuz2,zz)fi where m and z 4 = z[t\). The fourth order joint probability density function is required to evaluate the expected values and is given by 1 1 p(z1,z2,z+n+ivdt+v+rc r $ t « v t c (r + s + t)\(r + u + v)l(s + u + w)\(t + v + w)\ r\s\t\u\v\w\ (KZi Z 2 ) ^ KZl Z 3 ) ' {Kti Zt Y (KZ2Z3 )" (-ff*a zt Y (KZ3 Zi V The final part of this Appendix shows that the sum of the correlation functions Rgig{r) + R9ZI{T) is identically zero if the function g{z(t)) can be expanded in an Hermite polynomial series in z(t). In detail this sum is Rz'g{r) + RgA*) =< + T)> + < g{t)z\t + r) > d i M A ) = 5; + R9z(r)} Hermite Polynomials / A4.0 134 where the stationarity property described in Appendix A3 has been used. Assuming an Hermite polynomial series expansion for g{z(t)) as *(*(0) = 2>*ff**(*(0) and recalling that z = Hex(z) then Rzg(r) =< Hex{z{t))Y,ahHek{z{t + T)) > k Rgz{r) - < X ) HHe^zit^HeMt + r)) > k Both of these expected values are easily seen to be equal to ax by using the orthogonality property of the Hermite polynomials. Therefore equation A4A is indeed equal to zero as stated at the outset. APPENDIX A5 Convolution Spectrum Calculations It is instructive to illustrate how convolution spectra of any order can be obtained from a given spectral representation. The procedure will be illustrated for a simple harmonic wave for which convolution spectra of orders 2 and 3 will be computed by direct integration. For more complex spectral shapes such as the P - M spectrum the same procedure can be used but the Fourier transforms are efficiently computed by the fast Fourier transform technique. The two sided spectral representation of an harmonic wave with frequency Cl0 and mean square value a 2 is given by 2 2 Sr}(Q) = j6(Q-ClQ) + j5(n + Cl0) (A5.1) where 6( ) is the Dirac delta function. The procedure to find the convolution spectrum [