Botteldooren, Dick

Ding, Lei

Van Renterghem, Timothy

English

Ghent University Academic Bibliography

  An efficient approach to evaluate multiple scattering by foliage in a 3D-FDTD model  L. Dinga) T. Van Renterghemb) D. Botteldoorenc) Department of Information Technology, Ghent University Sint-Pietersnieuwstraat 41, Gent 9000, Belgium   In this paper, an approach is presented to study multiple scattering by periodically arranged obstacles in the 3D finite-difference time-domain (FDTD) method. The computational efficiency is improved by applying period boundary conditions; in this way multiple passages through the zone with scattering objects can be modeled. This method is exemplified by modeling sound propagation through foliage. Leaves are modeled as bending plates with damped vibrations. The separation of the incident wave and reflected wave is realized in frequency domain by using the least squares method. The energy attenuation obtained by propagation through this scattering medium can be calculated based on the cross power spectral density of pressure and velocity in the propagation direction. Finally, the method has been used to quantify the effect of different foliage properties on sound propagation through them.  1 INTRODUCTION  Studies on the interaction between sound and vegetation are mainly based on measurements. Often, many effects are observed simultaneously, making it difficult to gain knowledge on specific phenomena. Especially the importance of the interaction between leaves and sound waves remain a question. Scattering is expected to be the main effect, but also leaf vibrations1, 2 and (viscous) damping near the surface of leaves have been measured. Their relative importance is not known. Therefore, elaborating on numerical techniques is interesting. With numerical simulations, different effects can be more easily singled out. In order to numerically simulate the interaction between sound and plant leaves, 3D models are needed. Furthermore, since effects are expected to become important at wavelengths shorter than the dimensions of the leaves, fine numerical                                                  a) email:  lei.ding@intec.ugent.be b) email:  timothy.van.renterghem@intec.ugent.be c) email:  dick.botteldooren@intec.ugent.be  discretizations are mandatory. As a result, this leads to a high computational cost. Especially the memory requirement is a major bottleneck with the FDTD technique that will be used in this study. The use of periodic/cyclic boundaries is a possibility to limit the computational cost. Only a small vegetation volume is considered, where sound passes through a number of times. In this explicitly modeled volume, only a few leaves are placed, representing a realistic leaf area density. The total propagation time will then define the width of the vegetation zone. A drawback of this approach is that randomness in orientation and spacing of leaves cannot be considered. Wave reflection and scattering, and also energy dissipation caused by viscosity and structural damping, are captured by the FDTD model presented in this paper. In order to quantify transmission, the incident wave and reflected wave must be separated, which will be done in frequency domain. Finally, the transmitted amount of energy can be calculated, which could lead to improved engineering models.  2 NUMERICAL MODEL  2.1 Governing equations    In between the scattering objects (leafs) lossless isentropic sound propagation is assumed leading to the linearized continuity equation and momentum equation: 020 Vctp ,                                                                (1) ptV0  ,                                                                      (2) where p  is the pressure; 0  is the mass density of the medium; c  is the speed of sound in the medium; and V is the particle velocity vector. Close to the objects viscosity (and thermal conductivity) cannot be ignored since viscous energy decay in the boundary layer at the surface of the leaves is one of the mechanisms causing sound attenuation1. As in Ref.3 a time-domain approximation for a viscous boundary layer near an infinitely extended flat surface will be used. The viscosity adds an additional term (in frequency domain) to the linearized momentum equations in the directions which are parallel to the leaf surface plane:  Vdjkn , with 02 jkn ,                                                                (3) where velocity is averaged over a layer of thickness d ;  is the viscosity of air;  is the angular frequency and j is imaginary unit. The vibration of a leaf is modeled as a vibrating thin plate4. The viscoelastic damping accompanying the leaf vibration can be included by employing the generalized Maxwell model, which has been used by Antoine et al.5, 6. The leaf is approximated by a homogeneous plate having the shape of the leaf since taking into account the structural fine structure of the leave is beyond reach of the discretization. The bending wave can propagate in the two in-plane directions. Assuming that the plate is orthogonal to the x-direction, the velocity equation can be written as 2222202202zyzyhxpvRtvh xpLxpm ,                                      (4) where m  is the surface density of the plate material; h  is the thickness of the plate; xpv  is the plate velocity vector component in the x-direction; LR  denotes the viscous damping associated to the bending process; 0  denotes the bending and twisting moments per unit thickness;  denotes the viscoelastic damping during the bending of leaf and it can be described by the generalized Maxwell model (see Eqn. (6) ~ (7)). Assuming that the leaf is isotropic, the bending and twisting moments in Eqn. (4) can be formulated as 22220zwDywD ,                                                               (5) where D  is the bending stiffness; w  is the displacement component in x  direction and its time derivative is the velocity xpv . The damping term  in Eqn. (4) can be formulated as 1nnnDRt,                                                               (6) tntsn dsexwxwxwxwn0 22222222,                                   (7) where nR  and ns  are the viscoelastic damping parameters, which can be determined by data-fitting with measurement results.  2.2 Finite-difference Time-domain Method  In this paper, the finite-difference time-domain method (FDTD) is used to solve the set of equations presented in the previous section. The staggered grid organization, both in space and time, as suggested in Ref.7, is considered. A leap-frog scheme is used to update acoustic pressure and velocity components over time. For this specific scheme, the following notations are commonly used to represent the discrete pressures and velocity components in unbounded air ldtkdzjdyidxp ,, ,dtlkdzjdydxixv5.0,,5.0, dtlkdzdyjidxyv5.0,5.0,, dtldzkjdyidxzv5.05.0,,,                           (8) where dx , dy , and dz  are the spatial discretization steps in three directions and dt  denotes the time discretization step. The acoustic pressure is always updated at times ldt  and the velocity components at times dtl 5.0 .  The updating equation for velocity parallel to boundaries, v , is adapted to include the effect of the viscous boundary layer. The square root of  dependence in Eqn. (3) is hereby approximated by a ratio of polynomials of order m and n in frequency domain. Eventually, this leads to the adapted FDTD update equation: 112110112110021000022niiliimkklkkmkklklubbduaadtpaubdadt,   (9) where d  is the grid step in the direction orthogonal to the leaf plane;  denotes the directions parallel to the leaf surface; and  is the dynamic viscosity. For the simulations in this paper, m  and n are chosen equal to 2. The values for ka and ib  are the same as those used by Bockstael et al.8: 10a , 871.11a , 87213.02a , 02.3910b , 2.7691b , 2.3782b . Note that for the special case 10a , and other coefficients equal to 0, Eqn. (9) reduces to the discretized form of Eqn. (2). On the bending leaf, plate velocities follow a similar numerical discretization scheme as the particle velocities in air. The bending and twisting moments 0  are discretized in a grid collocated with the velocities in space and collocated with the pressures in time.  2.3 Cyclic Boundaries  Instead of simulating the whole vegetation volume, only a small cubic box is considered and multiple passages through this volume are modeled. For this, the concept of cyclic boundary condition is used: outgoing values on one boundary will be used as the ingoing values in the next time step at the other side of the simulation domain. For example, on the boundaries 1xx  and Nxx , the pressure equation and the velocity equation in x-direction can be written as kjlllNxlxlkjlkjdvvdxvvcdtpp,5.05.05.05.05.05.05.05.0120,,11,,1 ,                           (10) lkjNlkjlkjNxlkjNxppdxdtvv ,,,,105.0,,5.05.0,,5.0.                                       (11) This periodic extension introduces periodicity in the leaf placement which will lead to special effects of periodic structure but these will occur mainly below the frequencies of interest if the basic simulation cell is kept large enough.  Most applications of FDTD assume a spatially localized sound source. Because of the cyclic boundary, such a source would also be periodically extended which is not desired. Hence we opt for a plane wave as a starting field. To minimize numerical dispersion the initial plane wave is chosen to propagate along the diagonal of the cubic box. The initial values for pressures and velocities are Gaussian modulated sine waves both in space and time. The center plane of the Gaussian pulse must be chosen carefully in order to make the pressure and velocities on the edge of this cubic box match the pressure and velocities on the edge of its adjacent cubic boxes. The initial wave fields are therefore centered at three planes as shown in Fig. 1.  2.4 Wave Separation Approach  During the time-stepping, the amplitude of the plane wave propagating in the incident direction will be reduced because of scattering and dissipation. In order to quantify these effects this plane wave has to be separated out of the total field. For this, the velocities and the pressures on the three planes orthogonal to the diagonal propagation shown in Fig. 1 are recorded and the least-square-method9 is used to separate the wave fields. The waves in a 3D-problem can be written as destzyxs tinnnn ˆ,,, ,                                                 (12) where  is the angular frequency; tzyxs nnn ,,,  is the signals in time domain and sˆ  is its corresponding frequency spectrum, which can be locally approximated by  ikwikwikvikvikurefikuinc eCeCeBeBeAeA 2121 ,          (13)  where incA  denote the amplitude of the wave propagating in the original diagonal direction; and refA , 1B , 2B , 1C  and 2C  denote the amplitudes of scattered waves in three directions. One of the three directions, nnn zyxu ,,  is parallel to the diagonal direction and the other two directions, nnn zyxv ,,  and nnn zyxw ,, , are orthogonal to the diagonal direction. These coefficients can be found for each point on the planes orthogonal to the propagation direction shown in Fig. 1 based on the recorded data at m  neighboring measurement points (one of them lying outside the plane) by minimizing the quadratic error  221211ˆnikwikwikvikvmnikurefikuinc seCeCeBeBeAeA .   (14)  In order to find the minimum value, the derivatives to these coefficients are set equal to 0. As a result, a linear system of equations can be constructed which must be solved.  The procedure described above is applied to the pressure field and to the u -component of the velocity field uV –  where in the latter case the B and C terms can be assumed zero. The cross power spectrum density is used to denote the energy or intensity propagating in the original plane direction u , and it is given by incuincI VPS ,*,                                                          (15) where incP  and incuV ,  are the frequency spectrum of the pressure and velocity propagating forwards in the diagonal direction; and the * in the superscript denotes the complex conjugate. Finally, intensity is calculated on all points and averaged over the plane.  3 NUMERICAL CALCULATIONS  The size of the unit cubic box in the simulation has dimension 0.3m×0.3×m0.3m. The cell size is 0.01m×0.01m×0.01m and the time step is chosen to make the Courant number equal to 1. Three types of leaves, namely Prunus Laurocerasus, Tilia and Prunus Serrulata, are considered. The size (width×length) of the Prunus Laurocerasus leaf, Tilia and Prunus Serrulata is 8cm×16cm, 12.1cm×13.1cm and 7.5cm×12cm, respectively. Their corresponding leaf areas are approximately 1.01×10-2m2, 1.36×10-2m2 and 0.69×10-2m2; and the leaf surface density is 271 g/m2, 104 g/m2 and 185 g/m2. For each type of vegetation, a calculation was made with 1, 2 or 3 leaves placed in this cubic box. The simulation with 1 leaf corresponds to a leaf area density (LAD) of 0.374m-1 for Prunus Laurocerasus, 0.504m-1 for Tilia and 0.255 m-1 for Prunus Serrulata.  The bending stiffness of the leaf can be calculated by the following equation4  23112hED ,                                                               (16) where D  is the bending stiffness for a beam; E  is the Young's modulus; h  is the leaf thickness; and  is the Poisson's ratio. Published work by Takami Saito et al.10 suggests that the Young's modulus of the leaves from Quercus glauca and Quercus serrata plants are in the range of 200-800 MPa. Other work suggests that the Poisson's ratio of an isotropic leaf specimen can be taken close to 0.2511. The thickness of the leaf has the order of magnitude of 0.0005m. Based on these values, the bending stiffness can be estimated and in this paper it is assumed to be 0.0025N·m for all leaves. The damping terms during the leaf vibration have the same values as the spruce wood5, 6. In order to present the sound transmission loss through the foliage, the (dimensionless) approach presented by Aylor12 is used. The excess attenuation EAdivided by the square root of the product of LAD and the breadth/width of the vegetation L , is show as a function of ka , where k is the wave number and a  is the typical leaf width. This way of representing was shown to be independent of species when looking at sound transmission through reeds and corn12.  4 RESULTS AND DISCUSSION   Because of the presence of the leaves, scattering and dissipation will lead to transmission loss. Fig. 2~Fig. 4 show the decrease in intensity of the plane wave when it propagates through an area filled with Prunus Laurocerasus, Tilia, or Prunus Serrulata leaves. For the results at 4000Hz, the decrease in intensity is enhanced with propagated distance, when compared to the results in 2000Hz. 2 or 3 leaves attenuate more energy than 1 leaf because of stronger backscattering and higher energy dissipation. Initially, the attenuation is stronger than linear while after some propagation distance the attenuation is less than linear. The latter is due to multiple scattering resulting in part of the waves returning to the initial direction. Fig. 5~Fig. 7 show the attenuation normalized by leaf area density and the length of the propagation path. These curves show a similar behaviour as in the measurements performed by Aylor12. Furthermore, the values obtained have the same order of magnitude. The dimensionless transmission loss values still depend on the number of leaves (or LAD), in constrast to Ref.12. Note however that the LAD is much smaller in the current simulations. The current approach is expected to become invalid above ka=10~15. Besides, Fig. 7 shows that at low frequency the excess attenuation is negative, which can be attributed to the fact that the current wave separation approach is not valid at low frequencies. Further work is needed to validate the current approach, and parameter studies will be performed to reveal the relevant leaf parameters as for their noise reducing potential.  5 ACKNOWLEDGEMENTS  The research leading to these results has received funding from the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement No.234306, collaborative project HOSANNA.  6 REFERENCES  1. M. J. M. Martens and A. Michelson, “Absorption of acoustic energy by plant leaves”, J. Acoust. Soc. of Am., 69(1), 303-306, (1981).  2. S. H. Tang, P. P. Ong and H. S. Woon, “Monte Carlo simulation of sound propagation through leafy foliage using experimentally obtained leaf resonance parameters”, J. Acoust. Soc. of Am., 80(6), 1740-1744, (1986).  3. D. Botteldooren, “Vorticity and entropy boiundary conditions for acoustical finite-difference time-domain simulations”, J. Acoust. Soc. of Am., 102(1), 170-178, (1997).  4. F. Fahy, Foundation of Engineering Acoustics, Academic Press, (2001).  5. Antoine Chaigne and Christophe Lambourg, “Time-domain simulation of damped impacted plate. I. Theory and experiments”, J. Acoust. Soc. of Am., 109(4), 1422-1432, (2001).  6. Christophe Lambourg, Antoine Chaigne and Denis Matignon, “Time-domain simulation of damped impacted plate. II. Numerical model and results”, J. Acoust. Soc. of Am., 109(4), 1433-1447, (2001).  7. D. Botteldooren, “Acoustical finite-difference time-domain simulation in a quasi-Cartesian grid”, J. Acoust. Soc. of Am., 95(5), 2313-2319, (1994).  8. A. Bockstael, B. de Greve, T. Van Renterghem, D. Botteldooren, W. D’Haenens, H. Keppler, L. Maes, B. Philips, F. Swinnen and B. Vinck, “Verifying the attenuation of earplugs in situ: Method validation using artificial head and numerical simulations”, J. Acoust. Soc. of Am., 124(2), 973-981, (2008).  9. E. P. D. Mansard and E. R. Funke, “The measurement of incident and reflected spectra using a least square method”, Proc. 17th Coastal Engrg. Conf. 1, 154-172 (1980).  10. T. Saito, K. Soga, T. Hoson and I. Terashima, “The bulk elastic modulus and the reversible properties of cell walls in developing Quercus leaves”, Plant Cell Physiol., 47(6), 715-725 (2006).  11. K. J. Niklas, Plant Biomechanics: An Engineering Approach to Plant form and Function, The University of Chicago Press, Chicago, (1992).  12. D. Aylor, “Sound transmission through vegetation in relation to leaf area density, leaf width and breadth of canopy”, J. Acoust. Soc. of Am., 51(2), 411-414, (1972).    00.51 00.5100.51 Fig. 1 - Diagram showing the three planes on which the wave fields are centered.  Fig. 2 -  Plane sound wave attenuation in Prunus Laurocerasus. The left figure shows the results at 2000Hz and the right one at 4000Hz.  Fig. 3 -  Plane sound wave attenuation in Tilia. The left figure shows the results at 2000Hz and the right one at 4000Hz.  Fig. 4 -  Plane sound wave attenuation in Prunus Serrulata. The left figure shows the results at 2000Hz and the right one at 4000Hz. 0 2 4 6 8 100123kaExcess Attenuation/(F*L)1/2  1 Leaf2 Leaves3 Leaves Fig. 5 - Normalized Excess attenuation for the leaf of Prunus Laurocerasus.  0 5 10 150123kaExcess Attenuation/(F*L)1/2  1 Leaf2 Leaves3 Leaves Fig. 6 - Normalized Excess attenuation for the leaf of Tilia.  0 2 4 6 8 10-101234kaExcess Attenuation/(F*L)1/2  2 Leaf4 Leaves6 Leaves Fig. 7 - Normalized Excess attenuation for the leaf of Prunus Serrulata.   

An efficient approach to evaluate multiple scattering by foliage in a 3D-FDTD model

https://biblio.ugent.be/publication/3163759/file/3163762

An efficient approach to evaluate multiple scattering by foliage in a 3D-FDTD model

Abstract

Similar works

Full text

Available Versions

Ghent University Academic Bibliography