Semi-Analytical Integration of the Accelerogram and Equations of Motion of a System Subjected to a Paraseismic Shock

This contribution is an extension of a conference paper [1]. Motivation was an issue that my former student Kinga Zemła in her research of the master thesis. She compared results of velocity diagrams obtained from two numerical algorithms of accelerogram integration and they were different. The problem was related to the baseline correction described in the review paper [2] and cited in its source research papers.

Analysis of paraseismic shocks in mining activities areas is a difficult engineering task, and the reliability of the results is very important.

Usually, numerical integration algorithms are used to evaluate velocity and displacement functions [3, 4].

This paper presents an analytical approach to integration that is possible to apply thanks to the Mathematica system [5]. The procedure obtained is straightforward and surprisingly effective. It is presented as an alternative to numerical methods.

Moreover, the approach may be applied to integration of the equation of motion, and the crucial steps of this algorithm are also presented. The numerical approach to this problem is described in many books, for example [6, 7, 8].

Figure 1 shows an accelerogram from the research mentioned in the Introduction, as an example of the method used. The entire recorded signal lasted 5 seconds, but the largest acceleration amplitudes were recorded in the subinterval shown on the right side of the figure. The set of recorded data within this time period consists of 327680 elements. Thus, the time distance between them is 0.000152588 s. The analysed signal is highly oscillating, and we cannot say that measurement points are very dense; see Fig. 2.

Figure 1.

Accelerogram in the entire measurement range and in the subrange of high accelerations

Figure 2.

Interpolation of measurement data in the range of extremal accelerations

Numerical integration of this signal in Mathematica is possible, but the only effective algorithm with the use of a built-in function requires the division of this time interval into 327679 subintervals. Computation took an unreasonably long time.

In the Mathematica system, we can interpolate a set of measurement points with a function. The argument of this function is the order of the interpolation function. This corresponds roughly to the degree of polynomial fit between the measurement points of the set ts.

a[i_]:=

 a[i]:=

  Interpolation[ts, Interpolation → i]

The interpolation with interpolation order equal to 6 is presented in Fig. 2. From this moment we can deal with the interpolation function almost like with the analytical one. This “almost” and the fact that it is based on numerical data explains the word “semi-analytical” in this title of the article.

The function of acceleration is highly oscillating. As has already been mentioned, numerical integration fails. It turns out that the analytical integration of this function is very fast.

We compute the velocity (speed) function as an integral: v[ i− ][ t− ]:=v[i][t]=∫Evaluate[a[i][t]]dt

Figure 3 shows that the highest velocity values are reached around the 2 s of measurement, when the vibrations seem to stop. The speed value is about three times higher than at the moment of the greatest accelerations. Further vibrations with small accelerations cause the speed to drop to zero in the fifth second.

Figure 3.

Speed chart in the entire measurement range and in the subrange of high accelerations

Figure 4.

Displacement graph in the entire measurement range and in the subrange of high accelerations

By integrating the semi-analytical velocity function, we obtain the displacement function. δ[ i− ][ t− ]:=δ[i][t]=∣ Evaluate [v[i][t]]dt

The graphs shown in Fig. 4 may suggest that the displacement function is quite smooth. However, a “microscopic” close-up (Fig. 5) to the time interval as in Fig. 2 shows the actual course of the function after successive integrations.

Figure 5.

Velocity and displacement diagrams of displacements in the range as in Fig. 2

We can check the quality of integration by differentiation. So, if we differentiate the velocity function and subtract the acceleration function from it, Abs[ ∂{t}v[6][t]−a[6][t] ] we get what can be seen in Fig. 6, that the error is negligible since even at points where acceleration is equal to 3000 mm/s² the error is not bigger than 15 mm/s2 (0.5%). When we integrate the error function and divide it by an integral of an absolute acceleration value, we obtain an average error that is equal to 0.2%.

Figure 6.

Absolute error of interpolation

The above method can be used to integrate equations of motion [3, 4]. We will show it in the example of a dynamical system with one degree of freedom excited by the analysed acceleration function.

Mathematica can provide a general solution to this problem using a function to analytically solve differential equations. DSolve [{q''[t] + 2cq′[t] + ω2q[t]==a[t],q[θ]==θ,q′[θ]==θ},q[t],t]/⋅{K[i−]→τ} $$\eqalign{ & {\rm{DSolve }}\,[\{ {{\rm{q}}^{''}}[{\rm{t}}]{\rm{ }} + {\rm{ }}2{\rm{cq}}\prime [{\rm{t}}]{\rm{ }} + {\rm{ }}{\omega ^2}{\rm{q}}[{\rm{t}}] = = {\rm{a}}[{\rm{t}}], \cr & {\rm{q}}[\theta ] = = \theta ,{\rm{q}}\prime [\theta ] = = \theta \} ,{\rm{q}}[{\rm{t}}],{\rm{t}}] \cr & / \cdot \{ {\rm{K}}[{\rm{i}} - ] \to \tau \} \cr} $$

The result of such a generally formulated task is given by the program in the following form. { { q[t]→−et(−c+c2−ω2)∫10e2cτ+τ(−c−c2−ω2)c2−ω2a[τ]2(c−ω)(c+ω)dτ+et(−c+c2−ω2)∫1te2cτ+τ(−c−c2−ω2)c2−ω2a[τ]2(c−ω)(c+ω)dτ−et(−c−c2−ω2)∫10−e2cτ+τ(−c+c2−ω2)c2−ω2a[τ]2(c−ω)(c+ω)dτ+et(−c−c2−ω2)∫1t−e2cτ+τ(−c+c2−ω2)c2−ω2a[τ]2(c−ω)(c+ω)dτ } }

This expression – as beings endowed with real intelligence, not artificial, still superior to machines – we can simplify to the form: q[ t− ]:=e−t(c−c2−ω2)∫eτ(c−c2−ω2)a[τ]2c2−ω2dτ−e−t(c+c2−ω2)∫eτ(c+c2−ω2)a[τ]2c2−ω2dτ which in traditional notation can be denoted as: q(t):=e−t(c−c2−ω2)∫eτ(c−c2−ω2)a(τ)2c2−ω2 dτ−e−t(c+c2−ω2)∫eτ(c+c2−ω2)a(τ)2c2−ω2 dτ

This solution satisfies the initial conditions that, at time t = 0, both the displacement and the velocity are equal to zero. This is done taking into account the assumption that in the launching moment of motion registration the system do not move. It makes it possible to omit the discussion on constants of integration.

It is necessary to explain why definite integrals have been replaced by indefinite ones. It is possible since the definite integral is equal to the difference of antiderivatives of the integral at both ends of the integration interval and the antiderivative of the interpolated function at the moment of time equal to 0 is equal to 0, too.

In this case, the interpolation procedure requires an additional step. The integrands in the above formula are a mixture of analytical functions and interpolated acceleration. f1=τ↦eτ(c−c2−ω2)a[6][τ]2c2−ω2;f2=τ↦eτ(c+c2−ω2)a[6][τ]2c2−ω2;

Such mixtures cannot be integrated analytically.

Let us set the cyclic frequency and dumping parameter equal to ω=50;c=2

To make it possible to integrate these functions analytically, we have to reinterpolate them: cf1 = Interpolation [Table[ {t,f1[t]},{ t,0,5,1327679×2 } ], InterpolationOrder →6]; cf2 = Interpolation [Table[ {t,f2[t]},{ t,0,5,1327679×2 } ], InterpolationOrder →6];

Now, the calculations are straightforward. Input in the form. y[ t− ]=e−t(c−c2−ω2)∫cf1[t]dt−e−t(c+c2−ω2)∫cf2[t]dt produces the following output: −e−((2+8i39)t)

This result is presented in Fig. 7

Figure 7.

Result of integration of the equation of motion

It could be added that the general solution of a differential equation of motion can be presented in the equivalent form: q[ t − ]:= 1 e ct (cos[ ω 2 − c 2 t ] ∫ ​ − e cτ a[ τ ]sin[ ω 2 − c 2 τ] ω 2 − c 2 dτ+ sin[ ω 2 − c 2 t ] ∫ ​ e cτ a[ τ ]cos[ ω 2 − c 2 τ ] ω 2 − c 2 dτ)

in traditional notation: q(t):=1ect( −cos(ω2−c2t)∫ecτa(τ)sin(ω2−c2τ)ω2−c2 dτ +sin(ω2−c2t)∫ecτa(τ)cos(ω2−c2τ)ω2−c2 dτ )

This form seems to be more suitable when c < ω. In this case, all the terms in the expression above are real numbers. In the previous case, they are complex numbers, but the final result is a real function. It is not a problem for Mathematica to deal with such functions.

Both forms of equations exactly satisfy the differential equation and initial conditions. However, the solution was compared with numerical approaches. The numerical solution is less precise and requires significant dense steps, so the computation time is comparable. This comparison will be shown in the following contributions.

The numerical problem may occur with both formulas mentioned above if we are close to critical dumping c ≈ ω. In the case of critical dumping c = ω we have the following solution: q[ t− ]:=e−tω(t∫eτωa[τ]dτ−∫eτωτa[τ]dτ).

in traditional notation: q(t):=e−tω(∫eτωa(τ)dτ−∫eτωτa(τ)dτ)

Numerical problems connected with this special case can be overridden with higher precision of computation, since Mathematica can do calculations with an arbitrary precision.

The presented approach can be an alternative and verification tool for numerical algorithms embedded in accelerometer and system software for numerical calculations in the near future. This statement requires further comparison with contemporary professional numerical algorithms. Such comparisons will be made in further papers. The high speed of semi-analytical calculations and the verifiability of the results speak to their use in the analysis of para-seismic shocks in mining areas.