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'/<v a = ratio of current to harmonic velocity = V/(Umw) ratio of current to standard deviation of fluid velocity = V/crtti p = ratio of current to standard deviation of relative velocity = V / a y angle 1 = constant 6 = constant 5(fl) = Dirac delta function Sij- = Kronecker delta At = time increment A n = frequency increment s = nonlinearity parameter, perturbation parameter e = perturbation parameter £ = damping ratio TJ = water surface elevation fi , = mean value v = viscosity of water, index p = density of water ps = density of steel <j = standard deviation r = correlation time interval <j> = 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~ <PU* 1 n „ dU" fi + fD = p—c™-]pr + 2pDCd*!F rZ)2 d2x* 1 dx' dt* dx* (2.4) dt* Here the fluid force is assumed to be the sum of two independent flow situations, a far field flow and a near field flow. The force due to the far field flow is given by the first two terms on the r.h.s. of equation 2.4. These two terms correspond to the original form of the Morison equation. The force due to the near field flow, the last two terms on the r.h.s. of equation 2.4, are again given by the Morison equation, but now the structure's velocity dx*/dt* and acceleration d2x*/dt*2 are used. Again the wave force is a nonlinear function of the response. Laya et al. (1984) discuss the relative merits of each of these two formulations and qualitatively outline the conditions for which each formulation can be applied. For example, if the flow is characterized by a well defined wake then the flow can be considered to be quasi-steady and the relative velocity formulation will apply. This condition occurs for large values of K where the wave particle displacement is very much greater than the cylinder diameter. However, if the structure is rapidly oscillating in an otherwise slowly varying external flow then it may be more appropriate to consider two independent flow situations for which the independent flow fields formulation might apply. The first flow is the accelerating flow past a fixed structure, while the second is the rapid oscillation of a structure in an otherwise quiescent fluid. In general terms Laya et al. (1984) suggest that the relative velocity formulation applies for large K while the independent flow fields formulation applies for small K. It should be noted here that for K less than about 2 the effects of diffraction become important and must be considered. Determinbtic Wave Forcing / 2.1 15 A distinctive feature of the independent flow fields formulation is the presence of two inertia and two drag coefficients. The coefficients Cmv and Cjt are based on a Reynolds number found using the maximum wave velocity [dU'/dt*)mt and a Keulegan-Carpenter number found using the maximum wave particle displacement. However, the coefficients Cmx and Cjx are based on the Reynolds number now defined in terms of the structure's maximum velocity (dx*/dt*)m and on the Keulegan-Carpenter number now defined in terms of the maximum displacement of the structure. The advantage of having two more empirical coefficients in the independent flow fields formulation is that the formulation is more flexible. However, probably outweighing this is the necessity to try and determine these coefficients experimentally for a complete range of flow parameters. Note that in the limiting case of purely inertial flow, where the drag force terms vanish, an equivalence between equations 2.3 and 2.4 should occur with Cmv = Cm and Cmx = Cm — 1. This relationship gives some insight into the nature of these inertia coefficients but it is not known whether this relationship will still be valid for the general situation where drag forces are also important. It should also be noted that at present there is no theoretical basis for assuming that either formulation to modify the Morison equation is appropriate in any given situation. In addition, experimental verification is not presently available. Inclusion of the effects of wave-structure interaction by either formulation introduces a nonlinearity in the structure's response velocity. Physically, the principal effect of this interaction is to reduce the structure's response and therefore the phrase hydrodynamic damping is often used to describe this phenomenon. The magnitude of the overall damping, which is the sum of the structural damping and the hydrodynamic damping, can significantly influence the response and thus the accurate estimation of the damping is of paramount importance. Specifically, it is important not to overestimate the damping for reasons of safety and not to underestimate the damping for reasons of economy. The resulting equation using either formulation can be described as a non-linear ordinary differential equation with the the forcing specified by the wave particle Deterministic Wave Forcing / 2.1 16 kinematics. Even for the simplest forms for the forcing, no analytical solution in terms of elementary functions is known to exist, and either numerical or approximate methods must be used for a solution. Of the approximate methods, linearization and perturbation are often useful, especially when the effect of the nonlinearity is small. Numerical solutions are used to check the accuracy of the approximate methods. For the deterministic waves to be considered in this chapter the wave velocity is given by ^ = V-Umudnv? (2.5) where V is a constant current, Umu is the amplitude of the harmonic velocity term, u is the wave frequency and Um is the maximum wave particle displacement. The case of simple harmonic flow is obtained by setting V =0. Equations 2.1 through 2.5 are in dimensional form and their dimensionless counterparts can be obtained by using the substitutions m = m o + P—(Cm-l) tor-j^ *' rr V* 4 N OJ ( 2 ' 6 ) x = — U = — t = uNt* ft = — In addition it is useful to define the relative displacement between the fluid and the struc-ture as r = U — x. Substitution of equation 2.6 into equation 2.1 and using equation 2.3 results in a dimensionless equation for the relative velocity formulation which is x"+ 2;Nx' + x = iU"+ e{U'-x')\U'-x'\ (2.7) where the superscript' means differentiation with respect to t} 7 = (pD2/m)(ir/4)Cm and e = {pD2/m)C<i/2. In terms of the relative velocity r the dimensionless equation is r" + 2(Nr' + cr'\r'\ + r = (1 - i)U" + 2$NU' + U (2.8) The corresponding equation for the independent flow fields formulation is, using equation 2.4 x" + 2$Nx' + exx'\x'\ + x = inU" + e9U'\U'\ (2.9) Deterministic Wave Forcing / 2.2 17 where 7 , = {pD2/m){ir/4)Cmn, e% = {pD2/m)Cdj2 and eg = {pD2/m)Cdx/2. Note that in the independent flow fields formulation the added mass contains the inertia coefficient Cmx . The parameter e or ex represents the importance of the nonlinearity and will be used as the perturbation parameter in the perturbation solution. The wave particle displacement U, needed in the relative velocity formulation is given in dimensionless form as UmVl Um U = -^-at+-^-costlt (2.10) where a = V/(Umu) is the ratio of the constant current to the harmonic velocity am-plitude. The wave particle velocity and acceleration are easily found by differentiating equation 2.10. Simple harmonic flow is recovered when a = 0 . Note the similarity be-tween equations 2.8 and 2.9 in that the forcing term in each is a known function of the wave kinematics while the nonlinear term on the left side of each equation is of the same functional form. Therefore a solution method applicable to one equation should be equally applicable to the other. The solution of the nonlinear equations 2.7 - 2.9 can be obtained by numerical integration, linearization and perturbation methods. However, before pro-ceeding to describe the solution methods in detail it is instructive to look at the typical magnitudes of the parameters that appear in the governing equations. 2.2 P A R A M E T E R V A L U E S 2.2.1 Damping The estimation of damping coefficients has been and continues to be a very difficult task. There are many sources of damping in offshore structures; the main ones are structural damping due to internal friction, material damping, foundation damping and fluid or hydrodynamic damping. As will be seen later, the estimates of hydrodynamic damping depend on the solution method that is used. Typically all damping other than the hydrodynamic damping is lumped into an equivalent linear viscous dashpot with a damping ratio between 1% and 5% of critical damping. This will be assumed here with a Deterministic Wave Forcing / 2.2 18 damping ratio fjy = 0.02 being used in all cases. The addition of hydrodynamic damping will result in an effective damping ratio. 2.2.2 Inertia and drag coefficients A comprehensive review of the parameters which influence the inertia and drag coefficients in an oscillatory flow is given by Sarpkaya and Isaacson (1981). Repro-duced in Figures 2.2 and 2.3 are the values of Cm and Cd as a function of Reynolds number Re and Keulegan-Carpenter number K from the above reference. These have been found by extensive laboratory U-tube experiments. As can be seen these two coefficients vary from 0.6 to 2-0 over a wide range of Re and K. These curves are for smooth cylinders; the effects of roughness are not considered here but again reference to Sarpkaya and Isaacson (1981) may be made for more details. Note that for a given K there is an inverse relationship between the inertia and drag coefficient as a function of Re . That is, for constant K and low Re, Cm is relatively small while Cd is relatively large. For constant K and large Re, the opposite is true; Cm is relatively large and Cd is relatively small. At intermediate Re and constant K the values of inertia and drag coefficients are approximately the same. Blevins (1977) uses typical values of Cm = 1.25 and Cd = 1.45 and these values are used here. The above drag and inertia coefficients apply strictly to the case of harmonic flow only. With the presence of a constant current the inertia and drag coefficients no longer retain their harmonic values and Cm and Cd will depend on the velocity ratio a in addition to Re and K. When the constant current becomes larger than the harmonic velocity, a > 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<NU' + U-a Using the wave displacement given by equation 2.10 it is easy to show that the relative velocity is of the form r' = B{8 + s in(m + <f>)} 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 " = <N + E 37 Once the effective damping is known the solution can be found by the Fourier series expansion described in Appendix A l . The notation RVL4 will be used for the method of equivalent linearization. 2.3.3 Perturbation Method Often perturbation methods are useful to find approximate solutions to non-linear equations when the effect of the nonlinearity is small. Books by Nayfeh (1973) and Kevorkian and Cole (1981), for example, review many different types of perturbation methods that have been successful for such problems. For the steady state solution of harmonically forced oscillators, Lindstedt's method, which is a multiple scale pertubation method, is an effective solution technique. Hooft (1982) has adapted Lindstedt's method to consider nonlinearly damped harmonically forced oscillators and this method is described below for the relative velocity equation r" + 2(Nr' + er' \r'\ + r = c 0 cos ftt + dQ sin tit + e0t + g0 (2.15) Here r* |r ' | is the nonlinear damping function and the forcing terms are from equations 2.8 and 2.10 with the coefficients c 0 , d0> 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 + (<f2i + 2 f 2 a 0 n) 2 ] = o The desired value of ft is _ c 2 1 6 0 — rf^ ap ^ ~ 2n(ag + 6j) The effective damping to this order of approximation is now given by the root of ?jv = ft + ^i(Co) + *2ft(ft) where again the explicit dependence of ft and ft on ft is noted. Once ft has been found, the solution of equation 2.17c can be found in the standard way and is r2(t) = a2o + ^2 a2n cos nQt + b2n sin nClt n = l Note that in this case the sum starts from n = 1 as the fundamental frequency terms are not completely suppressed. The solution to 0(e2) is thus given by x(t) = u - r0(t) - erx{t) - e2r2{t) As will be seen in the next section an improved solution results by including terms of 0(e2). This solution will be applicable for small values of the perturbation parameter s, the exact limitations being assessed by comparison of the perturbation solution to the numerical integration solution. 2.4 RESULTS Now that numerous solution procedures have been described, the equation of motion can be solved for a range of parameter values by each solution method. The following parameter values are used: Cm = 1.25, Cj = 1.45, ftv = 0.02 and UmCl/D = 2.55. Four values of the nonlinearity parameter e = 0.05,0.15,0.25 and 0.50 are considered for the Deterministic Wave Forcing / 2.4 33 Table 2.1 Solution methods for deterministic forcing. N O T A T I O N S O L U T I O N D E S C R I P T I O N E X A C T Numerical integration RVL1 Zero drag damping RVL2 Time dependent drag damping R V L 3 Average drag damping RVL4 Equivalent linearization PERT1 Perturbation O(s) case of zero current, a = 0. The solutions are presented in terms of the maximum dynamic response divided by the static response as a function of the dimensionless frequency f). Table 2.1 is a useful reminder of the notation used to identify the different solution methods discussed above. Figures 2.4 through 2.7 show the solutions for the cases e = 0.05,0.15,0.25 and 0.50 respectively. The conclusions regarding these results are summarized in point form as follows: 1. Even for very small nonlinearity, the zero drag damping linearization R V L l , does not adequately represent the response at the resonance peaks at all. This method greatly overpredicts the response which is not surprising for a method with no added hydrodynamic damping. 2. The effect of the nonlinearity firstly results in an amplified response at a frequency of CI = 1/3 with a peak in the response curve occuring there. This effect is due to the quadratic nature of the drag force. Peaks in the response curves also occur at frequencies 0 = 1/5,1/7,. . . , but are insignificant due to the damping. The peak at f) = 1/5 does show up in the case of zero drag damping. These multiple peaks can be explained by the Fourier series ex-pansion of the nonlinearity which only has odd components when a = 0. As Deterministic Wave Forcing / 2.4 34 t — i r 0.0 02 0.4 0.6 0.8 1.0 12 1.4 1.6 1.8 2.0 DIMENSIONLESS FREQUENCY F I G U R E 2.4 Frequency response for a = 0 and e — 0.05. r 0.0 02 0.4 0.6 0.8 1.0 \2 1.4 1.6 1.8 DIMENSIONLESS FREQUENCY zo F I G U R E 2.5 Frequency response for a = 0 and e = 0.15. Deterministic Wave Forcing / 2.4 35 0.0 0.2 0.4 0.6 0.8 1.0 12 1.4 1.6 1.8 DIMENSIONLESS FREQUENCY F I G U R E 2.8 Frequency response for a = 0 and e = 0.25. EXACT RVL1 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 DIMENSIONLESS FREQUENCY F I G U R E 2.7 Frequency response for a = 0 and e = 0.50. Deterministic Wave Forcing / 2.4 36 will be seen later for the case of nonzero current additional resonant peaks will be present at the even fractions 1/2,1/4,. . . with the peak at = 1/2 being the most significant. The second effect of the nonlinearity is to distort the response curves near the primary resonance fi = 1.0. For larger nonlin-earity the distortion becomes more pronounced but due to the damping does not result in a multivalued response curve. Obviously increasing s increases the hydrodynamic damping resulting in a decrease in the dynamic/static re-sponse. The absolute response can be obtained from these figures using the static response computed at 0 = 0.05 which are given in Table 2.2. Table 2.2 Static response for various nonlinearity and a = 0. e Static response 0.05 0.327 0.15 0.979 0.25 1.627 0.50 3.217 3. There is only a small difference between the time dependent drag damping solution RVL2, and the average drag damping solution R V L 3 , the differ-ence increasing with increasing e. They both give accurate results for fre-quencies less than about ft = 0.6 and greater than about ft = 1.5 for all values of e except e = 0.5. Generally both methods underpredict the re-sponse at ft = 1.0 and the time dependent drag damping is slightly better than the more approximate average drag damping method. As previously stated, the approximation involved for these two methods is that the maxi-mum response velocity is much less than the maximum wave velocity, that is {dx/dt)max/(dU/dt)max <S 1 must be satisfied if either method is to yield ac-curate results. With UmCl/D = 2.55 this condition becomes X m f t / 2 .55 -C 1, Deterministic Wave Forcing / 2.4 37 where Xm is the maximum dynamic response amplitude. This condition is satisfied at low frequency and high frequency as stated above. Two examples will illustrate this relationship. Consider first the case where £ = 0.05 and ft = 1.0. The maximum dynamic response for the average drag damping case is Xm = (4.954)(0.327) =• 1.62 so that the above condition becomes (X m f t /2.55) = 0.64 meaning that the linearization condition is not satisfied. The solution from the numerical integration is 2.24 so that the above con-dition accounts for the average drag damping solution being 28% low. The second example is for e = 0.15 and ft = 1/3. The dynamic reponse using the average drag damping method is 1.17 and the linearization condition becomes X m f t /2 .55 = 0.14. The numerical solution is 1.09 showing that when the linearization condition is satisfied, a more accurate solution results. These examples show that care must be taken when using the linearization assump-tion of either RVL2 or RVL3 and that the condition (dx/dt)m/(dU/dt)m « : 1 should always be checked. Based on the computations performed here an er-ror of less than about 5% results if (dx/dt)m/(dU/dt)m < 0.1. 4. There is no visible difference in Figures 2.4 and 2.5 the method of equiva-lent linearization, RVL4, and the numerical solution for frequencies greater than about ft = 0.6. In addition, it can be seen from Figures 2.6 and 2.7 that the method of equivalent linearization gives very accurate results for frequencies greater than about ft = 0.6 regardless of the size of the nonlin-earity. This feature of the method of equivalent linearization has also been found in the random forcing CcLSCj 3»S will be seen later, making the method very powerful for the solution of nonlinear equations. However, due to the way the linearization is achieved, linearizing both the nonlinearity in dU/dt and in dx/dt, this method is unable to capture the nonlinear behavior at the secondary resonance peaks. Also, the method of equivalent linearization underestimates the response at the lower frequencies. Deterministic Wave Forcing / 2.4 38 5. The effect of increased damping with increasing e is analogous to that found in linear systems. The response curve near ft = 1.0 is broadened and de-creases in magnitude with increasing e. The maximum value occurs at ft = 1 for light damping but as the damping increases the maximum value shifts to lower frequencies. For example when e = 0.5 the maximum value is the static response. Recall again that e represents the magnitude of the nonlinearity and that increasing s corresponds to increasing the hydrodynamic damping. 6. Overall the perturbation method, here taken only to 0(e), gives the most accurate solution for all values of frequency and nonlinearity. Similar to the method of equivalent linearization, the perturbation method is indistinguish-able from the numerical solution for ft > 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! = <ru,{a + z) (3.1) where u' is a zero mean Gaussian random process, a = V/ani with being the standard deviation of u' and z = u'/aui is a zero mean unit variance Gaussian random process. A random sea state can be characterized by the wave height spectrum Sv(Cl) which indicates the frequency distribution of the wave energy. As linear wave theory is assumed here, a linear relationship exists between the surface elevation and the water particle kinematics. Linear spectral theory can then be used to obtain the velocity and acceleration spectra from the wave height spectrum. These spectra can then be used to find the spectrum of the force and ultimately the response spectrum. Appendix A2 reviews linear wave theory and derives the relationships between the wave height spectrum and the velocity and acceleration spectra. Specifically the relationship between the surface elevation n and the water particle displacement U at location s is cosh ks sinh kdP 49 Random Wave Forcing / 3.0 50 where k is the wave number, d is the water depth and s is the distance above the seabed. Based on this linear relationship the following spectral relationships can be derived which will be useful in what follows Su(Q) = \HHn(U)\2Sn(n) 5,,(o) = n 2 5 u ( n ) (3.2) 5a»(fi) = n25'u<(fi) = n 4 5 8 (f i) where the transfer function is i r r ,^s\2 cosh2ks For the location s = d, corresponding to the the mean water surface, and for large water depths for which tanh kd « 1, the wave height spectrum is equal to the wave displacement spectrum. Measurements of wave heights in the ocean have led to numerous analytical forms for the wave height spectrum; a popular one being the Pierson-Moskowitz (P-M) spectrum. The specific functional form of this spectrum is given in Appendix A2 along with some of its important properties. It should be mentioned at this point that the wave height spectrum can change in the presence of a current. The effect of the current is mainly to cause a scale change in the spectral amplitude but because consideration is given to arbitrary spectral amplitudes, the effects of current on the spectrum are not considered here. Further details of the effect of currents on the wave spectum can be found in Tung and Huang (1976). In what follows, equations 2.7 and 2.8 are solved for the random wave case using various linearization techniques, the perturbation method and numerical simulation. Two general linearization methods will be described, the first based on the linearized version of equation 2.7 which neglects the quadratic velocity term while the second is the method of equivalent linearization applied to equation 2.8. Both Eatock-Taylor and Rajagopalan (1982) and Gudmestad and Connor (1983) consider the linear version of Random Wave Forcing / 3.1 51 equation 2.7 while Spanos and Chen (1981) have considered the method of equivalent linearization applied to equation 2.8. Excepting numerical simulation, the general approach that is taken here is to replace the governing non-linear equation by one or more linear equations which are a suitable approximation to the original equation. With linear equations, the powerful spectral analysis methods can be used to obtain the response spectrum and the response mean square value which are of interest. A review of the methods used in linear spectral analysis is given in Appendix A3. Ideally the response probability density is the statistic of interest and some of the techniques associated with finding the probability density will also be discussed. For the case of deterministic waves attention was focused on the steady state solution, with the transient part assumed to be unimportant. A similar argument is taken for the random case and the sea state is assumed to be stationary and only stationary response statistics are computed. 3.1 TIME DEPENDENT DAMPING LINEARIZATION The first method to be discussed is based on the linearized form of equation 2.7 which neglects the the quadratic velocity term. As in the deterministic case this results in a time dependent drag damping as is evident from the form of the linearized equation x" + 2$Nx' + x = iU" + eU' \U'\ - 2e\V\x' where the last term on the r.h.s. represents the time dependent hydrodynamic drag damping. A perturbation solution of this equation is considered by Eatock-Taylor and Rajagopalan (1982) and is summarized here. Gudmestad and Connor (1983) use the ap-proximation of a constant drag damping which, as will be shown, represents the first term of this perturbation method. Random Wave Forcing / 3.1 52 The analysis starts by rewriting the hydrodynamic damping term as 2e\U'\ x' = 2e\U'\ x' - 2s (\U'\) x' + 2e (\U'\) x' where (U1) denotes the expected value of U'. This substitution results in x" + 2(^v + e (\U'\))x' + x = nU" + eU'\U'\ - e2s{\U'\ - (\U'\))x' where the last term on the r.h.s now represents the time dependent drag damping. The parameter € is attached to this term for identification purposes. At the end of the analysis setting e = 1 recovers the desired solution. The effective damping in this case is easily seen to be given by f«// = fcr + «<|tf , |> (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'\=<ju,\a + z\ = J2bkHek(z) (3.5) Jt=o using equation 3.1. Appendix A4 details several methods that can be used to evaluate the coefficients bk. It is appropriate to use Hermite polynomials to represent functions of zero mean unit variance Gaussian random variables because of the orthogonality property of Hermite polynomials and because of the resulting easy evaluation of the expected values that are of interest. Again the details can be found in Appendix A4. The coefficients bk are found to be h = - / / e j t - 2 ( _ a ) e " a 2 / 2 ° ' 1 and fc0 = °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'\ = <r2, (a + *) |a + z\ = £akHek(z) (3.8) with the coefficients ak given by H = ^=Hek_3(-a)e-a2'2 k # 0,1,2 «o = 4 ( ( « 2 + l)erf {a/y/2) + ^ | a<T"2/2} a x = 2a 2 ,{aerf(a/v /2) + e ~ ° 2 / 2 } a 2 = o-2,erf(ar/\/2) The expected value (17' is again the coefficient OQ which is identical to the expression given by Spanos and Chen (1981) who approached the problem in a slightly different fashion. Now with this background information complete, the perturbation method described by Eatock-Taylor and Rajagopalan (1982) can be discussed. Although they consider a multi-degree of freedom system, the method is easily adapted to the single Random Wave Forcing / 3.1 54 degree of freedom system under consideration here. The problematic time dependent damping term is identified by the parameter e and a regular perturbation series using this parameter is assumed x = x0 + exl + e2x2 H This expansion reduces the linear time dependent coefficient equation to a recursive set of more tractable constant coefficient equations of the form 4 + Kff4 + *o = lU" + eU'\U'\ = f0(t) (3.10a) A + Kfft'i + * i = - 2 * {\U'\ - (\U'\)) x'0 = hit) (3.106) A + MeffA + *2 = -2s (\U'\ - (\U'\)) x\ = f2{t) (3.10c) It can be seen that the first of these equations respresents the constant drag damping method of Gudmestad and Connor (1983) with the forcing given by the Morison equation and the effective damping given by equation 3.7. Of interest is the response spectrum which is given as SX(V) = SX0(Q) + e{SX0Xl(Q) + SXIX0(Q)} + e2{sX0X2in) + sX2X0in) + sxl(Q)} + • • • Due to the linear nature of the above equations the response spectrum for each equation is simply related to the forcing spectrum. This also applies to the cross spectral terms, for example (n) = |# x(n)| 2s / o / l(n) where 1^(0)1 is the transfer function for the above set of linear equations l#,(n)l2 = ( l - f i 2 ) 2 + ( 2 ? e / / n ) 2 Taking terms to 0(e 2 ) the response spectrum is 5x(fi) = \HX(Q)\2{S/0(Q) + e (S / o / l(fi) + Sflfo(Q)) + e2 (3.11) Random Wave Forcing / 3.1 55 In order to evaluate this expression it is necessary to find the relationship between the various forcing spectra and the wave height spectrum. For example /0(<) is seen to depend on the wave particle kinematics by equation 3.10a. Therefore it is possible to find the spectrum Sf0(fl) in terms of the spectrum of the wave particle velocity which is itself simply related to the input wave height spectrum by equation 3.2. The method that will be used to find Sf0(Q) is to first find the correlation function Rf0{f) and then take the Fourier transform. The details of t h i 3 procedure can be found in Appendix A3. As discussed by Eatock-Taylor and Rajagopalan (1982) the cor-relation functions # / 0 / 1 ( r ) » - # / i / 0 ( r ) » - f y , o / 2 ( r ) a n c * ^/2/o( T) a r e difficult to evaluate exactly and by using the cumulant discard hypothesis they are assumed to be negligible. This assumption is used here, simplifying the response spectrum to sx(n) = \Hx(Q)\2{sf0(Q) + sfl(n)} where e has been set equal to one as required. Here Sf (Q) represents the contribution of the time dependent damping term. Neglect of this term corresponds to the constant drag damping method as described by Gudmestad and Connor (1983). Some detail is shown here for the steps required to compute Sf0(Cl), as the same procedure is followed in future situations. From equations 3.1 and 3.10(a), fo(t) can be written as /o(') = 7<V*'(0 + + *(<)) \<* + *(0I where the Hermite polynomial expansion of equation 3.8 has been used. Now the correla-tion function is */o(') = (/o(0/o(* + r)> = 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 <S/0(O) is usually only taken to include terms up to r = 3 (see for example, Borgman 1969). This will be assumed to apply here. For the case of zero current a0 = 0 and all even coefficients in the above summation are also equal to zero. In order to find S/,(0) a similar procedure is used. That is the Fourier transform of R/^T) is taken to find the spectrum. Again Hermite polynomials can be used to write fx{t) as f^t) =-2eJ2bkHek{z(t))x'Q{t) k=i using equation 3.5 for the expansion of \U'\. Note that the sum in this case starts from A; = 1 as the expected value < \U'\ > 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))) <z'0(<)4(' + 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{<y + W + M { 3 l 7 ) In order to perform the expected value operations in equation 3.17, the probability density of r1 needs to be known. A Gaussian probability density is usually assumed for r \ This is due to the fact that for a linear equation and Gaussian forcing, the response is also Gaus-sian. Atalik and Utku (1976) prove a useful identity valid under quite general conditions provided that r1 is Gaussian which is (/(r')r') _ / 4 f ( Q \ (r'2) \ dr1 I This allows equation 3.17 to be simplified to C = 2e{\V + r'\) The expected value of a function of r' is given by <<7(r') >= 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 <rr» has replaced crai. Also note that because C depends on tjft, an iterative procedure is usually required in order to find the solution as was the case for deterministic equivalent linearization. Random Wave Forcing / 3.2 60 The resulting linear equation is r" + 2$effr' + r = (1 - 7)11" + 2 ^ u ' + u (3.18) where the effective damping is Uff = to + «v{0«rf(0/v5) + \jl^/2} (3.19) When the constant current is zero the effective damping reduces to the well known form (see for example, Malhotra and Penzien 1970a) $eff = to + e These expressions for the effective damping are the same as those given by the previous linearization method in equation 3.7 with crr» replacing aai. Equation 3.18 is linear and is easily solved by the spectral methods described in Appendix A3. The response velocity spectrum is given by \H«W\' where ( i - n ' ) a + ( 2 f c / / i i ) 2 \HV(Q)\ 2 1 (3-20) ~ ( i - ( i - 7 ) n 2 ) 2 + (2tofi)2 and again 5 B(fi) is simply related to the wave height spectrum using equation 3.2. The mean square relative velocity is simply aj = f SAW The nature of the iterative process required here is now evident. Assuming an initial value for <TTI specifies the effective damping, and then the spectrum Sy(fl) can be found and integrated to give a new estimate of ari. The procedure is repeated until convergence in cv is achieved. Random Wave Forcing- / 3.2 61 The offset r can be found by ensuring that equation 3.16 is satisfied on average. Since both r and u are zero mean processes the offset is given by which again is evaluated assuming that r1 is a zero mean Gaussian random process. The required integrations give rise to Equation 3.9 gives a similar expression for a 0 that was found using the Hermite polynomial expansions with ari replacing <ra. For zero current, /? = 0, and the offset will also be equal to zero as expected. So far only the spectrum of the relative process has been considered, but the absolute response x is of primary interest. The absolute response is related to the relative response by where u and r are zero mean processes. Obviously the absolute mean response wil l simply r =-e{(V + r')\V + r'\) (3.21) x = u — r — r be (i) = - f Again for zero current, (x) = 0. SX(Q) = f 25(fl) + SH(Q) + Sr(Q) - {S^Q) + Srn(U)} (3.22) As is shown in Appendix A3 this expression can be simplified to (3.23) Random Wave Forcing / 3.3 62 where ax = 2cV(/?erf(/?/v /2) + \j\^'2) The similarity between equation 3.23 and the spectrum given in equation 3.14 is striking. They are the same except that ari appears in equation 3.23 instead of aui in equation 3.14. The absolute mean square response is then simply the integral of the response spectrum which can be efficiently evaluated. As is expected, the mean square value contains the term r 2 due to the mean value of (x) — —f. 3.3 PERTURBATION METHOD The perturbation method described in the previous chapter is extended here to random forcing. The method starts with equation 3.16 for the relative velocity, repeated here for convenience r" + 2$Nr' + e(V + r')\V + r'\ + r + f = (1 - 7 ) u " + 2;Nu' + u (3.16) The constant part of the relative response is obtained by averaging this equation , f=-£({V + r')\V + r'\) where the fact that both u and r are zero mean variables is used. This indicates that F is a term of magnitude 0(s) which is helpful in properly placing it in the perturbation equations. Again the perturbation expansion to be used is given by equation 2.16 r = r 0 + erx + e 2 r 2 H (2.16) to = iTo + «fi + *zfe + • • • resulting in a recursive set of linear equations: r'o + 2$0r'0 + r 0 = (1 - 7 ) u " + 2to" ' + « = /o(<) (3.24a) r'l + 2s0r'x + rx = -(V + r'0)\V + r'0\- 2$xr'0 - f/e = fx(t) (3.246) A + 2<To'2 + r2 = -2r'x\V + r'0\ - 2$xr\ - 2 ^ = f2(t) (3.24c) Random Wave Forcing / 3.3 63 Now consider the first of these equations, 3.24a. Since u is a zero mean Gaussian random variable and since equation 3.24a is linear, then both r 0 and Tq will be zero mean Gaussian random variables. It is advantageous for the subsequent development to define the zero mean unit variance random variable to as Also, functions of to are to be expanded in Hermite polynomial series. For example, since w — Hei(w(t)), the forcing function fQ(t) can be written as Mt) = <rr>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<v o f lfTe 1(w(t)) " £ (3.26) Random Wave Forcing / 3.3 64 where again /? = V/v^ • The coefficients can be found using the procedure described in Appendix A4 and are identical to the coefficients given in equation 3.5 with <xri replacing crut. Also note that since o 0 is simply the expected value of the nonlinear function it can be shown that a 0 = - - = a 2 , {{P + w)\P + v>\) £ 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 <JUI. Due to the similarity between the expressions for the response spectrum from the linear constant drag damping method and the method of equivalent linearization and between the cubic constant drag damping method and the perturbation method, the results predicted by each method should also be similar. This will be shown in the results section. for random forcing is achieved by simulating a representative time history of the forc-ing and then numerically integrating the equation of motion to obtain a representative response time history. This method, often called Monte Carlo simulation, has two signif-icant disadvantages over any of the other methods described above. First, a very large computational effort is required, often two to three orders of magnitude greater than the method of equivalent linearization (see for example Spanos and Chen 1981). Second, and SX(Q) = f26(fl) + Sn(Cl) + Sr(V) - {Sur(Cl) + Sru(Q)} (3.28) 3.4 NUMERICAL SIMULATION The brute force method of obtaining a solution to the equation of motion Random Wave Forcing / 3.4 67 probably more important, is that the simulation solution only applies for the specific equa-tion that is under consideration. If one or more of the parameters of the problem were to change then a new simulation procedure would be required to find a solution to the new problem. Therefore a parametric study over a range of parameter values would involve prohibitive computational costs in contrast to the ease with which a parametric study could be achieved with one of the approximate methods described above. Despite these disadvantages, numerical simulation solutions often are the only means available to obtain a solution. Also they are usually the only means of assessing the accuracy of an approximate method. Some of the approximate methods suffer from the disadvantage that they are only able to compute the moments of the response. Ideally what is required is the probability structure of the response from which all statistics of interest can be found. Because a representative response history is obtained from the simulation procedure, the response probability can be found directly. The equation that is to be solved by the simulation process is equation 2.7, repeated here for convenience x" + 2$Nx' + x = qU" + e{U' - x')\U' - x'\ (2.7) where all that is known about the forcing is that the wave height spectrum is given by one of the commonly used spectra such as the P - M wave height spectrum. Using linear wave theory the wave velocity spectrum is given by 5..(0) = \Htt,„(n)\2sv(n) as given in Appendix A2. What is required is a time history of the velocity which is compatible with the above velocity spectrum. Two methods that can be used to simulate a velocity history, and subsequently an acceleration history, are described here. The first method assumes that the velocity is given by the summation of a finite number of harmonic terms, each with a random phase angle. That is M u'(t) = YfAksm(Qkt + <f>k) k=l Random Wave Forcings / 3.4 68 The random phase angles <pk are uniformly distributed over the interval [0,2ir] and the Ak are spectral amplitudes given by Ak = V2AnS„, ( f i f c ) where Aft is a frequency increment and Qk are the frequencies at which the spectrum is sampled Clk = kACl k = l,2,...,M The velocity u' becomes Gaussian as M —> 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<ru'U' Random Wave Forcing / 3.5 85 while the equation for the method of equivalent linearization can be written in terms of the absolute response as x" + 2$effx' + x = ^U" + £\j^crr,U' One difference between these two equations is that the values of effective damping predicted by each method are different as given in Table 3.3. Additionally, the last term on the r.h.s of the first equation depends on <7„», while the corresponding term of the second equation depends on ari. A s r = u - i , the standard deviation of u will be greater than that of r and the forcing term in the first equation above is greater than the forcing term of the second equation. This can be seen from Figure 3.9 which shows the relative velocity spectrum and the wave velocity spectrum for Vw = 31.6 and s = 0.25. o o E-° s E-•—• 0 O «5 O M O O CM O d WAVE VELOCITY RELATIVE VELOCITY i 1 1 r— 0.0 02 0.4 0.8 0.8 F I G U R E 3.9 1.0 12 1.4 1.6 1.8 ^0 DIMENSIONLESS FREQUENCY Relative velocity and wave velocity spectra for =31.6 and s = 0.25. As the variance is the integral of the spectrum it is obvious that oTt < o~ui. The increased forcing results in an increased response which offsets the decrease in response Random Wave Forcing / 3.5 Table 3.4 Standard Deviation of Response for s = 0.05. 86 a 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 0.0 0.248 0.251 0.250 0.253 0.253 1.0 0.270 0.273 0.272 0.276 0.278 2.0 0.329 0.331 0.330 0.333 0.337 3.0 0.391 0.392 0.392 0.394 0.398 Vw = 21.8 a 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 0.0 0.395 0.414 0.399 0.419 0.419 1.0 0.490 0.530 0.494 0.534 0.547 2.0 0.688 0.715 0.690 0.721 0.745 3.0 0.890 0.902 0.890 0.909 0.933 1^ = 31.6 a 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 0.0 0.573 0.629 0.579 0.638 0.638 1.0 0.762 0.855 0.769 0.865 0.913 2.0 1.160 1.221 1.163 1.235 1.304 3.0 1.576 1.606 1.577 1.619 1.676 due to the larger effective damping. A similar comparison between the cubic constant drag damping case and the perturbation method can be made to illustrate the reason for the larger response from the cubic drag damping method. Now consider the effect of a constant current. Table 3.4 gives the standard deviation of the response for all three P - M spectra with e = 0.05 and four values of velocity ratio: a = 0.0, 1.0, 2.0 and 3.0. Observe that as the current ratio increases, the Random Wave Forcing- / 3.5 87 response increases as expected. For the case of a non zero current the time dependent drag damping method now gives a slightly larger value than does the cubic constant drag damping method. These are followed in order by the perturbation method, the linear constant drag damping method and the method of equivalent linearization. The largest difference between these methods occurs for a = 1.0 with the difference increasing with current ratio. Figures 3.10 and 3.11 show the response spectra predicted by each method for Vw = 21.8 and = 31.6 respectively. In both cases e = 0.05 and a = 1.0. The differences in the spectra are most evident in Figure 3.11 for Vw = 31.6. Again the perturbation method and the cubic damping method are very similar as are the two linear methods. The time dependent drag damping predicts slightly larger spectral amplitudes than does any of the other methods. More evident in Figure 3.11 is the low frequency contribution to the spectrum as given by the cubic and time dependent drag damping methods and the perturbation method. This is a result of the contribution of the *2 convolution spectrum, which has a maximum value at Q, = 0, and also from the *3 convolution spectrum, which is nonzero near fi = 0. No simulations were performed for the case of non zero current. However, Spanos and Chen (1981) have reported that there is " quite good agreement " between their simulation results and the method of equivalent linearization. As the perturbation method is a consistent improvement to the method of equivalent linearization, it is felt that the perturbation method again gives the best estimate of the true solution among all of the approximate methods described above. Finally, Table 3.5 shows the values of effective damping for the nonzero cur-rent examples. The effective damping increases as the current ratio increases as expected, with identical values of damping being predicted by the two different methods. Random Wave Forcing / 3.5 88 0.0 02 0.4 0.8 0.8 1.0 12 1.4 1.8 DIMENSIONLESS FREQUENCY F I G U R E 3.10 Response spectra for Vw = 21.8 and a = 1.0. 0.0 02 0.4 0.6 0.8 1.0 12 1.4 1.6 1.8 DIMENSIONLESS FREQUENCY F I G U R E 3.11 Response spectra for Vw = 31.6 and a = 1.0. Random Wave Forcing / 3.6 89 Table 3.5 Effective Damping for e = 0.05. Vw = 11.6 a E Q L I N + P E R T L I N E A R + C U B I C 0.0 0.048 0.050 1.0 0.062 0.064 2.0 0.095 0.095 3.0 0.132 0.132 Vw = 21.8 a E Q L I N + P E R T L I N E A R + C U B I C 0.0 0.075 0.076 1.0 0.101 0.102 2.0 0.162 0.162 3.0 0.231 0.231 Vw = 31.6 a E Q L I N + P E R T L I N E A R + C U B I C 0.0 0.102 0.104 1.0 0.141 0.142 2.0 0.231 0.231 3.0 0.334 0.334 3.6 PROBABILITY DENSITY OF T H E RESPONSE The preceding sections have focused on approximate methods that can be used to determine the response spectrum and response variance for the problem of wave forces on a single degree of freedom model of an offshore structure. Specifically, the non-linear equation of motion was approximated by one or more linear equations and then spectral analysis was employed to obtain a solution of the second order response statis-Random Wave Forcing' / 3.6 90 tics. As stated previously the response spectrum is sufficient to completely characterize the response probability density for a linear system subjected to Gaussian forcing. The method of equivalent linearization reduces the original nonlinear equation to one with Gaussian forcing and a complete solution can be found. For the constant drag damping method and the perturbation method, at least one of the equations contains non-Gaussian forcing, meaning the response will also be non-Gaussian. Since any probability density is completely specified by it's moments, computation of all of the response moments will completely characterize the response probability density. For the case of an equation with a small nonlinearity the response should be only slightly non-Gaussian and an Edgeworth series expansion (see for example Kendall and Stuart 1977) can be used. The Edgeworth series of a random variable x is where Pg{x) 1 5 the Gaussian density, m 3 and m 4 are the third and fourth moments of the probability density, <r2 is the second moment or variance and Hen is an Hermite polynomial of order n. For practical use the Edgeworth series must be truncated at some point. Often only the terms shown in the above equation are used. The definitions of skewness A 3 = m 3 /cr 3 and kurtosis A 4 = ( m 4 / a 4 ) — 3 are used to measure the departure from a Gaussian distribution. Obviously skewness and kurtosis as defined above are zero for a Gaussian distribution. Non-zero skewness is a measure of the nonsymmetry of the distribution while non-zero kurtosis is a measure of central peakedness with a positive kurtosis being more peaked than a Gaussian distribution. A negative kurtosis has a flatter peak than a Gaussian distribution. A difficulty in using a truncated Edgeworth series occurs for certain ranges of values of skewness and kurtosis because equation 3.30 can predict unrealistic negative probabilities and spurious peaks in the probability density, especially in the tails of the distribution (see for example Draper and Tierney 1972). Random Wave Forcing J 3.6 91 As the higher order response moments characterize a non-Gaussian proba-bility density, a method is needed for their computation. It is useful at this point to recall the relationship between the response spectrum and the response variance. That is, the variance or second moment is the area under the spectral density al = Rx(0)= f Sx(fl)dCl In order to find higher order moments, higher order spectra can be defined and integrated in a similar fashion. For example the bispectrum Bx(Qx,Vl2) and a second order correlation function i2 I ( r 1 , r 2 ) can be defined by the two dimensional Fourier transform pair B x (n l l n 2 ) = (^52/ Jr Rx(Ti,r2)^p{-i(fiiTi + n 2 T 2 ) } < M r 2 R x { h , T 2 ) = f f B^MexpiiiCl^+^T^dCl^ where Rx{rUT2) =< x(t)x(t + Tl)x(t + r 2 ) > 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<Tf 2aj then the response probability density is P(x) = - ^ e x p ( - ^ ^ ) V2wax 2aZ That is, to obtain the response probability density function the variable x has replaced the variable / and the moments of x have replaced the moments of / . The same functional form is retained for the response probability. For the case of a non-Gaussian process a Random Wave Forcing / 3.6 93 similar result should also be approximately true. Tickell (1977) has hypothesized this for the case of the response probability density for the constant drag case and has verified this hypothesis based on experimental results. The hypothesis that the response probability density is of the same form as the forcing probability density with response moments being used instead of the forcing moments is assessed here using numerical simulation. The use of the above method requires the complete probability structure of the forcing. The constant drag damping method therefore is ideal for analysis by this technique because the forcing is given by the Morison equation and many authors have derived the appropriate probability density functions. Also the equation of motion is linear in the response variable. Due to the nonlinear dependence of the drag force on the velocity in the Morison equation, the forcing probability density is non-Gaussian. Therefore the response probability density will also be non-Gaussian but of unknown form. What is assessed in the next section is the hypothesis that for a linear equation, the probability density of the response is of the same form as the probability density of the forcing, even for the case of non-Gaussian forcing. 3.6.1 Constant Drag Damping In what follows only the case of zero current is considered. As such the mean response value is zero. In addition, the odd response moments are also zero, implying that the probability density of the response is symmetric with skewness equal to zero. Therefore the departure from a Gaussian distribution is indicated by a non zero kurtosis in this case. The equation for the fluid force based on the method of constant drag damp-ing is given by the Morison equation f = lU" + sU'\U'\ (3.32a) where U' and U" at any given time are independent and Gaussian so that the joint prob-ability density p(U', U") is known explicitly. By defining an auxiliary variable g = U' (3.326) Random Wave Forcing- / 3.6 94 the joint probability density p{f,g) can be found by the technique of transformation of variables as p{f,g) = p(u',u")\j\ where J is the Jacobian of the transformation of equations 3.32a and 3.326. The desired probability density p(f) is found by integrating out the auxiliary variable g. Pierson and Holmes (1965) have derived p(/) by the above method which results in an integral form for the probability density which must be evaluated numerically. Borgman (1967) has derived an alternative expression for the probability of the force in terms of the parabolic cylinder function. Using the defininitions ipx = ~\cru» and 0 2 = ecrl'> 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<J\< a s before. Recall that ^1 and ip2 a r e related to the second and fourth force moments by equation 3.33. Using the same hypothesis as outlined above the probability density of the response maxima is assumed to be of the form where xm is the response maxima and V*i and rp2 are given by equation 3.33 with the moments of x replacing the moments of / . The hypothesis that the probability density of the response maxima is given by the above equation can be checked using the results of the Monte Carlo numerical simulations. Figure 3.16 shows the probability density of the response maxima for e = 0.15, Vw = 31.6 and a = 0. A total of 106 time steps were used in the Monte Carlo simulation procedure. A larger response record is required in this case in order to generate a suffi-cient number of peaks for an accurate representation of the peak probability density. The asymptotic nature of the response peak probability is evident in Figure 3.16 as it does not match the simulation for small values of x but is reasonably accurate for larger values of x. Also, the approximate peak probability overestimates the simulated probability density for large response values. This is likely a result of the approximate peak probability density being based on p(x, x') rather than p(x,x',x"). In addition the approximate peak proba-bility is based on the cubic constant drag damping equation but the simulation procedure numerically integrates the original equation of motion, equation 2.7. Once p(xm) is known it is possible to determine the probability density of the largest of the response maxima Random Wave Forcing / 3.6 0.0 2.0 APPROXIMATE o SIMULATION 4.0 X 6.0 8.0 F I G U R E 3.16 Response peak probability density for a = 0,VW = 31.6 and e = 0.15. out of a large number of maxima by using order statistics (see for example Ochi 1978). It is the largest of the response maxima that is used in the design process. Finally, as was mentioned in the introduction, a totally different approach for finding the response probability density is to use a Markov vector approach. This involves solving a Fokker-Planck-Kolmogorov equation ( F P K E ) for the stationary joint probability density p(x,x'). For example, if the equation of motion is x" + g(x, x') = w(t) where g(x, x') is a nonlinear function and w(t) represents white noise, then the F P K E for Random Wave Forcing / 3.6 101 the stationary response is (see for example Caughey 1963) Here $ is the constant spectral density of the white noise. Known solutions of this F P K E exist for only a small class of nonlinear functions and approximate or numerical solutions must be used in general. Nigam (1983) reviews some of the methods that can he used in this regard. Obviously the use of the Markov vector approach is not directly applicable to the problem of wave forces on offshore structures because the wave spectra used as input forcing representations are not white noise processes. However the Markov vector approach can be used if an enlarged system of equations is considered. For example consider an equation of motion that is represented by £{*{*)) = M) where t is the operator representation of the equation of motion, and f{t) is the non-white noise forcing function. Now consider f(t) to be the response of a second equation with white noise forcing 0 {/(')> = 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 <f> is the velocity potential. The boundary conditions that apply in linear wave theory are *L = o z * = _<r {A2.2a) - J L - J L = ° **=0 (A2.26) d<f> + 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 <f>{x*,z*,t*) = <f>(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<f>, 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) =<x{t-r)x{t) >=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<re//r' + r=/(0 = ( l - 7 ) u " + 2^ vu' + u Review of Linear Spectral Theory / A3.0 127 where the time domain solutions for r and u are r(0= f °M£i ) / ( * -£ i ) ^ i J —oo /oo -oo Therefore the cross correlation can be written RRU(T) =< r(t)u(t + T) > /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<zl)g{z2)h{zz)dzldz2dzz (,44.3) - 0 0 J - c o J—00 Hermite Polynomials / A4.0 133 Sustituting equation A4.2 into equation A4.3 and using the orthogonality relationship repeatedly gives, after some manipulation it p (t * M - ^ W „ h , (r + s)!(r + Q!(3 + <)! r » t Carrying this analysis one step further results in the correlation function Rfghqihihih***) defined as Rf9hq(h,h,h,h) =< f(z(tl))9(z(h))h(z(t3))q(z(U)) > 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<hzt) = j^2exp{--{z2l+z2,+zl + zl)} YYYYYY—-— ^ rlsltlulvlwl r » t n v IB ( i r z l 3 2 r ( ^ 1 2 s ) ' ( ^ 2 Z J ) ' ( ^ 2 , 3 ) « ( ^ 2 Z 4 r ( i r 2 s , 4 ) « ' Following the same procedure as above results in the correlation to be given by Rfghqihtht hi U) = X ) X X X ) X X ar+«+tbr+tt+vCl>+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 [<S',J(fl)]*r starts by finding the correlation function RV{T) which is the Fourier transform of the spectrum. That is J —CO = a 2 cosft 0 r after equation A5.1 is substituted into the integral. Now raise the correlation function to the power r and take the Fourier transform to obtain the desired convolution spectrum. This procedure will be shown for r = 2 and r = 3 below. Therefore a 4 o 4 JRJ(T) = o 4 cos2 n 0 r = — + — cos 2 f i 0 r 135 Convolution Spectrum Calculations / A5.0 136 and , « , 3o 6 a 6 R*(T) = a 0 cos J U0T = — - c o s n 0 r H cos3ft 0r i 4 4 Note the nonlinear effect in the above calculation. Originally the correlation function contained only one frequency QQ but in raising the correlation to some power the correlation function and hence the energy content is distributed to other frequencies. Taking the Fourier transform then gives the convolution spectrum which in general is {5,(n)}< r = ±- I^R'n{T)e-^ dr and for the specific examples being considered here are {S„(n)Y2 = y6{Q) + j6(tl + 2fi 0) + ~6(Q - 2fi 0 ) and The mean square values of the convolved spectra are easily found by integration and for the case of r = 2 the mean square value is a 4 , while for the case of r = 3 the mean square value is a 6 . Figure A5.1 shows the relationship between the original spectrum and the *2 and *3 convolution spectra graphically. For the case of a more complex initial spectral shape, such as the Pierson-Moskowitz spectrum, the convolution spectra are obtained in an analogous fashion. How-ever, now the Fourier transforms must be obtained numerically and are efficiently evaluated using the fast Fourier transform algorithm (see for example Newlands 1975). Convolution Spectrum Calculations / A5.0 137 s,"(n) 5*3(n) n 1 ' Q -4n 0 - 2 n 0 o 2 n 0 4 n o F I G U R E A5.1 Convolution spectra, for harmonic wave.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Nonlinear response of structures in regular and random...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Nonlinear response of structures in regular and random waves Lipsett, Arthur William 1985
pdf
Page Metadata
Item Metadata
Title | Nonlinear response of structures in regular and random waves |
Creator |
Lipsett, Arthur William |
Publisher | University of British Columbia |
Date Issued | 1985 |
Description | The problem of the dynamics of a flexible offshore structure in either a regular 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 description. 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. |
Subject |
Offshore structures -- Hydrodynamics Waves |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-06-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0062575 |
URI | http://hdl.handle.net/2429/25826 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1985_A1 L56_7.pdf [ 6.73MB ]
- Metadata
- JSON: 831-1.0062575.json
- JSON-LD: 831-1.0062575-ld.json
- RDF/XML (Pretty): 831-1.0062575-rdf.xml
- RDF/JSON: 831-1.0062575-rdf.json
- Turtle: 831-1.0062575-turtle.txt
- N-Triples: 831-1.0062575-rdf-ntriples.txt
- Original Record: 831-1.0062575-source.json
- Full Text
- 831-1.0062575-fulltext.txt
- Citation
- 831-1.0062575.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062575/manifest