An Efficient Local Artificial Boundary Condition for Infinite Long Finite Element Euler–Bernoulli Beam

. To solve the wave propagation problems of the Euler–Bernoulli beam in an unbounded domain efectively and efciently, a new local artifcial boundary condition technology is proposed. It replaces the residual right-hand side of the truncated discrete equation with an equivalent linear algebraic system. First, the equivalent Schrodinger equation is discussed. Its artifcial boundary condition is obtained by frst rationalizing the Dirichlet-to-Neumann condition in the frequency domain with a Pade approximation and then inverse transforming each Pade term back into the time domain by introducing auxiliary degrees of freedom. Frequency shifting is employed such that it performs better near a prescribed frequency. Ten, the artifcial boundary condition of the fnite element Euler–Bernoulli beam is obtained by simple algebraic manipulations on that of the corresponding Schrodinger equation. Tis method only makes local changes to the original truncated discrete dynamic system and thus is very efcient and easy to use. Te accuracy of the proposed method can be improved by using more Pade terms and a proper shift frequency. Te numerical example shows, with only a few additional degrees of freedom, the proposed artifcial boundary condition efectively eliminates the spurious refection. Te idea of the proposed method can also be used in other dispersive wave systems.

Besides widely discussed second-order systems, higher order systems were also discussed.In this study, we focus on the fexible wave propagation on an infnitely long Euler-Bernoulli beam, which is meaningful in multiscale computation and civil engineering [22][23][24].Tang proposed an accurate ABC for the semidiscrete Euler beam system, which expresses the defection at the truncated boundary as the time convolution of its neighboring nodes' defection [25], i.e., Tis method is time consuming due to the convolution operation.Feng proposed a matching boundary condition (MBC) by looking for an approximate linear relation among the m + 1 nodes near the truncated boundary [19]: which is local in space and time and thus very efcient but less accurate.Based on Tang's work [25], Zheng expressed the convolution kernel in equation (1) as the summation of exponentials, i.e., where a j , b j are constant coefcients and thus simplify the computation by using the recursive relation of the discrete convolution of exponential function [26].Te method is accurate and efcient, but its implementation is tedious.
On the other hand, to the knowledge of the authors, the previous ABCs for beams consider the fnite diference method (FDM) only.Te ABC for the widely used fnite element method (FEM) beam model is seldom discussed.Te biggest diference between FEM and FDM beams is that the former involves rotational degrees of freedom (DOF) (e.g., [27,28]).In this work, we propose a new ABC for FEM Euler-Bernoulli beam, which is local, efcient, accurate, and easy to implement.It is not a variant of existing ABCs in the literature such as [19,25,26].Te idea is to attach a carefully designed mass-damper super element to the truncation end, and ensuring it recovers the dynamics response of the removed half infnite segment.A similar idea can be found in literature [8,9], in which the continued fraction technique is employed to recover the dynamic response of a half infnitely long waveguide.Te diference is that, in [8,9], the auxiliary DOFs are connected in series, but in this study, they are connected in parallel.
Te rest of the paper is organized as follows: For simplicity, we frst propose an ABC for the equivalent Schrodinger equation and discussed its performance in Section 2. Ten, the ABC for the fnite element beam is derived by the transform of the Schrodinger equation in Section 3. In Section 4, a numerical example is given to test its performance.Finally, there is a conclusion in Section 5.

ABC for Schrodinger Equation
where φ 0 (x), φ 1 (x) are the initial defection and initial velocity, respectively, and q(x, t) is the distributed transversal force per length.In this study, we always assume φ 0 (x), φ 1 (x) and q(x, t) are compactly supported in the spatial domain.
It is well-known that the Euler-Bernoulli equation is closely related to the Schrodinger equation (25): In fact, the Schrodinger equation ( 5) is related to the Euler-Bernoulli system equation ( 4) in the real domain.
Here, the source terms and the initial values of the two equations should satisfy If φ 1 (x) and q(x, t) are compactly supported, then by employing proper integral constants in equation (6), p(x, t), ϕ(x) are also compactly supported.
Tus, theoretically, the wave propagation can be simulated by solving equation (5) rather than equation (4).In this study, we frst derive the ABC of equation ( 5) and then that of equation ( 4) by using some simple algebraic transformation.
To numerically solve equation ( 5), the unbounded domain is truncated.Te bounded computation zone should cover the support domains of φ 0 (x), φ 1 (x), and q(x, t).Without loss of generality, only the treatment of the upper boundary is considered.Although it is possible to build the proposed ABC for the continuous equation ( 5), here we consider the semidiscrete form such that we can see its numerical implementation more clearly.By using a uniform mesh size h, the last few rows/columns of the central difference approximation of truncated equation ( 5) can be expressed as 2 Shock and Vibration where N is the last preserved node, which is assumed far away from the supported domains of φ 0 (x), φ 1 (x) and q(x, t).It can be seen equation ( 7) is not closed because the right-hand side contains an unknown component (u N+1 − u N )/h 2 + iz t u N /2.Tis entry refects the efect of the removed part in the form of an additional source term.To close equation (7), we still need a proper DtN condition that represents the spatial diference z x u N as a function of u N .

Pade Approximation and New ABC.
In the frequency domain, it is not difcult to fnd such a DtN condition.By considering the monochromatic right-going incident wave u(x, t) � u * exp(ikx − iωt), where u * is the complex amplitude, one obtains in which k is the wavenumber and ω is the angular frequency.Te analytical dispersion relation of the Euler-Bernoulli system ω � k 2 has been used to get equation (8).Te literature recommend that the dispersion relation should be modifed for semidiscretized systems [19], i.e., ω � 4 sin 2 (kh/2)/h 2 .But that will make the following progress more difcult.
To implement equation (8) in the time domain, we rationalize it by employing the following Pade approximation [16]: Here, M is the number of the Pade approximation terms, and the right-hand side converges to the left when M approaches infnity.By letting s � ω/ω 0 − 1 in equation ( 9), we can express equation (8) as in which Figure 1 demonstrates how the approximation performs with diferent M and ω 0 .It can be seen that with a small M, the approximation works good near ω � ω 0 .Tus, one can choose a proper constant ω 0 according to the frequency spectrum of a specifed problem.Now, we treat each term of the right-hand side of equation (10) separately.Te f 0 term does not contain ω and thus can be transformed back into the time domain easily.As for the f j (j > 0) term, we introduce an auxiliary variable v j such that

Shock and Vibration
It is easy to verify that equation ( 12) equals to equation ( 11) by eliminating v j in equation (12).Obviously, the time domain form of equation ( 12) is By replacing equation (11) with equation ( 13) in equation (10) and substituting the result into equation (7), one can have the fnal dynamic equation with ABC implemented.
Te auxiliary DOF v j is introduced as additional unknowns of the equation.From the second equation of equation ( 12), we can see the additional DOFs are oscillating rather than divergent for a given monochromatic incident wave.
As an example, the last few rows/columns of the fnal algebraic equation for M � 2 are given as follows:

Shock and Vibration
Note that in the last few rows of equation ( 14) the source term is zero, because we assumed node N is far away from the support domain of q(x, t).But the right-hand side can be nonzero in other rows if q(x, t) exists.
It can be seen the ABC implementation equals to parallelly attaching several algebraic systems to the last node of the truncated system.Tis kind of matrix manipulation is quite common in numerical methods and thus can be easily implemented.Finally, equation ( 14) can be solved by using the central fnite diference method in the time domain.

Refection Coefcient Analysis.
For an ideal ABC, no refection should occur at the truncation boundary.Tus, the ratio between the artifcial refection wave amplitude and incident wave amplitude can be used to assess the performance of an ABC.After going through the proposed method, it can be seen there are two error sources which may cause artifcial refection for the proposed ABC.One is the discretization error (equation ( 8)), which decreases with decreasing mesh size and is not discussed here.Te other is the Pade approximation (equation ( 10)), whose error estimation is given by [16]: Obviously, it approaches zero when M approaches infnity for any ω ≥ 0 and vanishes automatically at ω � ω 0 .
Let the superposition of the incident and refection waves near the ABC be where the absolute value of R is known as the refection coefcient.By substituting equation ( 16) into equation (10), one obtains that ik exp(ikx) Substituting the dispersion relation and equation ( 15) into equation (17), it follows that It can be seen the refection coefcient approaches zero when M increases for ω > 0. Besides, it converges fast near ω � ω 0 because c vanishes there.
Figure 2 depicts the refection coefcient.It shows that for frequencies near ω � ω 0 , |R| is rather small even though a small M is used.Te coefcient by matching boundary condition (MBC) [19] is also shown for comparison.It can be seen that by introducing the same number of additional DOFs, the proposed method is more favorable when a single wide frequency band is of interest, while the MBC performs better for multiple narrow frequency bands.However, the conclusion assumes that h is adequately small.If the step size h gets larger, the MBC may perform better because it is based on the semidiscrete beam model.

ABC for FEM Euler Beams
3.1.FEM Discretization.As discussed before, the wave propagation on an infnitely long Euler beam can be simulated by solving the Schrodinger equation (5) with FDM.A complex value sparse matrix solver is necessary for this purpose.However, in mechanics, we prefer to solve the problem by FEM with a real-value solver.Based on previous discussions, here we consider the ABC of FEM beam which contains both transversal and rotational DOFs.When employing the classic 2-node Hermite beam element and the consistent mass matrix, the element equation of equilibrium for problem equation ( 4) is as follows [27,28]: Here, θ, F, m are the slope angle, the transversal force, and the concentrated couple, respectively.Te sign conventions of the variables are shown in Figure 3.
After assembling equation (19) for the truncated Euler beam, the last few rows/columns of the total equation are   Shock and Vibration Physically, F N , m N on the right-hand side are the timedependent force and the moment applied on the truncated beam by the removed part.

ABC Transformation.
Similar to the previous discussion, the point is to fnd the relation between the additional source terms F N , m N and DOFs on boundary u N , θ N .Recalling that the real part of the solution to the Schrodinger equation ( 5) is related to the Euler-Bernoulli system equation ( 4), this can be done by separating the real and imaginary parts of equation (13).Noticing in equation ( 12) k j is pure imaginary, it follows that Here, the subscript re/im denotes the real or imaginary part.From equation ( 5) and the mechanics of material, we have By taking time derivative of equation ( 21) and substituting equation (23), it follows that Equation ( 24) gives the relation between m N and θ N,re .In the same way, the time derivative of equation ( 22) can be written as

Shock and Vibration
Equation (25) gives the relation between F N and u N,re .It can be seen there are totally 4M + 2 equations in equations ( 24) and ( 25) and (4M) + 4 variables, i.e., 2M auxiliary variables (v j,re , v j,im ), 2M intermediate variables (f j z t u N,re , f j z t u N,im ), and 4 variables which already exist in Eq (20) (u N,re , θ N,re , F N , m N ).Note that u N,im at the right-hand side of (25) has no contribution to the equation and thus does not count.Tus, with the help of ( 24) and ( 25), equation ( 20) is closed by introducing 2M additional variables (v j,re , v j,im , j � 1, 2..., M).As an example, the matrix form of equations ( 24) and ( 25) for M � 2 is given as follows: It can be easily assembled into the total equation of equilibrium Eq. ( 20) by considering the matrices as the damper/mass matrices of a supper element which attached to the last node of the truncated FEM beam.All entries of the matrices are real, and thus, a real-value solver can be used for computation.

Dispersion Relation.
Although the previous discussion about the ABC does not really involve the mass matrix, here for simulating the wave propagation on unbounded FEM Euler beams, we recommend consistent mass matrix rather than lumped mass matrix, because the former leads to a more accurate dispersion relation.

8
Shock and Vibration By assembling equations ( 19) of two successive beam elements which have no external loads, one obtains Te third and fourth equations in equation ( 27) are Te frst equation in equation (28) represents the shear force equilibrium, and the second represents the bending moment equilibrium.Unlike in FDM, each of the equations in equation ( 28) will give an independent dispersion relation, i.e., However, the two relations are compatible in an approximate manner.Te frst few terms of the Taylor series of equation ( 29) near kh ⟶ 0 are It can be seen the error between the relation in equation ( 30) and the analytical one (ω � k 2 ) is O((kh) 6 ), which is good enough.Te lumped mass matrix will lead to a less accurate relation.On the other hand, for the discrete Schrodinger equation (7) and the FDM beam [19], the diference between resulted dispersion relation ω � 4 sin 2 (kh/2)/h 2 and the analytical one is only O((kh) 4 ).As a result, the group velocity of the wave on the FEM beam is closer to that of the analytical beam, and it is larger than that of the FDM model/Schrodinger model.

Numerical Example
To validate the proposed method, we consider the example from Feng, 2021 [19].
Example 1.In equation ( 4), let the source term be q(x, t) ≡ 0 and the initial values be Te computation time domain is from 0 to 1.5.
Two proposed methods, i.e., the equivalent Schrodinger equation way and the FEM way, are employed to save the problem, with M � 3, �� � ω 0 √ � 10.Te ratio between the time step and mesh size is fxed, i.e., τ � h/500 throughout the example.First, we use the same setup as in literature [19], i.e., mesh size h � 0.05 and truncation boundary |x| � 20.Te Newmark method (see e.g., [27,28]) is employed for time integration.Te MBC5 result from [19] and the fast convolution result from [26] are given for comparison.
Two reference solutions are also obtained by FDM.Reference solution A employs a much larger domain but the same step size (|x| < 500, h � 1/20) such that no wave reaches the boundary, while reference solution B further employs a smaller step size (|x| < 500, h � 1/320) and a higher internal precision.Obviously, reference B is more convincing than  Shock and Vibration reference A, because it not only eliminates the boundary efect as reference A does but also reduces the mesh size efect.
Te solutions at diferent time points are shown in Figure 4. From Figure 4(a), it can be seen the solutions fall into two groups.Te waves of the FEM model and reference B are similar, and they go faster than those of the other four FDM models.Te MBC5 solution, the fast convolution solution, and the Schrodinger solution are close to reference A. As discussed before, this is because the FEM model has a more accurate dispersion relation than FDM models and thus a more accurate wave velocity when using the same step size.Te same can be seen in Figure 4(b), except that the Schrodinger model falls even behind.
In Figure 4(c), when the wave should have almost left the computation domain, the wave amplitudes given by reference A, reference B, the fast convolution method, and the FEM model are extremely small.Tis shows the FEM model efectively suppresses the boundary refection, even though only a few additional DOFs are introduced (M � 3).In this example, the FEM model performs no worse than the fast convolution method, but it is more efcient because the latter uses a global ABC.However, for the same M, the Schrodinger model leads to a signifcant artifcial refection.Te MBC5 model performs better than the Schrodinger model but worse than the others.
It is interesting that the FEM model performs much better than the Schrodinger model for this example,

Shock and Vibration
considering the former ABC is derived from that of the latter.Te reason is that they are both based on the analytical dispersion relation (see equation ( 8)), and the FEM model is more honest to the relation (see equation ( 30)).Since the dispersion relation is accurate only when kh is infnitesimal, to improve the performance of the Schrodinger model, one should use a smaller h.
Diferent M and h are thus employed for the Schrodinger model, and the solutions at t � 1.5 are given in Figure 5.It can be seen that the result gradually converges to reference B. With a smaller h, the left-going refection wave goes faster, and its amplitude is smaller.But a greater M does not seem to make any diference for this example.
Te performances of the models are further quantifed by the L2 norm of their diferences to reference B at t � 1.5.Te errors are plotted against the mesh size h in Figure 6.Te Schrodinger model shows a second-order convergence, for both M � 3 and M � 30.Te errors of FEM models are much smaller, and again, the infuence of M is not signifcant.With decreasing mesh size, the errors frst drop dramatically and then no longer change.We suggest that it is because the error in each cycle has reached the limit of the double-precision computation for small mesh sizes.Tus, we further employ a high-precision computation (32 digits), as denoted by FEM HP in Figure 6.In this condition when M � 3, the curve is similar to that of the double-precision solution.However, when a greater M (�30) is employed, a fourth-order convergence can be observed, which is very encouraging.Tis is because the time step is adequately small in this example, and the cubic Hermite element is employed in FEM.If the time step were large, a second-order convergence would be noticed because the time integration scheme has a secondorder accuracy.In summary, Figure 6 shows the proposed method can reduce the error to the machine precision level with a small M, but for high-precision computation, a greater M leads to an even better result.
To show the long-term performance of the proposed method, we extend the simulation time to t � 10 and calculate the total energy evolution in the computation zone (|x| < 20).Obviously, the energy should gradually drop to zero as the wave propagates out of the computation zone.In FEM, the total energy is given by where K tot , M tot , u are the total stifness matrix, the total mass matrix, and the displacement vector of the preserved beam part (not including auxiliary DOFs), respectively.Te result for h � 0.05 is shown in Figure 7.
It can be seen at the early stage, the energy drops by the proposed method match the reference one quite well.Actually, before t � 3, the proposed method performs best.After that, the energy level has become very low (<1E-8), the energy of the proposed model is higher than that of others.It implies the artifcial refection wave is not eliminated perfectly after a long-time simulation.Tis may be because the wave components whose frequencies are far away from ω 0 cannot be efectively treated by our method (see Figure 2).Using a greater M does not help much.In contrast, the MBC5 method brings more refection energy when the main waves frst reach the boundary, but after the wave goes forth and back for several times in the interval, its fnal energy level is lower.Tat is because of its good capability to absorb waves with a frequency far away from ω 0 (Figure 2).Te fast convolution method leads to the lowest residual energy after a long simulation time, but it performs worse than the proposed method for the frst collision with the artifcial boundary.12 Shock and Vibration

Conclusion
In this study, we proposed a new ABC for an infnitely long FEM Euler-Bernoulli beam which has rotational DOFs.In the ABC, the infuence of the removed half infnite beam segment is represented by the residual right-hand side, i.e., the unbalanced transversal force and moment at the boundary.To close the equation, by Pade approximation, the right-hand side is rewritten as the linear combination of the DOFs at the truncation point, some introduced auxiliary DOFs, and their time derivatives.Tis treatment equals to parallelly attaching several dynamic systems to the truncated system and thus is easy to be implemented in the frame of FEM.Te proposed method is local in both space and time domains.Te modifcation to the original algebraic system occurs only near the truncation points and has a quite low rank comparing with the matrix size.Te time integration scheme has no diference with that used in common dynamic systems, and no convolution is needed.Te attached systems will automatically absorb the incoming waves.As a result, the method almost costs no extra computation time.
Te error of the proposed ABC is controlled by both mesh size and the number of terms used in Pade approximation.If the goal accuracy is not very high, a few Pade approximation terms and a proper mesh size are usually enough.A higher accuracy can also be achieved by employing a fner mesh and more approximation terms.Even for the high accuracy condition, it is still more efcient than the global ABCs.
An ABC for the Schrodinger equation is also proposed as a byproduct.It is less appealing because it requires a very fne mesh to efectively suppress the artifcial refection.However, it shows a second-order convergence as the mesh size decreases.

Figure 3 :
Figure 3: Sign conventions of the beam element variables.
2.1.FDM Discretization.Te normalized Euler-Bernoulli equation in an unbounded domain can be written as