Derivation of the equations of motion and boundary conditions of a thin plate via the variational method

Abstracts

Small deflections in both a thin rectangular plate and a thin circular plate are studied via the variational method. In order to apply Hamilton’s principle to this system, the potential energy is expressed in terms of strain and stress tensors. Quantities such as the gradient displacement tensor and the traction vector are reviewed. It is showed the advantage of the variational method as a technique which allows to obtain the equations of motion and the boundary conditions simultaneously.

Keywords
Stress; strain; thin plate; Hamilton Principle

Pequenas deflexões em uma placa retangular fina e uma placa circular fina são estudadas mediante o método variacional. Para aplicar o princípio de Hamilton neste sistema, a energia potencial é expressa em termos dos tensores de deformação e tensão. Quantidades como o tensor gradiente de deformaçao e o vetor de tração são revisadas. Mostra-se a vantagem do método variacional como técnica que permite obter as equações do movimento e as condições de contorno simultaneamente.

Palavras-chave
Tensão; deformação; lâminas; Princípio de Hamilton

1. Introduction

The mechanics of continuous systems is one of the branches of engineering and physics with most applications in the design of structures and tools that are stable under stress and deformation. Moreover, the difficulty of analytically solving the differential equations governing the behaviour of these continuous systems is well known in the literature. An example of these systems is a thin plate subjected to vibrations, this kind of system was studied throughout the 19th century.

In the late 18th century Ernst Florens Friedrich Chladni noticed that any glass or metal plate produced a variety of sounds whenever he held and stroked it at different positions. Inspired by the experiments of Lichtenberg, who had made the traces of electric discharges visible in insulators by sprinkling dust on the corresponding places, Chladni spread sand on a brass plate, stroked it with the bow of a violin, and the sand formed a star-shaped pattern with ten rays. In 1787 Chladni published his experiments in his first acoustic work [¹]. Chladni’s experiment attracted a great deal of attention in his time. Even Napoleon, through the mediation of Laplace, offered a reward of 3000 francs to anyone who could give an explanation of the phenomenon [²]. This award was given to Sophie Germain, who between 1811 and 1815 formulated the first mathematical model for the deformation of an elastic plate [³]. Although Sophie Germain’s work was an essential breakthrough, her explanation was incomplete until 1850, when Kirchhoff showed that the Chladni figures for a square plate correspond to eigenvalues of the biharmonic operator [⁴]. At the beginning of the 20th century, the expert in sound theory, John William Strutt, later Baron Rayleigh, summarised the situation in his treatise [⁵]. In the same document it is shown that the oscillating thin plate problem has an analytical expression involving the vibrations frequencies, only in the case of a circular plate. Meaning that vibration frequencies can be only computed by numerical means¹, i.e. the eigenfunctions are not computed. For more historical details, the references [³, ⁶] are suggested.

Current work on Chladni figures focuses on finding more accurate and efficient numerical solutions. There are for example solutions via the Q-R and Lanczos algorithms [³], the use of the Ritz-Galerkin method [⁷] and application of the finite quadrature method [⁸, ⁹]. Studies on different geometries and the study of their boundary conditions are also usual [⁸, ⁹, ¹⁰, ¹¹].

The aim of this paper is to present a modern version of the derivation of the equations of motion and boundary conditions of free vibrations of a continuous system via the variational principle. In that perspective, the present paper is organized as follows: In section ², a summary of the properties and interpretations of the stress and unitary deformations in Cartesian and polar coordinates, as well as the derivation of the energy associated with the system, is given. In section ³, the derivation of the equations of motion and boundary conditions is carried out. This technique is applied for rectangular and circular plates.

2. Continuum Dynamics

In this section, the deformations suffered by an elastic body are analysed. Such deformations, characterized by a mathematical entity known as tensor strain, are a consequence of external forces which will be addressed with the so-called tensor stress.

2.1. The strain matrix in Cartesian coordinates

An elastic solid is said to be deformed or strained when the relative displacements between points in the body are changed. This is in contrast to rigid-body motion where the distance between points remains the same. In order to quantify deformation, consider the general example shown in Figure 1 where a continuum region $V_0$ and a generic point $P_0(\mathbf{x})$, in such region are observed. After deformation, the new configuration of the body is denoted by V and the position of the generic point is denoted by $P(\mathbf{X})$.

The displaced position of $P_0$ can be related to displacement vector $\mathbf{u}$ by the relationship

(1) $$\mathbf{X}=\mathbf{x}+\mathbf{u}.$$

Using a Cartesian coordinate system, vectors $\mathbf{x},\mathbf{X}$ and $\mathbf{u}$ are expressed as

(2) $$\begin{array}{cccccccc}\hfill \mathbf{x}\hfill & \hfill =\hfill & \hfill [x\hfill & \hfill y\hfill & \hfill z]{}^T,\hfill & & & \\ \hfill \mathbf{X}\hfill & \hfill =\hfill & \hfill [X\hfill & \hfill Y\hfill & \hfill Z]{}^T,\hfill & & & \\ \hfill \mathbf{u}\hfill & \hfill =\hfill & \hfill [u\hfill & \hfill v\hfill & \hfill w]{}^T,\hfill & & & \end{array}$$

where ${[]}^T$ refers to the transpose.

Figure 1:
Deformation of the continuum region.

It is possible to express the equation (1) in terms of a matrix representing the transformation that the set of points that form the body undergoes when it is deformed, i.e. the Jacobian will be found. To achieve this, the vector $\mathbf{u}$ should be expressed as a Taylor expansion. As an example will be shown the Taylor expansion for the component u around zero, it is

(3) $$u=\frac{\partial u}{\partial x}x+\frac{\partial u}{\partial y}y+\frac{\partial u}{\partial z}z+\mathrm{\dots }$$

Since this work is concerned with small deformations, which means that linear elasticity is being developed, terms higher than the first order can be neglected. Then, from (1) and (3), the component X of the vector $\mathbf{X}$ is

(4) $$X=\left(1+\frac{\partial u}{\partial x}\right)x+\frac{\partial u}{\partial y}y+\frac{\partial u}{\partial z}z$$

Similarly, the Y and Z components of the vector $\mathbf{X}$ are obtained:

(5) $$Y={\displaystyle \frac{\partial v}{\partial x}}x+\left(1+{\displaystyle \frac{\partial v}{\partial y}}\right)y+{\displaystyle \frac{\partial v}{\partial z}}z$$

(6) $$Z={\displaystyle \frac{\partial w}{\partial x}}x+{\displaystyle \frac{\partial w}{\partial y}}y+\left(1+{\displaystyle \frac{\partial w}{\partial z}}\right)z$$

From equations (4), (5) and (6), the equation (1) in its matrix form is

(7) $$\mathbf{X}=F\mathbf{x}\text{where}F=\left(\begin{array}{ccc}\hfill \left(1+\frac{\partial u}{\partial x}\right)\hfill & \hfill \frac{\partial u}{\partial y}\hfill & \hfill \frac{\partial u}{\partial z}\hfill \\ \hfill \frac{\partial v}{\partial x}\hfill & \hfill \left(1+\frac{\partial v}{\partial y}\right)\hfill & \hfill \frac{\partial v}{\partial z}\hfill \\ \hfill \frac{\partial w}{\partial x}\hfill & \hfill \frac{\partial w}{\partial y}\hfill & \hfill \left(1+\frac{\partial w}{\partial z}\right)\hfill \end{array}\right).$$

As a consequence of the linearity, the Jacobian $|F|$ must be invertible and therefore different from zero. Furthermore, to be physically admissible it must also be positive [¹²].

The Jacobian, also known as deformation gradients matrix, can be written as $F=I+\widehat{F},$ where I is the identity matrix of order three and the matrix $\widehat{F}$, called displacements gradients matrix is

(8) $$\widehat{F}=\left(\begin{array}{ccc}\hfill \frac{\partial u}{\partial x}\hfill & \hfill \frac{\partial u}{\partial y}\hfill & \hfill \frac{\partial u}{\partial z}\hfill \\ \hfill \frac{\partial v}{\partial x}\hfill & \hfill \frac{\partial v}{\partial y}\hfill & \hfill \frac{\partial v}{\partial z}\hfill \\ \hfill \frac{\partial w}{\partial x}\hfill & \hfill \frac{\partial w}{\partial y}\hfill & \hfill \frac{\partial w}{\partial z}\hfill \end{array}\right),$$

Each element of the matrix (8) comes from the Taylor expansion of the displacement vector $\mathbf{u}$, as shown in (3), therefore this matrix describes the spatial change of the displacement field. In general, such spatial changes are the product of deformations and rotations of the analyzed element. Thus, by representing the matrix $\widehat{F}$ as a sum of an antisymmetric matrix and a symmetric matrix, that is $\widehat{F}=\omega +ϵ$, where

(9) $$\begin{array}{cc}\hfill \omega & =\frac{1}{2}\left[\widehat{F}-{\widehat{F}}^T\right]=\left(\begin{array}{ccc}\hfill 0\hfill & \hfill {\omega }_{xy}\hfill & \hfill {\omega }_{xz}\hfill \\ \hfill {\omega }_{yx}\hfill & \hfill 0\hfill & \hfill {\omega }_{yz}\hfill \\ \hfill {\omega }_{zx}\hfill & \hfill {\omega }_{zy}\hfill & \hfill 0\hfill \end{array}\right)=-{\omega }^T\hfill \\ \hfill ϵ& =\frac{1}{2}\left[\widehat{F}+{\widehat{F}}^T\right]=\left(\begin{array}{ccc}\hfill {ϵ}_{xx}\hfill & \hfill {ϵ}_{xy}\hfill & \hfill {ϵ}_{xz}\hfill \\ \hfill {ϵ}_{yx}\hfill & \hfill {ϵ}_{yy}\hfill & \hfill {ϵ}_{yz}\hfill \\ \hfill {ϵ}_{zx}\hfill & \hfill {ϵ}_{zy}\hfill & \hfill {ϵ}_{zz}\hfill \end{array}\right)={ϵ}^T.\hfill \end{array},$$

it is observed that the spatial change due to rotations is represented by the antisymmetric matrix w, while the spatial change due to deformations is represented by the symmetric matrix $ϵ$ [¹³].

Since, in this study is considered that the infinitesimal element analysed does not suffer rotation, but only strain, just the symmetric part, called strain matrix, is expressed explicitly as

(10) $$ϵ=\left(\begin{array}{ccc}\hfill \frac{\partial u}{\partial x}\hfill & \hfill \frac{1}{2}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\hfill & \hfill \frac{1}{2}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)\hfill \\ \hfill \frac{1}{2}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\hfill & \hfill \frac{\partial v}{\partial y}\hfill & \hfill \frac{1}{2}\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)\hfill \\ \hfill \frac{1}{2}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)\hfill & \hfill \frac{1}{2}\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)\hfill & \hfill \frac{\partial w}{\partial z}\hfill \end{array}\right).$$

The diagonal elements of the matrix are called normal or extensional strain components and represent the change in length per unit length. While the elements outside the diagonal are the shear strain components and are associated with the change in the angle between two originally orthogonal directions of the infinitesimal element analysed in the continuous material.

2.2. The strain matrix in Polar coordinates

The aim of this section is to obtain the strain matrix in polar coordinates from (10). It is necessary to use the transformation matrix between Cartesian coordinates $(x,y,z)$ and polar coordinates $(r,\theta ,z)$ to achieve such a goal. It is straightforward to demonstrate, for the displacement $\mathbf{u}$ in (1), the following transformation relationship

(11) $$\left(\begin{array}{c}\hfill u\hfill \\ \hfill v\hfill \\ \hfill w\hfill \end{array}\right)=\left(\begin{array}{ccc}\hfill \cos \theta \hfill & \hfill -\sin \theta \hfill & \hfill 0\hfill \\ \hfill \sin \theta \hfill & \hfill \cos \theta \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right)\left(\begin{array}{c}\hfill u_r\hfill \\ \hfill v_{\theta }\hfill \\ \hfill w\hfill \end{array}\right).$$

Thus, the strain matrix in Polar coordinates, denoted by ${ϵ}_P$ is obtained in as shown bellow

(12) $${ϵ}_P=\left(\begin{array}{c}\hfill {ϵ}_{rr}{ϵ}_{r\theta }{ϵ}_{rz}\hfill \\ \hfill {ϵ}_{\theta r}{ϵ}_{\theta \theta }{ϵ}_{\theta z}\hfill \\ \hfill {ϵ}_{zr}{ϵ}_{z\theta }{ϵ}_{zz}\hfill \end{array}\right)$$ $$=\left(\begin{array}{c}\hfill \cos \theta \sin \theta 0\hfill \\ \hfill -\sin \theta \cos \theta 0\hfill \\ \hfill 001\hfill \end{array}\right)\left(\begin{array}{c}\hfill {ϵ}_{xx}{ϵ}_{xy}{ϵ}_{xz}\hfill \\ \hfill {ϵ}_{yx}{ϵ}_{yy}{ϵ}_{yz}\hfill \\ \hfill {ϵ}_{zx}{ϵ}_{zy}{ϵ}_{zz}\hfill \end{array}\right)$$ $$\times \left(\begin{array}{c}\hfill \cos \theta -\sin \theta 0\hfill \\ \hfill \sin \theta \cos \theta 0\hfill \\ \hfill 001\hfill \end{array}\right).$$

To give an idea of how ${ϵ}_P$ is obtained, the term ${ϵ}_{rr}$ will be explicitly calculated. It is easy to obtain the following from equation (12):

(13) $$\begin{array}{cc}\hfill {ϵ}_{rr}& ={\cos }^2\theta \frac{\partial u}{\partial x}+{\sin }^2\theta \frac{\partial v}{\partial y}+\sin \theta \cos \theta \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right).\hfill \end{array}$$

By replacing the partial derivative of (13) with the transformed ones into polar coordinates, the term ${ϵ}_{rr}$ becomes

(14) $$\begin{array}{cc}\hfill {ϵ}_{rr}=& {\cos }^2\theta \left({\cos }^2\theta \frac{\partial u_r}{\partial r}-\cos \theta \sin \theta \frac{\partial v_{\theta }}{\partial r}+\frac{{\sin }^2\theta }{r}{u}_r-\frac{\sin \theta \cos \theta }{r}\frac{\partial u_r}{\partial \theta }+\frac{\sin \theta \cos \theta }{r}{v}_{\theta }+\frac{{\sin }^2\theta }{r}\frac{\partial v_{\theta }}{\partial \theta }\right)\hfill \\ & +{\sin }^2\theta \left({\sin }^2\theta \frac{\partial u_r}{\partial r}+\sin \theta \cos \theta \frac{\partial v_{\theta }}{\partial r}+\frac{{\cos }^2\theta }{r}{u}_r+\frac{\cos \theta \sin \theta }{r}\frac{\partial u_r}{\partial \theta }-\frac{\cos \theta \sin \theta }{r}{v}_{\theta }+\frac{{\cos }^2\theta }{r}\frac{\partial v_{\theta }}{\partial \theta }\right)\hfill \\ & +\sin \theta \cos \theta [2\sin \theta \cos \theta \frac{\partial u_r}{\partial r}+({\cos }^2\theta -{\sin }^2\theta )\frac{\partial v_{\theta }}{\partial r}-\frac{2\sin \theta \cos \theta }{r}{u}_r+\frac{\left({\cos }^2\theta -{\sin }^2\theta \right)}{r}\frac{\partial u_r}{\partial \theta }\hfill \\ & -\frac{\left({\cos }^2\theta -{\sin }^2\theta \right)}{r}{v}_{\theta }-\frac{2\sin \theta \cos \theta }{r}\frac{\partial v_{\theta }}{\partial \theta }].\hfill \end{array}$$

Finally, it is straightforward to reduce (14) into

(15) $${ϵ}_{rr}=\frac{\partial u_r}{\partial r}.$$

By using the same procedure, the other elements of the strain matrix in polar coordinates has the following form

(16) $${ϵ}_P=\left(\begin{array}{ccc}\hfill \frac{\partial u_r}{\partial r}\hfill & \hfill \frac{1}{2}\left(\frac{\partial v_{\theta }}{\partial r}-\frac{v_{\theta }}{r}+\frac{1}{r}\frac{\partial u_r}{\partial \theta }\right)\hfill & \hfill \frac{1}{2}\left(\frac{\partial u_r}{\partial z}+\frac{\partial w}{\partial r}\right)\hfill \\ \hfill \frac{1}{2}\left(\frac{\partial v_{\theta }}{\partial r}-\frac{v_{\theta }}{r}+\frac{1}{r}\frac{\partial u_r}{\partial \theta }\right)\hfill & \hfill \frac{u_r}{\partial r}+\frac{1}{r}\frac{\partial v_{\theta }}{\partial \theta }\hfill & \hfill \frac{1}{2}\left(\frac{\partial v_{\theta }}{\partial z}+\frac{1}{r}\frac{\partial w}{\partial \theta }\right)\hfill \\ \hfill \frac{1}{2}\left(\frac{\partial u_r}{\partial z}+\frac{\partial w}{\partial r}\right)\hfill & \hfill \frac{1}{2}\left(\frac{\partial v_{\theta }}{\partial z}+\frac{1}{r}\frac{\partial w}{\partial \theta }\right)\hfill & \hfill \frac{\partial w}{\partial z}\hfill \end{array}\right).$$

2.3. The stress matrix in Cartesian coordinates

The application of external forces to a region of a continuous and deformable-body will result in the development of stresses within it. The measurement of the stress was postulated by Augustin Louis Cauchy. Cauchy claimed that internal stresses developed within a deformable medium are similar, in character, to the stresses that can be applied externally to create the internal state of stress [¹⁴].

Figure 2:
Sectioned body subject to external loadings: concentrated (${\mathbf{P}}_1$ and ${\mathbf{P}}_2$ ) and distributed ($\mathbf{P}$).

In order to quantify the nature of the internal distribution of forces within a continuum solid, consider a general body subject to arbitrary (concentrated and distributed) external loadings, as shown in Figure 2. To investigate the internal forces, a section is made through the body as shown in Figure 2. On the section $S^{*}$, which divides the region V into two separate regions, consider a small area $\mathrm{\Delta }A$ with unit outward normal vector $\mathbf{n}$. The resultant surface force acting on $\mathrm{\Delta }A$ is defined by $\mathrm{\Delta }\mathbf{F}$, then the traction vector is defined by

(17) $$\mathbf{T}=\underset{\mathrm{\Delta }A\to 0}{lim}\left(\frac{\mathrm{\Delta }\mathbf{F}}{\mathrm{\Delta }A}\right).$$

In the equation (17), $\mathrm{\Delta }A$ is the current area of the element under consideration. Therefore, $\mathbf{T}$ depends on the point P and the orientation of $\mathrm{\Delta }A$. If the outward unit normal to the surface at P is denoted by $\mathbf{n}$, then Cauchy showed that the three components of $\mathbf{T}$ can be determined from the result

(18) $$T_i={\sigma }_{ik}{n}^k,$$

known as the Cauchy formula, where Einstein’s summation notation is used for the repeated index, and the ${\sigma }_{ik}$ are matrix elements of stress tensor in Cartesian coordinates as follow

(19) $$\sigma =\left(\begin{array}{ccc}\hfill {\sigma }_{xx}\hfill & \hfill {\sigma }_{xy}\hfill & \hfill {\sigma }_{xz}\hfill \\ \hfill {\sigma }_{yx}\hfill & \hfill {\sigma }_{yy}\hfill & \hfill {\sigma }_{yz}\hfill \\ \hfill {\sigma }_{zx}\hfill & \hfill {\sigma }_{zy}\hfill & \hfill {\sigma }_{zz}\hfill \end{array}\right).$$

Now, the symmetric property of $\sigma$ will be demonstrated. Assume that the body is in equilibrium when it is subjected to external tractions $\mathbf{T}$ acting on the surface S of a body which has forces per unit volume denoted by ${\mathbf{f}}_i$. From the Newton´s second law, it follows

(20) $${\oint }_S{T}_i\mathit{d}S+{\int }_V{f}_i\mathit{d}V=0,$$

where the first term of the left hand side of (20) is component i that comes from the definition of vector strain on (17). By making use of the divergence theorem, the surface integral can be transformed into a volume integral with the help of the equation (18). Then, from equation (20) is easy to obtain

(21) $$0={\partial }^k{\sigma }_{ik}+f_i.$$

On the other hand, the moment equilibrium about the origin O is expressed as

(22) $$0={\oint }_S\mathbf{x}\times \mathbf{T}\mathit{d}S+{\int }_V\mathbf{x}\times \mathbf{f}\mathit{d}V,$$

where $\mathbf{x}$ is the position vector of dV or dS. Next, the cross product is expressed in terms of the permutation symbol Levi Civita

(23) $$0={\oint }_S{ϵ}^{ijk}{x}_j{T}_k\mathit{d}S+{\int }_V{ϵ}^{ijk}{x}_j{f}_k\mathit{d}V.$$

Again, using the divergence theorem in the first term of the right part and expressing $T_k$ depending of ${\sigma }_{ik}$ as is shown in the equation (18), the surface integral in (23) becomes to a volume integral. By replacing $f_i$ from (21) into (23), is obtained

(24) $${\int }_V{ϵ}^{ijk}{\sigma }_{jk}\mathit{d}V=0.$$

Since ${ϵ}^{ijk}$ is antisymmetric in the subscripts jk and (24) holds to the whole body, ${\sigma }_{jk}$ must be symmetric, i.e. ${\sigma }_{ij}={\sigma }_{ji}$.

3. Equation of Motions and Boundary Conditions via Variational Methods

The equations of motion of a continuum can be obtained using Newton’s laws which requires a free body diagram of a volume element of the structure. This vector approach provides a direct way to derive the equations of motion. However, it is not always clear what kind of boundary conditions to use. Another way to get the equations of motion is through Hamilton’s principle, known as the energy approach. From this approach, it is taken into account that a dynamical system is characterized by two energy functions, kinetic energy, and potential energy [¹⁵].

3.1. The strain energy density

A body of volume V with surface area S is considered to be in static equilibrium under the action of traction T acting on the surface and body forces $\mathbf{f}$ acting on the volume. The resulting stress state in the body is given by $\sigma$, and the virtual displacements in the vicinity of a generic point are denoted by $\delta \mathbf{u}$.

It will be assumed that the strain energy U is equal to the work done W by the applied tractions $\mathbf{T}$ and body forces $\mathbf{f}$ in transforming the body from an undeformed to a deformed configuration [¹⁴, ¹⁶]. The work for a virtual displacement is given by

(25) $$\underset{V}{\int }\delta \overline{W}\mathit{d}V=\underset{V}{\int }{f}_i\delta u_i\mathit{d}V+{\oint }_S{T}_i\delta u_i\mathit{d}S,$$

where the over-bar denotes quantities per unit volume and the subscripts indicate the components of the vectors.

Taking into account that the displacement is virtual it follows that $\dot{\mathbf{u}}=0$, then from equation (21) it is obtained that $f_i=-{\partial }_k{\sigma }_{ik}$. In the second integral of the equation (25) the traction components $T_i$ can be replaced from (18) and via divergence theorem, the work for virtual displacement is

(26) $$\begin{array}{cc}\hfill \underset{V}{\int }\delta \overline{W}\mathit{d}V& =\underset{V}{\int }\left[{\partial }_k\left({\sigma }_{ik}\delta u_i\right)-{\partial }_k{\sigma }_{ik}\delta u_i\right]\mathit{d}V\hfill \\ & =\underset{V}{\int }{\sigma }_{ik}\delta \left({\partial }_k{u}_i\right)\mathit{d}V.\hfill \end{array}$$

In (26), the expression ${\partial }_k{u}_i$ will be written as the sum of a symmetric and a antisymmetric tensor². Therefore, due ${\sigma }_{ij}$ is symmetric, just the symmetric part of ${\partial }_k{u}_i$ survive. Furthermore, the symmetric part is in fact the elements of strain matrix written on (10). In conclusion, the result (26) is rewritten as

(27) $$\underset{V}{\int }\delta \overline{W}\mathit{d}V=\underset{V}{\int }{\sigma }_{ik}\delta {ϵ}_{ik}\mathit{d}V.$$

The work done on the body corresponds to a change in the potential energy U of the system, so that $\delta \overline{W}=\delta \overline{U}$. Thus, the strain potential energy density is defined as

(28) $$\delta \overline{U}={\sigma }_{ik}\delta {ϵ}_{ik}.$$

The relation above becomes an exact differential by assuming that the potential energy $\overline{U}$ is actually a function of the strain tensor ${ϵ}_{ik}$; in that case

(29) $${\sigma }_{ik}=\frac{\partial \overline{U}}{\partial {ϵ}_{ik}}.$$

The equation (29) is the fundamental relation between stress and strain. Also, since the function $\overline{U}$ depends on the terms ${ϵ}_{ij}$, then the differential of $\overline{U}$ is

(30) $$d\overline{U}=\frac{\partial \overline{U}}{\partial {ϵ}_{11}}d{ϵ}_{11}+\frac{\partial \overline{U}}{\partial {ϵ}_{22}}d{ϵ}_{22}+\frac{\partial \overline{U}}{\partial {ϵ}_{33}}d{ϵ}_{33}+\frac{\partial \overline{U}}{\partial {ϵ}_{12}}d{ϵ}_{12}$$ $$+\frac{\partial \overline{U}}{\partial {ϵ}_{13}}d{ϵ}_{13}+\frac{\partial \overline{U}}{\partial {ϵ}_{23}}d{ϵ}_{23}$$

In addition, the relation between $\sigma$ and $ϵ$ for a linear, isotropic and homogeneous material, in the indicial form is written as follows [¹⁴]:

(31) $${\sigma }_{ik}=\lambda {ϵ}_{jj}{\delta }_{ik}+2\mu {ϵ}_{ik}$$

where ${ϵ}_{jj}={ϵ}_{11}+{ϵ}_{22}+{ϵ}_{33}$ and $\lambda ,\mu$ are known as the Lamé constants³. Replacing (31) in (30), it is possible to obtain the form of the function $\overline{U}$ for a linear, isotropic and homogeneous material

(32) $$\overline{U}=\frac{1}{2}\lambda {ϵ}_{jj}{ϵ}_{ii}+\mu {ϵ}_{ik}{ϵ}_{ik}.$$

It is necessary to specify that the Lamé constants are related to the more common Young’s modulus E and Poisson’s ratio $\nu$ by

(33) $$E={\displaystyle \frac{\mu (3\lambda +2\mu )}{\lambda +\mu }},$$

(34) $$\nu ={\displaystyle \frac{1}{2}}{\displaystyle \frac{\lambda }{\lambda +\mu }}.$$

Now, multiplying (29) by ${ϵ}_{ik}$ and substituting (31), the result is as follows

(35) $${ϵ}_{ik}\frac{\partial \overline{U}}{\partial {ϵ}_{ik}}=\lambda {ϵ}_{jj}{ϵ}_{ii}+2\mu {ϵ}_{ik}{ϵ}_{ik},$$

comparing the right part in (35) with (32) the following relationship is obtained

(36) $${ϵ}_{ik}\frac{\partial \overline{U}}{\partial {ϵ}_{ik}}=2\overline{U}.$$

In conclusion, with the help of the equation (29) and (36), it is straightforward to obtain the density of stress potential energy as

(37) $$\overline{U}=\frac{1}{2}{\sigma }_{ik}{ϵ}_{ik}.$$

3.2. The Hamilton’s principle applied to a body in Cartesian coordinates

Hamilton’s principle states that of all the possible paths along which a particle could travel from its position at instant $t_0$ to its position at instant $t_1$, the real path, denoted by $\mathbf{u}=\mathbf{u}(x,t)$, will be the one for which the integral

(38) $${\int }_{t_0}^{t_1}(K-U+W_e)\mathit{d}t$$

is an extremum [¹⁵, ¹⁶]. Thus, the Hamilton’s principle is

(39) $$\delta {\int }_{t_0}^{t_1}(K-U+W_e)\mathit{d}t=0,$$

where K is the kinetic energy of the system, U is the potential energy, and $W_e$ is the work done by the external forces on the system. The symbol $\delta$ is intended in the sense of calculus of variations. Furthermore, in this approach, it is assumed that the varied path $\mathbf{u}+\delta \mathbf{u}$ differs from the real path $\mathbf{u}$ except at instants $t_0$ and $t_1$. In this way, an admissible variation $\delta \mathbf{u}$ satisfies the condition

(40) $$\delta \mathbf{u}(x,t_0)=\delta \mathbf{u}(x,t_1)=0,\forall x.$$

Now, the variation of K, U and $W_e$ will be computed separately in order to obtain an expression for (39). By considering that the body has a constant density, the kinetic energy density $\overline{K}$ is found to be

(41) $$\overline{K}=\frac{\rho }{2}{\dot{u}}_i{\dot{u}}_i,i=1,2,3;$$

where the subscript i denotes a component of the vector $\overrightarrow{u}$. Therefore, the variation of this quantity is

(42) $$\delta \overline{K}=\rho {\dot{u}}_i\delta {\dot{u}}_i.$$

Then, in order to obtain the kinetic component of the equation (39), $\delta \overline{K}$ must be integrated between instants $[t_0,t_1]$ and in the volume. This means that the following integral must be found

(43) $$\delta {\int }_{t_0}^{t_1}K\mathit{d}t=\delta \underset{V}{\int }{\int }_{t_0}^{t_1}\overline{K}\mathit{d}t\mathit{d}V=\underset{V}{\int }{\int }_{t_0}^{t_1}\rho \dot{u_i}\delta {\dot{u}}_i\mathit{d}t\mathit{d}V.$$

Integration by parts over time is done giving

(44) $$\delta {\int }_{t_0}^{t_1}K\mathit{d}t=-{\int }_V{\int }_{t_0}^{t_1}\rho {\ddot{u}}_i\delta u_i\mathit{d}t\mathit{d}V,$$

where the conditions (40) have been used.

For the potential energy, the variation $\delta U$ was founded in (28), and by using (10) for the factor $\delta {ϵ}_{ij}$, the following integral is obtained

(45) $$\delta {\int }_{t_0}^{t_1}{\int }_V\overline{U}\mathit{d}V\mathit{d}t=\frac{1}{2}{\int }_{t_0}^{t_1}{\int }_V\frac{\partial \overline{U}}{\partial {ϵ}_{ij}}\left(\delta {\partial }_j{u}_i+\delta {\partial }_i{u}_j\right)\mathit{d}V\mathit{d}t,$$

The latter expression will be integrated by parts in the volume via $a=\partial \overline{U}/\partial {ϵ}_{ij}$ and $db=(\delta {\partial }_j{u}_i+$ $\delta {\partial }_i{u}_j)dV$. For b, the Gauss theorem gives $b={\oint }_{\mathrm{\Omega }}\left(\delta u_i{n}_j+\delta u_j{n}_i\right)\mathit{d}\mathrm{\Omega }$, where $\mathrm{\Omega }$ is the body superficial surface. Therefore, the equation (45) becomes

(46) $$\delta {\int }_{t_0}^{t_1}\underset{V}{\int }\overline{U}\mathit{d}V\mathit{d}t$$ $$=\frac{1}{2}{\int }_{t_0}^{t_1}\{{\oint }_{\mathrm{\Omega }}(\frac{\partial \overline{U}}{\partial {ϵ}_{ij}}\delta u_i{n}_j+\frac{\partial \overline{U}}{\partial {ϵ}_{ij}}\delta u_j{n}_i)d\mathrm{\Omega }$$ $$-\underset{V}{\int }[{\partial }_j\left(\frac{\partial \overline{U}}{\partial {ϵ}_{ij}}\right)\delta u_i+{\partial }_i\left(\frac{\partial \overline{U}}{\partial {ϵ}_{ij}}\right)\delta u_j]dV\}dt.\mathit{}$$

In order to reduce (46), on the second terms of the surface and volume integrals, the indices i and j will be swapped by taking advantage the symmetric nature of ${\sigma }_{ij}$. Then, with the help of (29), we obtain

(47) $$\delta {\displaystyle {\int }_{t_0}^{t_1}}{\displaystyle {\int }_V}\overline{U}\mathit{d}V\mathit{d}t$$ $$={\displaystyle {\int }_{t_0}^{t_1}}\left({\displaystyle {\oint }_{\mathrm{\Omega }}}{\sigma }_{ij}\delta u_i{n}_j\mathit{d}\mathrm{\Omega }-{\displaystyle {\int }_V}{\partial }_j{\sigma }_{ij}\delta u_{i}dV\right)\mathit{d}t.$$

The last variation to calculate is the external work. This is done easily by considering the equation (25)

(48) $${\int }_{t_0}^{t_1}\underset{V}{\int }\delta {\overline{W}}_e\mathit{d}V\mathit{d}t={\int }_{t_0}^{t_1}\underset{V}{\int }{f}_i\delta u_i\mathit{d}V+{\oint }_{\mathrm{\Omega }}{t}_i\delta u_i\mathit{d}\mathrm{\Omega }\mathit{d}t.$$

Finally, by replacing (44), (47) and (48), into the Hamilton’s Principle given in the equation (39), the following relation is obtained

(49) $$\{{\displaystyle {\int }_{t_0}^{t_1}}{\displaystyle \underset{V}{\int }}[-\rho {\ddot{u}}_i+f_i+{\partial }_j{\sigma }_{ij}]dVdt$$ $$-{\displaystyle {\int }_{t_0}^{t_1}}{\displaystyle {\oint }_{\mathrm{\Omega }}}[{\sigma }_{ij}{n}_j-t_i]d\mathrm{\Omega }dt\}\delta u{}_i=0,$$

where it is evident that

(50) $$\rho {\ddot{u}}_i=f_i+{\partial }_j{\sigma }_{ij},$$

(51) $${\sigma }_{ij}{n}_j=t_i.$$

The equation (50), that comes from the volume integral is the equation of motion for the elements of the body; and the one which comes from the surface integral, equation (51), are the boundary conditions. Therefore, as mentioned before, the variational method allows to determine simultaneously the necessary equations to analyze the dynamic of a solid body.

4. Small Deflections of Thin Rectangular Plates

Usually a thin plate is defined as a continuous region that is delimited by two surfaces of little or no curvature and whose thickness is considerably less than its lateral dimensions. However, the definition of thin elastic plate can be approached by considering the various modes of energy stored in the plate during its deformation. In general, when a plate-shaped elastic structural member undergoes deformation, the stored energy is composed of flexural strain energy due to change in curvature, shear strain energy due to distortion and extensional strain energy due to stretching in the plane of the plate. In classical thin plate theory, the flexural energy is assumed to be the dominant component [¹⁴].

As an example of the previous section, the dynamic that consists of flexural⁴ vibrations, with no in-plane motion, of a rectangular plate will be determined. Consider a plate that occupies a portion of space V, composed of a rectangular surface S with thickness h, hence $V=[0,L_x]\times [0,L_y]\times [-h/2,h/2]$, with $h\ll L_x$, $h\ll L_y$, and a mid-surface, also known as neutral surface, on z=0 as shown in Figure 3.

Figure 3:
Thin plate with volumen V and its neutral surface in z=0.

Because the deflection is w, the kinetic energy is expressed easily as

(52) $$K=\frac{\rho h}{2}{\int }_0^{L_y}{\int }_0^{L_x}{(\dot{w})}^2\mathit{d}x\mathit{d}y.$$

Therefore, it follows that $\delta K\sim \delta \dot{w}$. This expression must be integrated over time in order to replace into (39). The integration by parts of the latter integral drop out the time derivative of $\dot{w}$ give the following result

(53) $$\delta {\int }_{t_1}^{t_2}K\mathit{d}t=-\rho h{\int }_0^{L_y}{\int }_0^{L_x}{\int }_{t_1}^{t_2}\ddot{w}\delta w\mathit{d}t\mathit{d}x\mathit{d}y.$$

Similar to the kinetic energy, is necessary to express $\overline{U}$ as a function of $w(x,y)$. To achieve this goal it is necessary to take into account the classical plate theory which is based on the Kirchhoff hypothesis assumptions⁵, and as consequence of such theory, results the equations in (54) [¹⁵]:

(54) $${ϵ}_{zz}=0,{ϵ}_{xz}=0,{ϵ}_{yz}=0$$

From (10) and (54) it is possible to obtain

(55) $$\begin{array}{cc}\hfill {ϵ}_{zx}& =0=\frac{1}{2}(u_{,z}+w_{,x}),\hfill \\ \hfill {ϵ}_{zy}& =0=\frac{1}{2}(v_{,z}+w_{,y}),\hfill \end{array}$$

where $u_{,z}=\partial u/\partial z$.

By solving the equations (55), the relations $u=-z{w}_{,x}$ and $v=-z{w}_{,y}$ are obtained. Again, from (10) it is straightforward to obtain ${ϵ}_{xx}$, ${ϵ}_{yy}$, and ${ϵ}_{xy}$

(56) $${ϵ}_{xx}=-z{w}_{,xx};{ϵ}_{yy}=-z{w}_{,yy};{ϵ}_{xy}=-z{w}_{,xy}$$

Other relationships necessary to progress in this section are the stress-strain relationships, which are obtained from the generalized Hooke’s law, such relations are given by

(57) $${\sigma }_{ik}={\displaystyle \frac{E}{1+\nu }}\left({ϵ}_{ik}+{\displaystyle \frac{\nu }{1-2\nu }}{ϵ}_{ll}{\delta }_{ik}\right),$$

(58) $${ϵ}_{ik}={\displaystyle \frac{1}{E}}[(1+\nu ){\sigma }_{ik}-\nu {\sigma }_{ll}{\delta }_{ik}],$$

where ${\sigma }_{ll}={\sigma }_{11}+{\sigma }_{22}+{\sigma }_{33}$, $\nu$ is the Poisson constant and E is the Young modulus [¹⁵].

Now, from ${ϵ}_{zz}$ in (54) and (58) it is true that

(59) $${\sigma }_{zz}=0$$

The strain component ${ϵ}_{zz}$ is derived easily by equating to zero the equation for ${\sigma }_{zz}$ in (57), and by using the strain components (56). Hence

(60) $${ϵ}_{zz}=\frac{\nu }{1-\nu }z(w_{,xx}+w_{,yy}).$$

Also, for the potential energy U, the summation (37) has null terms because of the Kirchhoff Plate Theory. In consequence

(61) $$\overline{U}=\frac{1}{2}({\sigma }_{xx}{ϵ}_{xx}+{\sigma }_{yy}{ϵ}_{yy}+{\sigma }_{xy}{ϵ}_{xy}+{\sigma }_{yx}{ϵ}_{yx}).$$

It is now possible to derive an expression for the strain potential energy density in terms of the displacement $w(x,y)$. First of all, the stress factors ${\sigma }_{ij}$ are replaced via (57). Then, the strain matrix elements will be replaced following the expressions (56) and (60). After that, the density deformation energy is

(62) $$\overline{U}=\frac{E}{1+\nu }{z}^2\left\{\frac{1}{2(1-\nu )}{(\mathrm{\Delta }w)}^2+w_{,xy}^2-w_{,xx}{w}_{,yy}\right\},$$

where $\mathrm{\Delta }$ is the Laplacian operator in two dimensions.

With $\overline{U}$ from (62), the potential energy is obtained by integrating over the volume element dV = dxdydz. By defining the rigidity constant $D\equiv E{h}^3/12(1-{\nu }^2)$, this is

(63) $$U=\frac{D}{2}{\int }_0^{L_y}{\int }_0^{L_x}\{{(\mathrm{\Delta }w)}^2+2(1-\nu )$$ $$\times [{(w_{,xy})}^2-w_{,xx}{w}_{,yy}]\}dxdy.$$

In order to replace into (39), the variation of U must be computed. As a first sight, the variation has the following form

(64) $$\delta U=D[{\displaystyle {\int }_0^{L_y}}{\displaystyle {\int }_0^{L_x}}{w}_{,xx}\delta w_{,xx}dxdy$$ $$+{\displaystyle {\int }_0^{L_y}}{\displaystyle {\int }_0^{L_x}}{w}_{,yy}\delta w_{,yy}\mathit{d}x\mathit{d}y$$ $$+\nu {\displaystyle {\int }_0^{L_y}}{\displaystyle {\int }_0^{L_x}}{w}_{,xx}\delta w_{,yy}\mathit{d}x\mathit{d}y$$ $$+\nu {\displaystyle {\int }_0^{L_y}}{\displaystyle {\int }_0^{L_x}}{w}_{,yy}\delta w_{,xx}\mathit{d}x\mathit{d}y$$ $$+2(1-\nu ){\displaystyle {\int }_0^{L_y}}{\displaystyle {\int }_0^{L_x}}{w}_{,xy}\delta w_{,xy}dxdy].$$

In the variation (64), all integrals are $\sim \delta w_{,xx}$. With the objective to obtain the equation of motion and the boundary condition, we must obtain $\delta U$ proportional to $\delta w$, $\delta w_{,x}$, and $\delta w_{,y}$. The latter is possible via integration by parts obtaining the following expression

(65) $$\delta U=D[{\displaystyle {\int }_0^{L_y}}{(w_{,xx}+\nu w_{,yy})}_0^{L_x}\delta w_{,x}dy$$ $$-{\displaystyle {\int }_0^{L_y}}{\left(w_{,xxx}+(2-\nu )w_{,xyy}\right)}_0^{L_x}\delta w\mathit{d}y$$ $$+{\displaystyle {\int }_0^{L_x}}{\left(w_{,yy}+\nu w_{,xx}\right)}_0^{L_y}\delta w_{,y}\mathit{d}x$$ $$-{\displaystyle {\int }_0^{L_x}}{\left(w_{,yyy}+(2-\nu )w_{,xxy}\right)}_0^{L_y}\delta w\mathit{d}x$$ $$+{\displaystyle {\int }_0^{L_x}}{\displaystyle {\int }_0^{L_y}}{\mathrm{\Delta }}^{2}w\delta w\mathit{d}x\mathit{d}y$$ $$+2(1-\nu )w_{,xy}\delta w{|}_{\mathrm{at}\mathrm{corners}}],$$

where ${\mathrm{\Delta }}^{2}w=w_{,xxxx}+w_{,yyyy}+2w_{,xxyy}$, and the differential operator ${\mathrm{\Delta }}^2$ is knowing as biharmonic operator.

Finally, by substituting (53) and (65) into (39) avoiding external forces, is obtained

(66) $${\displaystyle {\int }_{t_1}^{t_2}}{\displaystyle {\int }_0^{L_y}}{\displaystyle {\int }_0^{L_x}}(D\mathrm{\Delta }\mathrm{\Delta }w-\rho h\ddot{w})\delta w\mathit{d}x\mathit{d}y\mathit{d}t$$ $$+{\displaystyle {\int }_{t_1}^{t_2}}[{\displaystyle {\int }_0^{L_x}}{(w_{,yy}+\nu w_{,xx})}_0^{L_y}\delta w_{,y}dx$$ $$-{\displaystyle {\int }_0^{L_x}}{\left(w_{,yyy}+(2-\nu )w_{,xxy}\right)}_0^{L_y}\delta w\mathit{d}x$$ $$+{\displaystyle {\int }_0^{L_y}}{\left(w_{,xx}+\nu w_{,yy}\right)}_0^{L_x}\delta w_{,x}\mathit{d}y$$ $$-{\displaystyle {\int }_0^{L_y}}{\left(w_{,xxx}+(2-\nu )w_{,xyy}\right)}_0^{L_x}\delta w\mathit{d}y$$ $$+2(1-\nu )w_{,xy}\delta w{|}_{\mathrm{at}\mathrm{corners}}]=0.$$

Here, from Hamilton’s principle is followed that the integral is split between a surface integral and boundary integrals; in other words, part of the potential energy contributes to the inertial forces that appear in the equation of motion, and the other part has to be compensated by the boundary conditions, so that disappears along the contour [¹⁶]. The results are summarised as follows:

• Equation of motion

  (67) $$\rho h\ddot{w}=-D\mathrm{\Delta }\mathrm{\Delta }w$$

• Boundary conditions

  (68) $$w_{,nnn}+(2-\nu )w_{,ntt}=0$$
  (69) $$w_{,nn}+\nu w_{,tt}=0$$
  (70) $$w_{,nt}=0,$$

where $\{n,t\}=\{x,y\}$. The equations (68) and (69) are valid for the sides of the rectangular plate and the last one (68) is valid for the corners.

5. Small Deflections of Thin Circular Plates

As a second example the technique will be used to develop de equation of motion and the boundary condition for a circular plate. The first step is determine the kinetic energy and its variation in polar coordinates. This is straightforward by coordinates transformation of (53). After an integral by parts, the following expression is obtained

(71) $$\delta {\int }_{t_1}^{t_2}K\mathit{d}t=-\rho h{\int }_0^{2\pi }{\int }_0^R{\int }_{t_1}^{t_2}\ddot{w}\delta w\mathit{d}tr\mathit{d}r\mathit{d}\theta ,$$

where $w=w(r,\theta )$ and the vector deformation is denoted as $u=(u_r,u_{\theta },w)$.

For the computation of the variation of the potential energy, the expression (37) will be generalized for every geometry as

(72) $$\overline{U}=\frac{1}{2}\mathrm{Tr}(\sigma .ϵ),$$

where the stress and strain matrices could be denoted in polar coordinates as follow

(73) $$\sigma .ϵ=\left(\begin{array}{ccc}\hfill {\sigma }_{rr}\hfill & \hfill {\sigma }_{r\theta }\hfill & \hfill {\sigma }_{rz}\hfill \\ \hfill {\sigma }_{\theta r}\hfill & \hfill {\sigma }_{\theta \theta }\hfill & \hfill {\sigma }_{\theta z}\hfill \\ \hfill {\sigma }_{zr}\hfill & \hfill {\sigma }_{z\theta }\hfill & \hfill {\sigma }_{zz}\hfill \end{array}\right)\left(\begin{array}{ccc}\hfill {ϵ}_{rr}\hfill & \hfill {ϵ}_{r\theta }\hfill & \hfill {ϵ}_{rz}\hfill \\ \hfill {ϵ}_{\theta r}\hfill & \hfill {ϵ}_{\theta \theta }\hfill & \hfill {ϵ}_{\theta z}\hfill \\ \hfill {ϵ}_{zr}\hfill & \hfill {ϵ}_{z\theta }\hfill & \hfill {ϵ}_{zz}\hfill \end{array}\right).$$

Similar to the rectangular case, the Kirchhoff Plate Theory and the relations stress-strain, (57) and (58), implies the following

(74) $${\sigma }_{rz}={\sigma }_{\theta z}={\sigma }_{zz}=0.$$

Therefore, by replacing (74) into (72), the density deformation energy is obtained in terms of the components of the strain and stress matrices

(75) $$\overline{U}=\frac{1}{2}({\sigma }_{rr}{ϵ}_{rr}+2{\sigma }_{r\theta }{ϵ}_{\theta r}+{\sigma }_{\theta \theta }{ϵ}_{\theta \theta }).$$

The relations stress-strain given in equations (57) and (58), will be used in this section again. Thus, the terms ${ϵ}_{zr},{ϵ}_{z\theta }$, will be obtained from (58)

(76) $${ϵ}_{zr}={\displaystyle \frac{1}{E}}[(1+\nu ){\sigma }_{zr}]=0,$$

(77) $${ϵ}_{z\theta }={\displaystyle \frac{1}{E}}[(1+\nu ){\sigma }_{z\theta }]=0.$$

Now, the next step is to obtain the elements of the stress matrix in terms of the vertical deformation $w(r,\theta )$. The term $u_r$ will be obtained from the general form of ${\sigma }_{zr}$ in (16) by taking into account the thin condition (76). Similarly, the deformation component $u_{\theta }$ will be obtained by replacing the thin condition (77) into the general form of ${\sigma }_{z\theta }$. Therefore:

(78) $$u_r=-z{\displaystyle \frac{\partial w}{\partial r}},$$

(79) $$u_{\theta }=-{\displaystyle \frac{w}{r}}{\displaystyle \frac{\partial w}{\partial \theta }}.$$

With (78) and (79) it is possible to obtain all elements of strain matrix with the exception of ${ϵ}_{zz}$ which can be naively calculated null. In order to obtain it, the thin condition ${\sigma }_{zz}=0$ must be replaced into (57) which is coordinate invariant too. Consequently

(80) $${ϵ}_{zz}=-\frac{\nu }{1-\nu }\left(-z\frac{{\partial }^{2}w}{\partial r^2}-\frac{z}{r^2}\frac{{\partial }^{2}w}{\partial {\theta }^2}-\frac{z}{r}\frac{\partial w}{\partial r}\right).$$

With the equations (78)–(80), it is possible to compute all elements that we need to replace into (75) with the help of the stress strain relation. Then, integrating in the volume of the plate, the potential energy has the following form

(81) $$U={\displaystyle \frac{D}{2}}{\displaystyle {\int }_0^{2\pi }}{\displaystyle {\int }_0^R}\{{(\mathrm{\Delta }w)}^2+2(1-\nu )[{({\displaystyle \frac{1}{r^2}}{w}_{,\theta }-{\displaystyle \frac{1}{r}}{w}_{,r\theta })}^2$$ $$-w_{,rr}({\displaystyle \frac{1}{r^2}}{w}_{,\theta \theta }+{\displaystyle \frac{1}{r}}{w}_{,r})]\}rdrd\theta ,$$

where $\mathrm{\Delta }w$ is given in polar coordinate as

(82) $$\mathrm{\Delta }w=w_{,rr}+\frac{1}{r}{w}_{,r}+\frac{1}{r^2}{w}_{,\theta \theta }.$$

The variation $\delta U$ could be determined from (81), but it must be reduced in order to obtain it proportional to $\delta w$, $\delta w_{,r}$, and $\delta w_{,\theta }$. The latter could be made by many integration by parts, the final result has the following form

(83) $$\begin{array}{cc}\hfill \delta U=& {2(1-\nu ){\left(-\frac{1}{r^2}{w}_{,\theta }\delta w+\frac{1}{r}{w}_{,r\theta }\delta w\right)}_0^R|}_0^{2\pi }\hfill \\ & +{\int }_0^{2\pi }{\int }_0^R\underset{{\mathrm{\Delta }}^{2}w}{\underbrace{\left(w_{,rrrr}+\frac{2}{r}{w}_{,rrr}+\frac{1}{r^3}{w}_{,r}-\frac{1}{r^2}{w}_{,rr}+\frac{1}{r^4}{w}_{,\theta \theta \theta \theta }+\frac{2}{r^2}{w}_{,rr\theta \theta }+\frac{4}{r^4}{w}_{,\theta \theta }-\frac{2}{r^3}{w}_{,\theta \theta r}\right)}}\delta wr\mathit{d}r\mathit{d}\theta \hfill \\ & +{\int }_0^{2\pi }\left[{\left(r{w}_{,rr}+\nu w_{,r}+\frac{\nu }{r}{w}_{,\theta \theta }\right)}_0^R\delta w_{,r}+{\left(-r{w}_{,rrr}-w_{,rr}+\frac{3-\nu }{r^2}{w}_{,\theta \theta }+\frac{\nu -2}{r}{w}_{,r\theta \theta }+\frac{1}{r}{w}_{,r}\right)}_0^R\delta w\right]\mathit{d}\theta \hfill \\ & +{\int }_0^R\left[{\left(\frac{1}{r^3}{w}_{,\theta \theta }+\frac{\nu }{r}{w}_{,rr}+\frac{1}{r^2}{w}_{,r}\right)}_0^{2\pi }\delta w_{,\theta }+{\left(-\frac{1}{r^3}{w}_{,\theta \theta \theta }+\frac{2(1-\nu )}{r}{w}_{,rr\theta }+\frac{1-2\nu }{r^2}{w}_{,r\theta }+\frac{2\nu -2}{r^3}{w}_{,\theta }\right)}_0^{2\pi }\delta w\right]\mathit{d}r.\hfill \end{array}$$

Finally, the two variations, (71) and (83), must be replaced into (39) giving the following expression

(84) $$\begin{array}{cc}& {\int }_{t_1}^{t_2}\{{\int }_0^{2\pi }{\int }_0^R(D{\mathrm{\Delta }}^{2}w+\rho h\ddot{w})\delta wrdrd\theta \hfill \\ & +D{\int }_0^{2\pi }\left[{\left(r{w}_{,rr}+\nu w_{,r}+\frac{\nu }{r}{w}_{,\theta \theta }\right)}_0^R\delta w_{,r}+{\left(-r{w}_{,rrr}-w_{,rr}+\frac{3-\nu }{r^2}{w}_{,\theta \theta }+\frac{\nu -2}{r}{w}_{,r\theta \theta }+\frac{1}{r}{w}_{,r}\right)}_0^R\delta w\right]\mathit{d}\theta \hfill \\ & +D{\int }_0^R\left[{\left(\frac{1}{r^3}{w}_{,\theta \theta }+\frac{1}{r^2}{w}_{,r}\right)}_0^{2\pi }\delta w_{,\theta }+{\left(-\frac{1}{r^3}{w}_{,\theta \theta \theta }+\frac{2(1-\nu )}{r}{w}_{,rr\theta }+\frac{1-2\nu }{r^2}{w}_{,r\theta }+\frac{2\nu -2}{r^3}{w}_{,\theta }\right)}_0^{2\pi }\delta w\right]\mathit{d}r\hfill \\ & 2D(1-\nu ){(-\frac{1}{r^2}{w}_{,\theta }\delta w+\frac{1}{r}{w}_{,r\theta }\delta w)}_0^{R,2\pi }\}=0.\hfill \end{array}$$

The results are summarised as follows:

• Equation of motion

  (85) $$\rho h\ddot{w}=-D{\mathrm{\Delta }}^{2}w$$

• Boundary conditions

(86) $$-D{\left(w_{,rr}+\frac{\nu }{r}{w}_{,r}+\frac{\nu }{r^2}{w}_{,\theta \theta }\right)}_R=0$$

(87) $$-D{\left(w_{,rrr}+{\displaystyle \frac{1}{r}}{w}_{,rr}+{\displaystyle \frac{3-\nu }{r^3}}{w}_{,\theta \theta }+{\displaystyle \frac{\nu -2}{r^2}}{w}_{,r\theta \theta }+{\displaystyle \frac{1}{r^2}}{w}_{,r}\right)}_R=0$$

In the equation (86) is the radial bending moment equal to zero and the equation (87) is the effective radial shear force equal to zero in the border of the plate.

6. Discussion, Conclusions and Perspectives

In this paper, the equations of motion and boundary conditions for small oscillations in an elastic solid were reviewed by means of variational calculus. In particular, it has been applied, with a detailed calculation, to the case of transverse oscillations in a thin plate with rectangular and circular geometry, which is rarely found in the current literature.

In section (²), a review of continuous dynamics elements such as stress and strain was made. For the case of an infinitesimal element which only undergoes strain, but not rotations, it was shown that such a strain can be characterized by symmetric strain matrix. To quantify the distribution of internal forces, the traction vector was defined, which is related, by means of Cauchy formula, to a symmetric matrix called stress matrix. The Cauchy formula indicates that each traction component can be expressed as a linear combination of particular stress components. Furthermore, both traction and stress have the same units (force per unit area), however, they are fundamentally different, since traction is a vector and stress is a tensor.

The review was also focused on the interpretation of these quantities and their mathematical properties. This effort was useful to obtain the potential energy of the solid as well as for the calculation of virtual work. For the purpose of applications, their expressions in Cartesian and cylindrical coordinates were obtained explicitly.

The mechanism to obtain, simultaneously, the equations of motion and the boundary conditions for any three-dimensional solid were obtained in the section ³. The algorithm consists of determining the potential energy via the virtual work due to infinitesimal displacements. Then, Hamilton’s principle gives us the equations of motion via Gauss’ theorem to obtain the integrals proportional to $\delta u_i$. It should be noted that equations (50) and (51), are coupled with respect to the deformation displacements via the stress strain relations.

Finally, in sections ⁴ and ⁵-as a matter of example- it is shown the application of the algorithm to the two-dimensional problem of rectangular and circular thin plates. In both cases, the variations are reduced until factors proportional to $\sim \delta w_i$ and $\sim \delta w_i^{\prime }$ are obtained. It is interesting to note that in both cases, the vertical displacement w can be viewed as a field in the two-dimensional space of x and y coordinates. The latter can serve as a guide to obtain these equations from the general three-dimensional case represented in (50) and (51).

The solutions for the thin rectangular plate are not obtained analytically, but by numerical methods. In [¹⁷], the author assumes that the shape of the solution is similar to that of the beams, and calls this type of solution beam functions, uses Rayleigh’s method to derive an approximation for the frequencies for all modes of vibration.

However, it is possible to obtain analytically the eigenfunctions of the biharmonic equation for circular thin plates, as developed in [⁵, ¹⁸], here an analytical solution is found in terms of the first and second type Bessel functions and coefficients depending on the boundary conditions.

Acknowledgments

The authors thanks to Vicerrectorado de Investigación – Universidad Nacional de Ingeniería for supporting the project FC-PF-21-2022.

References

• ^[1]
  E.F.F. Chladni, Entdeckungen über die Theorie des Klanges (Weidmanns Erben und Reich, Leipzig, 1787).
• ^[2]
  H.J. Stöckmann, Eur. Phys. J. Special Topics 145, 15 (2007).
• ^[3]
  M.J. Gander and F. Kwok, SIAM Review 54, 573 (2012).
• ^[4]
  G. Kirchhoff, J. Reine Angew. Math. 40, 51 (1850).
• ^[5]
  J.W.S. Rayleigh, The Theory of Sound (Dover, New York, 1945), v. 1.
• ^[6]
  R.S. Santos, P.S. Camargo Filho and Z.F. Rocha, Rev. Bras. Ensino Fís. 40, e2602 (2018).
• ^[7]
  M.J. Gander and G. Wanner, SIAM Review 54, 627 (2012).
• ^[8]
  T.Y. Wu, Y.Y. Wang and G.R. Liu, Comput. Methods Appl. Mech. Eng. 191, 5365 (2002).
• ^[9]
  H.Z. Gu and X.W. Wang, Journal of Sound and Vibration 202, 452 (1997).
• ^[10]
  H.T. Saliba, Journal of Sound and Vibration 94, 381 (1984).
• ^[11]
  H.S. Yalcin, A. Arikoglu and I. Ozkol, Applied Mathematics and Computation 212, 377 (2009).
• ^[12]
  R.J. Atkin and N. Fox, An Introduction to the Theory of Elasticity (Dover, New York, 1980).
• ^[13]
  M.H. Saad, Elasticity: Theory, Aplications and Numerics (Academic Press, London, 2021).
• ^[14]
  A.P.S. Selvadurai, Partial Differential Equations in Mechanics 2. The Biharmonic Equation. Poisson’s Equation (Springer, Berlin, 2000).
• ^[15]
  J.N. Reddy, Theory and Analysis of Elastic Plates and Shells (CRC Press, Boca Raton, 2006).
• ^[16]
  M. Ducceschi, Nonlinear Vibrations of Thin Rectangular Plates. A Numerical Investigation with Application to Wave Turbulence and Sound Synthesis Doctoral Thesis, ENSTA ParisTech, Palaiseau (2014).
• ^[17]
  G. Warburton, Proceedings of the Institution of Mechanical Engineers 168, 371 (1954).
• ^[18]
  A.W. Leissa, Journal of Sound and Vibration 31, 257 (1973).

• 1
  Vibration frequencies are related to the zeros of Bessel functions which are computed in an approximated way.
• 2
  The term ${\partial }_k{u}_i$ is written as

  $$\frac{\partial u_i}{\partial x_k}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_k}+\frac{\partial u_k}{\partial x_i}\right)+\frac{1}{2}\left(\frac{\partial u_i}{\partial x_k}-\frac{\partial u_k}{\partial x_i}\right)$$
• 3
  In linear elasticity, the Lamé parameters are two elastic constants that completely characterize the linear elastic behaviour of an isotropic solid in small deformations, these two parameters are designated as: $\lambda$, known as the first Lamé parameter and $\mu$, known as the transverse modulus of elasticity
• 4
  In engineering, a flexure is the effect caused by loads external to the plate, which can be forces perpendicular to the plane of the plate, or moments contained in said plane.
• 5
  The Kirchhoff hypothesis consists of the following three parts

  1. Straight lines perpendicular to the mid-surface before deformation remain straight after deformation.

  2. The transverse normals are inextensible.

  3. The transverse normals rotate such that they remain perpendicular to the middle surface after deformation.

Publication Dates

• Publication in this collection
  09 Mar 2022
• Date of issue
  2022

History

• Received
  04 Nov 2021
• Reviewed
  13 Jan 2022
• Accepted
  14 Feb 2022