A Gaussian Beam Based Recursive Stiffness Matrix Model to Simulate Ultrasonic Array Signals from Multi-Layered Media

Ultrasonic testing using arrays is becoming widely used to test composite structures in the Aerospace industry. In recent years, the Full Matrix Capture (FMC) technique has been implemented to extract the signals for post-processing to form an image. The inherent anisotropy and the layering of the structure pose challenges for the interpretation of this FMC data. To overcome this challenge, modeling techniques are required that take into account the diffraction caused by finite-size transducers and the response of the structure to these bounded beams. Existing models either homogenize the entire structure, use computationally expensive finite difference time domain (FDTD) methods, or do not consider the shape of the bounded beam, which is used to test such structures. This paper proposes a modeling technique based on combining the Multi-Gaussian beam model with the recursive stiffness matrix method to simulate the FMC signals for layered anisotropic media. The paper provides the steps required for the modeling technique, the extraction of the system efficiency factor, and validation of the model with experimentally determined signals for aluminum as an isotropic material such as aluminum and Carbon Fiber Reinforced Plastic (CFRP) laminate as a layered material. The proposed method is computationally inexpensive, shows good agreement with the experimentally determined FMC data, and enables us to understand the effects of various transducer and material parameters on the extracted FMC signals.


Introduction
In recent years, the use of composite materials in aerospace structures has risen significantly [1,2]. These structures mainly consist of multi-layered composite laminates. This has led to an increase in ultrasonic phased array testing of such structures [3]. The testing of such multi-layered structures is complicated due to the multiple reflections at the layer boundaries and the direction-dependent velocity caused due to the inherent anisotropy of such structures [4]. There is a need for computational models to simulate, analyze, and study the interaction of finite beam transducers with such multi-layered composite materials, taking into account the effects of multiple reflections and direction-dependent velocity.
A variety of approaches to simulate the response from multilayer structures have been reported in literature. One of these approaches is applying ray methods to the homogenized multi-layer structure [5]. Though these have the advantage of being applicable to planar and curved geometries, These layers are assumed homogeneous and anisotropic. The anisotropic layers are assumed of infinite extent in the plane (x1-x2) normal to the thickness direction. The laminate is bounded by two semiinfinite bounding layers that are denoted by 0 and n + 1. The plane wave traveling from the upper bounding layer has an incident angle θ with respect to the x3 axis and its wave vector projection on the x1-x2 plane is denoted by φ. The plane wave displacement u in a layer is given by the Equation (1) exp( ( )) u i t  =− kx (1) where i is the imaginary number, k is the wavenumber vector, ω is the angular frequency, and t is the time. Due to the application of Snell's law [20], the wavenumber components in the plane of the interfaces should be equal throughout the laminate, i.e., k1 and k2 remain the same. The wavenumber component k3 can be calculated using the Christoffel equation as given below [15]: where cijkl is the stiffness tensor, ρ is the density of the material, δil is the Kronecker delta, dl is the polarization vector component for different wave modes , and i, j, k, l consist of values 1, 2, 3, corresponding to the three axes x1, x2, x3. Equation (2) can be solved to obtain the values of the wavenumber component k3. k3 will have two solutions for each propagating wave mode. One solution corresponds to the downward going wave in the layer and the other corresponds to the upward traveling wave. The downward traveling wave is denoted by '+' and the upward traveling wave in the layer by '-'. The wave modes are represented by p with values 1, 2, and 3. The quasi-longitudinal wave is represented by p = 3 and the quasi-shear waves are represented by p = 1,2. Hence, the displacement in the layer m is given below [16]: x is the local coordinate of the layer m. The relationship between the stress and displacement in the layer is given below: By substituting Equation (3) into Equation (4), a layer stiffness matrix Sm is defined that relates the stresses and displacements at the top and bottom of the layer. These layers are assumed homogeneous and anisotropic. The anisotropic layers are assumed of infinite extent in the plane (x 1 -x 2 ) normal to the thickness direction. The laminate is bounded by two semi-infinite bounding layers that are denoted by 0 and n + 1. The plane wave traveling from the upper bounding layer has an incident angle θ with respect to the x 3 axis and its wave vector projection on the x 1 -x 2 plane is denoted by ϕ. The plane wave displacement u in a layer is given by the Equation (1) where i is the imaginary number, k is the wavenumber vector, ω is the angular frequency, and t is the time. Due to the application of Snell's law [12], the wavenumber components in the plane of the interfaces should be equal throughout the laminate, i.e., k 1 and k 2 remain the same. The wavenumber component k 3 can be calculated using the Christoffel equation as given below [13]: where c ijkl is the stiffness tensor, ρ is the density of the material, δ il is the Kronecker delta, d l is the polarization vector component for different wave modes, and i, j, k, l consist of values 1, 2, 3, corresponding to the three axes x 1 , x 2 , x 3 . Equation (2) can be solved to obtain the values of the wavenumber component k 3 . k 3 will have two solutions for each propagating wave mode. One solution corresponds to the downward going wave in the layer and the other corresponds to the upward traveling wave. The downward traveling wave is denoted by '+' and the upward traveling wave in the layer by '-'. The wave modes are represented by p with values 1, 2, and 3. The quasi-longitudinal wave is represented by p = 3 and the quasi-shear waves are represented by p = 1, 2. Hence, the displacement in the layer m is given below [14]: where a m,p+/− are the wave amplitudes of the downward and upward traveling waves of mode p in the layer m, and d m,p+/− i is the ith component of the polarization vector of wave mode p in the layer m. The coordinates x m 3 is the local coordinate of the layer m. The relationship between the stress and displacement in the layer is given below: where σ m 1i (0), σ m 1i (h m ) are the stress components at the top and bottom of layer m, respectively, and u m i (0), u m i (h m ) are the displacement components at the top and bottom of layer m, respectively. In order to define the stiffness matrix for the entire structure, continuity of stress and displacement is applied at each interface. The equation relating the stress in the upper semi-infinite bounding layer and the lower semi-infinite bounding layer is given below: where S N is the combined stiffness matrix for the entire structure, σ 0 1i , σ n+1 1i are the stress components in the upper and lower bounding layers, respectively, and u 0 i (0), u n+1 i (h n ) are the displacement components in the upper and lower bounding layers, respectively. The above equation can be solved for the unknown reflection and transmission coefficients. Equation (6) leads to the calculation of 9 unknowns from 6 equations. If the bounding layers are considered to be water, then the above equation is simplified as we know the properties of water, the boundary conditions, and the wave modes supported in water. As water is the upper bounding layer, the incident wave can only be a longitudinal wave, hence we now know the stress and displacement on the top layer caused by the incident wave. Choosing water as the bounding layer reduces the number of unknowns, while also considering the no-slip boundary condition between the upper semi-infinite bounding layer and the first interface, and the lower semi-infinite bounding layer and the last interface, The reflection coefficient can be calculated using the below equation [18]:  (8) where S 33 mn is the (3,3) component of the constitutive matrices of S and where ρ f and V f are the density and velocity of sound in water, respectively. The next section shows the theoretical fundamentals of multi-Gaussian beams.

Modeling of the Transducer Gaussian Beams
The transducer response of phased array at a distance z from the face of the transducer can be modeled as a superposition of multi-Gaussian beams [7] as shown below: where X is the coordinates between the jth transmitting element and the receiving elements, c is the wave velocity, x 3 is the distance traveled along the x 3 axis in Figure 1, and d is the polarization vector.
In the above equations, k is the wavenumber and a 1 and a 2 are the width and length of the rectangular transducer, respectively. A n , A m , B n , and B m are the Wen and Breazeale coefficients [19]. Wen and Breazeale expressed the radiation from a circular transducer as a superposition of Gaussian beams with coefficients obtained by nonlinear optimization. These coefficients were expanded for a rectangular transducer by Ding et al. [20].
Hence, at the face of the transducer where x 3 = 0 the velocity distribution is given below: The velocity distribution in the wavenumber-frequency domain can be calculated as given below:

Angular Spectrum of Plane Waves
The method of angular spectrum of waves was first proposed by Goodman [15]. According to the method, finite beam from a transducer can be decomposed using Fourier decomposition into infinite number of plane waves with different angles of propagation in the spatial frequency domain. The angular spectrum method can be carried out by using the Discrete Fourier Transform to transform from wavenumber domain to the spatial domain. The angular decomposition is given by the below equation: where f (k x 1 , k x 2 ) is the acoustic wave field at the face of the transducer in the wavenumber domain and F(x 1 , x 2 ) is the acoustic wavefield in the spatial domain. k x 1 and k x 2 are the wavenumber components in the plane normal to the plane of propagation of the wave. For 2-D inspection, the wavenumber in the y direction can be considered as 0, and Equation (14) simplifies to the equation given below:

Modeling of Array Signals to Simulate FMC
For modeling the array signals from a Gaussian beam transducer, we combine the multi-Gaussian beam model for the transducer elements with the response of the layered material using the stiffness matrix approach as given below: (16) where β(ω) is the combined system function for a pair of transducer and receiver elements, R(k x 1 , 0) is the reflection coefficient of the entire structure calculated using Equation (8), and x 1 is the distance between jth element and the receiver position.
The system function for an array element can be calculated in the following way as proposed by Schmerr [21]: The backwall echo response F BWE for a transducer element from a known material such as aluminium, etc. is calculated experimentally.

2.
The backwall echo F BWA is then calculated analytically using a simple testing configuration.
It is assumed that the relationship between the experimental and analytical backwall echo is given by Equation (17) F Hence, the combined system response between a pair of elements is given below: The deconvolution process in Equation (18) is carried out by implementing a Weiner Filter to reduce the sensitivity to noise as given below [22]: where * refers to the complex conjugate, and ε is a small noise constant. It is assumed that the elements are linear, time-invariant, and identical in frequency response and directivity as those demonstrated by Huang [7]. Hence, the combined system response for just one pair of elements is required to characterize the other elements.

Simulation and Experimental Results
In this section, simulation and experimental results will be presented. The calculations are carried out using three transducer arrays of different center frequencies, array size, and number of elements.
The experimental FMC signals were acquired using the FI ToolBox from Diagnostic Sonar. The transducers used were phased array transducers supplied by Olympus© (Leiderdorp, Netherlands). The signals were captured using the Diagnostic Toolbox and were imported into MATLAB ® for plotting the data. The simulations were carried out using MATLAB 2017 ® . Figure 2 shows the transducer array configuration used for the experiments and simulations, where 1, 2 donate the array element number and n is the total number of array elements. The specifications of the transducers are shown in Table 1.  The specifications of the transducers are shown in Table 1. For simulation and experimental purposes, we consider an aluminum block 25 mm thick and a CFRP laminate, which is quasi-isotropic and 19 mm thick with (0/45/−45/90) layup. There are 169 layers of UniDirectional CFRP prepreg of 110 µm thickness in the laminate with layer of epoxy resin of thickness 5 µm between them. The properties of aluminum and unidirectional CFRP lamina [23] are given in Table 2. Table 2. Material properties [23].

Properties
Aluminum (GPa) Carbon/Epoxy >65% Fibre-Volume Fraction (GPa) By taking the complex material properties, we take into account the attenuation caused due to viscoelasticity in the CFRP lamina. The simulation consists of evaluating the Equations (8), (13) and (18) and then substituting the results in Equation (16).

Total Reflection Coefficient of the Materials under Inspection
The reflection coefficient was obtained by evaluating Equation (8) for the materials in Table 1 at different frequencies and for normal incidence. The reflection coefficient is given in Figure 3a Figure 3a,b shows the reflection coefficients for aluminium and CFRP, respectively. The plane wave reflection coefficient is frequency-dependent, which can be observed in Figure 3. The resonance is characterized by a reflection coefficient of 1, which is observed in Figure 3. At the resonant frequencies, the transmission coefficient is 0. The resonant frequencies can also be analytically calculated by the below equation:  characterized by a reflection coefficient of 1, which is observed in Figure 3. At the resonant frequencies, the transmission coefficient is 0. The resonant frequencies can also be analytically calculated by the below equation: where n = 1, 2, 3, . . . , c is the velocity of the ultrasonic wave in the material, and d is the total thickness of the material. The resonant frequencies are dependent on the thickness of the material system and velocity of the wave in it. The resonant frequencies calculated using Equation (19) are equal to the resonant frequencies observed by evaluation of Equation (8). Therefore, the reflection coefficients at various frequencies and wavenumbers can be calculated using the stiffness matrix method and substituted in Equation (16).

System Functions of the Transducer Arrays
The system function is calculated by using Equation (18). Figures 4 and 5 show the system functions for elements of center frequency 2.25 MHz and 5 MHz used for inspection of the aluminum block and the CFRP laminate.  Figure 3a,b shows the reflection coefficients for aluminium and CFRP, respectively. The plane wave reflection coefficient is frequency-dependent, which can be observed in Figure 3. The resonance is characterized by a reflection coefficient of 1, which is observed in Figure 3. At the resonant frequencies, the transmission coefficient is 0. The resonant frequencies can also be analytically calculated by the below equation: 2 nc RF d = (19) where n = 1, 2, 3 …, c is the velocity of the ultrasonic wave in the material, and d is the total thickness of the material. The resonant frequencies are dependent on the thickness of the material system and velocity of the wave in it. The resonant frequencies calculated using Equation (19) are equal to the resonant frequencies observed by evaluation of Equation (8). Therefore, the reflection coefficients at various frequencies and wavenumbers can be calculated using the stiffness matrix method and substituted in Equation (16).

System Functions of the Transducer Arrays
The system function is calculated by using Equation (18).   In Figures 4 and 5, it can be observed that the system function peak is at the center frequency of the transducer, and the width of the peak depends on the transducer bandwidth. For CFRP laminate in Figure 4b, it is observed that there is another peak closer to the central peak. This is attributed to the surface of the CFRP laminate under inspection, which can also affect the system function. The system function has to be calculated whenever there is a change in central frequency of the transducer or the material under inspection.

Comparison of Experimental and Simulated FMC Signals
For the purposes of this paper, in Figures 6-10, the 1st element, as shown in Figure 2, is the In Figures 4 and 5, it can be observed that the system function peak is at the center frequency of the transducer, and the width of the peak depends on the transducer bandwidth. For CFRP laminate in Figure 4b, it is observed that there is another peak closer to the central peak. This is attributed to the surface of the CFRP laminate under inspection, which can also affect the system function. The system function has to be calculated whenever there is a change in central frequency of the transducer or the material under inspection.

Comparison of Experimental and Simulated FMC Signals
For the purposes of this paper, in Figures 6-10, the 1st element, as shown in Figure 2, is the transmitting element while the others are receiving elements. Similar figures can be plotted for the other transmitting and receiving elements from the simulated FMC data.

Experimental and Simulated FMC Signals in Aluminum
Figures 6 and 7 present the experimental and simulated FMC signals for inspection of aluminum with Arrays 1 and 2, respectively. In Figures 6 and 7, the backwall echo can be clearly seen between 8 and 10 µs. The simulated signals agree with the experimentally determined signals. The front surface reflected signal can be observed at various elements. The slower quasi-shear waves can be seen between 12 and 16 µs. The second backwall echo is observed at 17 µs. The signals decrease in amplitude as we move away from the firing element due to material attenuation and the effects of diffraction.
In Figure 6a, a small signal is observed just before 8 µs. This signal is seen in the experimental result and is missing in the simulation. The signal is a relatively low amplitude signal that can be attributed to small inconsistencies in the experimental aluminum reference block provided by Olympus© (Leiderdorp, Netherlands). In Figure 6b, a low amplitude signal can be seen in the simulation results before the backwall echo for elements 56 and 61. These are attributed to the noise signals generated while synthesizing the simulated signals from a high sampling rate, which is done so as to correspond with the Nyquist frequency of 50 MHz of the experimental results.  Simulated FMC signals for aluminum with 5 MHz 16 element array. Table 3 presents the comparison of the first element backwall amplitude reduction between experimental and simulated results It can be observed from the table that the difference between the reduced amplitude between the simulated and the experimental results is less than 2 dB, showing good agreement between the experimental and simulated results. The percentage error between the experimental and simulated results for the 2.25 MHz and 5 MHz is calculated to be 1.89% and 3.63%, respectively, showing agreement between the simulation and experimental results. Figures 8 and 9 show the experimental and simulated signals for the CFRP laminate. The amplitude of the backwall echoes is reduced as compared to aluminium due to the increased attenuation of signals. The reduction of the signals is due to attenuation caused by viscoelasticity of the lamina, diffraction effects, and the scattering of the wave from the ply interfaces. The signals presented for the 64 and 128 elements array are up to element 30 as the elements beyond this do not receive the reflected signal owing to the losses as stated above.  Table 3 presents the comparison of the first element backwall amplitude reduction between experimental and simulated results It can be observed from the table that the difference between the reduced amplitude between the simulated and the experimental results is less than 2 dB, showing good agreement between the experimental and simulated results. The percentage error between the experimental and simulated results for the 2.25 MHz and 5 MHz is calculated to be 1.89% and 3.63%, respectively, showing agreement between the simulation and experimental results. Figures 8 and 9 show the experimental and simulated signals for the CFRP laminate. The amplitude of the backwall echoes is reduced as compared to aluminium due to the increased attenuation of signals. The reduction of the signals is due to attenuation caused by viscoelasticity of the lamina, diffraction effects, and the scattering of the wave from the ply interfaces. The signals presented for the 64 and 128 elements array are up to element 30 as the elements beyond this do not receive the reflected signal owing to the losses as stated above.

Experimental and Simulated Signals in CFRP
Ply resonances can also be observed in the signals. It is also seen that as the center frequency of the signal increases, the amplitude of the ply resonances and also the scattering from the interfaces increases as seen in Figures 8 and 9. The slower shear waves cannot be observed due to the increased attenuation of the laminate.
In Figure 8a for elements 1, 3, and 5, signals are observed that occur before the backwall echo. These are missing from the simulated results. This can be attributed to the fact that although the simulation takes the layer attenuation and reflections into consideration, the layers in the manufactured material are not of equal thicknesses and might have pockets of resin and other small defects that influence the signal.  Further in the simulated signals of Figure 8b and Figure 10b, low amplitude signals can be observed at 20 μs. These low amplitude signals arise as the sampling rate is high to avoid aliasing and hence leads to the computation of a large number of frequencies. Due to this, when the inverse Fourier transform is carried out, the time window is longer than 20 μs and contains noise generated due to large number of sampling frequencies. As the output signals are cut at 20 μs for comparison with the experimental results, some of these low amplitude noisy signals are seen.
As a comparison of the effect of the element pitch on the FMC signals, Array 3 is used to inspect the CFRP laminate. In Figure 10, the increased pitch and element width of Array 3 show less resonance from the plies as compared to Figure 9. This shows one of the ways the simulation model can be used to optimize the array parameters for different thickness and material properties . Further in the simulated signals of Figures 8b and 10b, low amplitude signals can be observed at 20 µs. These low amplitude signals arise as the sampling rate is high to avoid aliasing and hence leads to the computation of a large number of frequencies. Due to this, when the inverse Fourier transform is carried out, the time window is longer than 20 µs and contains noise generated due to large number of sampling frequencies. As the output signals are cut at 20 µs for comparison with the experimental results, some of these low amplitude noisy signals are seen.
As a comparison of the effect of the element pitch on the FMC signals, Array 3 is used to inspect the CFRP laminate. In Figure 10, the increased pitch and element width of Array 3 show less resonance from the plies as compared to Figure 9. This shows one of the ways the simulation model can be used to optimize the array parameters for different thickness and material properties.
Fourier transform is carried out, the time window is longer than 20 μs and contains noise generated due to large number of sampling frequencies. As the output signals are cut at 20 μs for comparison with the experimental results, some of these low amplitude noisy signals are seen.
As a comparison of the effect of the element pitch on the FMC signals, Array 3 is used to inspect the CFRP laminate. In Figure 10, the increased pitch and element width of Array 3 show less resonance from the plies as compared to Figure 9. This shows one of the ways the simulation model can be used to optimize the array parameters for different thickness and material properties . It can be observed that the experimental front wall echo consists of noise and hence is not suitable for comparison between the experimental and simulated results. Hence, Table 4 presents the difference of 1st and 3rd element backwall amplitude between experimental and simulated results. It can be observed from the table that the difference between the reduced amplitude between the simulated and the experimental results is less than 1 dB, showing good agreement between the results. The percentage error between the experimental and simulated results for the 2.2 5 MHz and 5 MHz is calculated to be 5.5% and 1.43%, respectively, showing agreement between the simulation and experimental results.

Discussions
As shown in Figure 3, the stiffness matrix method predicts the reflection coefficient and the resonant frequency of the system under inspection. This can help in determining the resonant frequency of the material under inspection and hence lead to a better choice of inspection as at resonant frequencies the wave undergoes total reflection, leading to no penetration of the material. Figures 6-9 show that proposed modeling techniques to simulate the received FMC signals are in good agreement to the experimental results as also shown in Tables 3 and 4 with percentage errors less than 6% between the results. Slight discrepancies are noted in the experimental and simulated results. As discussed, these discrepancies arise due to the different ultrasonic velocities used in the simulation and in the experimental acquisition owing to slightly different material properties of the experimental and simulated samples. The low amplitude noise is also present in the simulated results which is due to the signal being synthesized from a large number of frequencies.
The thickness of the couplant used and its properties also influence the time of flight as in the simulation the laminate is surrounded by water bounding layers of infinite thickness , whereas in the experimental scenario the couplant gel layer has a finite thickness and slightly different wave velocity. As discussed in Section 2, water bounding layers are assumed due to the fact that only longitudinal waves can travel through water, hence only incident longitudinal waves need to be It can be observed that the experimental front wall echo consists of noise and hence is not suitable for comparison between the experimental and simulated results. Hence, Table 4 presents the difference of 1st and 3rd element backwall amplitude between experimental and simulated results. It can be observed from the table that the difference between the reduced amplitude between the simulated and the experimental results is less than 1 dB, showing good agreement between the results. The percentage error between the experimental and simulated results for the 2.25 MHz and 5 MHz is calculated to be 5.5% and 1.43%, respectively, showing agreement between the simulation and experimental results.

Discussions
As shown in Figure 3, the stiffness matrix method predicts the reflection coefficient and the resonant frequency of the system under inspection. This can help in determining the resonant frequency of the material under inspection and hence lead to a better choice of inspection as at resonant frequencies the wave undergoes total reflection, leading to no penetration of the material. Figures 6-9 show that proposed modeling techniques to simulate the received FMC signals are in good agreement to the experimental results as also shown in Tables 3 and 4 with percentage errors less than 6% between the results. Slight discrepancies are noted in the experimental and simulated results. As discussed, these discrepancies arise due to the different ultrasonic velocities used in the simulation and in the experimental acquisition owing to slightly different material properties of the experimental and simulated samples. The low amplitude noise is also present in the simulated results which is due to the signal being synthesized from a large number of frequencies.
The thickness of the couplant used and its properties also influence the time of flight as in the simulation the laminate is surrounded by water bounding layers of infinite thickness, whereas in the experimental scenario the couplant gel layer has a finite thickness and slightly different wave velocity. As discussed in Section 2, water bounding layers are assumed due to the fact that only longitudinal waves can travel through water, hence only incident longitudinal waves need to be considered entering the material under inspection. Due to this, the number of unknowns reduces in Equation (6), enabling us to solve 6 equations for 6 unknowns.
The extensive electrical response of array elements is missing in the simulated results, but the computation of the system response by using the backwall method captures the other saliant frequency response of the array elements. The frequency response is observed to be not precisely Gaussian in shape due to different varying factors such as the electronic components of the setup, the top surface of the material inspected, thickness of couplant, etc.
The trend of the diminishing backwall echos as we move further away from the transmitting element is identical for both the experimental and the simulated results. In Figures 8-10, we can observe the reverberations from the layers in the simulated results, which agree with the experimental results. These reverberations tend to contribute to the noise of the signal and hence the ability of the modeling technique to simulate these for the material under inspection can help in optimizing the array parameters without extensive experimental analysis. The slight discrepancy in the results of the CFRP laminate can be attributed to the slightly varying thickness of layers and resin in the manufactured laminate as in the simulation it is assumed that the layers are of constant thickness and parallel to each other. Figures 9 and 10 also provide a comparison between the array signals received due to different element sizes and pitch, which influence the received FMC signals.

Conclusions
The paper proposes a Gaussian beam and recursive stiffness matrix-based modeling techniques to model the FMC signals from layered CFRP laminates. The simulated signals have good agreement, to within 2 dB and a percentage error of less than 6% with the experimental signals, and are able to simulate the different components of the experimental signal, which include ply resonances, front wall reflection, backwall echos, and also the backwall echos from the slower shear waves. The proposed model takes into account the diffraction effects caused by Gaussian beams and mimics the real-world scenario where the transducer emits Gaussian-shaped beams. It is also shown how the model can be used to optimize various parameters of the inspection process. The model can be used for both isotropic and anisotropic layered media, where the anisotropic group velocity is taken into account.