# Matching Boundary Conditions for the Euler–Bernoulli Beam

Authors: Yaoqi Feng; Xianming Wang
Journal: Shock and Vibration
 (2021)

---

## Abstract

Artificial boundary conditions play a crucial role in the dynamic simulation of infinite Euler–Bernoulli beams. In this paper, a class of artificial boundary conditions, matching boundary conditions (MBCs), is presented to provide effective absorption of incident waves in numerical simulations of the Euler–Bernoulli beam. First, matching boundary conditions are proposed based on the space central difference scheme of the Euler–Bernoulli beam, and then, the specific coefficients of MBCs are determined by matching the dispersion relation. Moreover, reflection coefficient study and numerical tests are carried out to analyze the effectiveness of the proposed MBCs, indicating a remarkable agreement. Taken together, the proposed boundary conditions herein can absorb dispersive waves efficiently and are more compact than previous artificial boundary conditions, particularly suitable for real-time simulation.


---

## Body

## 1. Introduction

The Euler–Bernoulli beam is one of the most important models in the study of one-dimensional solid structures, which is widely applied in engineering, such as railway and bridge engineering [1, 2]. The beam always has a finite length in most real practical applications. However, the infinite beam models still exist in some special problems, for instance, the development of multiscale algorithms for beams [3–5] and the study of beams’ response subjected to moving load [6, 7]. With the development of computer technology, numerical simulation is increasingly used in the study of wave propagation in infinite media. Since numerical calculation could only be performed in a finite domain, it is necessary to introduce some methods, such as artificial boundary conditions, to deal with the boundary of domain, so that the simulation of the finite domain is close to the infinite domain. The artificial boundary conditions need to be treated carefully for simulations of an infinite beam model. Otherwise, large spurious reflections may severely pollute the solutions and, sometimes, even cause numerical instability.Over the last few decades, artificial boundary conditions have received considerable attention by several researchers. In the early stage, the artificial boundary conditions proposed by Engquist, Majda [8–10], and Higdon [11, 12] mainly dealt with the acoustic and elastic wave. The wave propagation speed in these problems is constant. Among these boundary conditions, the most common one-dimensional artificial boundary condition is the single-wave equation, which captures the characteristics of constant wave speed. However, the absorbing performance of these boundary conditions is severely impaired for dispersive media. Then, the perfectly matched layer (PML), introduced by Berenger et al. [13], provided a new idea for artificial boundary processing, which means to set a suitable damping layer outside the calculation area. The method was improved and extended successfully to solve dispersive wave problems by Chew and Weedon [14] and showed good performance.However, in the case of high-order spatial derivatives, such as the Euler–Bernoulli beam equation, the realization of the PML needs some special operations [15] which may lead to long-term instability [16] and drastically increase the computational effort. Besides, these boundary conditions are based on continuous governing equations. As an alternative, Tang et al. [4] regarded the semidiscretized Euler–Bernoulli beam equation as an atomic lattice with a three-atom potential and then designed two artificial boundary conditions. The first one is an exact boundary condition via using the time history treatment. Nevertheless, it is sensitive to numerical implementations. Thereupon, the Alex (Almost Exact) boundary condition was improved based on this boundary condition, which is numerically more effective but complex. In addition, this boundary condition is not suitable in actual simulation, which will be described in detail later. Another type of absorbing boundary condition, introduced by Ritzberger [17], was the absorbing boundary control method. The main idea of this method is combining a dispersion-relation-based prediction of the unbounded wave solution. Although this boundary condition has a good effect on the moving tensile Euler–Bernoulli beam, it is very complicated to use because the coefficient of the boundary condition needs to be optimized according to different problems. Besides, the Euler–Bernoulli beam problem in this article is more dispersive.In this paper, we design a class of matching boundary conditions for the Euler–Bernoulli beam by matching the dispersion relation, which was firstly proposed for molecular dynamics simulation [18–21]. The proposed boundary conditions are expressed as a linear constraint for the displacement and velocity at nodes near the artificial boundary. These nodes are derived from the semidiscrete scheme of the Euler–Bernoulli beam equation; hence, the MBCs are convenient in simulation. The combination coefficients are determined by matching the dispersion relation with a minimal number of nodes involved. The goal of this method is to optimize the spurious reflection suppression and the locality in space and time. That is to improve the accuracy and reduce computational costs. Several nodes near the boundary usually suffice for accuracy requirements in most applications.The outline of this paper is as follows. InSection 2, we firstly give a brief description of the dispersion relation for both the continuous and semidiscretized Euler–Bernoulli beam equation, and then, we construct the MBCs for the Euler–Bernoulli beam. After that, we present the reflection coefficient analysis to verify the effectiveness of the proposed boundary conditions in Section 3. Numerical studies in Section 4 further demonstrate the effectiveness and stability of MBCs. Finally, we make some concluding remarks and future directions in Section 5.


## 2. Design of Matching Boundary Conditions

In this section, we design a class of matching boundary conditions, which is convenient and effective for flexural wave propagation problems in the infinite Euler–Bernoulli beam.


### 2.1. Dispersion Relation of the Euler–Bernoulli Beam Equation and Its Difference Scheme

The following shows the Euler–Bernoulli beam vibration equation without loads:(1)EIuxxxx+ρAutt=0.Here,ρ, A, E, and I represent the material density, the area of the beam cross section, Young’s modulus, and the bending moment of inertia, respectively. For a uniform beam, all these parameters are constant, and the equation can be rewritten as follows:(2)utt+EI/ρAuxxxx=0.One common spatially discrete model of the beam equation in application is the central finite-difference scheme,(3)d2ujdt2+EIρAuj−2−4uj−1+6uj−4uj+1+uj+2Δx4=0,which can be rewritten as(4)d2ujdt2=−β2Δx4uj−2−4uj−1+6uj−4uj+1+uj+2,with β=EI/ρA.We consider the wave form solution(5)u=eiknΔx+ωt.Substituting it into the semidiscrete scheme (3), we obtain a relation between the wave number k and the frequency ω, which gives the dispersion relation.(6)ω=±4βΔx2sin2kΔx2.The dispersion relation of the continuous Euler–Bernoulli beam (2) reads(7)ω=±βk2.According to the abovementioned dispersion relations, we can get the group velocitycg and the phase velocity cp, which represent the wave propagation law of the Euler–Bernoulli beam. The wave velocities are shown as follows:(8)cp=βk,cg=2βk.The wave propagating speed is obviously dependent on the wave number, which indicates the Euler–Bernoulli beam is a dispersive medium. To further discuss the difference of two dispersion relations, we plot the wave velocity versus wave number of each model in Figure1, with β = 1, Δx = 0.05. Also, the wave velocities of the discrete model are shown as follows:(9)cp=4βkΔx2sin2kΔx2,cg=2βΔxsinkΔx.Figure 1
Wave velocity in the Euler–Bernoulli beam (β = 1, Δx = 0.05): (a) phase velocity; (b) group velocity.
(a)(b)We notice that the deviation between the dispersion relation between the discrete model and the continuous model increases with the wave number. Hence, the long wave is a primary concern in designing matching boundary conditions. The boundary condition mentioned by Tang et al. [4] has a good absorption effect for short waves. For long waves such as k = 5 or longer with β = 1, Δx = 0.05, the absorption effect of their boundary condition will become worse. In addition, this boundary condition cannot be adjusted according to the wave number; that is, the absorption of this boundary condition for long waves cannot be improved. Therefore, it is necessary to propose a class of boundary conditions that can absorb long waves and can be adjusted for different wave numbers.


### 2.2. Matching Boundary Conditions

We construct boundary conditions by matching the dispersion relation of Euler–Bernoulli beams. To make the boundary conditions more convenient, we will consider the discrete scheme. Since the difference equation (3) uses a fourth-order finite-difference scheme in the spatial domain, it is necessary to design two-layer boundary conditions, as shown in Figure 2.Figure 2
Boundary condition design diagram: (a) schematic diagram of the surface layer boundary condition; (b) schematic diagram of the internal layer boundary condition at the left end.
(a)(b)It is naturally to take the same form for both the two artificial boundary layers. Besides, the same form makes the boundary conditions more compact. Equation (10) is the boundary condition of the surface layer in the left side, and equation (11) is the internal layer’s boundary condition. Therefore, we just need to determine the coefficients in the outer boundary layer conditions.(10)∑j=0Najdujdt+∑j=0Nbjuj=0,(11)∑j=1N+1ajdujdt+∑j=1N+1bjuj=0.Among them, the linear combination coefficientsajs and bjs are determined by matching the dispersion relation (7). Specifically, we use the following linear algebraic equations to calculate the linear combination coefficients. After matching the special wave numbers, the proposed boundary conditions (10) and (11) are named as MBCN (k1, ..., kN).(12)a0=1,Δ0=0,Δkq=0,q=1,2,...,N,where Δk is called the frequency domain residual function, which is defined by(13)Δk=iωk∑j=0NajeijkΔx+∑j=0NbjeijkΔx,where ω=βk2; β=EI/ρA.The first row is to give a fixed value for the homogeneous form of the boundary condition. The second row is the minimum “consistent” requirement. The third row means matching the dispersion relation at the selected wave number(s). Splitting it into the real and imaginary parts, we obtain the following two real-valued equations:(14)−ωkq∑j=0NajsinjkqΔx+∑j=0NbjcosjkqΔx=0,ωkq∑j=0NajcosjkqΔx+∑j=0NbjsinjkqΔx=0.For the sake of brevity, we transform equation (12) to the matrix form,(15)e1TzTzTe2TWRAWRBWIAWIBAB=e1z,where the column vector is as follows:(16)A=a0a1...aN,B=b0b1...bN,e1=10...0N+1×1,e2=11...1N+1×1,z=00...0N+1×1.The elements of the remaining four matrices are formed by(17)WRAqp=−ωkqsinpkqΔx,(18)WIAqp=ωkqcospkqΔx,WIBqp=sinpkqΔx,q=1,…,N,p=0,…,N.For example, if we choose the wave numberk to calculate MBC1 (k). according to equation (15), the linear algebraic equations for the MBC1 (k) are as follows:(19)100000110−ωksinkΔx1coskΔxωkωkcoskΔx0sinkΔxa0a1b0b1=1000,where ωk=βk2.Then, the MBC1 (k) is found to be(20)du0dt+du1dt+k2EIρAcotkΔx2u0−k2EIρAcotkΔx2u1=0.The boundary condition MBC1 (k) can perfectly absorb the outgoing monochromatic wave ux,t=eikx+ωkt. Meanwhile, as N increases, more waves can be absorbed by boundary condition. Considering the abovementioned discussion of dispersion relation, we choose some special wave numbers to construct the boundary conditions. The coefficients of these boundary conditions are shown in Table 1. In addition, there is some symmetry in the boundary condition coefficients, which is shown in the table.Besides, whenN is odd, the stability of boundary conditions is good, while in other cases, the boundary conditions are generally unstable. This phenomenon is associated with the dispersion relation of the Euler–Bernoulli beam, and its reason needs to be further studies.


## 2.1. Dispersion Relation of the Euler–Bernoulli Beam Equation and Its Difference Scheme

The following shows the Euler–Bernoulli beam vibration equation without loads:(1)EIuxxxx+ρAutt=0.Here,ρ, A, E, and I represent the material density, the area of the beam cross section, Young’s modulus, and the bending moment of inertia, respectively. For a uniform beam, all these parameters are constant, and the equation can be rewritten as follows:(2)utt+EI/ρAuxxxx=0.One common spatially discrete model of the beam equation in application is the central finite-difference scheme,(3)d2ujdt2+EIρAuj−2−4uj−1+6uj−4uj+1+uj+2Δx4=0,which can be rewritten as(4)d2ujdt2=−β2Δx4uj−2−4uj−1+6uj−4uj+1+uj+2,with β=EI/ρA.We consider the wave form solution(5)u=eiknΔx+ωt.Substituting it into the semidiscrete scheme (3), we obtain a relation between the wave number k and the frequency ω, which gives the dispersion relation.(6)ω=±4βΔx2sin2kΔx2.The dispersion relation of the continuous Euler–Bernoulli beam (2) reads(7)ω=±βk2.According to the abovementioned dispersion relations, we can get the group velocitycg and the phase velocity cp, which represent the wave propagation law of the Euler–Bernoulli beam. The wave velocities are shown as follows:(8)cp=βk,cg=2βk.The wave propagating speed is obviously dependent on the wave number, which indicates the Euler–Bernoulli beam is a dispersive medium. To further discuss the difference of two dispersion relations, we plot the wave velocity versus wave number of each model in Figure1, with β = 1, Δx = 0.05. Also, the wave velocities of the discrete model are shown as follows:(9)cp=4βkΔx2sin2kΔx2,cg=2βΔxsinkΔx.Figure 1
Wave velocity in the Euler–Bernoulli beam (β = 1, Δx = 0.05): (a) phase velocity; (b) group velocity.
(a)(b)We notice that the deviation between the dispersion relation between the discrete model and the continuous model increases with the wave number. Hence, the long wave is a primary concern in designing matching boundary conditions. The boundary condition mentioned by Tang et al. [4] has a good absorption effect for short waves. For long waves such as k = 5 or longer with β = 1, Δx = 0.05, the absorption effect of their boundary condition will become worse. In addition, this boundary condition cannot be adjusted according to the wave number; that is, the absorption of this boundary condition for long waves cannot be improved. Therefore, it is necessary to propose a class of boundary conditions that can absorb long waves and can be adjusted for different wave numbers.


## 2.2. Matching Boundary Conditions

We construct boundary conditions by matching the dispersion relation of Euler–Bernoulli beams. To make the boundary conditions more convenient, we will consider the discrete scheme. Since the difference equation (3) uses a fourth-order finite-difference scheme in the spatial domain, it is necessary to design two-layer boundary conditions, as shown in Figure 2.Figure 2
Boundary condition design diagram: (a) schematic diagram of the surface layer boundary condition; (b) schematic diagram of the internal layer boundary condition at the left end.
(a)(b)It is naturally to take the same form for both the two artificial boundary layers. Besides, the same form makes the boundary conditions more compact. Equation (10) is the boundary condition of the surface layer in the left side, and equation (11) is the internal layer’s boundary condition. Therefore, we just need to determine the coefficients in the outer boundary layer conditions.(10)∑j=0Najdujdt+∑j=0Nbjuj=0,(11)∑j=1N+1ajdujdt+∑j=1N+1bjuj=0.Among them, the linear combination coefficientsajs and bjs are determined by matching the dispersion relation (7). Specifically, we use the following linear algebraic equations to calculate the linear combination coefficients. After matching the special wave numbers, the proposed boundary conditions (10) and (11) are named as MBCN (k1, ..., kN).(12)a0=1,Δ0=0,Δkq=0,q=1,2,...,N,where Δk is called the frequency domain residual function, which is defined by(13)Δk=iωk∑j=0NajeijkΔx+∑j=0NbjeijkΔx,where ω=βk2; β=EI/ρA.The first row is to give a fixed value for the homogeneous form of the boundary condition. The second row is the minimum “consistent” requirement. The third row means matching the dispersion relation at the selected wave number(s). Splitting it into the real and imaginary parts, we obtain the following two real-valued equations:(14)−ωkq∑j=0NajsinjkqΔx+∑j=0NbjcosjkqΔx=0,ωkq∑j=0NajcosjkqΔx+∑j=0NbjsinjkqΔx=0.For the sake of brevity, we transform equation (12) to the matrix form,(15)e1TzTzTe2TWRAWRBWIAWIBAB=e1z,where the column vector is as follows:(16)A=a0a1...aN,B=b0b1...bN,e1=10...0N+1×1,e2=11...1N+1×1,z=00...0N+1×1.The elements of the remaining four matrices are formed by(17)WRAqp=−ωkqsinpkqΔx,(18)WIAqp=ωkqcospkqΔx,WIBqp=sinpkqΔx,q=1,…,N,p=0,…,N.For example, if we choose the wave numberk to calculate MBC1 (k). according to equation (15), the linear algebraic equations for the MBC1 (k) are as follows:(19)100000110−ωksinkΔx1coskΔxωkωkcoskΔx0sinkΔxa0a1b0b1=1000,where ωk=βk2.Then, the MBC1 (k) is found to be(20)du0dt+du1dt+k2EIρAcotkΔx2u0−k2EIρAcotkΔx2u1=0.The boundary condition MBC1 (k) can perfectly absorb the outgoing monochromatic wave ux,t=eikx+ωkt. Meanwhile, as N increases, more waves can be absorbed by boundary condition. Considering the abovementioned discussion of dispersion relation, we choose some special wave numbers to construct the boundary conditions. The coefficients of these boundary conditions are shown in Table 1. In addition, there is some symmetry in the boundary condition coefficients, which is shown in the table.Besides, whenN is odd, the stability of boundary conditions is good, while in other cases, the boundary conditions are generally unstable. This phenomenon is associated with the dispersion relation of the Euler–Bernoulli beam, and its reason needs to be further studies.


## 3. Effectiveness of the MBC Reflection Coefficient Analysis

Reflection coefficient analysis is a very common method to evaluate the absorption effectiveness of artificial boundary conditions. In this section, we perform reflection coefficient analysis in the following way.Firstly, we consider a left-going harmonic wave and its reflected waves in the following form. Due to the order of the Euler–Bernoulli beam equation, there are two reflected waves generated at the boundary, including a propagating wave and an evanescent wave.(21)ujt=eiωt+kxj+Reiωt−kxj+R2eiωt−kxj,where R measures the amplitude of the reflected propagating wave. Also, R2 measures the amplitude of the reflected evanescent wave.Plugging (21) into the two-layers’ MBCs (8, 9), we obtain(22)Δk+RΔ−k+R2Δ2−k=0,eikΔxΔk+Re−ikΔxΔ−k+R2e−kΔxΔ2−k=0,where Δ2k=iωk∑j=0NajejkΔx+∑j=0NbjejkΔx.After some calculations,R and R2 are found to be(23)R=−ΔkΔ−ke1+ikΔx−1e1−ikΔx−1,R2=−ΔkΔ2−k−2ekΔxisinkΔxe1−ikΔx−1.For the general case, since what we usually care about is propagating waves, we defineR as the reflection coefficient of the boundary conditions. Noticing(24)e1+ikΔx−1e1−ikΔx−1=1,we may find the modulus of the reflection coefficient R.(25)R=ΔkΔ−k.In Figures3 and 4, we depict the variation of the reflection coefficients modulus of the MBCs listed in Table 1. According to the previous discussion, the incident harmonic wave with wave number kp can be perfectly absorbed by the proposed boundary conditions. Also, the downward spikes in Figures 3(b) and 4(b) correspond to this property. Additionally, as N increases, the reflection coefficient reduces for the wide range of wave numbers. The modulus of the reflection coefficient is 1. This phenomenon is a natural limitation because of the zero group velocity. However, since the group velocity is zero at k = 0, the reflection at the boundary does not have great influences on the numerical simulation.Figure 3
Reflection coefficient modulus for MBCs with a single wave number matched: (a) linear axis for reflection coefficients modulus; (b) logarithmic axis for reflection coefficients modulus.
(a)(b)Figure 4
Reflection coefficient modulus for MBCs with multiple wave numbers matched: (a) linear axis for reflection coefficients modulus; (b) logarithmic axis for reflection coefficients modulus.
(a)(b)Table 1
Coefficients in matching boundary conditions MBCN (with Δx = 0.05).
a0, aNa1, aN-1a2, aN-2b0, −bNb1, −bN-1b2, −bN-2MBC1 (10)1391.631MBC3 (5, 10, 15)1−1.605921096.36−3353.83MBC5 (2, 5, 10, 15, 18)1−4.614473.713471569.14−8330.8617151.09MBC2 (5, 10)1−1.5835117.90110MBC4 (2, 5, 10, 15)1−2.650173.09359193.181−388.8410The reflection coefficient study clearly reifies the effectiveness of the proposed boundary conditions can effectively suppress wave reflection. We remark that MBC3 and MBC5 are competent for general applications.


## 4. Effectiveness of the MBC Wave Packet Tests

To further verify the absorption effectiveness of the proposed boundary conditions, we perform two kinds of numerical tests, with either single or multiple main wave numbers. Also, for convenience, we takeEI/ρA = 1, Δx = 0.05 in all the following examples. Moreover, we use the velocity-Verlet method for time integration and let the time step Δt be 1 × 10−4.


### 4.1. Wave Packet Tests with Single Main Wave Numbers

In this section, we take the following form of wave packet as the initial value. The main wave number can be adjusted by the value ofk. Also, in the numerical tests of this section, k will take the values of 10 and 12.(26)fx=coskxcosπx5+1/2,x≤5,0,x>5.To make comparison, the reference solution is obtained by simulating a Euler–Bernoulli beam long enough such that the wave cannot reach the boundary in computational time. Furthermore, we adopt a new boundary condition MBC0 (k) to better reflect the effectiveness of the proposed boundary conditions. The left form of these boundary conditions is shown below, and the index of the first node on the left is 0. Essentially, MBC0 (k) is the discrete scheme of the single-wave equation.(27)du0dt+EIρAkΔxu1−u0=0,du1dt+EIρAkΔxu2−u1=0.The numerical solutions ofk = 10 with MBC1 (5), MBC0 (10), and MBC1 (1) are displayed in Figure 5. Since the numerical solutions have symmetry, we only show the right half of the wave profile evolves in several stages.Figure 5
Time evolution of the wave packet with the major wave number (k) = 10: (a) t = 0, (b) t = 0.5, (c) t = 1.0, and (d) t = 1.5.
(a)(b)(c)(d)Fort < 0.75, the initial wave packet is separated into two wave packets, propagating in two opposite directions, respectively. Due to dispersion, the length of the wave packet will increase with time. For t > 0.75, the wave packet frontiers arrive at the boundaries and reduce gradually. As the beam dispersion relation dictates, there is a very slow wave with a small wave number. Hence, the attenuation process is very slow.With all the three boundary conditions for the segment computations, we observe that the numerical solutions agree very well att = 0.5 from Figure 5(b). Later, as the wave packet gradually passes through the boundary, the reflection occurs at the boundary, as shown in Figures 5(c) and 5(d). To further demonstrate the effectiveness of the proposed boundary conditions, we depict the displacement error at t = 1 and t = 1.5 in Figure 6. Figure 6(a) at t = 1 shows observable deviations, more readable for the displacement near the boundaries. After this stage, further reflections occur and propagate back to the other end of the segment. In Figure 6(b), at t = 1.5, the deviations grow to an observable level. MBC1 (10) induces reflection error smaller than MBC1 (5) and MBC0 (10).Figure 6
Error of displacement with the major wave number (k) = 10: (a) t = 1.0; (b) t = 1.5.
(a)(b)The numerical solutions ofk = 12 with MBC1 (10), MBC3 (5, 10, 15), and MBC5 (2, 5, 10, 15, 18) are displayed in Figure 7. As N increases, we involve more notes in an MBC. Similarly, we depict the displacement error at t = 1 and t = 1.5 in Figure 8. Although the wave number of k = 12 is not matched, the error is gradually reduced with the increase of N.Figure 7
Time evolution of wave packet with the major wave numberk = 12: (a) t = 0, (b) t = 0.5, (c) t = 1.0, and (d) t = 1.5.
(a)(b)(c)(d)Figure 8
Error of displacement with the major wave numberk = 12: (a) t = 1.0; (b) t = 1.5.
(a)(b)Besides the displacement, the energy of the solutions is also shown in Figure9. To initially investigate the stability of the boundary conditions, numerical simulations are performed for a long term. The total energy is calculated without the complex boundary and normalized by the initial energy, specifically by(28)Energyt=12Energyinitial∑i=6Nd−5duidt2+∑i=7Nd−6ui+1−2ui+ui−1Δx22,where Nd is total number of nodes.Figure 9
Energy of wave packet with the major wave number: (a)k = 10; (b) k = 12.
(a)(b)Before the wave packet frontiers arrive at the boundaries, the total energy conserves. For a matched wave, even MBC1 captures the energy of magnitude order 10−4 at t = 1.5. Also, for an unmatched wave, MBC3 captures the energy of magnitude order 10−5 at t = 1. In addition, there is no divergence in these results.


### 4.2. Wave Packet Tests with Multiple Main Wave Numbers

We simulate two cases with multiple wave numbers to further demonstrate the effectiveness of the proposed boundary condition in a wide range of wave numbers. To better display the effectiveness of the boundary conditions, we magnify the computational domain since the wave propagates in the beam with strong dispersion. Furthermore, suitable initial velocity is used in each case so that the initial wave packets can propagate substantially only to the right.The wave packet of first case contains three main wave numbers,k = 5, 10, and 15. This case is treated with different boundary conditions and is used to verify that the proposed boundary condition can simultaneously treat the waves which contain more than one matched wave numbers.Figure10 shows the waveforms of the examples at different stages. For t < 1.1, the initial wave packet is split into three wave packets due to the strong dispersion, propagating to the right. Then, the wave packets reach the boundary and gradually disappear. Due to the dispersion relation and the wave number, the propagating speed of the fastest wave packet is almost three times as large as the slowest. Thus, when the slowest wave packet basically enters the right boundary, the reflected wave of the fastest wave packet is nearly to arrive at the left boundary.Figure 10
Time evolution of the waveform with three major wave numbersk = 5, 10, and 15: (a) (t) = 0, (b) t = 1.0, and (c) t = 6.0.
(a)(b)(c)With all the three boundary conditions for the segment computations, we observe that the numerical solutions agree very well att = 1 from Figure 10(b). Figure 10(c) shows observable deviations at t = 6, where MBC1 (10) has larger reflection of the wave packets of k = 5 and 15 and MBC1 (15) has larger reflection of the wave packets of k = 5 and 10. In contrast, MBC3 (5, 10, 15) can absorb the three wave packets simultaneously. This phenomenon is basically consistent with the reflection coefficient of the boundary conditions.Like Section4.1, we performed a long-term numerical simulation. Figure 11 shows the energy variation of the numerical testing in t < 20, with no significant divergence.Figure 11
Energy of the wave packet with three major wave numbersk = 5, 10, and 15.The second case contains two main wave numbers,k = 4 and 12 and will be treated with different boundary conditions. This case is used to verify that as the number of matching wave numbers N increases, the absorption effectiveness over a larger range of wave number will also increase.The results of this example at different stages are shown in Figure12. For t < 1.4, the initial wave packet is split into two wave packets, propagating to the right. Then, the wave packet reaches the boundary and gradually disappears. At t = 7, there are some observable deviations, which decrease as N increases, shown in Figure 12(c). This phenomenon is also consistent with the reflection coefficient of the boundary conditions. Figure 13 shows that the stability of the proposed conditions is good.Figure 12
Time evolution of the waveform with two major wave numbersk = 4 and 12: (a) t = 0, (b) t = 1.0, and (c) t = 6.0.
(a)(b)(c)Figure 13
Energy of the wave packet with two major wave numbersk = 4 and 12.In summary, we demonstrate the effectiveness of the proposed boundary conditions with the numerical tests, consistent with the reflection coefficient analysis.


## 4.1. Wave Packet Tests with Single Main Wave Numbers

In this section, we take the following form of wave packet as the initial value. The main wave number can be adjusted by the value ofk. Also, in the numerical tests of this section, k will take the values of 10 and 12.(26)fx=coskxcosπx5+1/2,x≤5,0,x>5.To make comparison, the reference solution is obtained by simulating a Euler–Bernoulli beam long enough such that the wave cannot reach the boundary in computational time. Furthermore, we adopt a new boundary condition MBC0 (k) to better reflect the effectiveness of the proposed boundary conditions. The left form of these boundary conditions is shown below, and the index of the first node on the left is 0. Essentially, MBC0 (k) is the discrete scheme of the single-wave equation.(27)du0dt+EIρAkΔxu1−u0=0,du1dt+EIρAkΔxu2−u1=0.The numerical solutions ofk = 10 with MBC1 (5), MBC0 (10), and MBC1 (1) are displayed in Figure 5. Since the numerical solutions have symmetry, we only show the right half of the wave profile evolves in several stages.Figure 5
Time evolution of the wave packet with the major wave number (k) = 10: (a) t = 0, (b) t = 0.5, (c) t = 1.0, and (d) t = 1.5.
(a)(b)(c)(d)Fort < 0.75, the initial wave packet is separated into two wave packets, propagating in two opposite directions, respectively. Due to dispersion, the length of the wave packet will increase with time. For t > 0.75, the wave packet frontiers arrive at the boundaries and reduce gradually. As the beam dispersion relation dictates, there is a very slow wave with a small wave number. Hence, the attenuation process is very slow.With all the three boundary conditions for the segment computations, we observe that the numerical solutions agree very well att = 0.5 from Figure 5(b). Later, as the wave packet gradually passes through the boundary, the reflection occurs at the boundary, as shown in Figures 5(c) and 5(d). To further demonstrate the effectiveness of the proposed boundary conditions, we depict the displacement error at t = 1 and t = 1.5 in Figure 6. Figure 6(a) at t = 1 shows observable deviations, more readable for the displacement near the boundaries. After this stage, further reflections occur and propagate back to the other end of the segment. In Figure 6(b), at t = 1.5, the deviations grow to an observable level. MBC1 (10) induces reflection error smaller than MBC1 (5) and MBC0 (10).Figure 6
Error of displacement with the major wave number (k) = 10: (a) t = 1.0; (b) t = 1.5.
(a)(b)The numerical solutions ofk = 12 with MBC1 (10), MBC3 (5, 10, 15), and MBC5 (2, 5, 10, 15, 18) are displayed in Figure 7. As N increases, we involve more notes in an MBC. Similarly, we depict the displacement error at t = 1 and t = 1.5 in Figure 8. Although the wave number of k = 12 is not matched, the error is gradually reduced with the increase of N.Figure 7
Time evolution of wave packet with the major wave numberk = 12: (a) t = 0, (b) t = 0.5, (c) t = 1.0, and (d) t = 1.5.
(a)(b)(c)(d)Figure 8
Error of displacement with the major wave numberk = 12: (a) t = 1.0; (b) t = 1.5.
(a)(b)Besides the displacement, the energy of the solutions is also shown in Figure9. To initially investigate the stability of the boundary conditions, numerical simulations are performed for a long term. The total energy is calculated without the complex boundary and normalized by the initial energy, specifically by(28)Energyt=12Energyinitial∑i=6Nd−5duidt2+∑i=7Nd−6ui+1−2ui+ui−1Δx22,where Nd is total number of nodes.Figure 9
Energy of wave packet with the major wave number: (a)k = 10; (b) k = 12.
(a)(b)Before the wave packet frontiers arrive at the boundaries, the total energy conserves. For a matched wave, even MBC1 captures the energy of magnitude order 10−4 at t = 1.5. Also, for an unmatched wave, MBC3 captures the energy of magnitude order 10−5 at t = 1. In addition, there is no divergence in these results.


## 4.2. Wave Packet Tests with Multiple Main Wave Numbers

We simulate two cases with multiple wave numbers to further demonstrate the effectiveness of the proposed boundary condition in a wide range of wave numbers. To better display the effectiveness of the boundary conditions, we magnify the computational domain since the wave propagates in the beam with strong dispersion. Furthermore, suitable initial velocity is used in each case so that the initial wave packets can propagate substantially only to the right.The wave packet of first case contains three main wave numbers,k = 5, 10, and 15. This case is treated with different boundary conditions and is used to verify that the proposed boundary condition can simultaneously treat the waves which contain more than one matched wave numbers.Figure10 shows the waveforms of the examples at different stages. For t < 1.1, the initial wave packet is split into three wave packets due to the strong dispersion, propagating to the right. Then, the wave packets reach the boundary and gradually disappear. Due to the dispersion relation and the wave number, the propagating speed of the fastest wave packet is almost three times as large as the slowest. Thus, when the slowest wave packet basically enters the right boundary, the reflected wave of the fastest wave packet is nearly to arrive at the left boundary.Figure 10
Time evolution of the waveform with three major wave numbersk = 5, 10, and 15: (a) (t) = 0, (b) t = 1.0, and (c) t = 6.0.
(a)(b)(c)With all the three boundary conditions for the segment computations, we observe that the numerical solutions agree very well att = 1 from Figure 10(b). Figure 10(c) shows observable deviations at t = 6, where MBC1 (10) has larger reflection of the wave packets of k = 5 and 15 and MBC1 (15) has larger reflection of the wave packets of k = 5 and 10. In contrast, MBC3 (5, 10, 15) can absorb the three wave packets simultaneously. This phenomenon is basically consistent with the reflection coefficient of the boundary conditions.Like Section4.1, we performed a long-term numerical simulation. Figure 11 shows the energy variation of the numerical testing in t < 20, with no significant divergence.Figure 11
Energy of the wave packet with three major wave numbersk = 5, 10, and 15.The second case contains two main wave numbers,k = 4 and 12 and will be treated with different boundary conditions. This case is used to verify that as the number of matching wave numbers N increases, the absorption effectiveness over a larger range of wave number will also increase.The results of this example at different stages are shown in Figure12. For t < 1.4, the initial wave packet is split into two wave packets, propagating to the right. Then, the wave packet reaches the boundary and gradually disappears. At t = 7, there are some observable deviations, which decrease as N increases, shown in Figure 12(c). This phenomenon is also consistent with the reflection coefficient of the boundary conditions. Figure 13 shows that the stability of the proposed conditions is good.Figure 12
Time evolution of the waveform with two major wave numbersk = 4 and 12: (a) t = 0, (b) t = 1.0, and (c) t = 6.0.
(a)(b)(c)Figure 13
Energy of the wave packet with two major wave numbersk = 4 and 12.In summary, we demonstrate the effectiveness of the proposed boundary conditions with the numerical tests, consistent with the reflection coefficient analysis.


## 5. Conclusions and Discussion

In this paper, we propose a class of matching boundary conditions for Euler–Bernoulli beams. Based on the dispersion relation of the Euler–Bernoulli beam, a class of effective and explicit artificial boundary conditions is designed, which can be conveniently used in actual research. The presented boundary conditions take the form of a linear combination of displacement and velocity at nodes near the artificial boundary, and its coefficients is determined by matching the dispersion relation. Also, the calculation process of the boundary condition coefficient is very concise; thus, only a simple adjustment is needed to obtain the boundary conditions for different problems. Furthermore, users of boundary conditions can choose some suitable wave numbers for the specific simulations, and perfect absorption may be reached for the outgoing waves at selected wave numbers. Reflection suppression is enhanced progressively by matching more wave numbers. The reflection coefficient analysis and numerical tests demonstrate the effectiveness of the proposed MBCs. Besides, the stability of artificial boundary conditions is often difficult to fully prove, and in this work, the stability of the odd-numbered boundary conditions is preliminarily verified by long-time numerical simulation, while the even-numbered boundary conditions need to be further studied. In sum, with a simpler form, greater computational efficiency, as well as better absorbing performance, MBCs could be used in real-time simulation of infinite Euler–Bernoulli beams. With some extensions, the construction method of the matching boundary condition in this article can handle some more complex structures (e.g., Timoshenko beam, Kirchhoff board, and Mindlin board).



---

*Source: 6685852-2021-03-15.xml*