Branched flows of flexural waves in non-uniform elastic plates

Flexural elastic waves and sound in solids are of great interest in wide-ranging contexts such as ultrasound in plates, geophysics, ocean engineering, aerospace and automotive structures, and musical acoustics. Despite bending waves being the most important elastic waves for such surface structures, their propagation in the presence of the inevitable non-uniformity is poorly understood. Here we show the branching and focusing behaviour of highly dispersive flexural waves travelling in elastic plates of non-uniform thickness. The thickness profile has isotropically correlated spatial randomness. The correlation length is much larger than the wavelength. The location of wave focusing shows a scaling relationship with randomness, which is consistent with those previously reported in other random media. We show this analytically and numerically. This suggests a universality in the scaling between the location of wave focusing with randomness and the correlation length, regardless of the physics of the waves in question. Despite the relevance of flexural elastic waves to many diverse physical systems, modelling their propagation is non-trivial. Here, branching and focusing of highly dispersive flexural waves in heterogeneous elastic plates is described numerically and analytically.

W aves in random media often exhibit spatial patterns of branched flow 1 . Spatial focusing and branching of waves in random media have been reported for electron waves 2 , microwaves 3 , tsunami waves 4,5 , light 6,7 , and ocean hydroacoustics 8 . The statistics of the spatio-temporal features of propagation, as induced by the randomness of the medium, have been of great current interest in numerous theoretical investigations [9][10][11] . The phenomenon has been studied for nonlinear waves too 12 . While branching flow of waves has been reported in seemingly unrelated fields, the underlying governing equation in most cases is the classical non-dispersive wave equation ∂ tt η ¼ c 2 ∇ 2 η, where ηðx; y; tÞ is the field variable that depends on the position x; y and time t; c is the wave speed, which is independent of the wavelength. By contrast, bending waves in elastic plates are dispersive, as shorter waves travel faster. Such structures are frequently encountered in the physics of geological plates, glaciology of ice sheets, thin films, nanomechanics of 2D materials such as graphene, aerospace skin structures, etc. Waves in thin elastic plates are governed by the biharmonic Kirchhoff-Love equation 13 Here, η is the transverse displacement, α 2 ¼ EH 2 =12ρð1 À ν 2 Þ combines the thickness of the plate H, Young's modulus E, Poisson's ratio ν and the material density ρ. The group velocity c g is proportional to the wave number k, i.e. c g $ αk, as evident from the fourth spatial derivative.
Consider plates of thickness much smaller than the wavelength of interest λ, i.e. H ( λ, so that higher order effects such as shear through the thickness 14 can be ignored. Likewise, rotary inertia [15][16][17] and kinematic non-linearity due to large rotations associated with the so called von Karman strains 18 are neglected. We report the existence of branched flows of flexural waves in such plates when their thickness has modest spatially correlated randomness.
Branched flows lead to regions of localised high amplitude called caustics or focusing events. In this work, we demonstrate the existence of an elegant scaling between the location of the first caustic with a suitably non-dimensionalised measure of the variance in the thickness of the elastic plate, and with the spatial correlation length of the randomness of the thickness. We do so using analytical approaches and two different numerical methods. The scaling we report here is consistent with similar scaling observed in other different physical contexts 1 .

Results and discussion
Consider plates with small random non-uniformity in thickness Hðx; yÞ, described by the variable h (see Fig. 1) that denotes small fluctuations over the mean thickness H 0 , i.e. Hðx; yÞ ¼ H 0 ð1 À hðx; yÞÞ; the mean of h is zero. While the results are developed here for a plate with non uniform thickness, the results can readily be generalised for plates with other non-uniform properties as long as the non-uniformity is small.
A pulse with carrier frequency ω that is modulated by a Gaussian envelope is launched at the left edge (see Methods section). The excitation thus provided resembles the motion of a rug transversely shaken at the left end. A plane wavefront with a predominant wavelength λ % 2π ffiffiffi α p ω À1=2 thus propagates from left to right. Figure 2 shows the spatio-temporal evolution, obtained using finite element elastodynamics simulations. A plane wave front is launched at the left end of a plate of non-uniform thickness. The four panels on the left show the amplitude of transverse displacement at four successive instants of time (top to bottom; x-y axes are normalised, see caption). Regions of highest displacement at each successive time instant are indicated by square boxes labelled A, B, C and D. Widening of the wavefront due to the dispersive nature of the medium is apparent. This can also be seen from the energy density plots (shown in light red above each panel in Fig. 2). We observe the emergence of distinct branches as the wave propagates. Regions of high amplitude at each instant of time shown are zoomed into, and reproduced as 3D plots on the right (with labels A, B, C and D, indicating the correspondence with the four figures on the left).
In a homogeneous dispersive medium of average properties, one would observe a monotonic reduction in the highest amplitudes of the wave front as it propagates (see Fig. S1, supplementary information). This follows from energy conservation since the widening of the wavefront must accompany reduction in amplitude. The maximum amplitude seen at the same time instant in a homogeneous plate is shown in the plots on the right with transparent planes. The levels of these planes decrease monotonically with time, as anticipated. However, the highest amplitudes do not decrease for the non-uniform plate. This is due to the emergence of branches that is attributed to scattering of the plane wave by the non-uniformity. These branches interfere and form locations of amplitudes much higher than what one may expect from the dispersive wave dynamics in an equivalent homogeneous medium. Regions with amplitudes that are lower than those for the homogeneous case also appear, as expected from energy conservation. See Supplementary Movie 1 for an animation comparing flexural wave transmission in uniform and non-uniform plates.
Scaling of the location of the first caustic with randomness. We are now interested in the distance l f of the first high amplitude peak from the left edge, as this quantitatively characterises the emergence of the branching flow of waves in random elastic plates. If the correlation length of randomness field h is L c , then the three length scales in the problem have the relationship H 0 ( λ ( L c . The first part of the inequality enables us to ignore shear effects through the thickness, whereas the latter justifies the use of the eikonal equation 19 , which reduces the plate flexure wave dynamics to one of ray dynamics. Then, using arguments from stochastic dynamics (see Methods section), one obtains where the angular brackets denote mean of the respective quantities; L c is the correlation length of non-uniformity; l f is the distance of the focusing point from the left end. Rays form envelopes known as caustics. The first focusing point in the branched structures of caustics will be referred to as the first caustic. In the ray limit, the only length scale in the problem is the correlation length. Hence, the linear scaling hl f i $ L c follows immediately from dimensional arguments.
Scaling law from numerical integration of ray equations. The scalings hl f i $ L c and hl f i $ hh 2 i À1=3 can be obtained numerically by integrating ray equations using a forward Euler method. As rays propagate, those that were initially parallel to the propagation direction begin to scatter due to inhomogeneities in the thickness of the plate. Some rays eventually cross over causing focusing events as seen in Fig. 3a. These focusing points can be numerically detected as the location where the curve of position vs wavenumber in the y-direction becomes multi-valued 20 . In the instance shown in Fig. 3, the location of the first such focusing point has been indicated with a circular marker. By running multiple instances of these ray simulations for different values of hh 2 i, we obtain the dependence of hl f i on hh 2 i which, is in excellent agreement with Eq. 2 (see, Fig. 3).
The linear scaling hl f i $ L c is also observed in the ray simulations. In Fig. 3b, taking the L c ¼ 0:1 m simulations as a reference, we denote the trend predicted from Eq. 2 using dotted lines. The ray simulations for L c ¼ 0:07 m and L c ¼ 0:15 m are in excellent agreement with the predicted trend. Ray simulations also confirm that hl f i is independent of wavelength (Fig. S2,  supplementary information). The expected location of the first caustic remains largely unchanged even with a 100% change in the wavelength. Nominal thickness H 0 ¼ 1 mm is used for all finite element as well as ray dynamics simulations.
Scaling law from elastodynamic simulations. Let us return to the full wave dynamics using finite element simulations, as the true dispersive character of the flexural waves propagating in elastic plates is best captured there. To detect focusing events, we use the time integrated intensity defined as Iðx; yÞ ¼ R T 0 η 2 ðx; y; tÞdt where T is, approximately, the time when the wavefront reaches the right edge. By definition, the local peaks in Iðx; yÞ correspond to focusing events. In Fig. 4b, a typical instance of Iðx; yÞ is shown. While plotting, the time integrated intensity is averaged along the y axis. This makes the local peaks in intensity visually pronounced. At each point on the x-axis one can calculate the scintillation index, σðxÞ ¼ ðhI 2 i y =hIi 2 y Þ À 1 (see Fig. 4a). The first prominent peak in σðxÞ provides the location of the first caustic. By varying hh 2 i, L c , and λ, we can, as in the case of ray integration, bring out the effects of the variation of the severity of randomness in the thickness, correlation length and wavelength. Figures 4c, d show that the finite element simulations also largely agree with the trend predicted by our analysis. The independence of hl f i with wavelength is also confirmed via our finite element simulations (Fig. S3, supplementary information). There is some deviation from predicted scaling relationships at higher values of hh 2 i, but this is expected, as the weak scattering assumption in the analysis begins to break down at higher value of hh 2 i. In finite element simulations (Fig. S3, supplementary material) we see that even at lower values of hh 2 i, the deviation for the shortest waves considered ( e λ ¼ 0:07) seems to be slightly more pronounced than for longer wavelengths. This may point to the minor role of through-thickness shear which is expected to get more severe as the wavelength decreases. Similar deviation is not seen in ray dynamics simulations (Fig. S2, supplementary material), which is consistent with the formulation of ray equations that ignores through thickness shear.

Conclusions
The emergence of focusing of bending waves in random thin elastic plates was demonstrated and quantitatively characterised. The scaling of the location of the first focusing event with the mean square randomness of the plate is determined when a pulse is launched from one end of plate. The ray description of the problem, arising from the eikonal equation for plate flexural waves, was reduced to a system of equations consistent with previously known scalings in other physical contexts, viz. hl f i $ L c and hl f i $ hh 2 i À1=3 . The ray equations were then integrated numerically and found to agree with these scaling relationships. The location of the first caustic was found to be independent of the wavelength and the correlation function that describes the disorder (Methods section). The scaling relationships between the location of focusing and disorder, the correlation length, were also obtained by elastodynamic simulations using the finite element method.
This work focused on thin plates where thickness shear effects are insignificant. The existence, or absence, of branched flows in thick plates where the wavelength is comparable to the thickness and, therefore, thickness shear is significant, remains an open question.
The formation of channels suggests the potential of tailoring and shaping wave paths within such elastic waveguides, especially when combined with acoustic metamaterials 21 . Since particle separation and manipulation 22 involve steering particles to nodal lines, we anticipate that this work may potentially inform the design of such technologies. While the analysis and computations presented here are restricted to thickness variations, the phenomenon reported here would, in principle, be observed for spatial heterogeneity in material properties such as the modulus of elasticity or density. Results presented here could also have interesting implications to metrology of elastic plates, geophysics and seismology.

Methods
Analytical proof of scaling from ray dynamics. The equation of motion for flexural waves in a thin elastic plate with spatially non-uniform flexural stiffness is given by 19 where, M x ¼ Dð∂ xx η þ ν∂ yy ηÞ, and M y ¼ Dðν∂ xx η þ ∂ yy ηÞ, are bending moments per unit edge length and M xy ¼ ð1 À νÞD∂ xy η is the twisting moment per unit edge length; ρ is the density of the material, Hðx; yÞ is the thickness of the plate, E is the Young's modulus, ν is the Poisson ratio, D ¼ EH 3 12ð1Àν 2 Þ is the bending rigidity of the plate. The above governing equation is valid when D; ρ; E; ν and H are slowly varying functions of x; y. We can write the eikonal approximation of a wave in an inhomogenous thin plates as follows 19 : where, ω is the angular frequency associated with the wavelength of interest, τ ¼ ωt and k ¼ ½k x k y T is the two dimensional wave vector. In order to demonstrate the scaling relationship associated with branching flow, we consider a plate with spatially uniform density, Poisson's ratio and Young's modulus but slowly changing thickness. Thickness can be expressed as Hðx; yÞ ¼ H 0 ð1 À hðx; yÞÞ where H 0 is the nominal thickness and jhðx; yÞj ( 1. Then we can introduce the constant γ 2 ¼ EH 2 0 12ð1Àν 2 Þρ to rewrite D=ρH ¼ γ 2 ð1 À hðx; yÞÞ 2 . Assume that hðx; yÞ is a random field with zero mean and isotropic spatial correlation characterised by correlation length L c ) λ. Then eikonal equations can be rewritten as: Barring some constants, the structure of these equations is identical, for monochromatic waves, to those in 4,5 . Following the arguments presented therein we can arrive at the scaling relationship Finite element simulation of elastodynamics. To numerically explore branching flows in elastic plates with randomly varying properties, we consider a thin rectangular plate of non-uniform thickness and dimensions L x L y modelled using shell finite elements in within ANSYS, a commercial finite element code. The main propagation direction is along x-axis. The length in this direction is L x ð) L y Þ. The thickness of the plate is given by H 0 ð1 À hðx; yÞÞ where hðx; yÞ is a random field with zero mean and isotropic spatial correlation length L c . A MATLAB code was used to generate multiple instances of the random field with desired standard deviation and correlation length. Simulations were carried out within the ANSYS APDL environment using SHELL181 finite elements. SHELL181 finite elements have both membrane and flexural strain energy taken into account. The degrees-of-freedom for each node are 3 displacements and 3 rotations. Rotational degrees-of-freedom are required in plate/shell finite elements, unlike many 3D elasticity formulations because of the intrinsic simplification through the thickness within the theory of plates and shells. During the simulations, the two in-plane degrees-of-freedom were locked and so was the rotation about the axis normal to the plate. Thus the remaining active degrees-offreedom for each finite element node are the transverse displacement and two rotations. The element allows specification of different thickness values at each of the nodes and also interpolates thickness within an element. This is crucial as we have a spatially varying flexural rigidity field Dðx; yÞ. The approach uses variational minimisation of the integral of the Lagrangian of the elastic continuum L ¼ T À U, where T is the kinetic energy and U the potential energy. Local basis functions, known as shape functions are used to describe the elastic field in conjunction with unknowns that are treated as the generalised coordinates of the problem. Setting δLdt ¼ 0 within two arbitrary instances of time, where δ is "the first variation of" 23 leads to a set of coupled ordinary differential equations, the spatial variable having been eliminated in the process. These coupled ordinary differential equations were solved using a Newmark-β scheme 24 . The details are omitted as the implementation is a part of the commercial code used. Although the element is capable of simulating dynamics of surface structures with intrinsic curvature, such as those in elastic shells, we have simulated a structure with a plate-like flat equilibrium shape.
The right edge is fixed and the two long edges running x-wise are free. The left edge is excited transversely with a Gaussian pulse of dominant frequency ω to give rise to a pulse composed predominantly of wavelength λ. Using the dispersion relation of a homogeneous plate whose thickness is the same as nominal thickness of the plate under consideration we can obtain ω ¼ γð2π=λÞ 2 . The temporal profile of the excitation is shown in Fig. 5.
Simulations were carried out on the normalised domain ðe x; e yÞ 2 ½0; 40 ½0; 8 except for the lowest three standard deviation which were done over ðe x; e yÞ 2 ½0; 80 ½0; 8. Note that even for the lowest amount of randomness studied h e l f i < 40, but it was ensured for all simulations the propagation distance was at least twice as much as the expectation value of the first focusing event to avoid biasing the sample. This means that for the highest 7 standard deviations studied, a field of normalised length 40 was sufficient but for that latter a field of normalised length 80 was used. This was done due to the computational load associated with running these simulation. It is apparent that this has not had any significant bias on the sample since there is no discontinuity in the trend of the expectation value of the first focusing events. This can also be seen by looking at the spread of data in the main text of the paper. More details about simulations can be found in Supplementary note 1, Supplementary information.
Energy densities shown in Figs. 2 and S1 are obtained from the finite element elastodynamic simulations by taking temporal and spatial derivatives, numerically, to find kinetic and potential energies at each time instant shown.
Other correlation functions. A Gaussian correlation function was used to construct the thickness fluctuation field hðx; yÞ for all the results reported in the main text. We also used a power law correlation function of the form Aðx; yÞ ¼ ð1 þ ðx 2 þ y 2 Þ=L 2 c Þ À2 . The scaling remains unaffected (Fig. 6). This is consistent with literature on branching phenomenon in other physical contexts that the scaling of hl f i remains unchanged if other correlation functions are used 4 . The only requirement is that the integral of this function over R 2 be well defined. Excitation applied to the left edge of the plate as a function of time in finite element based elastodynamic simulations. The zoomed out inset shows the full applied excitation and it can be seen that the excitation is a short initial pulse. Fig. 6 Scaling of location of caustics for power law correlation function. Finite element elastodynamic simulations confirm that using a power law correlation function of the form Aðx; yÞ ¼ ð1 þ ðx 2 þ y 2 Þ=L 2 c Þ À2 instead of a Gaussian correlation to construct hðx; yÞ does not affect the predicted scaling of the location of the first caustic. Note that this a log-log plot. Markers indicate the mean. A 2-standard-deviation width centred at the mean is indicated by a vertical line.