Abstract
To solve the wave propagation problems of the Euler–Bernoulli beam in an unbounded domain effectively and efficiently, a new local artificial 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 artificial boundary condition is obtained by first 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. Then, the artificial boundary condition of the finite element Euler–Bernoulli beam is obtained by simple algebraic manipulations on that of the corresponding Schrodinger equation. This method only makes local changes to the original truncated discrete dynamic system and thus is very efficient and easy to use. The accuracy of the proposed method can be improved by using more Pade terms and a proper shift frequency. The numerical example shows, with only a few additional degrees of freedom, the proposed artificial boundary condition effectively eliminates the spurious reflection. The idea of the proposed method can also be used in other dispersive wave systems.
1. Introduction
To numerically solve partial differential equations in an infinitely large domain, it is common to truncate the domain and introduce the artificial boundary conditions (ABC) at the truncation boundaries. There are rather extensive studies on ABCs for PDEs such as the heat equation [1, 2], classical Schrodinger equation [2–5], wave equation [6–9], and Klein–Gordon equation [10–13]. Generally speaking, when dealing with Cauchy problems of partial differential equations, the exact ABCs are nonlocal in time, involving temporal convolutions in their formulations. The convolutions often require a significant amount of computational resources in application. Thus, modification schemes were proposed to accelerate the convolution computation, and many of them involve expressing the kernel as the sum of exponential functions [14–17]. On the other hand, the approximate local ABCs are more efficient. Typical local ABC techniques include Higdon-like condition [7, 12], perfectly matched layer (PML) [6, 11], viscoelastic boundary [8, 9, 18], and matching boundary condition [19–21]. In general, the local ABCs are less accurate, and/or they include empirical parameters.
Besides widely discussed second-order systems, higher order systems were also discussed. In this study, we focus on the flexible wave propagation on an infinitely long Euler–Bernoulli beam, which is meaningful in multiscale computation and civil engineering [22–24]. Tang proposed an accurate ABC for the semidiscrete Euler beam system, which expresses the deflection at the truncated boundary as the time convolution of its neighboring nodes’ deflection [25], i.e.,
This 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 efficient 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 are constant coefficients and thus simplify the computation by using the recursive relation of the discrete convolution of exponential function [26]. The method is accurate and efficient, but its implementation is tedious.
On the other hand, to the knowledge of the authors, the previous ABCs for beams consider the finite difference method (FDM) only. The ABC for the widely used finite element method (FEM) beam model is seldom discussed. The biggest difference 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, efficient, accurate, and easy to implement. It is not a variant of existing ABCs in the literature such as [19, 25, 26]. The 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 infinite 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 infinitely long waveguide. The difference is that, in [8, 9], the auxiliary DOFs are connected in series, but in this study, they are connected in parallel.
The rest of the paper is organized as follows: For simplicity, we first propose an ABC for the equivalent Schrodinger equation and discussed its performance in Section 2. Then, the ABC for the finite 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.
2. ABC for Schrodinger Equation
2.1. FDM Discretization
The normalized Euler–Bernoulli equation in an unbounded domain can be written aswhere are the initial deflection and initial velocity, respectively, and is the distributed transversal force per length. In this study, we always assume and 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 and are compactly supported, then by employing proper integral constants in equation (6), are also compactly supported.
Thus, theoretically, the wave propagation can be simulated by solving equation (5) rather than equation (4). In this study, we first 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. The bounded computation zone should cover the support domains of , and . 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 aswhere N is the last preserved node, which is assumed far away from the supported domains of and . It can be seen equation (7) is not closed because the right-hand side contains an unknown component . This entry reflects the effect 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 difference as a function of .
2.2. Pade Approximation and New ABC
In the frequency domain, it is not difficult to find such a DtN condition. By considering the monochromatic right-going incident wave , where is the complex amplitude, one obtainsin which k is the wavenumber and is the angular frequency. The analytical dispersion relation of the Euler–Bernoulli system has been used to get equation (8). The literature recommend that the dispersion relation should be modified for semidiscretized systems [19], i.e., . But that will make the following progress more difficult.
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 infinity. By letting in equation (9), we can express equation (8) asin which
Figure 1 demonstrates how the approximation performs with different M and . It can be seen that with a small M, the approximation works good near . Thus, one can choose a proper constant according to the frequency spectrum of a specified problem.
(a)
(b)
Now, we treat each term of the right-hand side of equation (10) separately. The term does not contain and thus can be transformed back into the time domain easily. As for the term, we introduce an auxiliary variable such that
It is easy to verify that equation (12) equals to equation (11) by eliminating 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 final dynamic equation with ABC implemented. The auxiliary DOF 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 final algebraic equation for are given as follows:
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. This 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 finite difference method in the time domain.
2.3. Reflection Coefficient Analysis
For an ideal ABC, no reflection should occur at the truncation boundary. Thus, the ratio between the artificial reflection 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 artificial reflection for the proposed ABC. One is the discretization error (equation (8)), which decreases with decreasing mesh size and is not discussed here. The other is the Pade approximation (equation (10)), whose error estimation is given by [16]:
Obviously, it approaches zero when M approaches infinity for any and vanishes automatically at .
Let the superposition of the incident and reflection waves near the ABC bewhere the absolute value of R is known as the reflection coefficient. By substituting equation (16) into equation (10), one obtains that
Substituting the dispersion relation and equation (15) into equation (17), it follows that
It can be seen the reflection coefficient approaches zero when M increases for . Besides, it converges fast near because vanishes there.
Figure 2 depicts the reflection coefficient. It shows that for frequencies near , |R| is rather small even though a small M is used. The coefficient 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.
3. ABC for FEM Euler Beams
3.1. FEM Discretization
As discussed before, the wave propagation on an infinitely 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, are the slope angle, the transversal force, and the concentrated couple, respectively. The 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
Physically, on the right-hand side are the time-dependent force and the moment applied on the truncated beam by the removed part.
3.2. ABC Transformation
Similar to the previous discussion, the point is to find the relation between the additional source terms and DOFs on boundary . 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) 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 and . In the same way, the time derivative of equation (22) can be written as
Equation (25) gives the relation between and .
It can be seen there are totally 4M + 2 equations in equations (24) and (25) and (4M) + 4 variables, i.e., 2M auxiliary variables (), 2M intermediate variables (), and 4 variables which already exist in Eq (20) (). Note that at the right-hand side of (25) has no contribution to the equation and thus does not count. Thus, with the help of (24) and (25), equation (20) is closed by introducing 2M additional variables (). 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.
3.3. 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.
By assembling equations (19) of two successive beam elements which have no external loads, one obtains
The third and fourth equations in equation (27) are
The first 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. The first few terms of the Taylor series of equation (29) near are
It can be seen the error between the relation in equation (30) and the analytical one () is , which is good enough. The 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 difference between resulted dispersion relation and the analytical one is only . 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.
4. 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 and the initial values beThe 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 . The ratio between the time step and mesh size is fixed, i.e., throughout the example. First, we use the same setup as in literature [19], i.e., mesh size and truncation boundary . The Newmark method (see e.g., [27, 28]) is employed for time integration. The 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 reference A, because it not only eliminates the boundary effect as reference A does but also reduces the mesh size effect.
The solutions at different time points are shown in Figure 4. From Figure 4(a), it can be seen the solutions fall into two groups. The waves of the FEM model and reference B are similar, and they go faster than those of the other four FDM models. The 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. The 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. This shows the FEM model effectively suppresses the boundary reflection, 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 efficient because the latter uses a global ABC. However, for the same M, the Schrodinger model leads to a significant artificial reflection. The 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, considering the former ABC is derived from that of the latter. The 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 infinitesimal, to improve the performance of the Schrodinger model, one should use a smaller h.
Different 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 reflection wave goes faster, and its amplitude is smaller. But a greater M does not seem to make any difference for this example.
The performances of the models are further quantified by the L2 norm of their differences to reference B at t = 1.5. The errors are plotted against the mesh size h in Figure 6. The Schrodinger model shows a second-order convergence, for both M = 3 and M = 30. The errors of FEM models are much smaller, and again, the influence of M is not significant. With decreasing mesh size, the errors first 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. Thus, 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. This 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 second-order 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 bywhere are the total stiffness matrix, the total mass matrix, and the displacement vector of the preserved beam part (not including auxiliary DOFs), respectively. The 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 artificial reflection wave is not eliminated perfectly after a long-time simulation. This may be because the wave components whose frequencies are far away from cannot be effectively treated by our method (see Figure 2). Using a greater M does not help much. In contrast, the MBC5 method brings more reflection energy when the main waves first reach the boundary, but after the wave goes forth and back for several times in the interval, its final energy level is lower. That is because of its good capability to absorb waves with a frequency far away from (Figure 2). The fast convolution method leads to the lowest residual energy after a long simulation time, but it performs worse than the proposed method for the first collision with the artificial boundary.
(a)
(b)
(c)
(a)
(b)
5. Conclusion
In this study, we proposed a new ABC for an infinitely long FEM Euler–Bernoulli beam which has rotational DOFs. In the ABC, the influence of the removed half infinite 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. This treatment equals to parallelly attaching several dynamic systems to the truncated system and thus is easy to be implemented in the frame of FEM.
The proposed method is local in both space and time domains. The modification to the original algebraic system occurs only near the truncation points and has a quite low rank comparing with the matrix size. The time integration scheme has no difference with that used in common dynamic systems, and no convolution is needed. The attached systems will automatically absorb the incoming waves. As a result, the method almost costs no extra computation time.
The 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 finer mesh and more approximation terms. Even for the high accuracy condition, it is still more efficient 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 fine mesh to effectively suppress the artificial reflection. However, it shows a second-order convergence as the mesh size decreases.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the Science and Technology Research Program of Chongqing Municipal Education Commission (Grant no. KJQN202301130) and NSFC (Grant no. 11702046).