# Oscillating Nonlinear Acoustic Waves in a Mooney-Rivlin Rod
## 3 Numerical Analysis  
### Finite Element Model  
Consider a semi-infinite one-dimensional rod with harmonic force loading applied at the left end of the rod, Figure 1.  
The considered problem was solved by applying the explicit Lax-Wendroff [39] numerical scheme modified in [40] (Sec. 7.16.2.3) for a two-step predictor-corrector time-integration coupled with the finite element model for spatial discretization; the corresponding increments are denoted as \(\Delta x\) and \(\Delta t\). To achieve conditional stability, the Courant-Friedrichs-Lewi (CFL) condition was applied for choosing a numerically stable time-increment [41]:  
\[\Delta t\leq\frac{\Delta x}{\max(\varepsilon)}, \tag{29}\]  
where \(c(\varepsilon)\) is defined by Equation (18). Direct verification shows that in a range of admissible negative strain \(\varepsilon\in(-1;\ 0)\), the \(c(\varepsilon)\) value is unlimited tending to infinity at \(\varepsilon\to-1\) regardless of (positive) values of \(C_{01},\ C_{10}\).  
Another remark concerns the Lax-Wendroff energy preserving scheme; in a two-step predictor-corrector form [40] its first step may be written as:  
\[\begin{array}{rl}u\left(x_{j+1/2};t_{n+1/2}\right)=&\frac{1}{2}\left(u\left(x_{ j+1};t_{n}\right)+u\left(x_{j};t_{n}\right)\right)-\\ &\frac{\Delta t}{2\Delta x}\left(F\left(u\left(x_{j+1};t_{n}\right)\right)-F \left(u\left(x_{j};t_{n}\right)\right)\right),\end{array} \tag{30}\]  
and the corresponding second step written as:  
\[u\left(x_{j};t_{n+1}\right)=u\left(x_{j};t_{n}\right)-\frac{\Delta t}{\Delta x }\left(F\left(u\left(x_{j+1/2};t_{n+1/2}\right)\right)-F\left(u\left(x_{j-1/2}; t_{n+1/2}\right)\right)\right). \tag{31}\]  
Herein, \(x_{j}\) are the spatial grid nodes and \(t_{n}\) are the nodes along the time-axis. It should also be noted that Lax-Wendroff scheme, along with the Godunov scheme, is recommended for solving nonlinear hyperbolic differential equations with discontinuities in coefficients [40] and evaluating discontinuities in stress and strain fields [40]; see also [6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23] where this scheme was used in constructing solutions with the appearance and propagation of shock wave fronts.  
For computations, the two-node linear truss elements without bending and torsional stiffness were used [42]. The overall element number varied in a range \(6200\leq N\leq 16\),400. The mesh convergence test revealed almost indistinguishable results for any mesh with the overall element number \(N\geq 12\),800; the latter value was chosen for the main computations. The following median filtration was used to reduce the non-physical oscillations.  
\[f*(t_{i})=\frac{1}{2m+1}\sum_{k=1}^{m}f(t_{i\pm k}), \tag{32}\]  
where \(2m+1\) is the order of the filter, \(f\) is the unfiltered function, and \(f^{\prime}\) is filtered one. The median filters of orders \(5\div 9\) were used in the main computations. It was observed that even a small number of elements (actually the smallest number \(N=6200\) we have used in computations) gave vertical straight lines associated with the appearance of shock wave fronts when faster pulses overtake slower ones. However, at \(N=6200\), the plots for the displacement field contained noticeable saw-tooth jerks, which practically disappeared at \(N\geq 12\),800. Meanwhile, the field variables containing time and spatial derivatives, actually the velocity and strain which were used in computing kinetic and strain energy, needed filtering for all \(N\).  
In computations, the rod length was chosen to be large enough to avoid the arrival of reflected waves coming from the right edge of the rod to the points of observation; similarly, the overall time was chosen such that the reflected waves could not reach the points of observation. And finally, no parallelization algorithms were used for solving the considered nonlinear partial differential equation.