Branched flows of flexural waves in non-uniform elastic plates

11:31 21 junio in Artículos por Website by Consultora CIS

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¹⁹

$$\rho H{\partial }_{{tt}}\eta +{\partial }_{{xx}}{M}_{x}+{\partial }_{{yy}}{M}_{y}+2{\partial }_{{xy}}{M}_{{xy}}=0,$$

(3)

where, ${M}_{x}=D({\partial }_{{xx}}\eta +\nu {\partial }_{{yy}}\eta )$, and ${M}_{y}=D(\nu {\partial }_{{xx}}\eta +{\partial }_{{yy}}\eta )$, are bending moments per unit edge length and ${M}_{{xy}}=(1-\nu )D{\partial }_{{xy}}\eta$ is the twisting moment per unit edge length; $\rho$ is the density of the material, $H(x,y)$ is the thickness of the plate, $E$ is the Young’s modulus, $\nu$ is the Poisson ratio, $D=\frac{E{H}^{3}}{12(1-{\nu }^{2})}$ is the bending rigidity of the plate. The above governing equation is valid when $D,\,\rho ,\,E,\,\nu$ 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¹⁹:

$${\partial }_{\tau }{{{{{\bf{x}}}}}}={(D/\rho H)}^{1/2}{{{{{\bf{k}}}}}}/\omega ,$$

(4)

$${\partial }_{\tau }{{{{{\bf{k}}}}}}=-\frac{1}{4}(\rho H/D)\nabla (D/\rho H),$$

(5)

where, $\omega$ is the angular frequency associated with the wavelength of interest, $\tau =\omega t$ and ${{{{{\bf{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 ${{{{{\rm{|}}}}}}h(x,y){{{{{\rm{|}}}}}}\ll 1$. Then we can introduce the constant ${\gamma }^{2}=\frac{E{H}_{0}^{2}}{12(1-{\nu }^{2})\rho }$ to rewrite $D/\rho H={\gamma}^{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}_{{{{{{\rm{c}}}}}}}\gg \lambda$. Then eikonal equations can be rewritten as:

$${\partial }_{t}{{{{{\bf{x}}}}}}=\gamma (1-h){{{{{\bf{k}}}}}}$$

(6)

$${\partial }_{t}{{{{{\bf{k}}}}}}=\frac{\omega }{2}\frac{\nabla h}{(1-h)}.$$

(7)

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

$${{\langle }}{l}_{{{{{{\rm{f}}}}}}}{{\rangle}}\propto {L}_{{{{{{\rm{c}}}}}}}{{\langle }}{h}^{2}{{{\rangle}}}^{-1/3}.$$

(8)

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}\times {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}(\gg {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}_{{{{{{\rm{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-of-freedom 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 $\delta {{{{{\rm{\int }}}}}}L{{{{{\rm{d}}}}}}t=0$ within two arbitrary instances of time, where $\delta$ is “the first variation of”²³ 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-$\beta$ scheme²⁴. 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 $\omega$ to give rise to a pulse composed predominantly of wavelength $\lambda$. Using the dispersion relation of a homogeneous plate whose thickness is the same as nominal thickness of the plate under consideration we can obtain $\omega =\gamma {(2\pi /\lambda )}^{2}$. The temporal profile of the excitation is shown in Fig. 5.

**Fig. 5: Displacement forcing applied in elastodynamic simulations.**

Simulations were carried out on the normalised domain $(\widetilde{x},\widetilde{y})\in [0,\,40]\,\times\, [0,\,8]$ except for the lowest three standard deviation which were done over $(\widetilde{x},\widetilde{y})\in [0,\,80]\,\times\, [0,\,8]$. Note that even for the lowest amount of randomness studied ${{\langle }}{\widetilde{l}}_{{{{{{\rm{f}}}}}}}{{\rangle}} \, < \, 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}_{{{{{{\rm{c}}}}}}}^{2})^{-2}$. The scaling remains unaffected (Fig. 6). This is consistent with literature on branching phenomenon in other physical contexts that the scaling of ${{\langle }}{l}_{f}{{\rangle}}$ remains unchanged if other correlation functions are used⁴. The only requirement is that the integral of this function over ${{\mathbb{R}}}^{2}$ be well defined.

**Fig. 6: Scaling of location of caustics for power law correlation function.**

Source link

Blog