Development of Post Processing for Wave Propagation Problem: Response Filtering Method

Abstract

This study develops a new response filtering approach for recovering dynamic mechanical stresses under impact loading. For structural safety, it is important to consider the propagation of transient mechanical stresses inside structures under impact loads. Commonly, mechanical stress waves can be obtained by solving Newton’s second law using explicit or implicit finite element procedures. Regardless of the numerical approach, large discrepancies called the Gibb’s phenomenon are observed between the numerical solution and the analytical solution. To reduce these discrepancies and enhance the accuracy of the numerical solution, this study develops a response filtering method (RFM). The RFM averages the transient responses within split time domains. By solving several benchmark problems and analyzing the stresses in the frequency domain, it was possible to verify that the RFM can provide an improved solution that converges toward the analytical solution. A mathematical theory is also presented to correlate the relationship between the filtering length and the frequency components of the filtered stress values.

1. Introduction

This study develops a new response filtering approach for recovering dynamic mechanical stress under impact loads. From a structural engineering point of view, a mechanical stress wave, due to structural impacts such as an impact force, explosion, or collision, that propagates through a structural medium is important. When this wave reaches weak parts of the structure, it can cause unexpected catastrophic failure or fracturing of structural components. Thus, it is important to compute and consider these stress waves in mechanical or civil engineering [1,2,3,4]. With the use of state-of-the-art computational theory and computer hardware, numerical solutions can be regarded as acceptably close to the analytical solution. However, from a theoretical perspective, some discrepancies and errors remain, and verification of the convergence to the true solution is still important. Some efforts have been undertaken to reduce the differences and the numerical errors between the numerical solutions and analytical solutions [3,5,6]. To this end, this study develops a new post-processing approach for transient stress waves named the response filtering method (RFM) with an average operator.

Numerical methods, e.g., the finite element (FE) method, finite volume method, or finite difference method, have been employed for the transient analysis of complex structures. The numerical solutions obtained using proper computational theory usually converge to the analytical solutions with a refined time step and refined mesh [3,5,6]. However, it has been observed that it is $\tilde{f}=\frac{f_s}{FL}\cdot m\pm f$ difficult for a numerical solution to describe an analytical solution that has significant discontinuity, even when proper meshes and time steps are adopted. In finite element (FE) numerical solutions, overshooting, undershooting, and oscillations occur around the discontinuity point known as the Gibbs phenomenon [1,7,8,9]. Various studies have been conducted aiming to reduce these types of discrepancies between the analytical and numerical solutions. Some relevant research can be found to improve the numerical accuracy by introducing the numerical damping [10]. For instance, the discontinuous Galerkin method and variational integrator have been proposed to accurately describe discontinuities in the numerical solutions [11,12,13,14,15]. Other approaches have been proposed to reduce the dispersive and dissipative errors in the numerical solution of partial difference equations [3,4,16,17]. Klaus and Bathe make a significant contribution to the finite element method for wave propagation problem. They present higher-order spatially discretized finite elements with an extra computational cost for the wave propagation problem [3,16] and also present new implicit and explicit integration methods with Noh [4,17]. Linear combinations of consistent and lumped mass matrices [18,19,20] and modified spatial integration rules for mass and stiffness matrices [21,22] have been considered to ensure high-fidelity solutions. Comprehensive and considerable benchmark problems for wave propagation can be found in the work of Idesman et al. [23]. However, in numerical solutions, errors due to Gibb’s phenomenon cannot be eliminated completely, and some improvements are required to the oscillation of the numerical solution [9]. A few post-processing methods have been developed to reduce these kinds of numerical errors [21,23,24].

The objective of this study is not to develop a new finite element or integration method but to develop a simple post-processing approach, i.e., the RFM, to reduce discrepancies between the numerical and analytical solutions. In this study, only one-dimensional problems are considered to calculate exact solutions (analytical solution) as references. The transient stress values in the discrete time domain are simply divided into several subset domains, and each subset domain is averaged. The subsets have a specific number of stress values, which, in this paper, is called the filtering length. The averaged stress values are set as the filtered stress values, and the original stress values are replaced with the filtered values in the corresponding time domain, as shown in Figure 1. The RFM not only efficiently solves the oscillation problem in the numerical solution by adjusting the filtering length, but it also creates discontinuous solutions that converge to the analytical solutions. Alternatively, it can be considered that the RFM can create new frequency components of the responses (especially higher frequency components). In this study, the relationship between the RFM and frequency components is mathematically verified, and a formulation is introduced to predict the generated frequency components using RFM.

This remainder of this paper is organized as follows. In Section 2, the concept, application, and limitations of finite element formulations for the mechanical wave problem and numerical integration scheme, i.e., Newmark’s scheme, are reviewed. In Section 3, to resolve the oscillation issue, i.e., Gibb’s phenomenon in the numerical solution, the RFM is presented. In Section 4, several numerical examples are presented to demonstrate that the present RFM can be utilized to improve the accuracy of numerical solutions. Finally, Section 5 provides conclusions and directions for future research.

2. Finite Element Analysis for the Transient Stress Response

2.1. Transient Finite Element Simulation

The transient response of a mechanical structure can be analyzed using the finite element procedure as follows [3,5,6]:

$${\ddot{U}}_i=M^{-1}\left[F_i-C{\dot{U}}_i-K{U}_i\right]$$

(1)

where the mass, damping, and stiffness matrices are denoted by $M$, $C$, and $K$, respectively, and the force vector at the i^th time step is denoted by $F_i$, and the displacement, velocity, and acceleration vectors at the i^th time step are denoted by $U_i$, ${\dot{U}}_i$, and ${\ddot{U}}_i$, respectively. Many numerical schemes can be applied to solve the above Newton’s law. This study implements both explicit and implicit integration method to investigate performance of the present RFM. Central-difference method and Newmark’s method are adopted as representation of explicit and implicit integration methods.

Central-difference method approximates velocity and acceleration by

$${\dot{U}}_i=\frac{1}{2\Delta t}\left(U_{i+1}-U_{i-1}\right)$$

(2)

$${\ddot{U}}_i=\frac{1}{\Delta t^2}\left(U_{i+1}-2U_i+U_{i-1}\right)$$

(3)

and $U_{i+1}$ and $U_{i-1}$ can be obtained by Taylor series about time iΔt:

$$U_{i+1}=U_i+\frac{\Delta t}{1!}{\dot{U}}_i+\frac{\Delta t^2}{2!}{\ddot{U}}_i+\frac{\Delta t^3}{3!}{\stackrel{⃛}{U}}_i+\cdots$$

(4)

$$U_{i-1}=U_i-\frac{\Delta t}{1!}{\dot{U}}_i+\frac{\Delta t^2}{2!}{\ddot{U}}_i-\frac{\Delta t^3}{3!}{\stackrel{⃛}{U}}_i+\cdots$$

(5)

Combining Equations (2)–(5) with the second-order accuracy provides the following equation:

$$U_{i+1}={\left[\frac{1}{\Delta t^2}M+\frac{1}{2\Delta t}C\right]}^{-1}\left\{F_i-K{U}_i+\frac{1}{\Delta t^2}M(2U_i-U_{i-1})+\frac{1}{2\Delta t}C{U}_{i-1}\right\}$$

(6)

In Newmark’s method [3,5,6], displacement vector of i + 1^th time step can be obtained as follows:

$${\dot{U}}_{i+1}={\dot{U}}_i+\left(\Delta t\right)\left[\left(1-\gamma \right){\ddot{U}}_i+\gamma {\ddot{U}}_{i+1}\right]$$

(7)

$$U_{i+1}=U_i+\left(\Delta t\right){\dot{U}}_i+{\left(\Delta t\right)}^2\left[\left(\frac{1}{2}-\beta \right){\ddot{U}}_i+\beta {\ddot{U}}_{i+1}\right]$$

(8)

where Δt is the time step, and the two parameters of Newmark’s scheme are β and γ [6]. With the above approximations of Newmark’s method, the following discretized equations in the time domain can be obtained:

$$K^{\prime }{U}_{i+1}={F^{\prime }}_{i+1},$$

(9)

$$K^{\prime }=\frac{M}{\Delta t^2\beta }+\frac{\gamma }{\Delta t\beta }C+K,$$

(10)

$$\begin{array}{cc}\hfill {F^{\prime }}_{i+1}& =F_{i+1}+\left[\frac{M+\Delta t\gamma C}{\Delta t^2\beta }\right]U_i+\left[\frac{M+\Delta t\gamma C}{\Delta t\beta }-C\right]{\dot{U}}_i\hfill \\ & +\left[\left(\frac{1}{2\beta }-1\right)\left(M+\Delta t\gamma C\right)-\Delta t\left(1-\gamma \right)C\right]{\ddot{U}}_i.\hfill \end{array}$$

(11)

In this study, the coefficients of Newmark’s method, $\beta$ and $\gamma$, are set to $\frac{1}{4}$ and $\frac{1}{2}$, respectively for unconditional stability.

The dynamic stress at the ith time step can be calculated from the displacements at the ith time step as follows:

$${\sigma }_i={\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right]}_i=DB{U}_i^{}.$$

(12)

$$D=\frac{E}{1-{\nu }^2}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right]$$

(13)

$$2\mathrm{D}\mathrm{plane}\mathrm{stress}:{\sigma }_i=\sqrt{{\sigma }_x^2-{\sigma }_x{\sigma }_y+{\sigma }_y^2+3{\tau }_{xy}}.$$

(14)

The plane stress (von Mises stress), stress vector, displacement vector, constitutive matrix, and strain–displacement matrix at the ith time step are denoted by ${\sigma }_i$, ${\sigma }_i$, $U_i$, $D$, and $B$, respectively. The Young’s modulus and Poisson’s ratio are denoted by $E$ and $\nu$, respectively.

2.2. Wave Propagation Analysis and Gibb’s Phenomenon

Facilitated by the developments in hardware and software, finite element simulations of transient systems have been applied in many application areas [3]. Despite their successful applications, one of the remaining issues with transient finite element simulations is the Gibb’s phenomenon, with its associated overshooting and undershooting phenomena [5]. This issue becomes a serious problem in structural analyses. To illustrate this problem, let us consider the stress analysis in Figure 2 as an example [5]. The right side of a 20-in-long rectangular bar ($L_{\mathrm{T}}$) with a cross-sectional area (A) of 1 in² is fixed, and a step load of 100 lb is applied in the x-direction on the left side of the bar at $t=0$. The Young’s modulus (E), density (ρ), and Poisson’s ratio ($\nu$) of the bar are $30\times {10}^6$ psi, $7.4\times {10}^{-4}$ lb-s²/in⁴, and zero (for a one-dimensional problem), respectively. In this study, the analysis domain is discretized using a $40\times 1$ mesh of 4-node rectangular elements (Q4). The mesh size of this example is chosen by following the parameters and setting of the reference [5]. In addition, it is worth to note that the issues of the stability and the accuracy of the numerical integration schemes should be considered, and the time of this example is chosen to consider the stabilities of the numerical integration schemes. The time step (sampling time), $\Delta t$, is set to $4.0\times {10}^{-7}\mathrm{s}$ smaller than the critical time ($\Delta t_{cr}=L_c/\mathrm{v}$$(=2.483\times {10}^{-6})$, $\mathrm{v}:\sqrt{E/\rho }$, $L_c$: characteristic length) for the stability of the central-difference method.

The transient finite element simulation determines the transient solution with application of a sudden load to the cross section. The stress of the cross section in the steady-state solution becomes 100 psi. Figure 3 shows the mechanical stress value calculated at x = 9.75 in (20th element). The solid and the dotted lines represent the analytical and numerical solutions, respectively. With mechanical properties we can calculate speed of wave propagation ($\mathrm{v}=\sqrt{E/\rho }$$=2.012\times {10}^5\mathrm{in}/\mathrm{s}$) and time stress curve. Before the arrival time ($t_1$$=\frac{9.75}{2.012\times {10}^5}$$=4.7214\times {10}^{-5}\mathrm{s}$), the stress values at x = 9.75 in are zero because the stress wave has not yet arrived. After the first arrival time and before the next arrival time ($t_2$$=\frac{9.75+2\times (20-9.75)}{2.012\times {10}^5}$$=1.5158\times {10}^{-4}\mathrm{s}$), the analytical stress value becomes 100 psi. After $1.5148\times {10}^{-4}\mathrm{s}$, the analytical stress value increases to 200 psi owing to the arrival of the reflected wave at the right side. The steady-state stress is 100 psi. The stress values obtained with the explicit and implicit method are plotted with exact solution in Figure 3. The overall tendency is similar to that of the analytical solution, but some spurious oscillations naturally occur. It is observed that larger oscillations occur at the areas of sharp transition ($t_1=4.7214\times {10}^{-5}\mathrm{s}$ and $t_2=1.5158\times {10}^{-4}\mathrm{s}$). Comparing the analytical and numerical solutions, it seems that the numerical solution obtained using the transient finite element simulation contains additional frequencies. However, later we will demonstrate that the numerical solution actually cannot carry valuable information in the frequency domain compared with the analytical solution. From a mathematical perspective, this phenomenon is referred to as Gibb’s phenomenon and has been widely studied (see [4,5] and the references therein). In finite element theory, considering this phenomenon, the stability issue has been researched [4,5,9,14,16,20]. The present study revisits the Gibb’s phenomenon and presents a heuristic approach for improving the accuracy of finite element simulations.

3. Response Filtering Method for Improved Transient Stress

This section presents a new engineering approach, i.e., the RFM, to improve the transient stress calculated using the finite element approach. In static analysis, spatial-averaging approaches for the stress have been proposed, but to the best of our knowledge, an averaging scheme in the time domain has not been previously proposed.

3.1. Response Filtering Method

As its name implies, the RFM filters the transient stress values in the time domain. A standard FE simulation is conducted to compute the transient stress values. As illustrated in the previous section, the Gibb’s phenomenon hinders our ability to obtain accurate solutions using this method. By investigating the responses and their characteristics, this study proposes the concept of averaging the erroneous transient stress values using the RFM. Through averaging, or filtering, it is expected that the oscillations can be cancelled out, and thus the accuracy of the computed transient stress values can be improved. This concept has been applied for the stress analysis of static structure systems and is implemented in various commercial software [25,26]. As it is a relatively simple concept, it is easy to implement in the post-processing stage of any FE package. The present study is confined to averaging (or filtering) in the time domain in order to illustrate the benefit of the proposed RFM.

Figure 1 illustrates a schematic diagram of the proposed RFM. First, the stress and displacement responses are computed in the time domain using a conventional numerical approach, i.e., Newmark’s method. The resulting transient solutions are approximations of the analytical solution, and oscillations inevitably occur due to the Gibb’s phenomenon. The displacement or stress response at each time step in the numerical approach can be described by (12), and without the loss of generality, the stress values are employed to illustrate the RFM.

$$\sigma =\left[{\sigma }_1,{\sigma }_2,\cdots {\sigma }_{N-1},{\sigma }_N\right](N:\mathrm{number}\mathrm{of}\mathrm{time}\mathrm{steps}),$$

(15)

where the (von Mises) stress of target element is denoted by σ_i, and the subscript i denotes the time step and set of stresses with respect to time step is σ. The transient stress values can be divided into NR regions and renumbered for convenience as follows:

$$\begin{array}{cc}\hfill \sigma & =[\underset{\mathrm{the}1\mathrm{st}\mathrm{region}}{\underbrace{{\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_{FL}}},\underset{\mathrm{the}2\mathrm{nd}\mathrm{region}}{\underbrace{{\sigma }_{FL+1},\cdots ,{\sigma }_{2FL}}},\cdots ,\underset{\mathrm{the}NR-\mathrm{th}\mathrm{region}}{\underbrace{{\sigma }_{N-FL+1},\cdots ,{\sigma }_N}}]\hfill \\ & =[\underset{\mathrm{the}1\mathrm{st}\mathrm{region}}{\underbrace{{\sigma }_1^1,{\sigma }_2^1,\cdots ,{\sigma }_{FL}^1}},\underset{\mathrm{the}2\mathrm{nd}\mathrm{region}}{\underbrace{{\sigma }_1^2,\cdots ,{\sigma }_{FL}^2}},\cdots ,\underset{\mathrm{the}NR-\mathrm{th}\mathrm{region}}{\underbrace{{\sigma }_1^{NR},\cdots ,{\sigma }_{FL}^{NR}}}].\hfill \end{array}$$

(16)

The “region” defined in this research refers the time step region in which the solutions are averaged.

To avoid confusion in notation, the ith stress at the jth region is denoted by ${\sigma }_i^j$. The filter length, which is equal to the number of stress values in each region, is denoted by FL. The filtered stresses or responses are obtained by averaging with a weight factor as follows:

$$\tilde{\sigma }=[\underset{\begin{array}{c}\mathrm{the}1\mathrm{st}\\ \mathrm{region}\end{array}}{\underbrace{{\tilde{\sigma }}_{}^1}},\cdots ,\underset{\begin{array}{c}\mathrm{the}k-\mathrm{th}\\ \mathrm{region}\end{array}}{\underbrace{{\tilde{\sigma }}_{}^k}},\cdots ,\underset{\begin{array}{c}\mathrm{the}NR-\mathrm{th}\\ \mathrm{region}\end{array}}{\underbrace{{\tilde{\sigma }}_{}^{NR}}}],$$

(17)

$${\tilde{\sigma }}_{}^k=[{\tilde{\sigma }}_1^k,\cdots ,{\tilde{\sigma }}_j^k,\cdots ,{\tilde{\sigma }}_{FL}^k],({\tilde{\sigma }}_1^k={\tilde{\sigma }}_j^k={\tilde{\sigma }}_{FL}^k),$$

(18)

$${\tilde{\sigma }}_j^k=\frac{{\displaystyle \sum _{i=1}^{FL}{w}_i^k{\sigma }_i^k}}{{\displaystyle \sum _{i=1}^{FL}{w}_i^k}},k=1,2,\cdots ,NR,j=1,2,\cdots ,FL(N=FL\times NR,FR=\frac{FL}{N}\times 100),$$

(19)

where the set of original and filtered stresses (responses) at all time steps are denoted by σ and $\tilde{\sigma }$, respectively. The ith filtered stress and the weighting factor in the kth domain are denoted by ${\sigma }_j^k$ and $w_j^k$, respectively. The filtering length, filtering ratio and number of the filtered regions are FL, FR, and NR, respectively. In the present study, the weight factors are set as equal to one. For example, the filtered stress in (17) with an FL value of two becomes the following:

$$\tilde{\sigma }=\left[\underset{{\tilde{\sigma }}_1}{\underbrace{\frac{{\sigma }_1+{\sigma }_2}{2},\frac{{\sigma }_1+{\sigma }_2}{2}}},\underset{{\tilde{\sigma }}_2}{\underbrace{\frac{{\sigma }_3+{\sigma }_4}{2},\frac{{\sigma }_3+{\sigma }_4}{2}}},\cdots ,\underset{{\tilde{\sigma }}_{FL}}{\underbrace{\frac{{\sigma }_{N-1}+{\sigma }_N}{2},\frac{{\sigma }_{N-1}+{\sigma }_N}{2}}}\right]$$

(20)

This simple filtering approach, which is similar to the stress average in static FE analysis, is carried out in the time domain rather than in the spatial domain. To validate the performance of the proposed method, the benchmark problem in Section 2 can be revisited. The transient stresses of the benchmark problem in Figure 2 are filtered using the proposed RFM by varying the ratio of filter data, FR, in Figure 4. Moreover, the following norm-based absolute error for the stress (response) is employed to evaluate the performance of the present method for transient analysis:

$$\xi ={\displaystyle \int \left|\overline{\sigma }(t)-\sigma (t)\right|}dt\mathrm{or}{\displaystyle \int \left|\overline{\sigma }(t)-\tilde{\sigma }(t)\right|dt},$$

(21)

where the exact (analytical), numerical, and filtered numerical solutions at time t are denoted by $\overline{\sigma }(t)$ and $\tilde{\sigma }(t)$, respectively. Figure 4 presents the tendency of the accuracy improvement of the RFM with varying filtering ratio, FR.

Table 1. (a) Accuracy improvement of the RFM with different numerical integration method (${\mathrm{at}20}^{\mathrm{th}}\mathrm{element},\Delta t=\Delta t_1=4.0\times {10}^{-7}\mathrm{s}$). (b) Accuracy improvement of the response filtering method with different numerical integration method (${\mathrm{at}20}^{\mathrm{th}}\mathrm{element},\Delta t=\Delta t_1\times {10}^{-1}\mathrm{s}$).

(a)
Central-difference method
Filtering ratio (%)	w/o RFM	0.5	1	1.5	2	2.5	3	3.5
Error indicator $\xi (\times {10}^{-3})$	1.8472	1.8283	1.7097	1.5967	1.3808	1.2124	1.2241	1.1895
Number of over (under) shoots	67	65	58	45	28	21	19	19
Newmark’s method
Filtering ratio (%)	w/o RFM	0.5	1	1.5	2	2.5	3	3.5
Error indicator $\xi (\times {10}^{-3})$	1.7223	1.7053	1.5941	1.4387	1.3123	1.1205	1.1937	1.1340
Number of over (under) shoots	65	65	65	49	27	20	20	19
(b)
Central-difference method
Filtering ratio (%)	w/o RFM	0.5	1	1.5	2	2.5	3	3.5
Error indicator $\xi (\times {10}^{-3})$	1.7883	1.7539	1.6543	1.4968	1.3046	1.3640	1.2136	1.0409
Number of over (under) shoots	65	64	59	41	26	19	19	17
Newmark’s method
Filtering ratio (%)	w/o RFM	0.5	1	1.5	2	2.5	3	3.5
Error indicator $\xi (\times {10}^{-3})$	1.7876	1.7515	1.6510	1.4962	1.3064	1.3618	1.2165	1.0495
Number of over (under) shoots	65	64	59	41	26	19	19	17

shows the performance of RFM with respect to filtering ratio (FR), the error valuer in Equation (21) is reduced by increasing FR by around 2% and also the number of oscillations is reduced as shown in

Table 2. Accuracy improvement of the RFM with different numerical integration method (${\mathrm{at}1}^{\mathrm{st}}\mathrm{element},\Delta t=4.0\times {10}^{-7}\mathrm{s}$)

Central-difference method
Filtering ratio (%)	w/o RFM	0.5	1	1.5	2	2.5	3	3.5
Error indicator $\xi (\times {10}^{-3})$	1.6710	1.6400	1.5039	1.3723	1.0918	1.1012	0.9666	0.9809
Number of over (under) shoots	73	73	73	47	19	15	23	23
Newmark’s method
Filtering ratio (%)	w/o RFM	0.5	1	1.5	2	2.5	3	3.5
Error indicator $\xi (\times {10}^{-3})$	1.6025	1.5811	1.4269	1.3339	1.1746	0.9752	0.9827	0.9688
Number of over (under) shoots	71	71	71	51	21	13	21	23

. The table shows that the RFM could show similar performance regardless of integration methods. Approximately 2% of FR may show significant improvement in accuracy under different numerical conditions.

Because the proposed method yields an artificially generated piecewise constant solution, it has to be validated not only in time but also in space. Therefore, the benchmark problem in Figure 2 was conducted once more with different space condition. Both exact and numerical stress responses at x = 0.25 in (1st element) were calculated by the same procedure as the previous example. Lots of spurious oscillation was found around at second arrival time in the finite element method ($t=\frac{20+19.75}{2.012\times {10}^5}$) in Figure 5. The detail performance of RFM with different filtering ratio, FR, is written in Table 2. Similar to the previous example, the FR (filtering ratio) value around 2% shows the improvement in time and space.

3.2. Analysis of the Accuracy Improvement with the Response Filtering Method

As shown in the previous section, the simple averaging procedure with the RFM results in significant improvements in the prediction compared to the analytical solution by generating piecewise constant solutions. For the sake of conducting an accuracy analysis of the proposed method, this subsection investigates the solutions using the Fourier transformation.

Without the loss of generality, a simple signal, $y$, with a 50 Hz frequency component was generated for 0.1 s with 100 samplings in Figure 6. In other words, $y(t)=\sin (2\pi 50t)$, and the sampling time, $T_s=0.001\mathrm{s}$. Figure 6a (left) and (right) shows the original signal and its Fourier transformation, respectively, with the peak at 50 Hz. Figure 6b–d shows the filtered signals with different filtering lengths and their transformations. These figures illustrate that the Fourier transformations of the filtered signals contain some additional frequency components. For example, with FL = 2, the filtered signal bears additional frequency components at 450 Hz. With FL = 4, frequencies around 200, 300, and 450 Hz are added, while with FL = 5, frequencies around 150, 250, 350, and 450 Hz are added. Furthermore, the signals in the time domain become rectangular in shape owing to these added frequency components. These analyses indicate that the filtering procedures used for the time signals actually generate additional frequencies. In addition, it is found that there is a rule for the emerging frequencies. For example, Figure 6a shows that the original signal has a single frequency component of 50 Hz, whereas the filtered signal (FL = 2) has an additional frequency component at 450 Hz in addition to the frequency component (50 Hz) of the original signal, as shown in Figure 6b.

Prior to determining the relationship between the filtering length and the added frequency components, it should be noted that the discrete Fourier transform (DFT) defines the following relationship between the analysis time and frequency:

$$T_s=\frac{T_0}{N},f_0=\frac{f_s}{N},f_0=\frac{1}{T_0},f_s=\frac{1}{T_s},$$

(22)

where the number of time step, total analysis time and the sampling time ($\Delta t$) are denoted by $N$, $T_0$, and $T_s$, respectively, and the frequency step ($\Delta f$) and the sampling frequency are denoted by $f_0$ and $f_s$, respectively. Using the above notations, this study formulates the following theory based on the DFT.

Theorem 1.

Consider the signal, $\tilde{y}$, filtered from original signal, $y$, which has only a single frequency component, $f$, with a filtering length FL. With sampling frequency, $f_s$, the Fourier transformation of the filtered signal has the following frequency component:

$$\tilde{f}=\frac{f_s}{FL}\cdot m\pm f\mathrm{and}0\le \tilde{f}\le \frac{f_s}{2},$$

(23)

where m and $\tilde{f}$ denote an arbitrary integer and the frequency components generated by the frequency filtering method (i.e., the RFM), respectively. The frequency components generated by the filtering scheme in DFT can be mathematically found through the following proof.

Proof. 

The N-point DFT of discrete signal $y_k$ is obtained as follows:

$$Y_n=\frac{1}{N}{\displaystyle \sum _{k=0}^{N-1}{y}_k{e}^{-j\frac{2\pi n}{N}\cdot k}}.$$

(24)

Assuming $y_k=\exp (j\frac{2\pi f^{*}}{N}k)$ (a signal with a single frequency component, $f$, where $f^{*}=\frac{f}{f_0}$), the above DFT can be summarized as follows:

$$Y_n=\frac{1}{N}{\displaystyle \sum _{k=0}^{N-1}{e}^{-j\frac{2\pi }{N}(n-f^{*})\cdot k}={\delta }_f(n)},{\delta }_f\left(n\right)=\{\begin{array}{c}1(f^{*}=n)\\ 0(f^{*}\ne n)\end{array},$$

(25)

where ${\delta }_f\left(n\right)$ is the Kronecker delta function, and Equation (25) is used for the proof.

$${\displaystyle \sum _{k=0}^{N-1}a{r}^k}=\{\begin{array}{cc}\hfill Na& (r=1)\hfill \\ \hfill \frac{a\left(1-r^N\right)}{1-r}& (r\ne 1)\hfill \end{array},$$

(26)

$$r=\left(e^{-j\frac{2\pi n}{N}}\right)=1,(\frac{n}{N}=m,m\mathrm{is}\mathrm{arbitrary}\mathrm{integer}),$$

(27)

$${\left(e^{-j\frac{2\pi n}{N}}\right)}^N=1.$$

(28)

The condition of Equation (27) can be satisfied only if $m=0$ because $N$ is always larger than $n$ in Equation (25).

The DFT of the filtered signal, $\tilde{y}$, obtained from the RFM with FL can be formulated as follows:

$$\begin{array}{cc}\hfill {\tilde{Y}}_n& =\frac{1}{N}{\displaystyle \sum _{l=0}^{N/FL-1}{\tilde{y}}_{FL\cdot l}\cdot \left[e^{-j\frac{2\pi n}{N}FL\cdot l}+e^{-j\frac{2\pi n}{N}\left(FL\cdot l+1\right)}+\cdot \cdot \cdot +e^{-j\frac{2\pi n}{N}\left\{FL\cdot l+FL-1\right\}}\right]}\hfill \\ & =\frac{1}{N}{\displaystyle \sum _{l=0}^{N/FL-1}\frac{\left[y_{FL\cdot l}+y_{FL\cdot l+1}+\cdot \cdot \cdot +y_{FL\cdot l+FL-1}\right]}{FL}\cdot \left[e^{-j\frac{2\pi n}{N}FL\cdot l}+e^{-j\frac{2\pi n}{N}\left(FL\cdot l+1\right)}+\cdot \cdot \cdot +e^{-j\frac{2\pi n}{N}\left\{FL\cdot l+FL-1\right\}}\right]}\hfill \\ & =\frac{1}{N\cdot FL}{\displaystyle \sum _{l=0}^{N/FL-1}{y}_{FL\cdot l}(y_0+\cdot \cdot \cdot +y_{FL-1})\cdot e^{-j\frac{2\pi n}{N}FL\cdot l}(e^0+\cdot \cdot \cdot +e^{-j\frac{2\pi n}{N}\left\{FL-1\right\}})}\hfill \\ & =\frac{1}{N\cdot FL}\left({\displaystyle \sum _{i=0}^{FL-1}{e}^{j\frac{2\pi f^{*}}{N}\cdot i}}\cdot {\displaystyle \sum _{i=0}^{FL-1}{e}^{-j\frac{2\pi n}{N}\cdot i}}\right)\cdot {\displaystyle \sum _{l=0}^{N/FL-1}{e}^{j\frac{2\pi FL{f}^{*}}{N}\cdot l}{e}^{-j\frac{2\pi FLn}{N}\cdot l}}\hfill \\ & =\frac{1}{N\cdot FL}\cdot \frac{1-e^{j\frac{2\pi FL}{N}{f}^{*}}}{1-e^{j\frac{2\pi }{N}{f}^{*}}}\cdot \frac{1-e^{-j\frac{2\pi FL}{N}n}}{1-e^{-j\frac{2\pi }{N}n}}\cdot {\displaystyle \sum _{l=0}^{N/FL-1}{e}^{-j\frac{2\pi FL(n-f^{*})}{N}\cdot l}}\hfill \\ & =\frac{1}{N\cdot FL}\cdot \frac{1-e^{j\frac{2\pi FL}{N}{f}^{*}}}{1-e^{j\frac{2\pi }{N}{f}^{*}}}\cdot \frac{1-e^{-j\frac{2\pi FL}{N}n}}{1-e^{-j\frac{2\pi }{N}n}}\cdot {\displaystyle \sum _{l=0}^{N/FL-1}{r}^l},(r=e^{-j\frac{2\pi FL(n-f^{*})}{N}})\hfill \end{array}$$

(29)

Therefore, the Fourier transform of the filtered signal is given by Equations (25)–(28):

$${\tilde{Y}}_n=\{\begin{array}{cc}\hfill \frac{1}{F{L}^2}\frac{1-e^{j\frac{2\pi FL}{N}{f}^{*}}}{1-e^{j\frac{2\pi }{N}{f}^{*}}}\frac{1-e^{-j\frac{2\pi FL}{N}n}}{1-e^{-j\frac{2\pi }{N}n}}& (r=1)\hfill \\ \hfill 0\hfill & (r\ne 1)\hfill \end{array}.$$

(30)

The condition for $r=e^{-j\frac{2\pi FL(n-f^{*})}{N}}=1$ can be simplified using Equation (27) as follows:

$$\frac{FL}{N}(n-f^{*})=m,$$

(31)

$$n=\frac{N}{FL}m+f^{*},$$

(32)

where m represents an arbitrary integer. The integer domain can be transformed to the discrete frequency domain by multiplying the frequency step, $f_0$, as follows:

$$\begin{array}{cc}\hfill f_0\cdot n& =f_0\cdot \frac{N}{FL}m+f_0\cdot f^{*}\hfill \\ & =\frac{f_s}{FL}m+f\hfill \end{array}$$

(33)

The frequency components, $f_0\cdot n$, generated by the RFM can be simply denoted by $\tilde{f}$, as in (23). With $y_k=\sin \left(k\frac{2\pi f^{*}}{N}\right)=\sin \left(k\frac{2\pi f}{f_s}\right)=\frac{e^{jk\frac{2\pi f}{f_s}}-e^{-jk\frac{2\pi f}{f_s}}}{2j},(f^{*}=f/f_0,N=f_s/f_0)$, the DFT for this example can be expressed as follows:

$$Y_n=\{\begin{array}{cc}\hfill \frac{1}{2j}\frac{1}{F{L}^2}\frac{1-e^{jFL\frac{2\pi f}{f_s}}}{1-e^{j\frac{2\pi f}{f_s}}}\frac{1-e^{-jFL\frac{2\pi \tilde{f}}{f_s}}}{1-e^{-j\frac{2\pi }{f_s}}}& (\tilde{f}=\frac{f_s}{FL}m+f)\hfill \\ \hfill \frac{1}{2j}\frac{1}{F{L}^2}\frac{1-e^{jFL\frac{2\pi f}{f_s}}}{1-e^{j\frac{2\pi f}{f_s}}}\frac{1-e^{-jFL\frac{2\pi \tilde{f}}{f_s}}}{1-e^{-j\frac{2\pi }{f_s}}}& (\tilde{f}=\frac{f_s}{FL}m-f)\hfill \\ \hfill 0\hfill & \left(\mathrm{otherwise}\right)\hfill \end{array}$$

(34)

□

The range of (33) corresponds exactly to Theorem 1. Therefore, Theorem 1 is mathematically validated.

The locations of the newly created frequency components can be predicted using the sampling frequency (or sampling time), filtering length, and original frequency component with (23). For instance, a sampling frequency of 1000 Hz is used for the example in Figure 6. Therefore, the generated frequency components with respect to the filtering length can be calculated as summarized in

Table 3. Generated frequencies predicted using Theorem 1: $0\le \tilde{f}\le \frac{f_s}{2}$ ($N=100$, $T_0=0.1\mathrm{s}$, $f_s=\frac{N}{T_0}=1000\mathrm{Hz}$).

$\mathbf{Filtering}\mathbf{Length}\left(\mathit{F}\mathit{L}\right)$	$\mathbf{Generated}\mathbf{Frequencies}(\tilde{\mathit{f}})$
1	50 Hz
2	50 Hz, 450 Hz
4	50 Hz, 200 Hz, 300 Hz, 450 Hz
5	50 Hz, 150 Hz, 250 Hz, 350 Hz, 450 Hz

. It is verified that Theorem 1 predicts the newly generated frequency components with the RFM in Figure 6 and Table 3.

Based on the above consideration, we revisit the previous example in Figure 2 with Fourier transformation. The time domain mechanical stress response at x = 9.75 in is transformed to frequency domain by DFT (the time domain response is calculated by Newmark’s method with $\Delta t=4.0\times {10}^{-7}\mathrm{s}$). As shown in Figure 7, when the RFM is not applied, only the low-frequency components can be obtained. However, with the current RFM, several high-frequency components, which only can be obtained by piecewise constant components (red line) are generated. In addition, RFM, with a filtering ratio of 2% showing significant performance in previous time domain analyses, creates additional new frequency components that can cover all frequency ranges.

4. Numerical Examples

This section illustrates examples to demonstrate the benefits of the RFM in terms of accuracy as well as its limitations.

Example 1. 

Bar with three different segments

For the first example, a one-dimensional bar with different cross-section areas is considered. The problem geometry is composed of three different segments, as shown in Figure 8a, and each segment has a different cross-sectional area:$0.5\times 1$, $1.5\times 1$, and $0.5\times 1{\mathrm{in}}^2$, respectively. A step axial load of 100 psi in the x-direction is applied on the left side of the bar, and the example structure is discretized using $0.5\mathrm{in}\times 0.5\mathrm{in}$ Q4 element meshes. The detailed geometry, boundary conditions, and material properties are illustrated in Figure 8.

Different from the previous benchmark example, the impedance mismatches at $x=8$ and 11 in cause the occurrence of reflection and transmission waves owing to the change in the cross-sectional area. It should be noted that when a wave that propagates through a medium encounters the boundary of another medium having a different impedance, it produces reflection and transmission waves, as shown in Figure 9 [27].

The wave propagation impedance, Z, transmission coefficient, T, and reflection coefficient, R, can be defined by the following equations:

$$Z=\frac{p}{v},$$

(35)

$$T=\frac{2Z_2}{Z_2+Z_1}=\frac{2}{1+Z_1/Z_2},$$

(36)

$$R=\frac{Z_2-Z_1}{Z_2+Z_1}=\frac{1-Z_1/Z_2}{1+Z_1/Z_2},$$

(37)

where $p$ and $v$ denote the pressure (stress) and velocity vectors of the wave propagation, respectively. The transmission coefficient, T, and reflection coefficient, R, have the following relationship:

$$1+R=T.$$

(38)

Based on Equations (35)–(38), the values for each coefficient can be obtained in this example. First, the impedance of Domain 1 and Domain 3 is three times greater than that of Domain 2. The transmission coefficients of the $(+)xdir$ wave are 1/2 and 3/2 at $x=8\mathrm{in}$ and $x=11\mathrm{in}$, respectively, and the reflection coefficient can be obtained from (37) accordingly ($R_{x=8\mathrm{in}}=-1/2$, $R_{x=11\mathrm{in}}=1/2$). In the same way, the transmission and reflection coefficients for the $(-)xdir$ wave can also be obtained. The analytical stress value at $x=9.75\mathrm{in}$ can be obtained through a step-by-step procedure, as shown in Figure 10.

Figure 11a shows the differences between the numerical solution and the analytical solution before applying the RFM, while Figure 11b–d shows the shape of the filtered numerical solution as the filtering length of the RFM changes. The shape of the filtered numerical solution changes according to the size of the filtering length of the RFM. The filtered numerical solution gradually becomes increasingly similar to the analytical solution as the filtering length increases. In addition, it can be confirmed that, in this example, the oscillation of the numerical solution decreases or disappears with the application of the RFM. However, the larger filtering length does not always guarantee better accuracy improvement (see

Table 4. Accuracy improvement of the RFM with varying filtering length in Example 1.

Newmark’s Method
Filtering ratio (%)	w/o RFM	0.5	1	1.5	2	2.5	3	3.5
Error indicator $\xi (\times {10}^{-3})$	2.5622	2.5463	2.4712	2.3883	2.3330	2.0453	2.0096	2.0617
Number of over (under) shoots	40	40	40	34	28	22	22	22

). It should be noted that the RFM also plays an important role in the frequency domain. To explain this, the frequency components of the numerical solution without the RFM and analytical solution are presented in Figure 12.

Figure 12 shows that there is an enormous difference between the frequency components of the analytical solution and the numerical solution without the RFM. In particular, the analytical solution contains all of the frequency components, whereas the numerical solutions contain only the low-frequency components. Generally, a finer sampling time is required to consider higher frequencies, thereby increasing the computational cost. However, even finer sampling does not always guarantee inclusion of the higher-frequency components that appear in the analytical solution, as shown in Figure 12. In this study, however, it was discovered that a higher frequency range can be covered by using the RFM without modification of the sampling time.

As shown in Figure 13a, when the RFM is not applied, only the low-frequency components are present (black line). However, with the current RFM, several high-frequency components are generated (red line). In addition, increasing the filtering length creates additional new frequency components (see Figure 13b), and if the filtering length is greater than a certain length (in this example, greater than 8), the filtered numerical solutions can expand to cover all the frequency components (see Figure 13c). Note that the frequency components of the filtered numerical solution are changed without changing the sampling time. In general, the frequency components of the numerical solution are affected by the sampling time in a discrete system. In other words, with a chosen sampling time, the higher-frequency component is also determined. Applying the RFM can increase the higher frequency, which cannot be obtained by changing the sampling time.

Example 2. 

Complex geometric example

As the second example, Figure 14 considers a structure with many branches. Each bar has the same material properties as the previous problems: $E=30\times {10}^6$ psi, $\rho =7.4\times {10}^{-4}$ lb-s²/in⁴, $\nu =0$. A force is applied at the right arm, and the other arms are clamped. The stress at the center of the analysis domain is computed with and without the RFM.

Figure 14 shows the detailed shape of the complex beam structure with three types of rectangular beams. The example structure is discretized using $0.5\mathrm{in}\times 0.5\mathrm{in}$ Q4 elements. It is not clear whether exact solution to this problem can be calculated. At least, we can find the shape of the wave propagation by the numerical method in Figure 15. The numerical solution and the results of applying the RFM to the numerical solution ($\Delta t=4.0\times {10}^{-7}\mathrm{s}$) are shown in Figure 16. As in the previous examples, spurious oscillation of the numerical solution is observed.

This time, the effect of the RFM with different filtering length, FL (not FR), is describe in Figure 16; Figure 17 to investigate feature of the RFM with respect to FL. In the frequency domain, the increase in the frequency components can be predicted by Theorem 1. For example, the noticeable frequency component of numerical (FEM) solution is distributed from 0 to $4\times {10}^5$ Hz, and we can expect an additional frequency component by RFM (FL = 2) as $8\times {10}^5$ to $12\times {10}^5$ Hz by Theorem 1: $\tilde{f}=\frac{f_s}{FL}\cdot m\pm f$$=\frac{25\times {10}^5}{2}\cdot 1\pm 4\times {10}^5$. Without the RFM, the frequency components in the low-frequency region are dominant, whereas high-frequency components are generated with the RFM. This example shows that the RFM can work well even if the geometry of a structure is complex.

5. Conclusions

For future research, it is necessary to conduct a theoretical study to reveal a reasonable filtering ratio as well as a reasonable filtering length.

In numerical integration, the transient solutions of Newton’s second law obtained through either the implicit or explicit method inevitably exhibit undershooting or overshooting phenomena. Although both the stability and accuracy of the numerical methods are important, it is rare to investigate the accuracy of the time integration scheme compared with the stability issue. Thus, these numerical discrepancies are often neglected through the safety factor in mechanical or civil engineering. As a remedy for these discrepancies, this study develops an RFM for averaging the mechanical stress values in the time domain that improves the accuracy of the numerical solution. The RFM is a simple post-process averaging method to filter out the oscillations observed in the numerical solution and makes piecewise constant solutions of wave propagation problems. The performance of the present RFM was tested in both space and time with both the explicit and the implicit integration method. It was empirically found that a filtering ratio, FR, around 2% can improve the accuracy of numerical solutions. In addition, the frequency response functions for stress values before and after RFM application show that higher-frequency signals are added after filtering by generating piecewise constant responses. These higher frequency signals are essential to represent ramp phenomena in time domain. Though some benchmark problems, these higher-frequency components were investigated. One of the limitations of the present study is that it cannot provide the optimum filtering length for generalized wave propagation problems by a rigorous mathematical verification. In the same way, if the FR is too small it would not get rid of the oscillations, while if too large it would give big chunks of the solution into constant. By adopting a proper filtering ratio (around 2%) for the averaging operator in the RFM, it was observed that the averaged stress values converge toward the analytical stress values. For future research, it is necessary to conduct a theoretical study to reveal a reasonable filtering ratio as well as a reasonable filtering length.