Buckling Analysis of Laminated Anisotropic Kirchhoff’s Plates via The Boundary Element Method

Abstract

A new fundamental solution for laminated anisotropic Kirchhoff’s plates with out-of-plane and in-plane compressive loads is derived here. The multicompressed solution for both isotropic and anisotropic cases is obtained via the Radon Transform. Some fundamental kernels of the integral equations are described in detail. BEM results of displacements and critical buckling loads of several plates with different boundary conditions and geometries are presented. Comparisons with available analytical solutions and some published numerical results confirm the reliability and accuracy of the proposed formulation.

Keywords
Fundamental solution; plates; anisotropy; laminated composites; buckling

1 INTRODUCTION

Plates are important structural elements used in many engineering fields. For modern applications, such as in the aeronautical and automobile industries, there has been an increasing use of laminated anisotropic plates. Their usage in structural projects is highly advantageous since they are lighter and stronger than traditional materials like metals and metal alloys.

Many problems in anisotropic media cannot be solved analytically and so numerical methods are of paramount importance. The Boundary Element Method (BEM) is a well established alternative to the Finite Element Method (FEM). BEM discretizes the problem only on its boundary, which allows a reduction of the problem's dimension. When compared with FEM, BEM presents high precision on field variables.

Although BEM is highly efficient compared with FEM, it presents some difficulties such as the derivation of fundamental solutions and treatment of the domain integrals which appear e.g. when dealing with plate bending problems. On the other hand, unlike FEM, BEM does not suffer from numerical problems such as the shear locking when using the Shear Deformation Theory (SDT) of plates. Also BEM does not suffer from locking phenomenon due to high stiffness in preferred directions in anisotropic elasticity like FEM.

The bending problem of thin plates can be modelled by the Classical Plate Theory (CPT) (or Kirchhoff-Love theory) when the out-of-plane displacement is small and shear deformations are negligible. The first works using BEM to analyze the bending of thin plates with Kirchhoff-love's theory were presented in the 1970s. Some of those works can be found in Bézine (1978), Altiero and Sikarskie (1978), and Stern (1979). More recently, Dirgantara and Aliabadi (1999, 2006) investigated the bending and the large deflection of shells. Albuquerque et al. (2006) studied the bending of laminated composite Kirchhoff plates. Katsikadelis and Babouskos (2009) used BEM with the analog equation method to study bending and vibration of thick plates.

When there are compressive forces acting on the plate's midplane, an instability known as buckling may take place. Buckling may cause catastrophic failures to structures like plates. It is a highly non-linear phenomenon which is assessed with non-linear theories like the well known von Karman theory of plates. However, linearized theories are often used as means of obtaining critical loads for the pre-buckling stage. Such critical loads are extremely important as failure criteria for engineering applications. One of the first works applying BEM to plate buckling is Shi and Bézine (1990). Shi (1990) also studied the vibrational and buckling problem of thin plates with the direct boundary element method. Other works treating the buckling of plates and shells via the BEM can be seen in Syngellakis and Elzein (1994), and Baiz and Aliabadi (2007). Other numerical methods have been employed in the past for obtaining critical loads for plates within the realm of the linearized theory. For multi-compressed isotropic and/or anisotropic polygonal plates we may cite Liew and Wang (1995), Ferreira et al. (2011), and Shojaee et al. (2012).

In this work we present numerical results for the problem of anisotropic Kirchhoff plates with the action of compressive and shear in-plane forces combined with an out-of-plane perturbation. The critical loads for several isotropic/anisotropic plates are obtained using a concentrated out-of-plane force. Hence, we avoid the use of domain integrals using a pure Boundary Element Method. For multi-compressed isotropic and anisotropic plates, new fundamental solutions are derived here via the Radon Transform. The results are compared with existing analytical or FEM results, showing an excellent agreement.

2 Governing equations

According to Kirchhoff's theory the transverse displacement $w$ of multicompressed thin plates (Figure 1) can be described by the following differential equation (see Lekhnitskii (1968) and Timoshenko and Woinowsky-Krieger (1959))

$\begin{array}{l}{D}_{11}\frac{{\partial }^{4}w}{\partial x^4}+4D_{16}\frac{{\partial }^{4}w}{\partial x^3\partial y}+2(D_{12}+2D_{66})\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+4D_{26}\frac{{\partial }^{4}w}{\partial x\partial y^3}+D_{22}\frac{{\partial }^{4}w}{\partial y^4}+\\ N_x\frac{{\partial }^{2}w}{\partial x^2}+2N_{xy}\frac{{\partial }^{2}w}{\partial x\partial y}+N_y\frac{{\partial }^{2}w}{\partial y^2}=q.\end{array}$ (1)

Figure 1
Compressed plate.

With $D_{ij}$ denoting the flexural rigidities. Since we shall focus on symmetric laminated anisotropic plates, the final flexural rigidities can be obtained by adding the sequence of flexural rigidities of each rotated orthotropic ply.

$D_{ij}=\frac{1}{3}{\displaystyle \sum _{k=1}^l{({\overline{Q}}_{ij})}_k(z_k^3-z_{k-1}^3)}.$ (2)

With $z_k$ representing the distance from the middle plane to the ply's surfaces, k and k-1 (see Figure 2). ${({\overline{Q}}_{ij})}_k$ is the reduced rotated matrix in the plane stress state with the elastic constants of each ply (Eq. (3)). In Eq. (2) $i,j=\mathrm{1,2,6}$ and l is the number of plies. Notice that membrane forces are decoupled from bending forces when the laminate is symmetric around the midplane, allowing the use of Eq. (1) (see Christensen (1979)).

Figure 2
A symmetric laminated composite with N plies.

$\begin{array}{l}{\overline{Q}}_{11}^k=\frac{1}{8\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left[3E_l^k+3E_t^k+2E_t^k{\nu }_{lt}^k+4G_{lt}^k\left(1-{\left({\nu }_{lt}^k\right)}^2\right)\right]+\\ \frac{\cos \left(2{\theta }^k\right)}{2\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left(E_l^k-E_t^k\right)+\frac{\cos \left(4{\theta }^k\right)}{8\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left[E_l^k+E_t^k-2E_t^k{\nu }_{lt}^k-4G_{lt}^k\left(1-{\left({\nu }_{lt}^k\right)}^2\right)\right],\\ {\overline{Q}}_{12}^k=\frac{1}{8\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left[E_l^k+E_t^k+6E_t^k{\nu }_{lt}^k-4G_{lt}^k\left(1-{\left({\nu }_{lt}^k\right)}^2\right)\right]-\\ \frac{\cos \left(4{\theta }^k\right)}{8\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left[E_l^k+E_t^k-2E_t^k{\nu }_{lt}^k-4G_{lt}^k\left(1-{\left({\nu }_{lt}^k\right)}^2\right)\right],\\ {\overline{Q}}_{16}^k=-\frac{\sin \left(2{\theta }^k\right)}{4\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left(E_l^k-E_t^k\right)-\frac{\sin \left(4{\theta }^k\right)}{8\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left[E_l^k+E_t^k-2E_t^k{\nu }_{lt}^k-4G_{lt}^k\left(1-{\left({\nu }_{lt}^k\right)}^2\right)\right],\\ {\overline{Q}}_{22}^k=\frac{1}{8\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left[3E_l^k+3E_t^k+2E_t^k{\nu }_{lt}^k+4G_{lt}^k\left(1-{\left({\nu }_{lt}^k\right)}^2\right)\right]\\ -\frac{\cos \left(2{\theta }^k\right)}{2\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left(E_l^k-E_t^k\right)+\frac{\cos \left(4{\theta }^k\right)}{8\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left[E_l^k+E_t^k-2E_t^k{\nu }_{lt}^k-4G_{lt}^k\left(1-{\left({\nu }_{lt}^k\right)}^2\right)\right],\\ {\overline{Q}}_{26}^k=-\frac{\sin \left(2{\theta }^k\right)}{4\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left(E_l^k-E_t^k\right)+\frac{\sin \left(4{\theta }^k\right)}{8\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left[E_l^k+E_t^k-2E_t^k{\nu }_{lt}^k-4G_{lt}^k\left(1-{\left({\nu }_{lt}^k\right)}^2\right)\right],\\ {\overline{Q}}_{66}^k=\frac{1}{8\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left[E_l^k+E_t^k-2E_t^k{\nu }_{lt}^k-4G_{lt}^k\left(1-{\left({\nu }_{lt}^k\right)}^2\right)\right]\\ -\frac{\cos \left(4{\theta }^k\right)}{8\left(1-{\left({\nu }_{lt}^k\right)}^2\right)}\left[E_l^k+E_t^k-2E_t^k{\nu }_{lt}^k-4G_{lt}^k\left(1-{\left({\nu }_{lt}^k\right)}^2\right)\right].\end{array}$ (3)

With $E_l$representing the elastic modulus along the direction of the fiber, $E_t$representing the elastic modulus in the transversal direction to the fiber axis, ${\nu }_{lt}$ representing the Poisson ratio, $G_{lt}$representing the shear modulus and θk being the angle ply orientation. Eq. (3) was obtained based on Christensen (1979) and Herakovich (1998).

3 Boundary Integral Equation

Using the weighted residual method in Eq. (1), we obtain the singular Boundary Integral Equation (BIE) for compressed anisotropic Kirchhoff's plates (see Albuquerque et al. (2006)),

$\begin{array}{l}{c}_{i}w(x_i)+{\displaystyle \underset{\Gamma }{\oint }{V}_n^{*}(x,x_i)w(x)d\Gamma }-{\displaystyle \underset{\Gamma }{\oint }{M}_n^{*}(x,x_i)\frac{\partial w(x)}{\partial n}d\Gamma }+{\displaystyle \sum _{j=1}^{N_c}{R}_{c_j}^{*}(x,x_i)w_{c_j}(x)}={\displaystyle \underset{\Gamma }{\oint }{V}_n^{}(x)w^{*}(x,x_i)d\Gamma }\\ -{\displaystyle \underset{\Gamma }{\oint }{M}_n^{}(x)\frac{\partial w^{*}(x,x_i)}{\partial n}d\Gamma }+{\displaystyle \sum _{j=1}^{N_c}{R}_{c_j}^{}(x)w_{c_j}^{*}(x,x_i)+{\displaystyle \underset{\Omega }{\int }q(x)w^{*}(x,x_i)d\Omega .}}\end{array}$ (4)

With Γ and Ω representing the plate's boundary and domain, respectively; $x,x_i$ represent the field and source point, respectively; c_i is a free term depending on the smoothness of the plate’s boundary; $M_n$ represents the resultant moment (see Eq. (5)). $R_{c_j}$ is the corner reaction, described as $R_{c_j}=M_{ns}^{+}-M_{ns}^{-}$; $M_{ns}^{\pm }$ are the twisting moments after (+) and before (-) the corner (see Eq. (6)); $w_{c_j}$is the corner displacement and Nc is the total number of corners. $V_n$ represents here the resultant shear force for polygonal plates, which considers the effect of the out-of-plane and in-plane forces (see Eq. (8)). The starred terms are the known fundamental Kernels, which are obtained deriving the fundamental solution. $M_n^{*}$, $M_{ns}^{*}$, $Q_n^{*}$, $V_n^{*}$ are concisely expressed as

$M_n^{*}=n_x{}^2{M}_x^{*}+2n_x{n}_y{M}_{xy}^{*}+n_y{}^2{M}_y^{*}.$ (5)

$M_{ns}^{*}=n_x{n}_y\left(M_y^{*}-M_x^{*}\right)+\left(n_x^2-n_y^2\right)M_{xy}^{*}.$ (6)

With

$\begin{array}{l}{M}_x^{*}=-\left(D_{11}\frac{{\partial }^2{w}^{*}}{\partial x^2}+2D_{16}\frac{{\partial }^2{w}^{*}}{\partial x\partial y}+D_{12}\frac{{\partial }^2{w}^{*}}{\partial y^2}\right),\\ M_{xy}^{*}=-\left(D_{16}\frac{{\partial }^2{w}^{*}}{\partial x^2}+2D_{66}\frac{{\partial }^2{w}^{*}}{\partial x\partial y}+D_{26}\frac{{\partial }^2{w}^{*}}{\partial y^2}\right),\\ M_y^{*}=-\left(D_{12}\frac{{\partial }^2{w}^{*}}{\partial x^2}+2D_{26}\frac{{\partial }^2{w}^{*}}{\partial x\partial y}+D_{22}\frac{{\partial }^2{w}^{*}}{\partial y^2}\right).\end{array}$ (7)

With $n_x$ and $n_y$ denoting the direction cosines.

$V_n^{*}=Q_n{}^{*}+\frac{\partial M_{ns}^{*}}{\partial s}+N_n\frac{\partial w_{}^{*}}{\partial n}+N_{ns}\frac{\partial w_{}^{*}}{\partial s},$ (8)

$\begin{array}{l}{Q}_n{}^{*}=-\left[D_{11}\frac{{\partial }^3{w}^{*}}{\partial x^3}+3D_{16}\frac{{\partial }^3{w}^{*}}{\partial x^2\partial y^{}}+\left(D_{12}+2D_{66}\right)\frac{{\partial }^3{w}^{*}}{\partial x\partial y^2}+D_{16}\frac{{\partial }^3{w}^{*}}{\partial y^3}\right]n_x\\ -\left[D_{16}\frac{{\partial }^3{w}^{*}}{\partial x^3}+\left(D_{12}+2D_{66}\right)\frac{{\partial }^3{w}^{*}}{\partial x^2\partial y^{}}+3D_{26}\frac{{\partial }^3{w}^{*}}{\partial x\partial y^2}+D_{22}\frac{{\partial }^3{w}^{*}}{\partial y^3}\right]n_y,\end{array}$ (9)

$N_n=N_x{n}_x^2+2N_{xy}{n}_x{n}_y+N_y{n}_y^2,\begin{array}{cc}& N_{ns}=n_x{n}_y\left(N_y-N_x\right)+N_{xy}\left(n_x^2-n_y^2\right).\end{array}$ (10)

∂/∂n and ∂/∂s are the derivatives in the normal and tangential directions, respectively, given by

$\frac{\partial }{\partial n}=\frac{\partial }{\partial x}{n}_x+\frac{\partial }{\partial y}{n}_y,\begin{array}{cc}& \frac{\partial }{\partial s}=-\frac{\partial }{\partial x}{n}_y+\frac{\partial }{\partial y}{n}_x.\end{array}$ (11)

It is well known that due to the number of variables ($M_n$, $V_n$, $\partial w/\partial n$, $\partial w/\partial n$) a second equation is needed to fully solve the problem of plates. The other equation is obtained by deriving Eq. (4) in the direction of the outward vector, $n_i$, normal to the boundary surface at the source point. This is a hypersingular BIE since there are kernels, which present (1/r2) singularities. The latter require special analytical or numerical treatment since they are hypersingular in the sense of Hadamard.

$\begin{array}{l}{c}_i\frac{\partial w(x_i)}{\partial n_i}+{\displaystyle \underset{\Gamma }{\oint }\frac{\partial V_n^{*}(x,x_i)}{\partial n_i}w(x)d\Gamma }-{\displaystyle \underset{\Gamma }{\oint }\frac{\partial M_n^{*}(x,x_i)}{\partial n_i}\frac{\partial w(x)}{\partial n}d\Gamma }+{\displaystyle \sum _{j=1}^{N_c}\frac{\partial R_{c_j}^{*}(x,x_i)}{\partial n_i}{w}_{c_j}(x)}=\\ {\displaystyle \underset{\Gamma }{\oint }{V}_n^{}(x)\frac{\partial w^{*}(x,x_i)}{\partial n_i}d\Gamma }-{\displaystyle \underset{\Gamma }{\oint }{M}_n^{}(x)\frac{{\partial }^2{w}^{*}(x,x_i)}{\partial n\partial n_i}d\Gamma }+{\displaystyle \sum _{j=1}^{N_c}{R}_{c_j}^{}(x)\frac{\partial w_{c_j}^{*}(x,x_i)}{\partial n_i}+{\displaystyle \underset{\Omega }{\int }q(x)\frac{\partial w^{*}(x,x_i)}{\partial n_i}d\Omega .}}\end{array}$ (12)

We used Eq. (12) for the isotropic case and hence we will present on following the detailed derivation of each Kernel for the singular and hypersingular BIEs. For the anisotropic case, due to the high complexity of the numerical integrations, we used the strategy of outside points of the boundary plate as in Paiva and Venturini (1992), hence using only the singular BIE. Moreover, for a concentrated out of plane load (Dirac's Delta load) the domain integrations disappear. It is then only necessary to evaluate the fundamental solution at the point of application as $w^{*}(x_i,x=x_i)$.

4 Fundamental solution

4.1 Isotropic media

In this work, two fundamental solutions for isotropic plates are considered. One for the case with compressive forces only in one direction (unidirectionally compressed) and another for the case with the combined action of compressive forces in both direction with shear in-plane forces (multidirectionally compressed).

4.1.1 Unidirectionally compressed plates

For this case, the governing equation is

$D{\nabla }^{4}w+N_x\frac{{\partial }^{2}w}{\partial x^2}=q.$ (13)

Jahanshahi and Dundurs (1964) constructed a solution for the problem of plates compressed in one direction and a moving out-of-plane load. Adapting their solution to the case of a fixed out-of-plane load, we have

$w^{*}=\frac{1}{8D\kappa }{\displaystyle \underset{0}{\overset{|x-x_i|}{\int }}sin(\kappa \rho )}{\mathrm{Y}}_0\left(\kappa \sqrt{{\rho }^2+{(y-y_i)}^2}\right)d\rho +\frac{1}{8D}{\displaystyle \underset{0}{\overset{|y-y_i|}{\int }}(|y-y_i|-\rho )}{\mathrm{Y}}_0(\kappa \rho )d\rho ,$ (14)

Where $(x,y)$ and $(x_i,y_i)$ are the coordinates of the field and source points, respectively.$D$ is the flexural rigidity, given by $D=E{h}^3/(12(1-{\nu }^2))$; $\kappa =N_x/(4D)$; ${\mathrm{Y}}_0$ and ${\mathrm{Y}}_1$ are the Bessel functions of the second kind of order 0 and 1, respectively. Notice that $w^{*}$ and $\partial w^{*}/\partial n$ are not explicit, requiring further numerical integration. However, all other kernels are explicit. Deriving Eq. (14), we find

$\begin{array}{l}\frac{\partial w^{*}}{\partial x}=\frac{1}{8D\kappa }sin(\kappa (x-x_i)){\mathrm{Y}}_0\left(\kappa r\right),\\ \frac{\partial w^{*}}{\partial y}=\frac{1}{8D}\mathrm{sign}(y-y_i){\displaystyle \underset{0}{\overset{|y-y_i|}{\int }}{\mathrm{Y}}_0(\kappa \rho )d\rho }\\ -\frac{1}{8D}(y-y_i){\displaystyle \underset{0}{\overset{|x-x_i|}{\int }}\frac{sin(\kappa \rho )}{\sqrt{{\rho }^2+{(y-y_i)}^2}}}{\mathrm{Y}}_1(\kappa \sqrt{{\rho }^2+{(y-y_i)}^2})d\rho ,\\ \frac{{\partial }^2{w}^{*}}{\partial x^2}=\frac{1}{8D}\left[\cos (\kappa (x-x_i)){\mathrm{Y}}_0\left(\kappa r\right)-r_{,x}\sin (\kappa (x-x_i)){\mathrm{Y}}_1\left(\kappa r\right)\right],\\ \frac{{\partial }^2{w}^{*}}{\partial x\partial y}=-\frac{1}{8D}{r}_{,y}\sin (\kappa (x-x_i)){\mathrm{Y}}_1\left(\kappa r\right),\\ \frac{{\partial }^2{w}^{*}}{\partial y^2}=\frac{1}{8D}\left[\cos (\kappa (x-x_i)){\mathrm{Y}}_0\left(\kappa r\right)+r_{,x}\sin (\kappa (x-x_i)){\mathrm{Y}}_1\left(\kappa r\right)\right].\end{array}$ (15)

With $r=\sqrt{{(x-x_i)}^2+{(y-y_i)}^2}$, $r_{,x}=(x-x_i)/r$ and $r_{,y}=(y-y_i)/r$. Using the following useful identity

$\frac{\partial {\mathrm{Y}}_1(\kappa r)}{\partial r}=\kappa {\mathrm{Y}}_0(\kappa r)-\frac{{\mathrm{Y}}_1(\kappa r)}{r},$ (16)

we may express the third order derivatives as

$\begin{array}{l}\frac{{\partial }^3{w}^{*}}{\partial x^3}=\frac{1}{8D}[(r_{,y}{}^2-2)\kappa \sin (\kappa (x-x_i)){\mathrm{Y}}_0\left(\kappa r\right)-2r_{,x}\kappa \cos (\kappa (x-x_i)){\mathrm{Y}}_1\left(\kappa r\right)\\ +\frac{1-2r_{,y}{}^2}{r}\sin (\kappa (x-x_i)){\mathrm{Y}}_1\left(\kappa r\right)],\\ \frac{{\partial }^3{w}^{*}}{\partial x^2\partial y}=\frac{1}{8D}[-r_{,x}{r}_{,y}\kappa \sin (\kappa (x-x_i)){\mathrm{Y}}_0\left(\kappa r\right)-r_{,y}\kappa \cos (\kappa (x-x_i)){\mathrm{Y}}_1\left(\kappa r\right)\\ +\frac{2r_{,x}{r}_{,y}}{r}\sin (\kappa (x-x_i)){\mathrm{Y}}_1\left(\kappa r\right)],\\ \frac{{\partial }^3{w}^{*}}{\partial x\partial y^2}=-\frac{1}{8D}\left[r_{,y}{}^2\kappa \sin (\kappa (x-x_i)){\mathrm{Y}}_0\left(\kappa r\right)+\frac{1-2r_{,y}{}^2}{r}\sin (\kappa (x-x_i)){\mathrm{Y}}_1\left(\kappa r\right)\right],\\ \frac{{\partial }^3{w}^{*}}{\partial y^3}=\frac{1}{8D}[r_{,x}{r}_{,y}\kappa \sin (\kappa (x-x_i)){\mathrm{Y}}_0\left(\kappa r\right)-r_{,y}\kappa \cos (\kappa (x-x_i)){\mathrm{Y}}_1\left(\kappa r\right)\\ -\frac{2r_{,x}{r}_{,y}}{r}\sin (\kappa (x-x_i)){\mathrm{Y}}_1\left(\kappa r\right)].\end{array}$ (17)

Substituting the derivatives of Eq. (15) and Eq. (17) in Eqs. (5)-(11) we obtain the fundamental kernels of Eq. (4) and Eq. (12). For the sake of brevity we shall present here only one of the kernels, $\partial V_n^{*}/\partial n_i$ which has been derived, to the author's best knowledge, here for the first time.

$\begin{array}{l}\frac{\partial V_n^{*}}{\partial n_i}=D_1{\kappa }^2\cos \left[\kappa (x-x_i)\right]{\mathrm{Y}}_0\left(\kappa r\right)+D_2\frac{\kappa }{r}\cos \left[\kappa (x-x_i)\right]{\mathrm{Y}}_1\left(\kappa r\right)+D_3\frac{\kappa }{r}\sin \left[\kappa (x-x_i)\right]{\mathrm{Y}}_0\left(\kappa r\right)\\ +D_4{\kappa }^2\sin \left[\kappa (x-x_i)\right]{\mathrm{Y}}_1\left(\kappa r\right)+D_5\frac{1}{r^2}\sin \left[\kappa (x-x_i)\right]{\mathrm{Y}}_1\left(\kappa r\right).\end{array}$ (18)

Expressions for D₁ to D₅ can be found in ^{Appendix A}. The other kernels are presented in Monteiro Jr. and Daros (2017) and Monteiro Jr. (2017). For the sake of conciseness, we will only show here results for multidirectionally compressed plates.

4.1.2 Multidirectionally compressed plates

The governing equation in this case is

$D{\nabla }^{4}w+N_x\frac{{\partial }^{2}w}{\partial x^2}+2N_{xy}\frac{{\partial }^{2}w}{\partial x\partial y}+N_y\frac{{\partial }^{2}w}{\partial y^2}=q.$ (19)

We shall present here a new fundamental solution for this equation via the Radon Transform. The latter has been successfully used for several partial differential equations. For a complete review of the Radon Transform the reader is referred to Ludwig (1966) and Helgason (1980) where a thorough description of the Transform properties are given. The Radon Transform is defined as

$\mathcal{R}(f)=\widehat{f}(s,m)={\displaystyle {\int }_{x\cdot m=s}f(x)dS={\displaystyle \int f(x)\delta (s-x\cdot m)dx}}.$ (20)

Where $s=x\cdot m$is the transform hyperplane and $m$is the unit radius of the circle (Figure 3). Using the Radon Transform we find

$\begin{array}{l}\mathcal{R}\left(\delta (x,\xi )\right)=\delta (s-m\cdot \xi ),\\ \mathcal{R}(\mathcal{L}(\partial x)f(x))=\mathcal{L}(m)\frac{{\partial }^k}{\partial s^k}\widehat{f}(s,m).\end{array}$ (21)

Figure 3
Radon circle.

With $\xi$ denoting the source point and ℒ represents here a linear partial differential operator. Applying the Radon Transform to Eq. (19), we find

$\left[d^{2}s(d^{2}s+\beta )\right]w^{*}(s,m)=\frac{1}{D}\delta (s-m\cdot \xi ),$ (22)

With

$\beta =\frac{N_x}{D}{m}_1^2+2\frac{N_{xy}}{D}{m}_1{m}_2+\frac{N_y}{D}{m}_2^2.$ (23)

With $m_1=\cos (\theta )$ and $m_2=\sin (\theta )$. Based on Manolis et al. (2003), we propose, for $\beta >0$, Eq. (24) as a solution for Eq. (22),

$w^{*}(s,m)=\frac{1}{2D\beta }\left[|s-m\cdot \xi |+\frac{\text{i}}{\sqrt{\beta }}{e}^{\text{i}\sqrt{\beta }|s-m\cdot \xi |}\right].$ (24)

The inverse Radon transform is defined via the equation

$f(x)={\mathcal{R}}^{-1}(\widehat{f}(s,m))=\frac{1}{4{\pi }^2}{{\displaystyle \underset{\left|m\right|=1}{\int }\left({\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{{\partial }_{\sigma }f(p,m)}{s-\sigma }d\sigma }\right)}|}_{s=m\cdot x}dm.$ (25)

Applying the Inverse Radon Transform in Eq. (24), we find

$\begin{array}{l}{\mathcal{R}}^{-1}(|s-m\cdot \xi |)=\frac{1}{4{\pi }^2}{\displaystyle \underset{\left|m\right|=1}{\int }{2\log (|s-m\cdot \xi |)|}_{s=m\cdot x}}dm,\\ {\mathcal{R}}^{-1}\left(\frac{{\text{ie}}^{\text{i}k|s-m\cdot \xi |}}{2k}\right)=\frac{1}{4{\pi }^2}{\displaystyle \underset{\left|m\right|=1}{\int }\frac{1}{2}\{\text{i}\pi \cos \left(k|s-m\cdot \xi |\right)-}\\ {2\left[\text{ci}\left(k|s-m\cdot \xi |\right)\text{cos}\left(k|s-m\cdot \xi |\right)+\text{si}\left(k|s-m\cdot \xi |\right)\text{sin}\left(k|s-m\cdot \xi |\right)\right]\}|}_{{}_{s=m\cdot x}}dm,\end{array}$ (26)

With $k=\sqrt{\beta }$, $\mathrm{si}$ is the sine integral function and $\mathrm{ci}$ is the cosine integral function, defined in Eqs. (27.1) and (27.2).

$\text{ci}\left(x\right)=\gamma +\log \left(x\right)+{\displaystyle \underset{0}{\overset{x}{\int }}\frac{\cos \left(t\right)-1}{t}dt,}\begin{array}{cc}& \text{si}\left(x\right)={\displaystyle \underset{0}{\overset{x}{\int }}\frac{\sin \left(t\right)}{t}dt.}\end{array}$ (27)

With $\gamma$ representing the gamma-Euler constant. Applying Eq. (26) on Eq. (24), we have

$w^{*}(x,\xi )=\frac{1}{4{\pi }^{2}D}{\displaystyle \underset{0}{\overset{2\pi }{\int }}\frac{1}{k_{}^2}\left\{\log (z)+\frac{1}{2}\left[\text{i}\pi \cos \left(kz\right)-2(\text{ci}\left(kz\right)\text{cos}\left(kz\right)+\text{si}\left(kz\right)\text{sin}\left(kz\right))\right]\right\}d\theta ,}$ (28)

With $z=|s-m\cdot \xi |$. However, fundamental solutions of elliptic operators are not unique and we may try other solutions. As in Rangelov et al. (2005), we try a solution with only the real part of Eq. (28). Removing the complex term of Eq. (28), thus

$w^{*}(x,\xi )=\frac{1}{4{\pi }^{2}D}{\displaystyle \underset{0}{\overset{2\pi }{\int }}\frac{1}{k_{}^2}\left\{\log (z)-\left[\text{ci}\left(kz\right)\text{cos}\left(kz\right)+\text{si}\left(kz\right)\text{sin}\left(kz\right)\right]\right\}}d\theta .$ (29)

Substituting Eq. (29) in Eq. (19), we find

$D{\nabla }^{4}w+N_x\frac{{\partial }^{2}w}{\partial x^2}+2N_{xy}\frac{{\partial }^{2}w}{\partial x\partial y}+N_y\frac{{\partial }^{2}w}{\partial y^2}=\frac{1}{4{\pi }^2}{\displaystyle \underset{\left|m\right|=1}{\int }\frac{1}{{\left(m\cdot (x-\xi )\right)}^2}dm}.$ (30)

Which is the integral representation of Dirac's Delta generalized function (see Rangelov et al. (2005)).

A great advantage of using the Radon Transform is that the derivatives are straightforward. We present on following some of the kernels obtained for the singular and hypersingular BIEs.

$\frac{\partial w^{*}}{\partial n}=\frac{1}{4\pi {}^{2}D}{\displaystyle \underset{0}{\overset{2\pi }{\int }}\frac{1}{k}\left[\text{ci}\left(kz\right)\sin \left(kz\right)-\text{si}\left(kz\right)\cos \left(kz\right)\right]\text{sign}\left(m_1{r}_{,x}+m_2{r}_{,y}\right)\left(m_1{n}_x+m_2{n}_y\right)}d\theta .$ (31)

$M_n^{*}=-\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}a\left[\mathrm{ci}\left(kz\right)\cos \left(kz\right)+\mathrm{si}\left(kz\right)\sin \left(kz\right)\right]}d\theta ,$ (32)

$a=n_x^2\left(m_1^2+\nu m_2^2\right)+2n_x^{}{n}_y^{}\left(1-\nu \right)m_1^{}{m}_2^{}+n_y^2\left(m_2^2+\nu m_1^2\right).$ (33)

The resultant moment is a singular kernel in z, so it must be divided in two parts: one regular and other singular. Using the singularity subtraction method, the regular part is given by

$M_n^{R*}=-\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}a\left\{-\log \left(z\right)+\left[\mathrm{ci}\left(kz\right)\cos \left(kz\right)+\mathrm{si}\left(kz\right)\sin \left(kz\right)\right]\right\}}d\theta .$ (34)

Since $I_1$, the singular part can be rewritten as

$M_n^{S*}=\underset{I_1}{\underbrace{-\frac{\log \left(r\right)}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}ad\theta }}}-\underset{I_2}{\underbrace{\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}a\log \left(\left|m_1{r}_{,x}+m_2{r}_{,y}\right|\right)d\theta }}},$ (35)

The integral $I_1$ can be analytically evaluated and is given by $I_1=-(1-\nu )\log (r)/(4\pi )$. Integral $I_2$ can also be obtained in a closed form as

$I_2=\frac{1}{8\pi }\left\{\left(1-\nu \right)\left[2{\left(\nabla r\cdot n\right)}^2-1\right]+\left(1+\nu \right)\log \left(4\right)\right\}.$ (36)

So after treating the singularity, the resultant moment can be expressed as

$\begin{array}{l}{M}_n^{*}=-\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}-a\left\{-\log \left(z\right)+\left[\mathrm{ci}\left(kz\right)\cos \left(kz\right)+\mathrm{si}\left(kz\right)\sin \left(kz\right)\right]\right\}}d\theta +\\ \frac{1}{8\pi }\left\{\left(1-\nu \right)\left[2{\left(\nabla r\cdot n\right)}^2-1\right]+\left(1+\nu \right)\log \left(4\right)\right\}-\frac{\left(1+\nu \right)}{4\pi }\log \left(r\right).\end{array}$ (37)

Likewise, the total shear force is given as

$\begin{array}{l}{V}_n^{*}=-\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}\frac{d}{k}\left[\mathrm{ci}(kz)\sin (kz)-\mathrm{si}(kz)\cos (kz)\right]}\mathrm{sign}(m_1{r}_x+m_2{r}_y)d\theta \\ -\frac{1}{4\pi r}\left(\nabla r\cdot n\right)\left\{2+\left(1-\nu \right)\left[2{\left(\nabla r\cdot n\right)}^2-1\right]\right\}.\end{array}$ (38)

With $d=c+k^{2}b$ and

$\begin{array}{l}b=n_x\left[1+n_y^2\left(1-\nu \right)\right]m_1^3+n_y\left[1-2n_x^2\left(1-\nu \right)+n_y^2\left(1-\nu \right)\right]m_1^2{m}_2\\ +n_x\left[1-2n_y^2\left(1-\nu \right)+n_x^2\left(1-\nu \right)\right]m_1{m}_2^2+n_y\left[1+n_x^2\left(1-\nu \right)\right]m_2^3,\end{array}$ (39)

$c=\left[\frac{N_x}{D}\left(n_x^3+n_x^{}{n}_y^2\right)+\frac{N_{xy}}{D}\left(n_x^2{n}_y^{}+n_y^3\right)\right]m_1+\left[\frac{N_{xy}}{D}\left(n_x^3+n_x^{}{n}_y^2\right)+\frac{N_x}{D}\left(n_x^2{n}_y^{}+n_y^3\right)\right]m_2.$ (40)

The kernel representing the corner forces are given as

$R_{c_j}^{*}=-\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}f\left[\mathrm{ci}\left(kz\right)\cos \left(kz\right)+\mathrm{si}\left(kz\right)\sin \left(kz\right)\right]}d\theta ,$ (41)

With

$f=(1-\nu )\left\{\left(n_x^{+}{n}_y^{+}+n_x^{-}{n}_y^{-}\right)\left(m_2^2-m_1^2\right)+\left[\left(n_x^{+2}-n_y^{+2}\right)-\left(n_x^{-2}-n_y^{-2}\right)\right]m_1{m}_2\right\}.$ (42)

Doing a similar approach to the resultant moment, we obtain

$\begin{array}{l}{R}_{c_j}^{*}=-\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}f\left\{-\log \left(z\right)+\left[\mathrm{ci}\left(kz\right)\cos \left(kz\right)+\mathrm{si}\left(kz\right)\sin \left(kz\right)\right]\right\}}d\theta \\ -\frac{\left(1-\nu \right)}{4\pi }\left\{\left[\left(n_x^{+2}-n_y^{+2}\right)-\left(n_x^{-2}-n_y^{-2}\right)\right]r_{,x}{r}_{,y}+\left(n_x^{+}{n}_y^{+}-n_x^{-}{n}_y^{-}\right)\left(r_{,y}^2-r_{,x}^2\right)\right\}.\end{array}$ (43)

The other Kernels for the hypersingular BIE can be concisely expressed as

$\begin{array}{l}\frac{{\partial }^2{w}_{}^{*}}{\partial n\partial n_i}=-\frac{1}{4{\pi }^{2}D}{\displaystyle \underset{0}{\overset{2\pi }{\int }}g\left\{-\log \left(z\right)+\left[\mathrm{ci}\left(kz\right)\cos \left(kz\right)+\mathrm{si}\left(kz\right)\sin \left(kz\right)\right]\right\}}d\theta +\\ \frac{1}{8\pi D}\left\{n_x\left[n_{x_i}\left(2r_{,y}^2-1+\log (4)\right)-2n_{y_i}{r}_{,x}{r}_{,y}\right]-2n_y{r}_{,y}\left(\nabla r\cdot n_i\right)+n_y{n}_{y_i}\left(1+\log \left(4\right)\right)\right\}\\ -\frac{\left(n\cdot n_i\right)}{4\pi D}\log \left(r\right),\end{array}$ (44)

With

$g=n_x{n}_{x_i}{m}_1^2+\left(n_y{n}_{x_i}+n_x{n}_{y_i}\right)m_1{m}_2+n_y{n}_{y_i}{m}_2^2.$ (45)

$\begin{array}{l}\frac{\partial M_n^{*}}{\partial n_i}=-\frac{1}{4\pi {}^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}ak\left[\text{ci}\left(kz\right)\sin \left(kz\right)-\text{si}\left(kz\right)\cos \left(kz\right)\right]\text{sign}\left(m_1{r}_{,x}+m_2{r}_{,y}\right)\left(m_1{n}_{xi}+m_2{n}_{yi}\right)}d\theta \\ +\frac{1}{4\pi r}\left\{\left(1+\nu \right)\left(\nabla r\cdot n_i\right)-2\left(1-\nu \right)\left[{\left(\nabla r\cdot n\right)}^2\left(\nabla r\cdot n_i\right)-\left(\nabla r\cdot n\right)\left(n\cdot n_i\right)\right]\right\}.\end{array}$ (46)

$\begin{array}{l}\frac{\partial V_n^{*}}{\partial n_i}=\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}d\left\{\left[\mathrm{ci}\left(kz\right)\cos \left(kz\right)+\mathrm{si}\left(kz\right)\sin \left(kz\right)\right]\left(m_1{n}_{xi}+m_2{n}_{yi}\right)\right\}}\\ -\frac{\left(\nabla r\cdot n\right)\left(\nabla r\cdot n_i\right)}{4\pi r^2}\left\{4+2\left(1-\nu \right)\left[{\left(\nabla r\cdot n\right)}^2-1\right]+4\left(1-\nu \right){\left(\nabla r\cdot n\right)}^2\right\}\\ +\frac{\left(n\cdot n_i\right)}{4\pi r^2}\left\{2+\left(1-\nu \right)\left[{\left(\nabla r\cdot n\right)}^2-1\right]+4\left(1-\nu \right){\left(\nabla r\cdot n\right)}^2\right\}.\end{array}$ (47)

$\begin{array}{l}\frac{\partial R_{c_j}^{*}}{\partial n_i}=-\frac{1}{4\pi {}^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}fk\left[\text{ci}\left(kz\right)\sin \left(kz\right)-\text{si}\left(kz\right)\cos \left(kz\right)\right]\text{sign}\left(m_1{r}_{,x}+m_2{r}_{,y}\right)\left(m_1{n}_{xi}+m_2{n}_{yi}\right)}d\theta \\ -\frac{\left(1-\nu \right)}{4\pi r}\{\left[\left(n_x^{+2}-n_y^{+2}\right)-\left(n_x^{-2}-n_y^{-2}\right)\right]\left[\left(r_{,x}^2{r}_{,y}-r_{,y}^3\right)n_{x_i}+\left(r_{,y}^2{r}_{,x}-r_{,x}^3\right)n_{y_i}\right]+\\ 4\left(n_x^{+}{n}_y^{+}-n_x^{-}{n}_y^{-}\right)\left(r_{,x}{r}_{,y}^2{n}_{x_i}-r_{,x}^2{r}_{,y}{n}_{y_i}\right)\}.\end{array}$ (48)

Notice that Eq. (47) is hypersingular in r for straight boundary elements and must be treated in the sense of Hadamard.

4.2 Anisotropic media

Using a similar approach as in the isotropic case, we obtain Eq. (49) as a fundamental solution for the problem of multidirectionally compressed anisotropic plates

$w^{*}(x,\xi )=\frac{1}{4{\pi }^2{D}_{11}}{\displaystyle \underset{0}{\overset{2\pi }{\int }}\frac{1}{\alpha k_{}^2}\left\{\log (z)+\left[2(\text{ci}\left(kz\right)\text{cos}\left(kz\right)+\text{si}\left(kz\right)\text{sin}\left(kz\right)\right]\right\}}d\theta .$ (49)

Where $\text{k=}\sqrt{\beta /\alpha }$ and

$\alpha =m_1^4+\frac{4D_{16}}{D_{11}}{m}_1^3{m}_2^{}+\frac{2\left(D_{12}+2D_{66}\right)}{D_{11}}{m}_1^2{m}_2^2+\frac{4D_{26}}{D_{11}}{m}_1^{}{m}_2^3+\frac{D_{22}}{D_{11}}{m}_2^4.$ (50)

$\beta =\frac{N_x}{D_{11}}{m}_1^2+2\frac{N_{xy}}{D_{11}}{m}_1{m}_2+\frac{N_y}{D_{11}}{m}_2^2$ (51)

Some of the other kernels associated with the fundamental solution are given on following,

$M_n{}^{*}=-\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}\frac{a}{\alpha }\left\{-\log (z)+\left[\text{ci}\left(kz\right)\text{cos}\left(kz\right)+\text{si}\left(kz\right)\text{sin}\left(kz\right)\right]\right\}}d\theta ,$ (52)

$\begin{array}{l}a=n_x^2\left(m_1^2+2\frac{D_{16}}{D_{11}}{m}_1^{}{m}_2^{}+\frac{D_{12}}{D_{11}}{m}_2^2\right)+2n_x{n}_y\left(\frac{D_{16}}{D_{11}}{m}_1^2+2\frac{D_{66}}{D_{11}}{m}_1^{}{m}_2^{}+\frac{D_{26}}{D_{11}}{m}_2^2\right)+\\ n_y^2\left(\frac{D_{16}}{D_{11}}{m}_1^2+2\frac{D_{26}}{D_{11}}{m}_1^{}{m}_2^{}+\frac{D_{22}}{D_{11}}{m}_2^2\right),\end{array}$ (53)

$\begin{array}{l}{V}_n{}^{*}=-\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}\frac{d}{\alpha k}[\text{ci}\left(kz\right)\text{sin}\left(kz\right)-}\text{si}\left(kz\right)\text{cos}\left(kz\right)]\text{sign}\left(m_1{r}_{,x}+m_2{r}_{,y}\right)d\theta \\ -\frac{1}{4{\pi }^2}{\displaystyle \underset{0}{\overset{2\pi }{\int }}\frac{b}{\alpha }\frac{1}{\left(m_1{r}_{,x}+m_2{r}_{,y}\right)}}d\theta .\end{array}$ (54)

Like before $d=c+k^{2}b$, now $b$ and $c$ are described as

$\begin{array}{l}b=\left[n_x\left(1+n_y^2\right)+2\frac{D_{16}}{D_{11}}{n}_y^3-\frac{D_{12}}{D_{11}}{n}_x^{}{n}_y^2\right]m_1^3+\left[\frac{D_{22}}{D_{11}}{n}_y\left(1+n_x^2\right)+2\frac{D_{26}}{D_{11}}{n}_x^3-\frac{D_{12}}{D_{11}}{n}_x^2{n}_y^{}\right]m_2^3+\\ \left[4\frac{D_{16}}{D_{11}}{n}_x+\frac{D_{12}}{D_{11}}{n}_y\left(1+n_x^2\right)+4\frac{D_{66}}{D_{11}}{n}_y^3-n_x^2{n}_y^{}-2\frac{D_{26}}{D_{11}}{n}_x^{}{n}_y^2\right]m_1^2{m}_2^{}+\\ \left[4\frac{D_{26}}{D_{11}}{n}_y+\frac{D_{12}}{D_{11}}{n}_x\left(1+n_y^2\right)+4\frac{D_{66}}{D_{11}}{n}_x^3-\frac{D_{22}}{D_{11}}{n}_x^{}{n}_y^2-2\frac{D_{16}}{D_{11}}{n}_x^2{n}_y^{}\right]m_1^{}{m}_2^2,\end{array}$ (55)

$c=\left[\frac{N_x}{D_{11}}\left(n_x^3+n_x^{}{n}_y^2\right)+\frac{N_{xy}}{D_{11}}\left(n_x^2{n}_y^{}+n_y^3\right)\right]m_1+\left[\frac{N_{xy}}{D_{11}}\left(n_x^3+n_x^{}{n}_y^2\right)+\frac{N_x}{D_{11}}\left(n_x^2{n}_y^{}+n_y^3\right)\right]m_2.$ (56)

The other kernels can be found in Monteiro Jr. (2017).

5 Numerical Implementation

We used here quadratic isoparametric discontinuous elements with additional nodes on the corners. The matrix system is formed by the discretization of the plate's boundary in NE elements, given by

$Hu=Gq+p,$ (57)

with H representing the influence functions of the plate's forces, $V_n^{*}$, $M_n^{*}$, $R_{c_j}^{*}$, $\partial V_n^{*}/\partial n_i$, $\partial M_n^{*}/\partial n_i$, and $\partial R_{c_j}^{*}/\partial n_i$; G represents the influence functions of the plate's displacements, $w^{*}$, $\partial w^{*}/\partial n$, $w_{c_j}^{*}$, $\partial w_{}^{*}/\partial n_i$, ${\partial }^2{w}^{*}/\partial n\partial n_i$ and $\partial w_{c_j}^{*}/\partial n_i$; u is the vector representing the nodes’ displacements ($w^{}$,$\partial w/\partial n$, $w_{c_j}^{}$) along the boundary; q is the vector representing the nodes’ forces ($V_n$, $M_n$, $R_{c_j}^{}$,) along the boundary and p is the vector with the domain integrals terms of Eqs. (4) and (12).

The integral influence functions have been numerically integrated using a simple Gauss quadrature formula. Some of the influence functions associated with the fundamental solution present singularities. We have all types of singularities in the BIEs, weak (log(r)), strong (1/r) and hypersingular (1/r2). For the case of unidirectionally and multicompressed isotropic plates we were able to obtain analytical expression for the singular boundary elements. For the anisotropic case, singular boundary elements have been treated using special logarithmic Gauss quadrature and strong singular Gaussian quadrature formulae described in Kutt (1975a) and Kutt (1975b). The equations were implemented in a Matlab© code.

6 Numerical Results

In this section, some examples of polygonal plates with different boundary conditions and dimensions are presented. The validation of the proposed formulation was carried out with available analytical solutions and numerical results found in the literature.

The displacements and load factors have been defined as dimensionless quantities as $w_{nd}=wD/(q{a}^2)$ and $N_{nd}=N{a}^2/(4D)$ . The boundary conditions are arranged as Figure 4. For example, a SCFC plate means that edge 1 is simply supported, edges 2 and 4 are clamped, and edge 3 is free.

Figure 4
Plate’s boundary conditions

Two types of meshes used within this work are illustrated in Figure 5. The first one has 12 elements (3 elements per side), Figure 5a and the second has 20 elements (5 elements per side), Figure 5b.

Figure 5
Boundary Element Mesh

6.1 Isotropic simply supported multicompressed plate

The first problem analyzed was a simply supported square plate of side a = 1m and thickness h = 0.01m. The material properties are E = 200GPa and ν = 0.3. In order to verify the precision of the solution, the displacement of the plate’s center point has been compared with the analytical solution given by Timoshenko and Woinowsky-Krieger (1959). Applying $N_y=N_{cr}$ and$N_x=0.9N_{cr}$, the displacement using a 20 elements mesh is $w=1.1488797\times {10}^{-5}$m, while the analytical is $w=1.1295666\times {10}^{-5}$ m which represents an error of 1.71%. If a 28 elements mesh is used, the displacement is $w=1.1365988\times {10}^{-5}$m with an error of 0.62%. It can be seen that a good agreement is obtained with a 20 elements mesh and the error is reduced as the number of elements increase.

Figure 6 shows the displacement of the same plate with two meshes and a variable $N_x$. The asymptoptic behaviour of the plate’s displacement at the critical load can be clearly seen in Figure 6. Near the critical load the proposed solution has a poor performance with only 3 elements per side (12 elements mesh), but is in good agreement with the series solution of Timoshenko and Woinowsky-Krieger (1959) with the 20 elements mesh.

Figure 6
Critical load of a simply supported square plate $N_y=N_{cr}$ and $N_x$ variable.

6.2 Isotropic multicompressed plate - other boundary conditions

Liew and Wang (1995) used the Rayleigh-Ritz method to analyze the buckling of some multicompressed plates. They presented results for several cases of polygonal plates with constant linear distributed load perpendicular to the plate’s edges (see Figure 7). For the case of square plates, $N_x=N_y=N$.

Figure 7
Square plate compressed in x and y direction.

Table 1 presents the load factors obtained with the proposed solution with different boundary conditions. The following material properties are E = 200GPa, and ν = 0.3. The results are compared with the load factors of Liew and Wang (1995). With the 12 elements mesh most cases present a difference lower than 1%. The critical load in the formulation presented is calculated using the displacement. The case with only one edge clamped is the one that presents the most difficult convergence, but refining the mesh the difference decreases. Applying a 20 elements mesh the difference is approximately 2.29% and with a 28 elements mesh the difference is around 1.66%.

Table 1
Comparisons of buckling load factor

Boundary Condition	BEM Solution
	12 elements	Difference [%]
CCCC	13.088	0.10
CCSS	8.0106	0.2
FFFC	0.6089	3
FFCC	1.4368	0.65

6.3 Geometry influence on the critical buckling load on isotropic plates

A simply supported square plate has been bevelled to assess the influence of the geometry on the critical buckling load. The cases analyzed are illustrated in Figure 8.

Figure 8
Beveled plates

To analyze these cases we assume, $N_y=N_{xy}=0$ and $N_x$ varying from 50% to 150% of the critical buckling load. The material properties are the same as in the previous examples and a mesh with 3 elements per edge was used herein.

The effect of the bevels in the critical buckling load can be seen in Figure 9. As expected, the critical buckling load is increasing with the reduction of the square plate’s area when compared with the analytical solution of a simply supported square plate.

Figure 9
Critical buckling loads of the beveled isotropic plates

The critical load increase for case I is around 3.8%, which is less than the reduction of the area that is of 4.5%. Whilst cases II and III have the same area, both present different critical loads. Case II rises the critical load in approximately 6.4% and case III has an increase of approximately 8%, which is almost the reduced area for both cases, 9%. This difference is probably due to the symmetries present in case III. Case IV, similar to case I, presents an increase of approximately 12.6%, less than the area reduction (13.5%). The increase in the critical buckling load for case V (approximately 17.6%) is close to the area reduction (18%).

Refining the mesh, the critical load presents values of approximately 3.4% for case I (Figure 10) and 8% for case III (Figure 11), with a 7 elements mesh per side for both cases.

Figure 10
Convergence Study Case I

Figure 11
Convergence Study Case III

6.4 Compressive in-plane forces combined with shear in-plane forces for isotropic plates

In order to assess the effect of the shear in-plane forces applied together with compressive forces a clamped square plate has been analyzed. In this case $N_y=N_{cr}$, $N_{xy}=\gamma N_{cr}$ and $N_x$ variable, with $N_{cr}$ representing the critical buckling load of the bicompressed case. The plate’s dimension are a = b = 1 m, h = 0.01m and the material properties are E = 200GPa and $\nu =0.3$.

Figure 12 and Figure 13 show the effects of the shear in-plane force on the plate’s displacement field with $\gamma =\mathrm{0,0.3,0.5,0.8}$ and a 12 elements mesh. The plate’s deflection is increased with the presence of $N_{xy}$ and the displacement field is distorted (a rotated ellipsis). Figure 14 shows the effect in the critical buckling load. As expected, the critical load is reduced with the increase of $N_{xy}$. With $\gamma =0.3$, the reduction was approximately 3% and with $\gamma =0.5$ the critical load was reduced by approximately 8%.

Figure 12
Displacement isocurves of a clamped square plate.

Figure 13
Displacement fields of a clamped square plate.

Figure 14
Critical buckling load of a clamped square plate with combined compressive and shear in-plane forces.

6.5 Orthotropic plates

Displacement curves for the case where there are in-plane forces and out-plane forces acting on the plate are not available in the literature. Lekhnitskii (1968) developed a series solution for the case of orthotropic plates with a constant distributed out-of-plane force.

A simply supported square orthotropic plate of 1m dimension, thickness, h = 0.01m, material properties $E_l=200$GPa, $E_t=100$GPa, $G_{lt}=48$GPa, ${\nu }_{lt}=0.25$, and $q=1$N\m2 has been considered. The series solution at the center of the plate gives a displacement of $3.2624\times {10}^{-07}$m and the proposed BEM solution using a 20 elements mesh gives a displacement of $3.2736\times {10}^{-07}$m, which represents an error of 0.34%.

Lekhnitskii (1968) also developed solutions for the critical buckling load of some unidirectionally and bidirectionally compressed plates. Parameters for the three problems are the following: SSSS unidirectionally compressed and bicompressed. Figure 15 illustrates numerical BEM results for the buckling loads for three problems that have been analytically studied by Lekhnitskii. Using the 12 elements mesh for a square plate with 1m dimension, h = 0.01m, $E_l=20$GPa, $E_t=100$GPa, $G_{lt}=\mathrm{38,46}$GPa, ${\nu }_{lt}=0.3$ and a perturbation in the form of a concentrated force of $q=1$N.

Figure 15
Critical buckling loads for orthotropic plates

6.6 Symmetric laminated plates

Shojaee et al. (2012) applied the FEM to analyze several symmetric laminates compressed in the x direction. The material properties and geometry are: a = b =10m, h = 0.06m, El/Et = 2.45, Glt/Et = 0.48, ${\nu }_{lt}=0.23$. The buckling load factor is defined as $\beta =N{b}^2/({\pi }^2{D}_{\mathrm{0,1}})$, with $D_{\mathrm{0,1}}$ defined as the flexural rigidity along the fiber direction of a plate with a 0° angle ply. Table 2 shows the buckling load factors with a 12 BEM elements mesh for different configurations of ply orientation angles. The BEM proposed solution shows good agreement with FEM results.

Table 2
Buckling load factor of laminates

Laminated	Method	Boundary Conditions
		SSSS	CCCC	CSCS	SCCS
[0°,0°,0°]	Shojaee et al. (2012)	2.36	6.71	4.27	3.93
	BEM	2.36	6.727	4.27	3.93
[15°,-15°,15°]	Shojaee et al. (2012)	2.43	6.57	4.41	3.96
	BEM	2.426	6.57	4.399	3.94
[30°,-30°,30°]	Shojaee et al. (2012)	2.59	6.29	4.79	4.00
	BEM	2.577	6.258	4.742	3.98
[45°,-45°,45°]	Shojaee et al. (2012)	2.66	6.10	4.78	4.01
	BEM	2.646	6.10	4.828	3.989
[90°,0°,90°]	Shojaee et al. (2012)	2.36	6.70	4.36	3.93
	BEM	2.36	6.70	4.36	3.93

6.7 Influence of the angle ply on the critical load of an orthotropic plate

Figure 16 illustrates the critical load of a simply supported square plate with 1x1m, h = 0.01m, E_l = 200GPa, E_t = 100GPa, G_lt = 20GPa, ν_lt = 0.25 for different angle plies. The greatest increase in the critical load has been found for a 45° angle ply in the present example, around 35% increase. We see that the fibers' angle ply has an important effect on the critical load. This effect is highlighted in Figure 17, which illustrates the displacement behavior of the plate with a 0° ply-orientation angle compared to the plates with 30° and 60° ply-orientation angles. The here developed fundamental solution can be used as a tool for optimizing the stability project of laminated polygonal plates.

Figure 16
Critical buckling load of a simply supported square plate with different angle plies.

Figure 17
Critical buckling load of a simply supported square plate with ply-orientation angles of 0°, 30° and 60°.

6.8 Geometry influence on the critical buckling load of orthotropic plates

The same geometries of the isotropic case were studied for an orthotropic plate. The properties for this case are $E_l=200$. GPa, $E_t=20$GPa, $G_{lt}=7.69$GPa${\nu }_{lt}=0.3$. Figure 18 illustrates the plate’s displacement. The same effect of the isotropic case is observed in the orthotropic case. However, the increase in the critical buckling load for the orthotropic plate are greater than for the isotropic case. The orthotropic critical buckling loads have been increased approximately by 5%, 15%, 16%, 20%, and 31% for each case respectively with 3 elements per edge. This difference is probably due to the differences between the elastic modulus along the fiber and the elastic modulus transversal to the fiber. Refining the mesh, case I tends to decrease the critical load by approximately 3.5% (Figure 19) whilst case III decrease to approximately 15% (Figure 20).

Figure 18
Critical buckling load for beveled orthotropic plates

Figure 19
Case I convergence

Figure 20
Case III convergence

7 Conclusions

We have developed new fundamental solutions using the Radon transform to solve the buckling problem for isotropic and anisotropic plates under the combined action of in-plane compressive forces/shear forces and out-of-plane loads. Several numerical examples were shown for isotropic, orthotropic, symmetric laminated and monoclinic materials with different boundary conditions and geometries. The results obtained were compared with available analytical solutions and numerical results in the literature. The present formulation can be used without domain integrals for the case of concentrated loads. As long as $\beta >0$, this BEM solution can be used with the combined action of compressive and shear in-plane forces. Another interesting feature of the present formulation is the numerical evaluation of the actual displacement plate’s field $w(x,y)$ near the critical load. With the true displacement field it is possible to observe the numerical behavior of the plate’s displacement surface versus analytical expressions of $w(x,y)$, which are obtained with a great number of elements in a sum series. Hence, we can assess BEM’s numerical precision for the whole plate’s displacement surface $w(x,y)$. The formulation presented here is general, hence the plate’s displacement behavior can be observed for various kinds of loads, not just the point force used in this work. The obtained BEM numerical results show excellent agreement when compared with both analytical and previously published FEM results. The new fundamental solutions can be used for optimizing the stability project of polygonal, laminated anisotropic plates.

References

• Albuquerque, E. L., Sollero, P., Venturini, W. S., Aliabadi, M. H. (2006). Boundary element analysis of anisotropic Kirchhoff plates. International Journal of Solids and Structures. 43:4029-4046.
• Altiero, N. J., Sikarskie, D. L. (1978). A boundary integral method applied to plates of arbitrary plan form. Comp. & Struct. 9:163-168.
• Baiz, P. M., Aliabadi, M. H. (2007). Buckling analysis of shear deformable shallows shells by the boundary element method. Engineering Analysis with Boundary Elements. 31:361-372.
• Bézine, G. P. (1978). Boundary integral formulation for plate flexure with arbitrary boundary conditions. Mech. Res. Comm. 5:197-206.
• Christensen, R. M. (1979). Mechanics of Composite Materials. Dover Publications, Inc.
• Dirgantara, T., Aliabadi, M. H. (1999). A new boundary element formulation for shear deformable shells analysis. International Journal for Numerical Methods in Engineering. 45:1257-1275.
• Dirgantara, T., Aliabadi, M. H. (2006). A boundary element formulation for geometrically nonlinear analysis of shear deformable shells. Computer Methods in Applied Mechanics and Engineering. 195:4635-4654.
• Ferreira, A. J. M., Castro, L. M., Roque, C. M. C., Reddy, J. N., Bertoluzza, S. (2011). Buckling Analysis of laminated plates by wavelets. Computers and Structures. 89:626-630.
• Helgason, S. (1980). The Radon Transform. Birkhäuser.
• Herakovich, C. T. (1998). Mechanics of Fibrous Composites. John Wiley & Sons, Inc.
• Jahanshahi, A., Dundurs, J. (1964). Concentrated forces on compressed plates and related problems for moving loads. Journal of Applied Mechanics. 31:83-87.
• Katsikadelis, J. T., Babouskos, N. G. (2009). Nonlinear Flutter Instability of Thin Damped Plates. Journal of Mechanics Materials and Structures. 4 (7-8):1395-1414.
• Kutt, H. R. (1975a). Quadrature Formulae for Finite Part Integrals. The National Research Institute for Mathematical Sciences. Report WISK 178.
• Kutt, H. R. (1975b). On the Numerical Evaluation of Finite Part Integrals Involving an Algebraic Singularity. The National Research Institute for Mathematical Sciences. Report WISK 179.
• Lekhnitskii, S. G. (1968). Anisotropic Plates. Gordon Breach Science Publishers.
• Liew, K. M., Wang, C. M. (1995). Elastic buckling of regular polygonal plates. Thin-Walled Structures. 21:163-173.
• Ludwig, D. (1966). The Radon transform on Euclidian space*. Communications on Pure and Applied Mathematics. 19:49-81.
• Manolis, G. D., Rangelov, T. V., Shaw, R. P. (2003). The non-homogeneous biharmonic plate equation: Fundamental solutions. International Journal of Solids and Structures. 40:5753-5767.
• Monteiro Jr., J. (2017). Static Stability Analysis of Laminated Anisotropic Kirchhoff’s Plates via the Boundary Element Method. MS. Dissertation (in Portuguese), University of Campinas, Brazil.
• Monteiro Jr., J., Daros, C. H. (2017). Buckling Analysis of Thin Kirchhoff Plates Unidirectionally Compressed via the Boundary Element Method. In: Proceedings of International Symposium on Solids Mechanics - MecSol, 6., Joinville - Brazil: ABCM.
• Paiva, J. B., Venturini, W. S. (1992). Alternative Technique for the Solution of Plate Bending Problems Using the Boundary Element Method. Advances in Engineering Software. 14:265-271.
• Rangelov, T. V., Manolis, G. D., Dineva, P. S. (2005). Elastodynamic Fundamental Solutions for Certain Families of 2d Inhonogeneous Anisotropic Doamians: Basic Derivations. European Journal of Mechanics A/Solids. 24:820-836.
• Shi, G (1990). Flexural vibration and buckling analysis of orthotropic plates by the boundary element method. International Journal of Solids Structures. 26:1351-1370.
• Shi, G., Bézine, G. P. (1990). Buckling analysis of orthotropic plates by boundary element method. Mechanics Research Communications. 17:1-8.
• Shojaee, S., Valizadeh, N., Izadpanah, E., Bui, T., Vu, T. V. (2012). Free vibration and buckling analysis of laminated composite plates using NURBS-based isogeometric finit element method. Composite Strucutres. 94:1677-1693.
• Stern, M. A. (1979). A general boundary integral formulation for the numerical solution of plate bending problems. Int. J. Solids Struct. 15:769-782.
• Syngellakis, S., Elzein, A. (1994). Plate Buckling Loads by The Boundary Element Method. International Journal for Numerical Methods in Engineering. 37:1763-1778.
• Timoshenko, S. P., Woinowsky-Krieger, S. (1959). Theory of Plates and Shells. McGraw-Hill.

• Available Online: May 08, 2018

Appendix A - D Expressions.

$\begin{array}{l}{D}_1=-\frac{1}{8}\left\{\left(\nabla r\cdot n\right)\left[2+2n_y^2\left(1-\nu \right)\right]-\left(1-\nu \right)r_{,y}{n}_y\right\}\left(\nabla r\cdot n_{ì}\right)\\ -\frac{1}{8}\left\{2n_x+(1-\nu )\left[r_{,y}{}^2{n}_x{}^3+n_x{n}_y^2(2-3r_{,y}{}^2)-3r_{,x}{r}_{,y}{n}_x^2{n}_y+r_{,x}{r}_{,y}{n}_y{}^3-n_x^3\frac{N_x}{D{\kappa }^2}\right]\right\}{n}_{xi}\end{array}$ (A1)

$\begin{array}{l}{D}_2=\frac{1}{8}\left\{\left(-r_{,y}{}^2{n}_x+r_{,x}{r}_{,y}{n}_y\right)\left[2+2n_y^2\left(1-\nu \right)\right]-\left(1-\nu \right)r_{,x}{r}_{,y}{n}_y\right\}{n}_{xi}\\ +\frac{1}{8}\left\{-\left(1-\nu \right)\left[\left(r_{,x}{}^2-r_{,y}{}^2\right)\left(n_x^3-3n_x{n}_y^2\right)+2r_{,x}{r}_{,y}{n}_y\left(3n_x^2-n_y^2\right)\right]\right\}{n}_{xi}\\ +\frac{1}{8}\left\{\left(-r_{,x}{}^2{n}_y+r_{,x}{r}_{,y}{n}_x\right)\left[2+2n_y^2\left(1-\nu \right)\right]+\left(1-\nu \right)r_{,x}{}^2{n}_y\right\}{n}_{yi}\\ +\frac{1}{8}\left\{\left(\nabla r\cdot n\right)\left[2+2n_y^2\left(1-\nu \right)\right]-\left(1-\nu \right)r_{,y}{n}_y\right\}\left(\nabla r\cdot n_{ì}\right)\end{array}$ (A2)

$\begin{array}{l}{D}_3=\frac{(1-\nu )}{8}\{-\left[3n_x{n}_y{}^2\left(2r_{,y}{}^2-1\right)+2r_{,x}{r}_{,y}{n}_y\left(3n_x{}^2-n_y{}^2\right)+n_x{}^3\left(1-2r_{,y}{}^2\right)\right]\left(\nabla r\cdot n_{ì}\right)\\ +\left[2r_{,y}{}^2{r}_{,x}{n}_x\left(n_x{}^2-3n_y{}^2\right)+\left(r_{,y}{}^3-r_{,x}{}^2{r}_{,y}\right)\left(3n_x{}^2-n_y{}^2\right)n_y\right]n_{xi}\\ +\left[-2r_{,x}{}^2{r}_{,y}{n}_x\left(n_x{}^2-3n_y{}^2\right)+\left(r_{,x}{}^3-r_{,y}{}^2{r}_{,x}\right)\left(3n_x{}^2-n_y{}^2\right)n_y\right]n_{yi}\}\end{array}$ (A3)

$\begin{array}{l}{D}_4=\frac{1}{8}\left\{\left(\nabla r\cdot n\right)\left[2+2n_y^2\left(1-\nu \right)\right]-\left(1-\nu \right)r_{,y}{n}_y\right\}{n}_{xi}\\ +\frac{1}{8}\left\{2n_x+(1-\nu )\left[r_{,y}{}^2{n}_x{}^3+n_x{n}_y^2(2-3r_{,y}{}^2)-3r_{,x}{r}_{,y}{n}_x^2{n}_y+r_{,x}{r}_{,y}{n}_y{}^3-n_x^3\frac{N_x}{D{\kappa }^2}\right]\right\}\left(\nabla r\cdot n_{ì}\right)\end{array}$ (A4)

$\begin{array}{l}{D}_5=\frac{(1-\nu )}{8}\{2\left[3n_x{n}_y{}^2\left(2r_{,y}{}^2-1\right)+2r_{,x}{r}_{,y}{n}_y\left(3n_x{}^2-n_y{}^2\right)+n_x{}^3\left(1-2r_{,y}{}^2\right)\right]\left(\nabla r\cdot n_{ì}\right)\\ +\left[r_{,y}{}^2{r}_{,x}{n}_x\left(12n_y{}^2-4n_x{}^2\right)-2\left(r_{,y}{}^3-r_{,x}{}^2{r}_{,y}\right)\left(3n_x{}^2-n_y{}^2\right)n_y\right]n_{xi}\\ +\left[-r_{,x}{}^2{r}_{,y}{n}_x\left(12n_y{}^2-4n_x{}^2\right)-2\left(r_{,x}{}^3-r_{,y}{}^2{r}_{,x}\right)\left(3n_x{}^2-n_y{}^2\right)n_y\right]n_{yi}\}\end{array}$ (A5)

Publication Dates

• Publication in this collection
  2018

History

• Received
  31 July 2017
• Reviewed
  03 May 2018
• Accepted
  07 May 2018