Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 Porto, Portugal, 30 June - 2 July 2014 A. Cunha, E. Caetano, P. Ribeiro, G. Müller (eds.) ISSN: 2311-9020; ISBN: 978-972-752-165-4 Theoretical model for estimating the soil particular velocity due to a heavy structure fall Alain Pecker1, Charles Fernandez2*, Mickael Metzger3 Geodynamique et Structure, 157 rue des Blains, Bagneux, 92220, France 2 GDF SUEZ/CRIGEN, 361 av. du President Wilson, 93211, Saint-Denis, France 3 Immeuble BORA, 6 rue Raoul Nordling, 92270, France Email: [email protected] (*corresponding author), [email protected], [email protected] 1 ABSTRACT: Over the past few years, extensive construction of wind turbines has taken place all around the world in areas where many steel pipelines are already buried in the ground. The possible fall of these heavy machines may induce damageable vibrations to the pipeline. Therefore a theoretical predictive model has been developed in order to assess the particle velocity near a buried pipeline, may a wind-turbine fall nearby. A fall of a 1.5 t mass was performed and velocities were measured in order to update the model. An experimental soil characterization through MASW tests was performed in a representative soil in order to obtain relevant input parameters for a non linear FE model. The FE model allowed updating the theoretical model for very heavy structures as wind-turbines. The updated model is now part of the RAMCES software which has been developed for more than a decade at CRIGEN and is widely used in France by pipeline transportation operators. KEY WORDS: Modeling; Soil particle velocity estimation; Vibrations; Heavy structure fall. 1 INTRODUCTION This model presents a simplified analytical model for estimating the peak particle velocity (PPV) at a depth into the ground and at a distance from the impact of a heavy structure fall. That model is used by oil and gas operators to estimate safety distances from their buried pipelines and the implantation of heavy structures nearby. Over the past few years, different models have been used, from very simplified ones based on simple wave equations and considering that all the energy of the falling structure was transmitted to the wave near the buried pipeline. The wave equation was stated in simple terms of u(x,t) U0sin(ωt + kx) functions. These models were overconservative because they did not take into account the energy dispersion into the ground and due to the plastification of the soil with the craterization phenomenon. Then, empirical models were developed thanks to extensive tests campaigns (like the Mayne model presented in [1] or the Menard model from [2] (presented below) where the peak particle velocity is estimated from the fall of masses for dynamic compaction techniques). bound for the velocity under an elastic impact. An efficiency coefficient is then introduced to take into account the plastification phenomena. This coefficient is updated in a fourth part thanks to vibrations measurements and MASW tests used as input parameters in a 3D dynamic non linear finite element model. The terms in bold are vectors. 2 2.1 EVALUATION OF THE VIBRATIONS IN THE GROUND Problem and basic assumptions The problem considered in this paper is the fall of a tall and heavy structure in the gravity field g with an initial velocity V0 and an incident angle θ from a height h producing an impact at a distance d from the base of structure. The fall produces vibrations nearby a buried pipeline located at a length r from the impact and a depth z from the surface soil (see Figure 1). But once again, the Mayne PPV model was not adapted for all the field configurations met by a gas operator. Indeed, the masses at play in ground compaction are an order of magnitude less than the masses of the wind-turbines and the energy dissipation into the ground does not follow a linear law in function of the falling mass. In a first part, the problem is introduced. In a second part, the velocity due to the Rayleigh waves is obtained for some loading. The aim of this section is to show that Rayleigh waves are preponderant in the phenomenon of interest. In a third part, under elastic assumptions, the impact force is estimated. The aim of this section is to provide an upper Figure 1. Problem modeled. Velocity vz at the impact point located at a distance d like in Figure 1, is written as, 899 Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 where τ = t cT / . (1) If θ = π/2, vz is maximal but due to the configuration of the wind-turbine nacelle, this value is highly unlikely to appear. Moreover, a variation of the θ angle from 0° to 30° on the horizontal gives a very small variation of the impact velocity, thus the angle is considered to be null and therefore, the velocity at the impact point is written as, . 2.2 r w Q 0,75 P (2) , (3) where σ is the stress tensor, Q is the force amplitude and δ is the Kronecker symbol. The Heaviside function , H(t) is null for t<0 and equal to 1 for t>0. With the assumption that the soil Poisson’s coefficient ν = 0.25, the Lame’s coefficients, λ and μ, which characterize the elastic behavior of the half-space are equal. The vertical displacement w(t,r), at a distance r from the impact point and time t is written as, for R 0,25 0,00 The original formulation for the normal point-source load over a half-space was established by Lamb for a harmonic load. The soil surface displacements for a point-source impact load is given by Pekeris (see [3]). The displacement is reproduced below for a load with a Heaviside time dependency; the derivation of this solution with modern techniques in the complex plane is available in reference [4]. At the surface of the half-space, the boundary conditions are written as, , S 0,50 Point-source impact, general solution for is the scaled time with cT the celerity of the shear waves, the time of arrival of the Rayleigh wave which propagates with celerity cR. These expressions are plotted at Figure 2. ; , -0,25 -0,50 cT t r -0,75 0,0 0,5 1,0 1,5 2,0 2,5 Figure 2. Vertical surface displacement (Pekeris, 1955), from equations (4). Amplitude of the displacement in vertical axis, dimensionless time in horizontal axis. P, S and R mean arrival of P-waves, S-waves and Rayleigh waves respectively. 2.3 Vibrations induced by surface waves, point-source impact, Equations (4) show that the Rayleigh wave gives the larger perturbation at the soil surface; these expressions show moreover a geometric attenuation when r increases. The following assumption is then made, which is endorsed by seismic observations: the perturbations propagating with the Rayleigh waves give larger amplitudes of vibrations which are of interest in this paper. The problem is then considerably simplified because the displacement of interest can be obtained by solving the equation in the vicinity of the pole associated to the Rayleigh waves. That solution was calculated by Chao et al. in ref. [5] (see also [4]). These expressions are valid for times close to the arrival time of the Rayleigh waves and for depths z/r (see Figure 1) which are small. They are written as, horizontal displacement for , (5) (4) ; vertical displacement for , ; , with the following notations, 900 (6) Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 ; 1,5E-01 ; (7) 1,0E-01 , 5,0E-02 with , and are respectively the real part and the imaginary part of a complex number Z. These expressions are plotted at Figure 3. 0,0E+00 du/dt -5,0E-02 dw/dt 2,0E-02 -1,0E-01 -3,00 -2,00 -1,00 0,0E+00 1,00 2,00 3,00 Figure 4. Velocity due to the Rayleigh waves. Horizontal axis, the scaled time (cT t – τR r) / z ; vertical axis, the scaled radial -2,0E-02 -4,0E-02 -6,0E-02 u -8,0E-02 w -1,0E-01 -3,00 -2,00 -1,00 0,00 0,00 1,00 2,00 ( ) velocity components. This figure shows the amplitudes of the velocity components due to the Rayleigh waves. 3,00 2.4 Figure 3. Displacements due to the Rayleigh waves (Chao et al. 1961). Horizontal axis, the scaled time (cT t – τR r) / z ; vertical axis, the scaled radial ( ) and vertical ( ) and vertical ( ) displacement components. This figure shows the amplitudes of the displacement components due to the Rayleigh waves. Vibrations induced by surface waves, impulse load The case of an impulse load can be calculated with previous equations, considering at time t an impact of amplitude + Q and at time t + dt, where dt is finite, an impact of amplitude Q, if dt is small enough, the solution is obtained by differentiation and superposition of equations (8) and (9), which is written as, Expressions (5) to (7) show that the displacements depend only on the dimensionless time τ and that they decrease in function of 1/ compared to the other components of the displacement which decrease in function of 1/r. Therefore for high distances from the impact, z/r << 1, the Rayleigh wave contribution is predominant. The velocity is obtained from the displacement by differentiation, which is written as, (8) , (10) . (11) Considering equations (8) and (9), equations (10) and (11) are written as, ; . ; (12) . (13) (9) These expressions are plotted at Figure 4. These expressions are plotted at Figure 5. 901 Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 , 4,0E-01 3,0E-01 2,0E-01 1,0E-01 0,0E+00 -1,0E-01 -2,0E-01 -3,0E-01 -4,0E-01 -5,0E-01 -6,0E-01 -3,00 -2,00 -1,00 (18) . 3 3.1 du/dt dw/dt 0,00 1,00 2,00 3,00 Figure 5. Velocity due to the Rayleigh waves and an impulse force. Horizontal axis, the scaled time (cT t – τR r) / z ; vertical axis, the scaled radial ( ) and vertical ( ) velocity components. This figure shows the amplitudes of the velocity components due to the Rayleigh waves and an impulse load Qdt. The graphs presented in Figure 5 show that, for any time t, IMPACT FORCE DETERMINATION Circular foundation Under the assumption that the soil behavior remains elastic when impacted, the wind-turbine nacelle is assumed to be a circular foundation (mass m) animated with an initial velocity equal to the impact velocity given in equation (2). In that case, the foundation can be modeled by a spring and a damper (assumed to be frequency independent for the sake of simplicity). The values of the spring and the damper used to model the impedance of a circular foundation over a homogenous, isotropic half-space are given by Gazetas (see reference [6]). For a foundation of radius r0 over a half-space with behavior parameters μ and ν (=0.25, like previously), the expressions are written as, , (14) . 2.5 With Vibrations induced by surface waves, for any load Any given load can be modeled by a superposition of timeshifted impulse loads. The velocity due to that load and generated by the Rayleigh waves is given by: and , (19) . (20) , the response of the 1-DOF oscillator is immediate and the displacement and vertical velocity are written as, (15) , (21) , (22) ; with vz is the velocity of the missile (i.e. wind-turbine nacelle) at the impact and where the following notations are used, (16) , (23) . (24) , where T = N dt is the total duration of the impact force application, Φ and are functions not used below. When dt becomes close to 0, the last member of the equations (15) and (16) tends to the total impulse contact load, which is written as, The impact force can then be written as, = . (17) With equation (17), equations (14) to (16) provide the following equations, 902 . (25) As an example, Figure 6 depicts the variation of the missile displacement and velocity after the impact, and the developed contact force. This graph was plotted for a mass m = 50.103 kg Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 and r0 = 2 m, falling from 80 m over a soil with the following characteristics, cT= 300 m.s-1, μ = 162 .106 Pa. Force 600 (MN) 550 500 450 400 350 300 250 200 150 100 50 0 -50 -100 0 parameters above and the 3 fundamental terms). The following dimensionless parameters are chosen and written as, D. (mm) , V. (m.s-1) 0,02 t(s) 120 110 Force 100 90 Displacement 80 70 Velocity 60 50 40 30 20 10 0 -10 -20 0,04 0,06 , (26) . (27) and the relation is written as, The dimensionless parameters for i=1 to 3 are computed with the values shown in Table 1. Figure 7 and Figure 8 show the variation of as a function of and . 3,50 3,00 Figure 6. Impact characteristics. Horizontal axis, time t; left vertical axis, Force (106 N) in dotted line, right vertical axis, displacement (mm) continuous line and velocity (m.s-1) in discontinuous line. Force diagram can be seen as a triangular force diagram with maximum value at 390 MN, duration 0.013 s and total impulse 2.55 MN.s. The maximum force is obtained for t = 0 and equal to C vz. After 0.012 s, the force becomes negative, which means, taking into account the previous assumptions that there is a traction force from the missile to the soil; this phenomenon cannot be possibly physical and thus, after 0.012 s the mass lifts off and previous equations are not valid any more. 2,50 2,00 1,50 1,00 0,50 0,00 0 3,50 In this section, the impulse load (see eq. (18)) is computed numerically on the interval [0,T]; that corresponds to a positive force and the parameters used vary like in the following Table 1. 3,00 Table 1. Parameters variation for the computation of the total impulse load IQ. 1,50 Parameters Min. value Symbol Max. value Equivalent radius for 0.5 ≤ r0 ≤ 5m impact area (m) Nacelle mass 25 ≤ m≤ 70 t (t 103 kg) Fall height 50 ≤ h≤ 100 m (m) Shear waves celerity in the 150 ≤ cT ≤ 1000 m.s-1 soil (m.s-1) The total impulse IQ depends a priori on the following terms: m, vz, r0, ρ, cT, where ρ=μ/cT2 is the soil density. These parameters can be expressed with the fundamental terms L (length), T (time) and M (mass). The Vashy-Buckingham dimensional analysis theorem states that it exists a relation between 3 dimensionless parameters (relating the 6 40 60 Figure 7. Computation of parameters (vertical axis) versus (horizontal axis). Maximum velocity induced by an elastic impact 3.2 20 2,50 2,00 1,00 0,50 0,00 0 10 20 30 40 Figure 8. Computation of parameters (vertical axis) versus (horizontal axis). Figure 7 and Figure 8 show that the variation of versus and is relatively small. A simplified average value can be retained and is written as, . (28) For a short duration impact, the elastic forces linked to the stiffness K (see eq. (19)) become negligible in front of the damping force from C (see eq. (20)) and thus the impulse is simplified and its value becomes i.e. . In the following, that value will be taken as a reasonable approximation for the impact computation. The chosen upper 903 Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 bounds for the radial and vertical velocities induced by an elastic impact and propagated with Rayleigh waves are then written as, 3.4 , (29) . 3.3 Application to a wind-turbine fall case Under the assumptions of section 2, which are recalled below, Rayleigh waves predominance over the general movement; Linear elastic behavior for the soil at the impact point and in the medium, equations (29) give an evaluation for the particle velocity in the soil as a function of the soil characteristics ( ), the mass of the impact (m), the velocity at the impact ( ) and the distance to the impact point ( ). The depth ( ) is arbitrarily chosen to be equal to 5 m (the linear elastic soil behavior leads to infinite displacement and velocity at the soil surface). The velocity is estimated far enough (z/r << 1) from the impact point. Impulse: 1 MN.s 2 MN.s 3 MN.s 0.3 MN.s Experimental dynamic compaction 1,0E+01 1,0E+00 (29) with 1,0E-01 1,0E-02 is plotted so as to compare with the model of eq. = 0.3 MN.s (m ≈ 20 t and h ≈ 10 m). Discussion and efficiency coefficient When comparing the experimental results with the predicted model, it appears that the model assumptions lead to an overestimation of the peak particle velocity for dynamic compaction impacts. The overestimation is due to the fact that the impact from real masses, like in dynamic compaction, is far from being elastic. Moreover, the wave propagation takes place in a nonlinear medium. For these reasons, an efficiency coefficient e is introduced, which accounts for the plasticity and non linearities. Figure 9 shows that the efficiency coefficient e can be estimated, for the lengths of interest, to . (30) The assumptions for using the empirical Menard model (from [2]) are the following, the product mh should be included between 200 and 300 t.m (for example a fall of a 20 t-mass from 10 to 15 m). The assumptions from [1] are mh ≤ 1000 t.m (for example 50 t from 20 m). As a reminder, the aim of this paper is to estimate the in-ground velocity due to a heavy structure fall like wind-turbines. For wind-turbines, the product mh can be as large as 25000 t.m (a 250 t-nacelle falling from 100 m) and much more for concrete-mast windturbines (on-shore). The efficiency coefficient found for the dynamic compaction tests is thus not assumed to be valid for the fall of a windturbine. 1,0E-03 1,0E-04 4 1,0E-05 100 1000 4.1 Figure 9. Plot of the particle velocity (in m.s-1, vertical axis, log scale) in function of the distance (impact distance + 26 m) to the base of a wind-turbine (in m, horizontal axis, log scale) for 5 different impulse loads: big dashed line for 1 MN.s ; very big dashed line for 2 MN.s; continuous line for 3MN.s; dotted line for 0.3 MN.s; dashed line for experimental dynamic compaction estimation of the velocity. The results plotted in Figure 9 for a soil density of 1800 kg.m3 give the variation of the maximum amplitude for the velocity ( , computed with eq. (29)) as a function of distance to the foot of a wind-turbine mast which is equal to r -d (feedback studies state that the nacelle falls at most at d = 26 m from the foot of wind-turbine mast), where r and d are shown in Figure 1. The input parameters used are the following impact loads = 0.3, 1, 2 and 3 MN.s. An empirical relation based on tests (from [2]), where 904 UPDATED EFFICIENCY COEFFICIENT FOR WINDTURBINES Methodology of updating In the previous section, an efficiency coefficient has been updated for dynamic compaction phenomenon thanks to empirical data. In the case of this paper, a test with real masses of interest is not possible (because it is hard to find a field to instrument and make a wind-turbine fall). Therefore a 3D finite element computational model (Dynaflow, see [9]) has been used to model the fall of the masses of interest. This FE-model has been updated with soil parameters measured in field experiments performed by GDF SUEZ (see [7] and [8]) and the dynamic compaction parameters. The updated FE-model has then been used to compute a windturbine fall, which has provided an estimation of the efficiency coefficient for real wind-turbines. Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 4.2 Field tests The field tests used to update the FE-model have been extensively presented in ref. [7] and [8] and only main results are reproduced in this paper. 4.3 Figure 10. Position of the sensors for the 1.5 t-mass fall tests from 3, 5, 10 to 15 m high and measurement of the velocity from 3, 8, 15, 30, 50 to 80 m from the impact point. Frequency domain measured is in [0, 30 Hz]; Impact footprint is about 0.70 cm in diameter, the depth increased with the number of fall from 0.30 cm up to 1.00 m; Measured particle velocities show a non negligible scatter; Some sensors were occasionally saturated (150 mm.s-1). 3D Nonlinear finite element model for the soil The Dynaflow software has been used (see [9]). The impact simulation was used with the following parameters: 1.5 t mass, 0.35 m radius, 0.50 m depth of impact and initial velocity given in eq. (2). The Prevost constitutive model used in the computation takes into account the data gathered at Table 2 and some parameters were taken from literature because they were not measured (see ref. [10] and [11]). The mesh was adapted to the waves transmission up to 30 Hz. The criteria used was 10 elements by wavelength, the maximum mesh size was calculated, . Figure 11. Position of the sensors during the MASW tests. (31) The mesh dimensions presented in Table 2 are compatible with eq. (31). The non linearities were modeled up to 10.50 m from the impact point (it is considered that the medium is elastic beyond that radius, up to 150 m). The damping follows a Rayleigh law and is fixed at 2% for 0.5 Hz and 20 Hz. The FE-model has 14,523 elements, 14,886 nodes and 29,048 DOF. Nodes on the lateral boundary have only vertical DOF and nodes on the lower boundary are fixed (see Figure 12). Several heights of fall have been tested to quantify the influence of plastification at the impact point; for each fall, 5 to 6 impacts have been performed (see Figure 10). A scatter of about 2 has been observed on the results (up to twice the minimal value for a test). Figure 12. FE mesh. Geotechnical data have been measured at the same field with Multichannel Analysis Surface Wave method and were used in the FE-model (see Figure 11) and are gathered in Table 2. Table 2. Parameters measured from MASW test (thickness layer, celerities) and estimated from experience (density). Mesh dimensions for the FE-model. Layers 1 2 3 4 5 6 Thickness (m) 1.5 2.5 2.0 4.0 5.0 5.0 Density (t.m-3) 1.8 1.8 1.9 2.0 2.0 2.0 Shear waves celerity (m.s-1) Compression waves celerity (m.s-1) Mesh height (m) 150 220 220 330 450 450 400 400 800 800 1800 1800 0.25 0.24 0.35 0.50 0.50 0.50 Mesh width (m) The numerical integration scheme is detailed in ref. [12]; the non linear system of equations is solved with a modified Newton-Raphson algorithm. The integration step is taken equal to 10-3 s and total analysis duration is taken equal to 1 s. The first step of the computation, which is the initialization of the stresses in the soil, is statically realized; after stress initialization, displacements are then reset to 0 and the propagation is computed. 4.4 0.35 0.35 0.35 0.35 0.35 0.35 Results Figure 13 shows an example of result obtained with the FEmodel which gives conservative results compared to the field measurements. It can be said of the measurements: 905 Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 CONCLUSION 100 FE-model Field measurements 10 1 y = 141.95x-1.153 0,1 y = 71.966x-1.791 0,01 0 20 40 60 80 100 Figure 13. Radial velocity (m.s-1, log-scale, vertical axis) in function of distance to the impact point (m, horizontal axis). Continuous line: FE-model prediction with interpolation curve in dashed line; Triangular symbols: field measurements with interpolation curve in continuous line. All comparisons between FE-computations and field measurements have globally the same shape, which means that the relation between the velocity and the other parameters (see eq. (29)) is different from the relationship used with linear elasticity assumptions multiplied by a constant efficiency coefficient (except for distances far enough from the impact). In this paper, an analytical model for estimating the peak particle velocity near a buried pipeline due to the fall of a heavy structure like a wind turbine is presented. In a first step, it is assumed that the impact and then the propagation are elastic. It is also assumed that the Rayleigh waves are predominant. These assumptions lead to a model which can be updated with an efficiency coefficient when the impact is far enough from the buried pipeline, the efficiency coefficient is a simple constant scaled with empirical relations based on field experiments. However, an estimation of the velocity obtained with a 3D non linear plastic finite element model in time domain, updated with specific field test data showed that the approximation for distances close to the impact were not proportional to the linear elastic solution (it is still true with an efficiency coefficient, far enough from the impact), therefore a PPV relation based on empirical relations used in dynamic compaction was performed. The impact energy for windturbines falls is an order of magnitude larger than for dynamic compaction, then the PPV relation presented in this paper had to be updated numerically for the specific values of the windturbines parameters. A specific example is presented: fall of a 70 t wind-turbine nacelle from a height of 100 m. ACKNOWLEDGMENTS The authors thank the GRTgaz Engineering Center in Lille and Paris, who found the right place to perform the seismic tests, and also INNOGEO who performed the field tests. Therefore, a relation based on the Mayne model ( ) and the Menard model ( ) [1] is written as, , (32) where, PPV is the Peak Particular Velocity in mm.s-1; r is the distance to the impact point in m; m the impacting mass in kg; is the velocity at the impact in m.s-1; is the soil density in kg.m-3; and n depends on . Figure 14 shows the use of equation (32) for the fall of a 70 t nacelle from 100 m. 500,0 50,0 5,0 0 50 100 150 Figure 14. Peak particle velocity due to the fall of windturbine nacelle of 70 t from a height of 100 m. Vertical axis, PPV in mm.s-1(log-scale), horizontal axis, distance in m from the point of interest to the foot of the wind-turbine mast. 906 REFERENCES P.W. Mayne, J.S. Jones and J.C. Dumas, Ground response to dynamic compaction, Journal of Geotechnical Engineering (109), n°16, pp757774, ASCE, 1984. [2] Menard Sol Traitement. Vitesse particulaire resultante en function de la distance au point d’impact. Private conversation, 2001 [3] C. L. Pekeris, The seismic surface pulse, Proceedings of the National Academy of Sciences, (41), pp469-480, 1955. [4] J.D. Achenbach, Wave propagation in elastic solids, North Holland Publishing Company, 1973. [5] C.C. Chao, H.H. Bleich and J. Sackman, Surface waves in an elastic half space, Journal of Applied Mechanics (28), pp300-301, 1961. [6] G. Gazetas, Analysis of machine foundations vibrations – State of the Art, International Journal of Soil Dynamics and Earthquake Engineering (2), pp2-42, 1983. [7] Fernandez, C., Bourgouin, L. Riegert, F. & Pecker, A., Modelling of Windturbine fall-induced vibrations near buried steel transmission pipelines, an updated RAMCES software extension. Proceedings of the IPC 2012, Calgary, Canada. [8] Fernandez, C., Bourgouin, L. Riegert, F. & Pecker, A., The induction of vibrations in transmission pipelines by the fall of a heavy structure nearby: modelling the safety distances. The Journal of Pipeline Engineering (12) n°2, pp75-82, 2013. [9] J.H. Prevost, Dynaflow, A finite element analysis program for the static and transient response of linear and non linear two and threedimensional systems, Department of Civil Engineering – Princeton University – V2002, V02. [10] J.H. Prevost, Plasticity theory for soil stress-strain behavior, Journal of Engineering Mechanics, ASCE, (104), n°EM5, 1978. [11] J.H. Prevost, A simple plasticity theory for frictional cohesionless soils – Journal of Soil Dynamics and Earthquakes Engineering (4), n°1, pp917, 1985. [12] T.J.R. Hugues, Analysis of transient algorithms with particular reference to stability behavior, Journal of Computational Methods for Transient Analysis, Elsevier, 1983.
© Copyright 2024