## 3 Numerical Simulation of Wave Propagation in Glasses
a tetragonal grid or equivalently \(\Delta x=\Delta y\) with \(i=[1,n]\) and \(j=[1,n]\), the finite integration provides the following discretized forms for Equations (4) and (5):  
\[u_{x,ij}^{n}=u_{x,ij}^{n-1}+\frac{\Delta t}{\Delta x}B_{x}\Big{(}\tau_{xx,i+1,j }^{n-1/2}-\tau_{xx,i,j}^{n-1/2}+\tau_{xy,j,i}^{n-1/2}-\tau_{xy,j,i-1}^{n-1/2} \Big{)}+\Delta tB_{x}f_{x,i,j} \tag{6a}\] \[u_{y,ij}^{n}=u_{y,i,j}^{n-1}+\frac{\Delta t}{\Delta x}B_{y}\Big{(}\tau_{xx,i+1,j }^{n-1/2}-\tau_{xx,i,j}^{n-1/2}+\tau_{xy,j,i}^{n-1/2}-\tau_{xy,i,j-1}^{n-1/2} \Big{)}+\Delta tB_{y}f_{x,i,j}\] (6b) \[\tau_{xx,i,j}^{n+1/2}=\tau_{xx,i,j}^{n-1/2} +\frac{\Delta t}{\Delta x}\Big{\{}(\lambda+2\mu)\Big{(}u_{x,i,j}^{ n}-u_{x,i-1,j}^{n}\Big{)}+\lambda\Big{(}u_{y,i,j}^{n}-u_{y,i,j-1}^{n}\Big{)}\] (7a) \[+(\lambda^{\prime}+2\mu^{\prime})\Big{(}\Big{(}\frac{\partial u_{ x}}{\partial t}\Big{)}_{i,j}-\Big{(}\frac{\partial u_{x}}{\partial t}\Big{)}_{i -1,j}\Big{)}\] \[\rho\frac{\partial u_{y}}{\partial t}=\frac{\partial\tau_{xy}}{\partial x}+ \frac{\partial\tau_{yy}}{\partial y}+f_{y}\]0 \[\rho\frac{\partial u_{y}}{\partial t}=\frac{\partial\tau_{xy}}{\partial x}+ \frac{\partial\tau_{yy}}{\partial y}+f_{y}\]1 (7b) \[\rho\frac{\partial u_{y}}{\partial t}=\frac{\partial\tau_{xy}}{\partial x}+ \frac{\partial\tau_{yy}}{\partial y}+f_{y}\]2 \[\rho\frac{\partial u_{y}}{\partial t}=\frac{\partial\tau_{xy}}{\partial x}+ \frac{\partial\tau_{yy}}{\partial y}+f_{y}\]3 \[\rho\frac{\partial u_{y}}{\partial t}=\frac{\partial\tau_{xy}}{\partial x}+ \frac{\partial\tau_{yy}}{\partial y}+f_{y}\]4 (7c) \[\rho\frac{\partial u_{y}}{\partial t}=\frac{\partial\tau_{xy}}{\partial x}+ \frac{\partial\tau_{yy}}{\partial y}+f_{y}\]5 \[\rho\frac{\partial u_{y}}{\partial t}=\frac{\partial\tau_{xy}}{\partial x}+ \frac{\partial\tau_{yy}}{\partial y}+f_{y}\]6  
where \(\lambda^{\prime}\), \(\mu^{\prime}\) are the anelastic parameters; \(g_{i,j}\) represent the stress sources; \(B_{i}\) is the effective buoyancies; and \(\overline{\mu}_{ij}\); \(\overline{\mu}_{ij}^{\prime}\) are the visco-elastic effective rigidity constants, respectively, given by:  
\[\rho\frac{\partial u_{y}}{\partial t}=\frac{\partial\tau_{xy}}{\partial x}+ \frac{\partial\tau_{yy}}{\partial y}+f_{y}\]7  
\[\rho\frac{\partial u_{y}}{\partial t}=\frac{\partial\tau_{xy}}{\partial x}+ \frac{\partial\tau_{yy}}{\partial y}+f_{y}\]8  
\[\rho\frac{\partial u_{y}}{\partial t}=\frac{\partial\tau_{xy}}{\partial x}+ \frac{\partial\tau_{yy}}{\partial y}+f_{y}\]9  
\[\overline{\mu}_{ij}^{\prime}=\frac{4 Despite the fact that several \"defect\" acoustic phenomena including dispersion, scattering, reflection, and mode conversion phenomena, have been taken into consideration in the numerical simulation, there are minor differences between the experimental and the theoretically predicted spectrum. These differences are due to the non-perfect parallelism between the parallel phases of the glass samples. From the experimental point of view, the non-parallelism induces a dispersion of the acoustic beam that cannot be incorporated into the simulation procedure. An additional cause of the relative intensity differences is the shape of the initial pulse. In the numerical simulation, the shape of the pulse was raised cosine, while the experimental waveform is quite different since the probe is not ideal. In all pulse-echo simulations, we used 5 mm absorbing layers on the left and right side of the sample. The piezoelectric elements used to transmit and receive the longitudinal and shear ultrasonic waves had a constant diameter of 9 mm. The frequency of the ultrasonic wave was fixed at 10 MHz for all simulations. The shape of the pulse was raised cosine with initial amplitude set to 1 in arbitrary units and the total simulation time was 25 us for all simulations performed. This simulation duration was more than enough to receive a complete back wall echo train.