- Email: [email protected]

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

T. Weinberger ∗∗ A. Schirrer ∗∗ T. Weinberger ∗∗ A. Schirrer ∗∗ T. T. Weinberger Weinberger A. A. Schirrer Schirrer ∗ ∗ Institute of Mechanics and Mechatronics, Vienna University of Institute of Mechanics and Mechatronics, Vienna ∗ ∗ Technology, Institute of of Mechanics Mechanics and Mechatronics, Mechatronics, Vienna University University of of Austria, (e-mail: [email protected]) Institute and Vienna University of Technology, Austria, (e-mail: [email protected]) Technology, Technology, Austria, Austria, (e-mail: (e-mail: [email protected]) [email protected]) Abstract: Considering a linear one-dimensional wave propagation problem, a technique to Abstract: Considering a linear one-dimensional wave propagation problem, aa technique to Abstract: Considering linear one-dimensional wave problem, to predict wave transport a between computational subdomains is proposed. Abstract: Considering aacross lineargaps one-dimensional wave propagation propagation problem, a technique techniqueThe to predict wave transport across gaps between computational subdomains is proposed. The predict wave transport across gaps between computational subdomains is proposed. The observed waves at the “source” domain are disassembled into a modal description, and their predict wave transport across gaps between computational subdomains is proposed. The observed waves at the “source” domain are disassembled into aa modal description, and their observed waves the “source” are disassembled into description, their propagation intoat domain predicted using system’s dispersionand relation. observed waves atthe the“destination” “source” domain domain areis disassembled intothe a modal modal description, and their propagation into the “destination” domain is predicted using the system’s dispersion relation. propagation intoisthe the “destination” domain is predicted predicted using the the system’s system’s dispersion relation. This approach extended into a domain multi-modal wave prediction functionality which relation. enables propagation into “destination” is using dispersion This approach is extended into a multi-modal wave prediction functionality which enables This approach extended a multi-modal wave functionality which enables efficient domainis problems. Simulation results on This approach isdecomposition extended into intofor a time-critical multi-modal computation wave prediction prediction functionality which enables efficient domain decomposition for time-critical computation problems. Simulation results on domain decomposition for time-critical computation problems. Simulation results on aefficient railway catenary’s contact wire simulation highlights the attained high accuracy of the efficient domain decomposition for time-critical computation problems. Simulation results on a railway catenary’s contact wire simulation highlights the attained high accuracy of the a railway catenary’s contact wire simulation highlights the attained high accuracy of the prediction. This technique can serve as a cornerstone for the development of real-time control a railway catenary’s contact wire simulation highlights thedevelopment attained high accuracy control of the prediction. This technique can serve as a cornerstone for the of real-time prediction. This technique can serve serve as decomposition cornerstone for for the development development of real-time real-time control control tasks or anyThis applications where domain is applied to reduce computational costs. prediction. technique can as aa cornerstone the of tasks or any applications where domain decomposition is applied to reduce computational costs. tasks tasks or or any any applications applications where where domain domain decomposition decomposition is is applied applied to to reduce reduce computational computational costs. costs. © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Wave transport, Modal analysis, Domain decomposition, Dispersion relation Keywords: Wave Wave transport, Modal Keywords: Modal analysis, analysis, Domain Domain decomposition, decomposition, Dispersion Dispersion relation relation Keywords: Wave transport, transport, Modal analysis, Domain decomposition, Dispersion relation 1. INTRODUCTION tion for such problems. First, the solution (displacement 1. INTRODUCTION INTRODUCTION tion for for such such problems. problems. First, First, the the solution solution (displacement (displacement 1. tion field) within designatedFirst, “source” is disassem1. INTRODUCTION tion for suchaa problems. the subdomain solution (displacement field) within designated “source” subdomain is disassemdisassemfield) within aa designated subdomain is bled into a modal basis of“source” several harmonic functions that field) within designated “source” subdomain is disassembled into a modal basis of several harmonic functions that Wave propagation problems frequently occur in engineer- bled into aatravelling modal basis several harmonic functions that represent waveof components in a relevant range of bled into modal basis of several harmonic functions that Wave propagation problems frequently occur in engineerrepresent travelling wave components in a relevant range of Wave propagation problems frequently in ing, and solving them accurately andoccur efficiently, pos- represent travelling wave components in aa relevant range of frequencies. These modal contributions are being reassemWave propagation problems frequently occur in engineerengineerrepresent travelling wave components in relevant range of ing, and solving them accurately and efficiently, posfrequencies. These modal contributions are being reasseming, and solving them accurately and efficiently, possibly under real-time constraints, is challenging. In the frequencies. contributions are reassembled at the These distantmodal “destination” domain corresponding ing, and solving themconstraints, accurately isand efficiently,In posfrequencies. These modal contributions are being being reassemsibly under real-time challenging. the bled at the distant “destination” domain corresponding sibly under constraints, is the development of high-speed train catenaries and In pantoat distant “destination” domain corresponding to their specific propagation velocities, compare Fig. 1. sibly under real-time real-time constraints, is challenging. challenging. the bled bled at the the distant “destination” domain corresponding development of high-speed high-speed train catenaries catenaries and In pantoto their their specific propagation velocities, compare Fig. 1. 1. development of train and pantograph current collectors, accurate and efficient models of to specific propagation velocities, compare Fig. To the best of the authors’ knowledge such multi-modal development ofcollectors, high-speed train and catenaries and pantoto their specific propagation velocities, compare Fig. 1. graph current accurate efficient models of To the best of the authors’ knowledge such multi-modal graph current collectors, accurate and efficient models of the catenary crucial (Kim et al. (2003), the the knowledge multi-modal transportation characterised by such disassembling the graph current dynamics collectors, are accurate and efficient models of To To the best best of of method the authors’ authors’ knowledge multi-modal the catenary catenary dynamics are crucial (Kim et al. al. (2003), transportation method characterised by such disassembling the the are crucial (2003), Aschauer et al.dynamics (2016)). In the same(Kim field et control engi- transportation method characterised by disassembling the wave into modal parts, forward them separately over a the catenary dynamics are crucial (Kim et al. (2003), method characterised by separately disassembling the Aschauer et et al. al. (2016)). (2016)). In In the the same same field field control control engiengi- transportation wave into modal parts, forward them over a Aschauer neers developed predictive control algorithms to realistiwave into modal parts, forward them separately over gap in space and reassembling them for further usage is a Aschauer et al. (2016)). In the same field control engiwave into modal parts, forward them separately over a neers developed predictive control algorithms to realistigap in space and reassembling them for further usage is a neers developed predictive control algorithms to cally emulate complex catenary dynamics at a pantograph in space and them further usage is aa novel approach. Itreassembling fully captures thefor system’s dispersion neers developed predictive control algorithms to realistirealisti- gap gap in space and reassembling them for further usage is cally emulate complex catenary dynamics at a pantograph novel approach. approach. It It fully fully captures captures the the system’s system’s dispersion dispersion cally emulate complex catenary at aa pantograph test rig, so that pantographs ofdynamics high speed can be novel for efficient accurate domain cally emulate complex catenary at trains pantograph novel approach.and It allows fully captures theand system’s dispersion test rig, rig, so that that pantographs ofdynamics high speed speed trains can be be characteristics characteristics and allows for efficient and accurate domain test so pantographs of high trains can thoroughly tested in the lab (Schirrer et al. (2016), Schirrer characteristics and allows for efficient and accurate domain decomposition in the considered frequency range of intertest rig, so that pantographs of high speed trains can be and allows for efficient and accurate domain thoroughly tested tested in in the the lab lab (Schirrer (Schirrer et et al. al. (2016), (2016), Schirrer Schirrer characteristics decomposition in the considered frequency range of interthoroughly et al. (2017)). The controllers compute the displacement decomposition in the considered frequency range of interest. Moreover, due to the causal prediction architecture, thoroughly tested in the lab (Schirrer et al. (2016), Schirrer decomposition in the considered frequency range of interet al. (2017)). The controllers compute the displacement est. Moreover, due to the causal prediction architecture, et al. compute the of the(2017)). contact The wirecontrollers using a numeric of the est. Moreover, due to the prediction transient characteristics arecausal preserved well. architecture, et al. (2017)). The controllers computesimulation the displacement displacement est. Moreover, due to the causal prediction architecture, of the contact wire using a numeric simulation of the transient characteristics characteristics are are preserved preserved well. well. of the wire aa and numeric simulation the catenary, in which theusing carrier contact wires are of repreof the contact contact wire numeric simulation the transient transient characteristics are preserved catenary, in which which theusing carrier and contact wires are are of repreThe method proposed here serves as well. an initial building catenary, in the carrier and contact wires represented by axially pre-tensioned Euler-Bernoulli beams in a catenary, whichpre-tensioned the carrier and contact wiresbeams are repreThe method method proposed proposed here here serves serves as as an an initial initial building building sented by by in axially Euler-Bernoulli in aa The block to construct a full, bi-directional domain sented axially pre-tensioned Euler-Bernoulli beams in The method proposed here serves as initialdecompobuilding spatial subdomain surrounding the pantograph in moving sented by axially pre-tensioned Euler-Bernoulli beams in a block to construct a full, bi-directional an domain decompospatial subdomain surrounding the pantograph in moving block to construct aa full, bi-directional domain decomposition method for wave problems with relevant dispersion spatial subdomain the pantograph moving block to construct full, bi-directional domain decompocoordinates and insurrounding real-time (Ritzberger et al.in (2015)). spatial subdomain surrounding the pantograph in moving sition method method for for wave wave problems problems with with relevant relevant dispersion dispersion coordinates and and in in real-time real-time (Ritzberger (Ritzberger et et al. al. (2015)). sition However, the required (bi-directional coordinates sition for wave problemsextensions with relevant dispersion Thereby, absorbing layers were established for effects. coordinates and in boundary real-time (Ritzberger et al. (2015)). (2015)). effects.method However, the required required extensions (bi-directional Thereby, absorbing boundary layers were established for effects. However, the extensions (bi-directional modal disassembling, combination of absorption layers and Thereby, absorbing boundary layers were established for effects. However, the required extensions (bi-directional reflection-less absorption of travelling waves leaving the Thereby, absorbing boundary layers were established for modal disassembling, combination of absorption layers and reflection-less absorption of travelling travelling waves leaving the the modal disassembling, of absorption layers and boundary excitation) combination are not treated here and are subject reflection-less absorption of waves leaving modal disassembling, combination of absorption layers and computational subdomain (Ritzberger et al. (2016)). reflection-less absorption of travelling et waves leaving the boundary excitation) are not treated here and are subject computational subdomain (Ritzberger al. (2016)). boundary excitation) are not treated here and are subject of ongoing research. computational et al. (2016)). excitation) are not treated here and are subject However, long subdomain high-speed (Ritzberger trains may simultaneously use boundary computational subdomain (Ritzberger et al. (2016)). of ongoing research. However, long long high-speed high-speed trains trains may may simultaneously simultaneously use use of ongoing research. However, of ongoing research. several pantographs which trains are positioned far apart. These However, long high-speed may simultaneously use Section 2 defines the wave transport problem, Section 3 deseveral pantographs which are positioned far apart. These Section 2 defines the wave wave transport transport problem, problem, Section Section 33 dedeseveral pantographs which are positioned far apart. These long domains entirely spanning such train configurations several pantographs which are positioned far apart. These Section 2 defines the tails the wave transport Section long domains domains entirely entirely spanning spanning such such train train configurations configurations Section 2wave defines the waveprediction transportmethod, problem,and Section 3 de-44 tails the transport prediction method, and Section long currently cannot be computed in real-time with high aclong domains entirely spanninginsuch train configurations the wave method, and shows simulation resultsprediction to illustrate the method’s high44 currently cannot be computed computed real-time with high high ac- tails tails wave transport transport method, and Section Section showsthe simulation resultsprediction to illustrate illustrate the method’s method’s high currently cannot in with curacy. However, several non-overlapping of shows currently cannot be be computed in real-time real-timesubdomains with high acacsimulation results to the high accuracy, efficiency, and tuning possibilities. curacy. However, several non-overlapping subdomains of shows simulation results to illustrate the method’s high accuracy, efficiency, and tuning possibilities. curacy. However, several non-overlapping subdomains of computation around each pantograph could be modeled curacy. However, several non-overlapping subdomains of accuracy, efficiency, and tuning possibilities. computation around each pantograph could be modeled accuracy, efficiency, and tuning possibilities. computation each pantograph could be and computedaround in parallel meet real-time computation each to pantograph could computation be modeled modeled and computed computedaround in parallel parallel to meet real-time real-time computation and in to meet computation requirements at significantly higher fidelity. Then, an efand computedatinsignificantly parallel to meet real-time computation requirements higher fidelity. Then, an efefrequirements significantly higher fidelity. Then, an ficient way to at model the subdomain interaction is needed requirements at significantly higher fidelity. Then, an efficient way to model the subdomain interaction is needed ficient way the interaction is which motivates the approach presented in this work. ficient way to to model model the subdomain subdomain interaction is needed needed which motivates motivates the approach approach presented in this this work. work. which the presented in The main contribution of this work is a method to effiwhich motivates the approach presented in this work. The main main contribution contribution of of this this work work is is aa method method to to effieffiThe ciently formulate the transport of travelling waves across The main contribution of this work is a method to efficiently formulate the transport of travelling waves across ciently formulate transport of waves across spatial gaps in thethe considered problem domains, serving as ciently formulate the transport of travelling travelling waves across spatial gaps in the considered problem domains, serving as Fig. 1. Railway catenary; subdomains at pantographs spatial gaps in the considered problem domains, serving as key functionality for an envisaged full domain decomposispatial gaps in the considered problem domains, serving as Fig. 1. Railway catenary; subdomains at pantographs key functionality for an envisaged full domain decomposikey Fig. 1. 1. Railway Railway catenary; catenary; subdomains subdomains at at pantographs pantographs key functionality functionality for for an an envisaged envisaged full full domain domain decomposidecomposi- Fig. Copyright 2018 IFAC 1 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2018, IFAC (International Federation of Automatic Control) Copyright © 2018 IFAC 1 Peer review© of International Federation of Automatic Copyright 2018 IFAC 1 Copyright ©under 2018 responsibility IFAC 1 Control. 10.1016/j.ifacol.2018.03.095

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 countat 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