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

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.


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 andAliabadi (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.

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) Figure 1: Compressed plate.
With ij D 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.
With k z representing the distance from the middle plane to the ply's surfaces, k and k-1 (see Figure 2).( ) is the reduced rotated matrix in the plane stress state with the elastic constants of each ply (Eq.( 3)).In Eq. ( 2) , 1, 2, 6 i j  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) With l E representing the elastic modulus along the direction of the fiber, t E representing the elastic modulus in the transversal direction to the fiber axis, lt  representing the Poisson ratio, lt G representing the shear modulus and θ k being the angle ply orientation.Eq. ( 3) was obtained based on Christensen (1979) and Herakovich (1998).
With Γ and Ω representing the plate's boundary and domain, respectively; , i x x represent the field and source point, respectively; ci is a free term depending on the smoothness of the plate's boundary; n M represents the resultant moment (see Eq. ( 5)).
R is the corner reaction, described as M  are the twisting moments after (+) and before (-) the corner (see Eq. ( 6)); j c w is the corner displacement and Nc is the total number of corners.n V represents here the resultant shear force for polygonal plates, which considers the effect of the outof-plane and in-plane forces (see Eq. ( 8)).The starred terms are the known fundamental Kernels, which are obtained deriving the fundamental solution.
With x n and y n denoting the direction cosines.
∂/∂n and ∂/∂s are the derivatives in the normal and tangential directions, respectively, given by n n , n n .
x y y x x y x y It is well known that due to the number of variables ( n M , n V , w , / w 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, i n , normal to the boundary surface at the source point.This is a hypersingular BIE since there are kernels, 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 * ( , ) 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).

Unidirectionally compressed plates
For this case, the governing equation is 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 Where ( , ) x y and ( , ) i i x y are the coordinates of the field and source points, respectively.D is the flexural rigidity, given by   ; 0 Y and 1 Y are the Bessel functions of the second kind of order 0 and 1, respectively.Notice that * w and * / w  n are not explicit, requiring further numerical integration.However, all other kernels are explicit.Deriving Eq. ( 14), we find .Using the following useful identity we may express the third order derivatives as 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,  n which has been derived, to the author's best knowledge, here for the first time.
Expressions for D1 to D5 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.

Multidirectionally compressed plates
The governing equation in this case is 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 Where s   x m is the transform hyperplane and m is the unit radius of the circle (Figure 3).Using the Radon Transform we find With ξ denoting the source point and ℒ represents here a linear partial differential operator.Applying the Radon Transform to Eq. ( 19), we find  24) as a solution for Eq. ( 22), The inverse Radon transform is defined via the equation Applying the Inverse Radon Transform in Eq. ( 24), we find

 
With  representing the gamma-Euler constant.Applying Eq. ( 26) on Eq. ( 24), we have 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 Substituting Eq. ( 29) in Eq. ( 19), we find 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.
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 , the singular part can be rewritten as The integral 1 I can be analytically evaluated and is given by 1 (1 ) log( ) / (4 ) . Integral 2 I can also be obtained in a closed form as So after treating the singularity, the resultant moment can be expressed as Likewise, the total shear force is given as The kernel representing the corner forces are given as (1 ) .
x y x y Doing a similar approach to the resultant moment, we obtain The other Kernels for the hypersingular BIE can be concisely expressed as Notice that Eq. ( 47) is hypersingular in r for straight boundary elements and must be treated in the sense of Hadamard.

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 (51) Some of the other kernels associated with the fundamental solution are given on following, , now b and c are described as The other kernels can be found in Monteiro Jr. ( 2017).

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 with H representing the influence functions of the plate's forces,  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/r 2 ).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.

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 .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.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.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 x N .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.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%.

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.To analyze these cases we assume, 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.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 14: Critical buckling load of a clamped square plate with combined compressive and shear in-plane forces.

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 , and 1 q  N\m 2 has been considered.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,    16 illustrates the critical load of a simply supported square plate with 1x1m, h = 0.01m, El = 200GPa, Et = 100GPa, Glt = 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 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 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).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 outof-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 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.

Figure 2 :
Figure 2: A symmetric laminated composite with N plies.

,
si is the sine integral function and ci is the cosine integral function, defined in Eqs.(27.1) and (27.2).
functions of the plate's displacements, * w , along the boundary and p is the vector with the domain integrals terms of Eqs. (

Figure 6 :
Figure 6: Critical load of a simply supported square plate y

Figure 7 :
Figure 7: Square plate compressed in x and y direction.

N
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.

Figure 9 :
Figure 9: Critical buckling loads of the beveled isotropic plates

Figure 10 :
Figure 10: Convergence Study Case I

Figure 12 and
Figure12and Figure13show the effects of the shear in-plane force on the plate's displacement field with 0,0.3,0.5, 0.8 and a 12 elements mesh.The plate's deflection is increased with the presence of xy N and the displacement field is distorted (a rotated ellipsis).Figure14shows the effect in the critical buckling load.As expected, the critical load is reduced with the increase of x y N .With 0.3   , the reduction was approximately 3%

Figure 13 :
Figure 13: Displacement fields of a clamped square plate.

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

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

Figure
Figure 20: Case III convergence Buckling Analysis of Laminated Anisotropic Kirchhoff's Plates via The Boundary Element Method Latin American Journal of Solids and Structures, 2018, 15(10 Thematic Section), e793/23 2 ) singularities.The latter require special analytical or numerical treatment since they are hypersingular in the sense of Hadamard.
Buckling Analysis of Laminated Anisotropic Kirchhoff's Plates via The Boundary Element Method Latin American Journal of Solids and Structures, 2018, 15(10 Thematic Section), e79 5/23 which present (1/r Buckling Analysis of Laminated Anisotropic Kirchhoff's Plates via The Boundary Element Method Latin American Journal of Solids and Structures, 2018, 15(10 Thematic Section), Buckling Analysis of Laminated Anisotropic Kirchhoff's Plates via The Boundary Element Method Latin American Journal of Solids and Structures, 2018, 15(10 Thematic Section), e79 7/23 Buckling Analysis of Laminated Anisotropic Kirchhoff's Plates via The Boundary Element Method Latin American Journal of Solids and Structures, 2018, 15(10 Thematic Section), e79 8/23 Buckling Analysis of Laminated Anisotropic Kirchhoff's Plates via The Boundary Element Method Latin American Journal of Solids and Structures, 2018, 15(10 Thematic Section), e79 9/23 46)Buckling Analysis of Laminated Anisotropic Kirchhoff's Plates via The Boundary Element Method Latin American Journal of Solids and Structures, 2018, 15(10 Thematic Section), e79 10/23 Buckling Analysis of Laminated Anisotropic Kirchhoff's Plates via The Boundary Element Method Latin American Journal of Solids and Structures, 2018, 15(10 Thematic Section), e7911/23

Table 1 :
Comparisons of buckling load factor The series solution at the center of the plate gives a displacement of Analysis of Laminated Anisotropic Kirchhoff's Plates via The Boundary Element Method m, which represents an error of 0.34%.Buckling Latin American Journal of Solids and Structures, 2018, 15(10 Thematic Section), e79 17/23

Table 2 :
Table2shows 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.Buckling load factor of laminates Influence of the angle ply on the critical load of an orthotropic plate Figure