Absorbing boundary layer control for linear one-dimensional wave propagation problems

Abstract

In this work, a boundary layer control scheme for one-dimensional wave propagation problems is presented that provides reflection-less absorption of incident waves in numerical simulations. The desired absorption properties are formulated as an optimization problem. By using the dispersion relation to predict the unbounded wave propagation as a reference trajectory, a constant-gain state-feedback boundary layer controller is found. Since no additional auxiliary variables are introduced, this approach is highly computationally efficient, making it suitable for simulations under real-time requirements. The performance of the boundary layer controller is first evaluated and demonstrated on the scalar wave equation (vibrating string), for which reference absorbing boundary conditions are well established. Afterward, the moving Euler–Bernoulli beam under axial tension is considered, for which a good absorbing performance is achieved.

Keywords

Absorbing boundary control optimal control axially moving materials perfectly matched layer transverse vibrations

1. Introduction

One-dimensional wave propagation problems, described by partial differential equations are extensively used in practical engineering applications, for example, when modeling the vibrations of beams, railway tracks (Vostroukhov and Metrikine, 2003), overhead lines (Ritzberger et al., 2015), or bandsaws and conveyor belts (Tang et al., 2015). Approximate solutions can be obtained by applying numerical techniques, such as the finite-difference method. Thereby, the domain is truncated and suitable boundary conditions must be introduced. By applying absorbing boundary conditions, ideally waves are attenuated without generating reflections at the boundary; thus, the solution of the truncated domain approximates the unbounded wave solution (see Figure 1).

Figure 1.

When the wave is totally absorbed at the boundary, the solution inside the computational domain approximates the free-wave solution.

In early works (Engquist and Majda, 1977; Higdon, 1987), absorbing boundary conditions were derived under the assumption of constant wave propagation speed. These localized absorbing boundary conditions are easy to implement, but if dispersion is present, i.e. the wave propagation speed depends on the spatial frequency, the absorbing performance is severely impaired. Another type of absorbing boundary condition, introduced by Berenger (1994), is the perfectly matched layer (PML). Thereby, the computational domain is surrounded by a layer of finite thickness in which waves are efficiently attenuated. Chew and Weedon (1994) interpreted the PML approach in terms of using complex coordinate stretching, and extended its application to dispersive wave problems.

Although good performance can be achieved with the PML, its implementation relies on problem-specific transformations of the partial differential equation. Thereby, depending on the order of the spatial derivatives, the solution field is split or, in the case of high-order spatial derivatives, several auxiliary partial differential equation are introduced (Hesthaven, 1998). Doing this can lead to long-term instability (Abarbanel et al., 2002) and drastically increases the computational effort.

In this work, the construction of an absorbing boundary layer is perceived as an optimal boundary control problem. The main contribution is a novel method of combining a dispersion-relation-based prediction of the unbounded wave solution, which serves as a desired reference trajectory, and an associated optimal control problem, closely related to the model predictive control theory, with inputs distributed across the boundary layer. Effectively, this method realizes an impedance controller aiming at emulating the physical properties of the unbounded domain. By appropriately introducing the unbounded wave prediction in terms of linear operations, it is shown that online optimization of the control inputs, typical in model predictive control, is avoided and a time-invariant state-feedback controller is obtained.

The proposed control scheme provides the advantage that the partial differential equations under consideration do not need to undergo problem-specific transformations, and therefore no auxiliary fields are introduced. A simple state-feedback controller achieves the desired PML-like behavior instead, making this approach highly computationally efficient and suitable for real-time applications. Since the control design is carried out directly on the discretized system, dispersion due to discretization is inherently considered and does not decrease the absorbing performance. The controller design can be highly automated and is therefore easily adapted to different kinds of linear one-dimensional wave propagation problem.

The absorbing capabilities will be first investigated for the one-dimensional scalar wave equation, since absorbing boundary conditions are well established for this type of problem and the performance in comparison with these reference absorbing boundary conditions can be investigated through numerical simulations. Thereafter, the proposed control scheme is applied to the moving Euler–Bernoulli beam under axial load, for which, to the authors’ best knowledge, alternative absorbing boundary conditions have not yet been devised. In comparison, alternative absorbing boundary conditions for the simplified nonmoving Euler–Bernoulli beam without axial load, such as the classic PML approach derived by Arbabi and Farzanian (2014) or the procedure proposed by Montseny (2011), become highly problem-specific and several auxiliary partial differential equations or additional dynamic state variables are introduced.

2. Fundamentals

In this work, dynamic linear one-dimensional partial differential equations of the form

\sum_{k} a k \frac{\partial k w}{\partial t k} + \sum_{l} b l \frac{\partial l w}{\partial x l} + \sum_{m} \sum_{n} c m, n \frac{\partial m \partial n w}{\partial t m \partial x n} = 0

(1)

are considered for the implementation of a controlled damping layer. Thereby, w(x, t) is a scalar field and the coefficients a_k, b_l, c_m_, _n of the temporal, spatial, and mixed derivative terms are constant.

2.1. Finite-difference approximation

To obtain an approximate solution of a partial differential equation, various numerical techniques can be applied. One possibility is the finite-difference time-domain method. The partial derivatives in equation (1) are thereby approximated by their corresponding central finite differences (Strikwerda, 2007).

By defining a fixed step size Δt and Δx for time and space, the solution is only available at discrete points $w (j Δ t, n Δ x)$ , for which the notation $w_{n}^{j}$ is introduced. The discrete domain, including the controlled damping layer, is truncated to a finite length $N Δ x = L$ . Therefore, a set of N linear algebraic equations is obtained, which can be written as a first-order state-space system

x j + 1 = Ax j

(2a)

w j = Cx j

(2b)

where x ^j is the state vector, A is the system matrix, C is the output matrix, and w ^j is the discrete solution vector of equation (1) at time index j.

2.2. Boundary conditions

Since it is assumed that, owing to the controlled damping layer, all waves are sufficiently attenuated, simple clamped boundary conditions

w (x = 0) = w (x = L) = 0

(3a)

\frac{\partial w}{\partial x} | x = 0 = \frac{\partial w}{\partial x} | x = L = 0 :

(3b)

are applied to the outermost end of the damping layer, where the order of derivatives required at the boundary depends on the highest spatial derivative of the partial differential equation. Boundary conditions are discretized using one-sided finite differences.

2.3. Dispersion relation

For a linear wave-like partial differential equation with constant coefficients, the harmonic wave (as in Hagedorn and DasGupta, 2007)

w (x, t) = e i ω x x e i ω t

(4)

is a fundamental solution of the partial differential equation, where i is the imaginary unit, ω_t is the temporal frequency, and ω_x is the spatial frequency (or wave number). The discrete version of equation (4) can be expressed as

w_{n}^{j} = (e i ω x Δ x) n (e i ω t Δ t) j

(5)

Inserting equation (5) into the discretized partial differential equation yields the so-called discrete dispersion relation (Higdon, 1986), which is generally an implicit function of the form

F (ω t Δ t, ω x Δ x) = 0

(6)

and describes the dependency of ω_x and ω_t. The propagation speed of a fundamental wave (phase velocity) is defined as

v p = - \frac{ω x}{ω t}

(7)

The quantities

ω x Δ x

and

ω t Δ t

are confined to the interval

[- π, π]

. A wavenumber of

ω x Δ x = \pm π

corresponds to a sawtooth mode in space, and higher frequencies cannot be resolved by the grid.

3. Main results

In this section, the objective function that evaluates the absorbing performance of a controller and serves as the base of the optimization problem will be discussed. Afterward, two methods are presented that facilitate a useful solution: an analytical treatment of the optimization problem and numeric optimization using a genetic algorithm.

3.1. Optimization problem formulation

The state-space system (2) is to be controlled so that, based on the idea of complex coordinate stretching (see Appendix 1 for details), an exponential decay is added to the harmonic wave solution (5)

y_{n}^{j} = e - σ (Δ xn) (e i ω x Δ x) n (e i ω t Δ t) j

(8)

With additional control inputs at the nodes in the damping layer (see Figure 2), the state-space system is extended to

x j + 1 = Ax j + Bu j w j = Cx j

(9)

where u ^j is the control input and B the input matrix. A general solution is approximated as the superposition of a finite set of harmonic waves with spatial frequencies inside the discrete spectrum

ω x Δ x \in [0, π]

. The corresponding temporal frequencies are determined by the dispersion relation (6) and each harmonic wave is analytically continued in time according to equation (8), assuming that a suitable damping profile

σ (Δ xn)

is provided. To evaluate the absorbing performance of the state-space system with given inputs, an objective function is formulated. To do so, the damped harmonic wave (equation (8)) is evaluated for a predefined number of time steps, which serve as a reference trajectory. The state-space system is initialized with this harmonic wave and its solution over time is predicted. The error between the state-space system’s trajectory and the damped harmonic reference wave is then summed over all time instances N_p and sampled frequencies Ω

J (u) = \sum_{j = 1}^{N p} \sum_{i = 1}^{Ω} | y j (ω x, i) - w j (u j) |_{2}^{2}

(10)

An optimal control input sequence is sought to minimize equation (10). In the following, an analytic approach to derive and solve the optimization problem will be developed.

Figure 2.

Resulting control scheme for the analytic optimization, N_c = 6. Red nodes are pinned, realizing a clamped boundary condition. The complete state vector is used.

3.2. Analytic solution of the optimization problem

By applying a modal decomposition to the displacement vector w ^j and by using the information of the dispersion relation it is possible to construct a reference solution with the desired nonreflective damping behavior, which only depends on the current state information of time step j. If this reference trajectory can further be constructed by a linear operation, the optimization of equation (10) yields a time-invariant state-feedback control law.

Assumption 1

It is assumed that the displacement vector w^j at time instance j of the state-space system (9) can be approximated by a modal superposition of a finite number of harmonic waves with a spatial wave number $ω x, i, \forall i \in {1, \dots, Ω}$

w j = \sum_{i = 1}^{Ω} L i g_{i}^{j}

(11)

where

L i = [\cos (ω x, i Δ x) \sin (ω x, i Δ x) \cos (2 ω x, i Δ x) \sin (2 ω x, i Δ x) : : \cos (N ω x, i Δ x) \sin (N ω x, i Δ x)]

(12)

and

g_{i}^{j} = [G_{1, i}^{j} G_{2, i}^{j}]

(13)

The two modal amplitudes in the gain vector

g_{i}^{j}

thereby inherently contain the information of phase and amplitude of mode i at time instance j.

Lemma 1

Given Assumption 1, the modal decomposition of the current displacement vector w ^j allows the construction of a future desired reference state with augmented damping at time step $j + n p$ by

y_{ref}^{j + n p} = \sum_{i = 1}^{Ω} D T i, n p g_{i}^{j}

(14)

where

T i, n p = [\cos (ω x, i Δ x + n p ω t, i Δ t) \sin (ω x, i Δ x + n p ω t, i Δ t) \cos (2 ω x, i Δ x + n p ω t, i Δ t) \sin (2 ω x, i Δ x + n p ω t, i Δ t) : : \cos (N ω x, i Δ x + n p ω t, i Δ t) \sin (N ω x, i Δ x + n p ω t, i Δ t)]

(15)

and

D = [e - σ (0) e - σ (Δ x) e - σ (N Δ x)]

(16)

Proof

Applying a temporal shift $j \to j + n p$ to the damped harmonic wave solution (8) yields equations (15) and (16). To determine $T i, n p$ completely, the dispersion relation is evaluated for each specific spatial wave number $ω x, i$ of the finite frequency grid to obtain the corresponding temporal frequency $ω t, i$ .

Remark 1

In Lemma 1, it is assumed that the amplitudes of waves traveling in the reverse direction, from the damping layer into the computational domain, is zero. This is a justified assumption as long as there are no additional disturbances exciting the system inside the damping layer and if the desired attenuation of waves is sufficiently large that spurious reflections at the boundaries become negligible. Therefore, only waves traveling from the interior domain into the damping layer are considered in the prediction.

Theorem 1

Given a system of the form of equation (9), and using Lemma 1, an optimal linear state-feedback control law

u j = Kx j

(17)

can be found, which minimizes the quadratic objective function

J = (Y ref - Y) T (Y ref - Y) + U T RU

(18)

where

Y ref = [y_{ref}^{j + 1} y_{ref}^{j + 2} : y_{ref}^{j + N p}], Y = [w j + 1 w j + 2 : w j + N p], U = [u j u j + 1 : u j + N p - 1]

Thereby, Y denotes the stacked vector of the predicted system trajectory under the aggregated control inputs U and Y_ref contains the desired future reference trajectory with augmented damping at the boundary layer.

Proof

With the extended state-space system (9), the future states based on the current state information are predicted sequentially up to N_p time steps

w j + 1 = CAx j + CBu j w j + 2 = CA 2 x j + CABu j + CBu j + 1 : w j + N p = CA N p x j + CA N p - 1 Bu j + CA N p - 2 Bu j + 1 + \dots + CA 0 Bu j + N p - 1

(19)

Note that the superscript on the input vector u, state vector x, and displacement vector w denotes the time step, whereas for the matrix A it is an exponent. The system of equations is further consolidated into

Y = Fx j + Φ U

(20)

with

F = [CA CA 2 : CA N p]

(21)

Φ = [CB 00 \dots 0 CAB CB 0 \dots 0 CA 2 B CAB CB \dots 0 : : : : CA N p - 1 B CA N p - 2 B CA N p - 3 B \dots CA 0 B]

(22)

Inserting equation (20) into the objective function (18) and setting its derivative with respect to future control inputs U to zero leads to

U = (Φ T Φ + R) - 1 Φ T (Y ref - Fx j)

(23)

which is the solution for the optimal future control input according to the standard unconstrained model predictive control theory (see Rawlings, 2009).

The state vector x ^j can be expressed in terms of modal coordinates using equation (11)

x j = [w j w j - 1] = [L 00 L] [g j g j - 1]

(24)

L = [L 1 L 2 \dots L Ω]

(25)

g j = [g_{1}^{j} g_{2}^{j} : g_{Ω}^{j}]

(26)

and by sequentially constructing the desired reference trajectory according to equation (14)

Y ref = Tg j

(27)

T = [DT 1, 1 DT 2, 1 \dots DT Ω, 1 DT 1, 2 DT 2, 2 \dots DT Ω, 2 : DT 1, N p DT 2, N p \dots DT Ω, N p]

(28)

the optimal control input (equation (23)) becomes

u j = \underset{K gain}{\underset{︸}{S (Φ T Φ + R) - 1 Φ T (T [I 0] - F [L 00 L])}} [g j g j - 1]

(29)

which is a time-invariant feedback control law with respect to the modal gains. Matrix

S = [I N c \times N c 0]

selects the current control input u ^j from the stacked input vector U. Solving equation (24) for the stacked gain vector by a least-mean squares method yields

[g j g j - 1] = \underset{K mod}{\underset{︸}{([L 00 L] T [L 00 L] + α I) - 1 [L 00 L]}} x j

(30)

where α is a chosen small-valued constant that regularizes the least-mean squares solution. A time-invariant state-feedback control law is obtained as

u j = K gain K mod x j

(31)

and thus equation (17) holds.

Remark 2

In the control law (equation (17)) proposed in Theorem 1, the complete state vector of the system (9) is used. A schematic representation of the control scheme is shown in Figure 2.

3.3. Direct optimization under structural constraints

As an alternative, a direct optimization technique (e.g. genetic algorithm) can be used to optimize the objective function (10). This allows for additional structural assumptions about the controller to be made. An example structure that has proven to be useful, as it offers a good trade-off between absorbing performance and number of parameters to be optimized, is

K = [k 1 k N c + 1 k 2 k N c + 2 k N c k 2 N c]

(32)

The controlled input is then given by

u j = {Kx}_{DL}^{j}

(33)

where

x_{DL}^{j} = [w_{3}^{j} : w_{N c}^{j} w_{3}^{j - 1} : w_{N c}^{j - 1}]

(34)

is a reduced state vector that only consists of the nodes to be controlled inside the damping layer. Because of the diagonal sub-matrices of K, the control input at a node is only dependent on its collocated displacement at the current and previous time steps. The controller can therefore be interpreted as an optimized spring and damper for each node (elastic bedding).

A schematic illustration of the resulting control scheme is depicted in Figure 3. By contrast with the analytic optimization shown in Figure 2, a reduced state vector is fed back into the control gain matrix K, which only consists of the nodes inside the damping layer that are to be controlled.

Figure 3.

Schematic damping layer controller obtained by using a genetic algorithm, N_c = 6.

The goal of the genetic optimization is to find the set of parameters of the control gain matrix of equation (32). To achieve this goal, an initial set of solution candidates is created. In each iteration, these candidates are manipulated by selection, mutation, and recombination, and a new generation is formed. Ideally, the objective function is thereby minimized and the population converges after several generations. Owing to its random nature, the final generation differs for each optimization run and convergence to a useful solution is not guaranteed. By contrast, the analytic optimization approach is reproducible in the sense that with the same initial problem setup the same optimal controller is found. The advantage of using the genetic algorithm lies in the fact that the optimization framework can be set up fairly quickly, assuming that a suitable optimization software (e.g. MATLAB Optimization Toolbox™) is available. Furthermore, the structure of the controller can be assessed directly; this can be advantageous under certain circumstances, for example, if a sparse structure is desired to further accelerate the computation (at the cost of absorbing performance) or a certain physical structure, such as an elastic bedding, is to be used.

4. Numerical simulation

At first, the scalar one-dimensional wave equation will be considered, for which absorbing boundary conditions are well established.

Finally, the presented method will be applied to the moving Euler–Bernoulli beam under axial tension, for which no alternative absorbing boundary technique is available.

4.1. Scalar wave equation

The problem considered in this section is the one-dimensional wave equation given by

\frac{\partial 2 w}{\partial t 2} = c 2 \frac{\partial 2 w}{\partial x 2}

(35)

with c as a real-valued constant wave propagation speed.

4.1.1. Dispersion relation

Inserting the harmonic wave solution (4) into equation (35) yields the continuous dispersion relation

ω_{x}^{2} = c 2 ω_{t}^{2}

(36)

or, to write it another way, the phase speed

v p = \frac{ω t}{ω x} = \pm c

(37)

which is constant, equal to c, and symmetric for waves traveling in different directions. Discretizing equation (35) using central finite differences and inserting the discrete harmonic wave equation (5) leads to the discrete dispersion relation

\sin (\frac{ω t Δ t}{2}) 2 = c 2 \frac{Δ t 2}{Δ x 2} \sin (\frac{ω x Δ x}{2}) 2

(38)

For

Δ x, Δ t \to 0

, equation (38) approximates to equation (36). Furthermore, if the spatial and the temporal step size are set according to the condition

λ = c \frac{Δ t}{Δ x} = 1

(39)

where λ is known as the Courant number, the discrete dispersion relation equals the continuous dispersion relation. Figure 4 shows the phase velocity for different Courant numbers with the nominal wave speed c = 1. The phase velocity is constant if the computational grid can be conditioned such that equation (39) holds. For smaller Courant numbers, numerical dispersion is present, which leads to a deviation from the nominal wave speed especially for high spatial frequencies.

Figure 4.

Phase velocities for different Courant numbers λ.

4.1.2. Numerical setup and reference systems

The domain of interest is defined to have N = 200 nodes and absorption is desired at both boundaries. The nominal wave speed is set to c = 1, the time step size to Δt = 0.1 and the spatial step size to Δx = 1. Therefore, the Courant number is λ = 0.1, which causes significant numerical dispersion. To compare the performance of the controlled damping layer, the following reference systems are used:

An enlarged reference system with Dirichlet boundary conditions and N_ref = 20N spatial nodes. The reflection error is assessed with respect to this enlarged reference system, since spurious reflections do not reach the domain of interest during the simulation.

A reference system with first-order Engquist–Majda absorbing boundary conditions (Engquist and Majda, 1977), which are perfectly absorbing for the one-dimensional wave equation if λ = 1. Since λ = 0.1, a decrease of the absorbing performance for high spatial frequency waves is expected.

An unsplit PML formulation is obtained through complex coordinate stretching, which is described by Grote and Sim (2010), and uses a space-time staggered grid with explicit time stepping. Here, however, the same simple collocated grid scheme as in the other reference systems was used. As a result, numerical (long-term) instability was observed. As the unstable dynamics are only relevant at time scales of one order of magnitude greater than the desired simulation time, the comparison is deemed valid here.

4.1.3. Damping layer controller

With the dispersion relation (38), and the state-space system obtained by applying the central finite-difference approximation to equation (35) the analytic optimization can be carried out. The tunable parameters are set according to Table 1.

Table 1.

Optimization parameters.

Parameter	Symbol	Value
Control layer width	N _pml	10
Number of prediction steps	N _p	50
Number of sampling frequencies	Ω	250
Input weighting constant	r	1
Gain mapping weighting constant	α	0.1

The damping profile used in equation (16) is chosen as a quadratic function of the form

σ (n) = {\frac{σ m}{N_{c}^{2}} n 2 - 2 \frac{σ m}{N c} + σ m 1 \leq n \leq N c 0 n > N c

(40)

where σ_m denotes the maximum damping value at the outermost end of the damping layer. It should be chosen so that all waves are sufficiently attenuated and the Dirichlet boundary conditions do not evoke significant spurious reflections. In this simulation the damping parameter was set to σ_m = 20.

A useful property of the resulting control matrix is that the discrete system for which the optimization is carried out does not need to be of the same size as the system on which the controller is later applied. In this case, the system for the optimization consists of N_opt = 100 nodes. When applied to the system with N = 200 nodes only the states present at the optimization are fed back to the controller. Since the dispersion relation is symmetric for the scalar wave equation, the resulting controller can be applied on both boundaries. However, the elements of the controller have to be rearranged to be consistent with the node order of the displacement vector for the opposite boundary.

4.1.4. Simulation results for one-dimensional wave equation

The initial conditions for all system were set to be a pulse of the form

\frac{\partial u (x, 0)}{\partial t} = 0

(41b)

where x_c marks the center of the respective system and d was chosen such that the pulse-width is

d / Δ x = 5

nodes. In Figure 5, the simulation results of the system with controlled damping layer are compared with the solution of the enlarged reference system. The significant numerical dispersion leads to the separation of the high and low frequency content of the initial pulse. No reflections or deviations from the reference solution are visible.

Figure 5.

Wave equation with initial pulse. Comparison of simulation results between controlled damping layer and enlarged reference system.

4.1.5. Error evaluation

To quantify the error in comparison with the enlarged reference system, the following measure is used

E = \sqrt{\frac{1}{N d} \sum_{n = 1}^{N d} (w n - w ref, n) 2}

(42)

where d is the number of nodes inside the domain, w_n is the solution of the system under investigation and w_ref, _n the reference solution at node n.

Figure 6 shows the error over time for the controlled damping layer in comparison with the Engquist–Majda absorbing boundary condition, the PML obtained through complex coordinate stretching, and the clamped boundary. Since the finite-difference time-domain method used to obtain a numerical solution is explicit, the results for the clamped boundary, the absorbing boundary condition, and the PML do not deviate from the reference solution significantly until the wave front impacts the boundary (t ≈ 90 s). For the controlled damping layer, a deviation can be observed much earlier. This is because the inputs for the damping layer controller are the displacements of nodes that reach far into the computational domain. Therefore, slight control action occurs before the wave package reaches the interface. The majority of the wave package has left the computational domain at t ≈ 300 s. The Engquist–Majda absorbing boundary conditions suffer from the numerical dispersion caused by the Courant number λ = 0.1 and produce significant reflections, especially in the high spatial frequency range. The PML obtained through complex coordinate stretching shows a good performance at the impact of the wave package but, as mentioned earlier, the simple discretization scheme has a negative effect on long-term stability and the error starts to increase with time. The controlled damping layer, although deviating sooner from the reference solution, shows a good absorbing performance by having a smaller reflection error than the Engquist–Majda absorbing boundary condition and the PML.

Figure 6.

Error evaluated for the different absorbing boundary methods over time for the initial pulse simulation. ABC: absorbing boundary conditions; PML: perfectly matched layer.

4.2. Moving Euler–Bernoulli beam

In this test case, the controlled damping layer will be applied on the moving Euler–Bernoulli beam equation (Wickert and Mote, 1990)

ρ A \frac{\partial 2 w}{\partial t 2} = - EI \frac{\partial 4 w}{\partial x 4} + (T - ρ Av 2) \frac{\partial 2 w}{\partial x 2} + 2 v ρ A \frac{\partial 2 w}{\partial x \partial t}

(43)

4.2.1. Deriving the state-space system representation

After discretization, an algebraic equation for each grid point is obtained that describes the dependency to its surrounding nodes in space and time

a T [w_{n - 1}^{j + 1} w_{n}^{j + 1} w_{n + 1}^{j + 1}] = b T [w_{n - 2}^{j} w_{n - 1}^{j} w_{n}^{j} w_{n + 1}^{j} w_{n + 2}^{j}] + c T [w_{n - 1}^{j - 1} w_{n}^{j - 1} w_{n + 1}^{j - 1}]

(44)

Owing to the mixed derivative term being approximated by central finite differences, the resulting set of algebraic equations must be solved implicitly. The coefficient vectors are

a = [\frac{v ρ A}{2 Δ x Δ t} \frac{ρ A}{Δ t 2} + \frac{β}{2 Δ t} - \frac{v ρ A}{2 Δ x Δ t}] b = [- \frac{EI}{Δ x 4} 4 \frac{EI}{Δ x 4} + \frac{T - ρ Av 2}{Δ x 2} - 6 \frac{EI}{Δ x 4} - 2 \frac{T - ρ Av 2}{Δ x 2} + 2 \frac{ρ A}{Δ t 2} 4 \frac{EI}{Δ x 4} + \frac{T - ρ Av 2}{Δ x 2} - \frac{EI}{Δ x 4}] c = [\frac{v ρ A}{2 Δ x Δ t} \frac{ρ A}{Δ t 2} - \frac{β}{2 Δ t} - \frac{v ρ A}{2 Δ x Δ t}]

(45)

Aggregating equation (44) for every spatial node leads to a second-order system of difference equations

A so w j + 1 = B so w j + C so w j - 1

(46)

Matrices A_so and C_so are tridiagonal and B_so is pentadiagonal with the entries corresponding to equation (45). By introducing

x j = [w j, w j - 1] T

, Equation (46) can be cast into a first-order state-space system

x j + 1 = \underset{A}{\underset{︸}{[A_{so}^{- 1} B so A_{so}^{- 1} C so I O]}} x j w j = \underset{C}{\underset{︸}{[I O]}} x j

(47)

4.2.2. Dispersion relation

As with the numeric example of the scalar wave equation above, inserting equation (5) into equation (44) yields the so-called discrete dispersion. After simplification, the following implicit function for the moving Euler–Bernoulli beam is obtained

\sin (\frac{ω t Δ t}{2}) 2 = (\frac{T}{ρ A} - v 2) \frac{Δ t 2}{Δ x 2} \sin (\frac{ω x Δ x}{2}) 2 + \frac{EI Δ t 2}{ρ A Δ x 4} (4 \sin (\frac{ω x Δ x}{2}) 2 - \sin (ω x Δ x) 2) + \frac{v Δ t}{2 Δ x} \sin (ω t Δ t) \sin (ω x Δ x)

(48)

4.2.3. Simulation results for the moving Euler–Bernoulli beam

The beam parameters are set as shown in Table 2. The beam velocity v is significant compared with the wave propagation speed c

c = \sqrt{\frac{T}{ρ A}} \approx 3 v

(49)

causing a strong asymmetry in the dispersion relation (see Figure 7). For a spatial step size Δx = 0.6, the maximum stable time step size was found to be

Δ t max = 4.99 \times 10 - 3

by evaluating the eigenvalues of the discrete state-space system and applying a bisectional method to adapt the time step. The temporal step size used in the simulation is

Δ t = 5 \times 10 - 4

. The optimization parameters are set again as in Table 1, except with the prediction horizon increased to N_p = 200. Initial conditions are again defined by equation (41).

Table 2.

System parameters.

Parameter	Symbol	Value	Unit
Length	L	120	m
Tensile force	T	20 × 10³	N
Specific mass	ρA	1.35	kg/m
Bending stiffness	EI	150	N m²
Beam velocity	v	40	m/s

Figure 7.

Normalized phase velocity plot for the scalar wave (blue), nonmoving Euler–Bernoulli beam (EBB) (red) and moving Euler–Bernoulli beam (black).

Figure 8 gives a visual comparison between the system with controlled damping layer and an enlarged reference system. The direction-dependent propagation speed of waves is evident. A good agreement with no visible reflections at the interface can be observed. In Figure 9, the reflection errors according to equation (42) for the controlled damping layer obtained through analytic optimization and the control layer using a genetic algorithm are compared. The optimization parameters, namely the prediction horizon and the number of sampled spatial frequencies, are chosen to be equal. The analytic solution achieves a higher absorption performance than the structured genetic algorithm solution. The simple diagonal substructure (equation (32)) used in this example gives a reasonable trade-off between performance and robustness of the optimization. More complex structures can lead to better performance, as there are more degrees of freedom. However, with an increasing number of parameters to be optimized, convergence to a practical solution is more unlikely and time-consuming. Since for the numerically optimized controlled damping layer only the states inside the layer itself are fed into the gain matrix, a significant deviation from the reference system occurs when the wave package passes over the interface. By contrast, only a small deviation for the analytically optimized control layer is evident from the start.

Figure 8.

Moving Euler–Bernoulli beam with initial pulse. Comparison of the simulation results between controlled damping layer (black) and the enlarged reference system (green).

Figure 9.

Error with respect to the enlarged reference system for the controlled damping layer obtained through analytic optimization and genetic algorithm.

5. Conclusions

In this work, the construction of absorbing boundary conditions is interpreted as an optimal boundary control problem and solved using modal decomposition and modal prediction of the unbounded solution. As a result, a time-invariant state-feedback control law is obtained that emulates the absorbing properties of a PML. The scalar wave equation is investigated, for which reference absorbing boundary conditions are well known, and the effectiveness of the proposed method in comparison with standard absorbing boundary conditions is demonstrated in simulations. Additionally, extensive transformations of the original problem and the introduction of auxiliary partial differential equations are avoided. The same method is then applied to the moving Euler–Bernoulli beam equation under axial tension, for which alternative absorbing boundary conditions are not available. The proposed optimal control method again achieves high absorbing performance and computational efficiency. Since the control design uses the discretized system, numerical dispersion due to discretization is inherently considered and does not decrease the performance, as is the case for absorbing boundary conditions, which are designed in continuous time. The ease of adaptability, design, and low computational effort, combined with good absorbing performance, makes this approach highly valuable for time-critical simulations of unbounded one-dimensional wave propagation problems.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Appendix 1

References

Abarbanel

Gottlieb

Hesthaven

(2002) Long time behavior of the perfectly matched layer equations in computational electromagnetics. Journal of Scientific Computing 17(1–4): 405–422.

Arbabi F and Farzanian M (2014) Propagation of waves in infinite beams: PML approach. In: Proceedings: 11th world congress on computational mechanics, Barcelona, Spain, 20–25 July 2014. Barcelona: IACM.

Berenger

(1994) A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics 114(2): 185–200.

Chew

Weedon

(1994) A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates. Microwave and Optical Technology Letters 7(13): 599–604.

Chew

Jin

Michielssen

(1997) Complex coordinate stretching as a generalized absorbing boundary condition. Microwave and Optical Technology Letters 15(6): 363–369.

Engquist

Majda

(1977) Absorbing boundary conditions for numerical simulation of waves. Proceedings of the National Academy of Sciences 74(5): 1765–1766.

Grote

Sim

(2010) Efficient PML for the wave equation. arXiv preprint arXiv:1001.0319.

Hagedorn

DasGupta

(2007) Vibrations and Waves in Continuous Mechanical Systems, Hoboken, NJ: John Wiley & Sons.

Hesthaven

(1998) On the analysis and construction of perfectly matched layers for the linearized Euler equations. Journal of Computational Physics 142(1): 129–147.

10.

Higdon

(1986) Absorbing boundary conditions for difference approximations to the multidimensional wave equation. Mathematics of Computation 47(176): 437–459.

11.

Higdon

(1987) Numerical absorbing boundary conditions for the wave equation. Mathematics of Computation 49(179): 65–90.

12.

Hueppe

Kaltenbacher

(2012) Stable matched layer for the acoustic conservation equations in the time domain. Journal of Computational Acoustics 20(01): 1250004.

13.

Kaltenbacher

Sim

(2013) A modified and stable version of a perfectly matched layer technique for the 3-D second order wave equation in time domain with an application to aeroacoustics. Journal of Computational Physics 235: 407–422.

14.

Montseny

(2011) Diffusive wave-absorbing control: example of the boundary stabilization of a thin flexible beam. Journal of Vibration and Control 18: 1708–1721.

15.

Rawlings

(2009) Modal Predictive Control: Theory and Design, Madison, WI: Nob Hill Publishing.

16.

Ritzberger

Talic

Schirrer

(2015) Efficient simulation of railway pantograph/ catenary interaction using pantograph-fixed coordinates. Vienna International Conference on Mathematical Modelling 8: 61–66.

17.

Strikwerda

(2007) Finite Difference Schemes and Partial Differential Equations, Philadelphia, PA: SIAM.

18.

Tang

Zhang

Gao

(2015) Vibration characteristic analysis and numerical confirmation of an axially moving plate with viscous damping. Journal of Vibration and Control. Epub ahead of print 18 May. DOI: 10.1177/1077546315586311.

19.

Vostroukhov

Metrikine

(2003) Periodically supported beam on a visco-elastic layer as a model for dynamic analysis of a high-speed railway track. International Journal of Solids and Structures 40(21): 5723–5752.

20.

Wickert

Mote

(1990) Classical vibration analysis of axially moving continua. Journal of Applied Mechanics 57(3): 738–744.