# Multi-modal transport prediction for linear one-dimensional wave propagation problems

## Multi-modal transport prediction for linear one-dimensional wave propagation problems

Proceedings of the 9th Vienna International Conference on Proceedings ofModelling the 9th Vienna International Conference on Mathematical Proceedings ...
Proceedings of the 9th Vienna International Conference on Proceedings ofModelling the 9th Vienna International Conference on Mathematical Proceedings of of the the 9th 9th Vienna Vienna International International Conference Conference on on Proceedings Mathematical Modelling Vienna, Austria, February 21-23, 2018 Available online at www.sciencedirect.com Mathematical Modelling Mathematical Modelling Vienna, Austria, February 21-23, 2018 Vienna, Austria, Austria, February February 21-23, 21-23, 2018 2018 Vienna,

ScienceDirect

IFAC PapersOnLine 51-2 (2018) 565–570

Multi-modal transport prediction for linear Multi-modal transport prediction for linear Multi-modal transport prediction for linear one-dimensional wave propagation one-dimensional wave propagation one-dimensional wave propagation problems problems problems

Proceedings of the 9th MATHMOD 566 Vienna, Austria, February 21-23, 2018

T. Weinberger et al. / IFAC PapersOnLine 51-2 (2018) 565–570

2. PROBLEM STATEMENT

which describes the dependency of γ and ω. The phase velocity of travelling wave mode (γ, ω) is ω (8) vphase = − γ considered for real γ and ω, and generally changes with frequency (dispersion), compare Ritzberger et al. (2016).

2.1 Finite Element catenary model The catenary, here simplified to the contact wire’s dynamics, is modeled by an axially pre-tensioned Euler-Bernoulli beam. This leads to a Partial Differential Equation (PDE) of second order in time and fourth order in space, ∂4w ∂2w ∂2w ∂w + EI 4 − T =f (1) ρA 2 + β0 ∂t ∂t ∂x ∂x2 with homogeneous boundary conditions w(x = 0, t) = w(x = L, t) = 0, (2) ∂w ∂w (x = 0, t) = (x = L, t) = 0 (3) ∂x ∂x and initial conditions ∂w w(x, t = 0) = w0 (x), (x, t = 0) = v0 (x). (4) ∂t The total domain length is L, mass per length unit is ρA, the viscous damping coefficient β0 , bending stiffness EI, the axial pretension force T , the distributed vertical force per unit length f (x, t), and the vertical displacement field w(x, t). As in Aschauer et al. (2016) the solution will be approximated with the spatial Finite Element (FE) method, leading to Hermitian beam elements with thirdorder polynomial interpolation of the displacement field. The assembled FE model equations form a system of coupled 2nd-order ordinary differential equations in time, ¨ ˙ MQ(t) + CQ(t) + KQ(t) = F(t), (5)

2.3 Exemplary case The domain is decomposed into two subdomains of interest, A and B. Subdomain A, surrounding the point x = xA with displacement qA = w(x = xA , t) acts as the “source” domain where the signal will be analysed; subdomain B will be called the “destination” domain, into which the waves will be predicted. The reassembled modal wave components will be fed into B on the point x = xB as a time-variant boundary condition qB = w(x = xB , t), compare Fig. 1. 3. METHOD The method can be divided into three parts: (i) disassemble the wave into modal components in the “source” domain (ii) transport the modal components across a gap in space (iii) reassemble the parts in the “destination” domain In the following, “(t)” is equally used to “(n∆t)” to express the dependency on discrete time steps.

with the global degrees of freedom Q(t), the system’s mass, damping, and stiffness matrices M, C, and K, and F(t) as the vector of external nodal forces and moments. Equation (5) is solved by an efficient numerical Newmarktype integration scheme (Aschauer et al., 2016).

3.1 Modal disassembling The idea is first to approximate the solution in the “source” domain via a modal basis comprised of m travelling wave modes. These modes are chosen over a frequency grid over the range of frequencies of interest. Typically this range can be found through an offline simulation. The travelling wave solution of the displacement field in A will be approximated by m  (9) θi (t)κki ζin w(k∆x, ˆ n∆t) =

It is noted here that the considered Euler-Bernoulli beam dynamics under axial tension shows non-zero propagation velocity for zero frequency and unbounded propagation velocity for high-frequency components; however, its FE approximation is inherently limited in frequency, and also in typical applications (simulation, control), the frequency range of interest is bounded. Consequently, also the relevant travelling wave components have bounded propagation velocities, albeit showing significant physical dispersion. The proposed transport prediction technique is designed to correctly account for dispersion and capable of describing transient processes in real-time applications within a problem-specific, bounded frequency spectrum.

i=1

with

κi = ejγi ∆x and ζi = ejωi ∆t (10) constant, when ∆x and ∆t is fixed. Spatial and temporal frequencies of the i-th chosen harmonic function are denoted by γi and ωi . (11) θi (t) = αi (t) + jβi (t) is a complex coordinate of the i-th mode, with αi its real and βi its imaginary part, which leads to a modal representation linear in αi (t) and βi (t). To determine these, a spatially-weighted Least-Squares procedure is applied at each time step n, where the FE solution is represented w ˆ by the degrees of freedom, w, ˆ w ˆ˙ = ∂∂twˆ , w ˆ  = ∂∂x , and 2 ∂ w ˆ  ˙ w ˆ = ∂x∂t in subdomain A (0 ≤ k ≤ K). It is important to choose the “source” domain large enough, so that the lowest relevant frequency is observed for at least one full period.

2.2 Dispersion relation Travelling wave solutions of the undamped, constantcoefficient, linear PDE, discretised via spatial Finite Elements, can be expressed by: w(k∆x, n∆t) = (ejγ∆x )k (ejω∆t )n (6) with j the imaginary unit, γ ∈ R the spatial frequency, ω ∈ R the temporal frequency, ∆x, ∆t the discretisation step in space and time, and k, n as natural numbers. By inserting (6) into the autonomous discretised PDE one obtains the discrete dispersion relation of the form: F (γ, ω) = 0 (7) 2

Proceedings of the 9th MATHMOD Vienna, Austria, February 21-23, 2018

T. Weinberger et al. / IFAC PapersOnLine 51-2 (2018) 565–570

Eventually the extended, weighted Least-Squares problem for θ¯ reads:       λw DwFE λw DXw α1  λw˙ DXw˙   β1   λw˙ DvFE     LS  λφ DφFE   λφ DXw    .  =  . (15)  λ DX      ˙  .  λ w ˙  ˙ DφFE  φ˙   φ  λ I   αm    0 θ βm λdyn I λdyn θ¯∗ (n∆t) with ¯ − 1)∆t) θ¯∗ (n∆t) = Mθ((n (16)   ζ1 jζ1 0 .    −jζ1 ζ1 . .    M = Re  (17) .. ..  . .     ζm jζm 0 −jζm ζm The modal approximation of qA is built from the obtained ¯ The contribution of the i-th modal coordinate vector θ. mode to qA is

The linear system of equations to solve by Least-Squares is of the form: LS Xθ¯ = w (12)   α1     wFE Xw  β1    X   vFE  .  X =  w˙   θ¯ =  (13)  ..  w =  φFE  Xw α  m Xw˙  φ˙ FE βm   n n n jζ1n · · · ζm jζm ζ1 n n n n  κ1 ζ1 jκ1 ζ1 · · · κm ζm jκm ζm   , (14)  Xw = Re  .. .. ..   ... . . .

n K n K n K n κK 1 ζ1 jκ1 ζ1 · · · κm ζm jκm ζm and Xw˙ , Xw , and Xw˙  are built up analogously with the corresponding derivatives of the field variable. The nodal displacements wFE , velocities vFE , angles φFE , and angular velocities φ˙ FE have K + 1 entries each, corresponding to the values on the discrete FE nodes in the “source” domain.

K

qA,i (t) = θi (t)κi2 ζin , i = 1, ..., m. Finally, qA is formulated using (9): m m   K qA (t) = θi (t)κi2 ζin = qA,i (t)

The Least-Squares modal splitting is extended by weightings to improve performance: (i) Spatial weightings are introduced by a diagonal matrix D (for example, to implement triangular window function over subdomain A). D is multiplied by the regressor matrix X and the right-hand side vector w to achieve smooth data windowing and thereby an improvement of robustness and transient performance of the fitting task. Fig. 2 illustrates the spatial weighting with a triangular window function, applied to the FE data wFE .

i=1

(iv) To suppress high-frequency noise in the parameters’ time evolution during transients, their deviation from the predicted (rotated) parameter values are penalized and weighted by a cost factor λdyn . This prediction vector θ¯∗ (n∆t) is obtained from the previous ¯ parameters, θ((n − 1)∆t), by rotating each modal coordinate tuple by ζi .

0.02

0.5

0.01 0 -8

-6

-4

-2

0

2

4

6

8

3.3 Reassembling

weight value

displacement (m)

1

solution weighted solution spatial weighting D

(19)

i=1

With the distance between the two marked points of the subdomains, ∆L = xA − xB and the phase velocity of the i-th mode vphase,i , the time tlag,i required to bridge the gap for mode i equals ∆Lγi ∆L tlag,i = =− . (20) vphase,i ωi In the time-discrete form of the FE model, tlag,i can be expressed by: (21) tlag,i = ∆t(nlag,i + nroff,i ) where the integer time step count and the fractional remainder are:   tlag,i tlag,i nlag,i = round nroff,i = − nlag,i (22) ∆t ∆t The major part of discretised elapsed time is put into nlag,i which affects directly performance of transient processes in a positive way. Only a small number (−0.5 ≤ nroff,i < 0.5), describing the roundoff-error because of discrete time grid, is used to efficiently correct the phase angle, compare Fig. 3 where phase velocity and roundoff-error are interpreted on a two-dimensional space-time-grid. Using (20), the rounded time step countat the same phase angle of mode  ∆Lγi i at qB is nlag,i = round − ∆tωi .

(iii) A Tikhonov-type regularization is added by assigning cost λθ to the modal parameters directly to avoid overfitting.

subdomain A

(18)

3.2 Modal transportation

(ii) The regressor matrix X and the right-hand side vector w are additionally multiplied by weights λw , λw˙ , λφ , λφ˙ to adjust the relative importance of the involved field variables and hence to tune the fit over frequency.

0.03

567

To reassemble the signal at x = xB , it is noted that the individual modal contributions originate from different times at xA according to (20). The displacement qB,i (t) of the i-th mode arises from the corresponding earlier value of qA,i : (23) qB,i (t) = qA,i (t − tlag,i )

0

x (m)

Fig. 2. Displacement field with triangular weighting 3

Proceedings of the 9th MATHMOD 568 Vienna, Austria, February 21-23, 2018

T. Weinberger et al. / IFAC PapersOnLine 51-2 (2018) 565–570

Fig. 4. Principle to apply the multi-modal transportation method 4. SIMULATION The described method is demonstrated on a test problem in the following. A long axially pre-tensioned EulerBernoulli beam is formulated in a reference domain (0 ≤ x ≤ 1000m). It is discretised via the Finite Element method (spatial step size/element length ∆x = 0.5m) and integrated over time using a Newmark integration scheme (time step size ∆t = 0.001s). The axial pretension force is T = 20000N, the specific mass is ρA = 1.35kgm−1 , bending stiffness amounts to EI = 150Nm2 , and no gravity or damping β0 = 0 is considered. The beam is excited by an oscillating force F (t) at the point xC = L/2 = 500m of the reference domain. Subdomain A is defined at x ∈ [xA −7.5m, xA +7.5m] with xA = 400m, to observe and approximate the numerical solution with the proposed weighted Least-Squares approach at each time step. The target point of subdomain B is chosen at xB = 300m, see Fig. 5. The required phase velocities of each predefined mode are measured beforehand by simulating with this mode only and measuring the elapsed time for the wave to propagate a distance of 100m.

Fig. 3. Phase velocity and roundoff-error interpreted on a two-dimensional space-time-grid Using (21), the contributions at time t = n∆t are qB,i (n∆t) = qA,i ((n − (nlag,i + nroff,i ))∆t)

(24)

and, with (18), it follows that K

n−(nlag,i +nroff,i )

qB,i (n∆t) ≈ θi ((n − nlag,i )∆t)κi2 ζi ≈ qA,i ((n −

−n nlag,i )∆t)ζi roff,i

=: q˜B,i

(25) (26)

holds. The i-th mode of qB is built from previous modal contributions at qA , nlag,i time steps earlier and corrected −n by the roundoff phase increment ζi roff,i . The amplitude remains uncorrected with respect to nroff,i which is typically admissible: If θi is Lipschitz-continuous with Lipschitz-constant qˆ˙i , the absolute error is bounded by |qB,i − q˜B,i | ≤ qˆ˙i · ∆t/2. Finally qB is obtained by reassembling its modal components: qB (t) :=

m  i=1

q˜B,i (t) =

m  i=1

Fig. 5. Setup of simulation case study

−nroff,i

qA,i ((n − nlag,i )∆t)ζi

(27) The time shift operator z −nlag,i and the corrective roundoff −n term ζi roff,i can be computed in advance to establish a shift register and a corrective multiplication factor for each mode. The schematical application of the multimodal transport process, from the input parameters to the prediction of qB in a real-time-capable architecture is shown in Fig. 4.

4.1 Single-mode prediction First, the basic fitting and prediction capabilities are investigated for a single mode. A vertical force signal F (t) = 200N · sin(10s−1 · 2π · t) is applied at point xC . To approximate the displacement field in subdomain A, one mode with 10Hz and a static offset (0Hz) is used to form the regressor, compare Table 1. 4

Proceedings of the 9th MATHMOD Vienna, Austria, February 21-23, 2018

T. Weinberger et al. / IFAC PapersOnLine 51-2 (2018) 565–570

Table 1. Configuration single-mode prediction frequency 10Hz

q A at x A

phase velocity vphase,1 =

reference solution modal approximation

0.06

121.78 m/s 0.05

0Hz (static)

vphase,2 = vphase,1 0.04

displacement (m)

weightings triangular window λw = λw˙ = λφ = λφ˙ = 1 λθ = 20, λdyn = 0 without λθ , λdyn no imaginary part β2

569

The modal approximation of the displacement field in subdomain A shows a normalized RMS fit of 99.2% for the stationary wave. The mode is predicted to arrive nlag,1 = 821 time steps later at xB with a roundofferror of nroff,1 = 0.1604 time steps. There, the predicted displacement field in the surroundings of xB attains a spatial fit of 98.9%. These excellent fit values demonstrate that, as spatial steps and therefore roundoff-errors are small compared to the wave length, the stationary wave propagation is predicted highly accurately with an efficient modal description.

0.03

0.02

0.01

0

-0.01 0.5

0.7

0.9

1.1

1.3

reference solution predicted solution

0.06

Now, multiple modes are assessed. The beam is loaded with a point force F (t) = 800N · sin(27.4s−1 · 2π · t) · 10 (σ(t) − σ(t − 27.4 s)) + 200N · sin(11.3s−1 · 2π · t) · (σ(t) − 10 σ(t − 11.3 s)) at point xC , comprised of ten sine periods of 800N amplitude at 27.4Hz and ten sine periods of 200N amplitude at 11.3Hz. The parameters of the harmonic excitation were chosen arbitrarily, however deliberately off the frequency grid points to test real-world performance. The modal frequency grid is chosen from 10 to 30Hz spaced by 1Hz and a static offset (0Hz). The phase velocities at 10, 20, 30Hz are measured and interpolated quadratically to parametrize the method, and for simplicity, the static offset’s propagation velocity is set to the lowest phase velocity, see Table 2. In Fig. 6 (top) one can see the spatially fitted displacement qA (for each time step) and the reference solution wref (x = xA , t), plotted over time. When the main wave package propagates through subdomain A the algorithm provides excellent fitting. Also at the transient processes the fitted signal is able to follow the reference solution well, with only small deviations. Fig. 6 (bottom) shows the predicted signal qB versus the reference solution wref (x = xB , t), plotted over time. At some specific times, the prediction deviates from the solution visibly, notably at the initial transient of the wave package around t = 1.6s. Also, some high-frequency leakage artifacts (ringing) in the predicted data is seen around t = 1.3s, which can be suppressed by adjusting the weights λw˙ , λφ˙ to lower, or λdyn to higher values. Setting λθ to higher values also reduces ringing, but also leads to a significant performance loss in fitting the main wave package.

0.05

displacement (m)

0.04

0.03

0.02

0.01

0

-0.01 1.25

1.45

1.65

1.85

2.05

2.25

2.45

time (s)

Fig. 6. Multi-mode case: displacement qA at xA over time (top), displacement qB at xB over time (bottom) 4.3 Parameter tuning The ringing phenomenon at t = 1.3s (compare Fig. 6) motivates the following, deeper investigation. The observed ringing originates from an artifact oscillation of about 200Hz in the signal at qA around t = 0.5s. Because the Least-Squares approach tries to explain the signal with the highest frequency available in the regressor with additional high-frequency oscillations in the modal coordinates, the predicted signal shows significant errors. To counteract this phenomenon, the effects of the weighting parameters are demonstrated with the multi-modal case configuration, compare chapter 4.2. The chosen parameters for the method variants 0, 1, and 2 are listed in Table 3. In method 1 ringing is avoided by decreasing the costs of mis-fits in the velocities, which reduces the modal contributions of the higher frequencies. This effectively suppresses ringing, see detail “ringing” of Fig. 7, with the side effect that the high frequency part at the very beginning of the wave package is not reproduced as

Table 2. Configuration multi-modal prediction frequency 10Hz to 30Hz (1Hz spaced) 0Hz (static)

1.7

q B at x B

4.2 Multi-modal prediction

weightings triangular window λw = λw˙ = 1 λφ = λφ˙ = 1 λθ = 20, λdyn = 0 without λθ , λdyn no imaginary part β22

1.5

time (s)

phase velocity vphase,1 = 121.78 m/s vphase,11 = 122.21 m/s vphase,21 = 122.40 m/s (quadratic interpol.) vphase,22 = vphase,1

5

Proceedings of the 9th MATHMOD 570 Vienna, Austria, February 21-23, 2018

T. Weinberger et al. / IFAC PapersOnLine 51-2 (2018) 565–570

Table 3. Parameters of methods 0, 1, and 2 method 0 1 2

λθ 20 20 20

λ w = λφ = 1 1 1

λw˙ = λφ˙ = 1 0.1 1

Table 4. Computational efficiency ∆t (s) 0.001 0.001 0.001 0.001 0.0005 0.00025

λdyn 0 0 100

∆x (m) 0.0625 0.125 0.25 0.5 0.5 0.5

ref. sim. (s) 205.9 103.3 51.4 27.0 53.9 106.7

runtime (s) 41.5 23.0 12.9 7.9 15.1 29.5

ratio 5.0 4.5 4.0 3.4 3.6 3.6

5. CONCLUSION The achieved results show high accuracy also with transient processes provided that the basis of modes is chosen suitably and the corresponding phase velocities are accurate. The robustness of the method is satisfactory, if parameters, especially the time step size, are set correctly. The introduced weights can be adjusted to avoid ringing prediction errors and thus improve robustness and accuracy. This method is easy to implement and significantly reduces computational effort, which is highly important in real-time-capable systems. The possibility to store data in a shift-register-manner is a significant advantage in data storage when using this method. It would be of interest for prospective studies to correctly feed the wave predictions into the destination domain when an absorbing boundary layer is present. Also, directional modal decomposition is required to construct a full, bi-directional domain decomposition method for parallel real-time computation. REFERENCES Aschauer, G., Schirrer, A., and Jakubek, S. (2016). Realtime-capable finite element model of railway catenary dynamics in moving coordinates. In Control Applications (CCA), 2016 IEEE Conference on, 1078–1083. IEEE. Kim, W.M., Kim, J.T., Kim, J.S., and Lee, J.W. (2003). A numerical study on dynamic characteristics of a catenary. KSME international journal, 17(6), 860–869. Ritzberger, D., Schirrer, A., and Jakubek, S. (2016). Absorbing boundary layer control for linear onedimensional wave propagation problems. Journal of Vibration and Control, 1077546316679579. Ritzberger, D., Talic, E., and Schirrer, A. (2015). Efficient simulation of railway pantograph/catenary interaction using pantograph-fixed coordinates. IFACPapersOnLine, 48(1), 61–66. Schirrer, A., Aschauer, G., Talic, E., Kozek, M., and Jakubek, S. (2017). Catenary emulation for hardwarein-the-loop pantograph testing with a model predictive energy-conserving control algorithm. Mechatronics, 41, 17–28. Schirrer, A., Talic, E., Aschauer, G., Kozek, M., and Jakubek, S. (2016). Optimization based determination of highly absorbing boundary conditions for linear finite difference schemes. Journal of Sound and Vibration, 365, 45–69.

Fig. 7. Parameter tuning to avoid ringing good as before, compare detail “transient performance” of Fig. 7. Method 2 uses the dynamic weight λdyn to penalize too fast variations of modal coordinates over time. As one can see in Fig. 7, ringing is decreased significantly compared to method 0, though performance is not as high as in method 1. Nevertheless, the approach of method 2 is preferred, because the signal at the transient process is fitted even better than with the parameters of method 0, and method 1 performs worst in this aspect.

4.4 Computational efficiency To estimate computational efficiency, the runtime of the proposed approach (FE computation for subdomain A and wave transport prediction operations) is compared to the time needed to compute the corresponding FE solution over a continuous domain (FE only, subdomain x ∈ [xB , xA +7.5 m]). The proposed approach incorporates the multi-modal case configuration as in chapter 4.2, and the measured runtimes are listed for various choices of time and space gridding in Table 4 (column “runtime”). Therein, “ref.sim.” indicates the time to compute the covering continuous domain, and it is evident that the proposed concept significantly reduces computational effort already in the present first implementation. 6