Bending of elastic plates with micro-voids

Accepted Manuscript Bending of elastic plates with micro-voids M. Repka, V. Sladek, J. Sladek PII: DOI: Reference:

S0263-8223(18)30169-7 https://doi.org/10.1016/j.compstruct.2018.05.072 COST 9711

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

12 January 2018 7 May 2018 17 May 2018

Please cite this article as: Repka, M., Sladek, V., Sladek, J., Bending of elastic plates with micro-voids, Composite Structures (2018), doi: https://doi.org/10.1016/j.compstruct.2018.05.072

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Bending of elastic plates with micro-voids M. Repka1*, V. Sladek1, J. Sladek1 1

Institute of Construction and Architecture, Slovak Academy of Sciences, 84503 Bratislava, Slovakia

Abstract: In this paper, we present a systematic development of unified formulation for bending of elastic plates with micro-voids. The homogenization concept is applied with replacing real medium with micro-voids by a homogeneous medium with effective material coefficients. The effective material coefficients are calculated in advance by solving properly selected micro-structural 3D boundary value problems in representative volume elements with specified shape and volume content of micro-voids. The unified formulation for plate bending problems is developed from 3D elasticity with taking into account the assumptions of three plate bending theories such as Kirchhoff-Love theory for bending of thin plates, 1st and 3rd order shear deformation theories for thick plates. Switching among these theories is allowed by proper selection of two key factors introduced in the unified formulation. The influence of the porosity on bending of porous plates subjected to transversal loading is studied via the numerical simulations.

Keyword: empty voids, homogenization, RVE, anisotropy, plate bending, thin/thick plate

1. Introduction It is common practice in engineering and design process to utilize structural elements like beams shells and plates. The latter is a structural element with zero curvature (straight plate) and small thickness as compared with other dimensions [1]. The governing equations of the plate within plate theories can be derived from 3-D formulation assumed in semi-integral form across the plate thickness and due to this assumption the original 3-D problem is reduces to the 2-D problem. Since the end of 19th century several plate theories have been developed, for instance, the Kirchhoff-Love theory (KLT) which is the extension of the Euler-Bernoulli beam theory. The main shortcoming of this theory is the omission of the shear deformations which play an *Corresponding author, Tel.: +421 910 920 907 E-mail address: [email protected]

important role in bending of thick plates. In order to take in to account shear deformations, the generalized shear deformation theories have been proposed including the first order shear deformation theory (FSDT) [1] and higher-order shear deformation theory (HSDT) [2-5] . During recent decades many authors have contributed to developing and studying these theories (see e.g [6-11]), In the FSDT, the shear strains are considered as constants across the plate thickness and the theory requires the shear correction factor to compute transverse shear forces, while in the HSDT no shear correction factors are needed and the zero tangential traction boundary conditions are satisfied on the bottom and top surfaces of the plate. During the past decades the advances in material science and engineering technology have enabled to manufacture more complex and sophisticated materials which are often nonhomogeneous even because of technological microstructural inhomogeneity like pores or voids. The voids and pores in structures affect their mechanical behaviour, therefore it is necessary to capture these effects in engineering design. The design is mostly the matter of numerical simulations and the detailed modelling of micro-structural inhomogeneity can be computationally expensive, time consuming and practically inconvenient for description the behaviour of macrostructures. One of the possible ways how to circumstance this problem is to utilize a concept of homogenization [12] which can be applicable under certain conditions. Usually it is possible, if the material defect is much smaller than the characteristic length in structural design. Although the concept of homogenization is classical it is still an area of open research [13]. The idea of the homogenization is applicable also to microscopically inhomogeneous materials as long as these materials are macroscopically or statistically homogeneous in considered macrostructure

[14], [15]. When the concept of homogenization is employed in analysis of

bending problems for plates with microstructure, it is necessary to evaluate the effective coefficients by solving 3D boundary value problems in representative volume elements, since the 2D formulation for plate bending is derived from the 3D elasticity formulation under certain plate bending assumptions. Also Dong et al [16] utilized a 3D finite element discretization to extract the homogenized material properties from unit cell analysis, which were then used to compute the global homogenized plate response.

The main objective of this paper is to study a bending response of thin and/or thick plates with empty micro-voids within three plate bending theories (KLT, FSDT, TSDT). The material inhomogeneity is taken into account by using the concept of homogenization in which the effective material parameters are calculated with modelling 3D representative volume element (RVE). Note that it is insufficient to calculate the effective material coefficients from 2D representative volume element and replace the material coefficients in classical plate bending formulation by such effective values. The correct 2D formulation for plate bending must be derived from the 3D elasticity of homogenized medium with effective material elasticity coefficients. For this purpose, several properly selected microstructural 3D boundary value problems have to be solved in advance. To this end the finite element method has been applied [17], [18]. In general, the matrix of effective material coefficients is fully occupied matrix corresponding to anisotropic elastic medium. In the developed of unified theory for bending of elastic plates with micro-voids, both the governing equations and boundary conditions are derived with using the variational principle. The influence of porosity on bending response is investigated in numerical simulations for a square plate with empty micro-voids, clamped edges and subjected to uniform transversal loading. The paper is organized as follows: section 2 is devoted to calculation of the effective elastic material coefficients, while section 3 details the derivation of the unified plate bending formulation in case of general elastic continua. In Section 4, the numerical examples and numerical results are presented and discussed. Conclusions are summarized in section 5.

2. Evaluation of effective material coefficients In order to eliminate the necessity to model micro-voids and facilitate the numerical study of macrostructures with micro-voids, the homogenization concept is acceptable under certain conditions [14]. Note that the material properties of effective elastic medium can differ from those of the elastic skeleton. Since the 2D plate bending formulations are derived from 3D elasticity, we need to know and calculate the effective material coefficients in 3D elasticity in advance. Around each pore a cube cell is considered and the representative volume element (RVE) can be composed of one or more cells (Fig. 1).

Fig.1 Scheme of the RVE composed of one cell Note that B is regular sub-region of the macroscopically homogeneous body (homogenized body) but also occupied by a RVE of the microscopically inhomogeneous body composed of the homogeneous skeleton   B and empty voids B   with B   ,     0 and  0 being the boundary of the RVE, skeleton and voids in the RVE, respectively. The volume average or the mean value of the field variable f is defined as f

B



1  f ( x) d  . B B

The constitutive relationships in the linear elastic homogenized continuum are given as eff eff  ijeff (x)  cijkl  kl (x) ,

eff  const . Hence,  ijeff where cijkl

B

eff  cijkl  kleff

B

.

However, these relationships cannot be used for evaluation of the stiffness coefficients in the homogenized effective continuum, since the mean values of the strains and stresses in effective continuum are not known (we cannot solve any direct boundary value problem in effective medium because of absence of material coefficients). Therefore, we should find a physically meaningful correlation between   ijeff 

B

eff ,  kl

 and   ms ,  ms  , because the latter kl  ij B B 

B 

can be found from the solution of certain micro-structural boundary value problems (BVP) in the RVE [14]. For this purpose, there are crucial the Hill’s theorem [12] and certain approaches, e.g. such as [15]: (i) average stress equivalence (ASE), (ii) boundary density equivalence (BDE).

Specific boundary conditions are required in order to satisfy the Hill’s conditions, e.g. the Kinematic Uniform Boundary Conditions (KUBC). Let us consider the KUBC for solution of microstructural BVP on the outer boundary of the RVE with micro-voids in elastic skeleton, while the surface of empty micro-voids is assumed to be traction free uims   ik xk





, ti

0

  ikms nk

0

.

Assuming the same boundary conditions on the outer boundary of the RVE also for the effective homogenized medium, we can find the average value of effective strains

 ijeff 

B

 lk

 

2B

 u

1 1  ijeff d    B B 2B



eff i, j

 u eff j ,i  d  

B

x n j   lj xk ni  d  

li k



 lk 2B

 

x

li k , j

1 2B

 u

eff i

n j  u eff j ni  d  



  lj xk ,i  d  

 ij   ji

B

2

(1)

  ij

The average value of the density of deformation energy in effective medium is calculated as

2U   ijeff  ijeff 

B

1 1 1  ijeff  ijeff d     ijeff uieff, j d    B B B B B



 

eff ij

uieff  d   ,j

B

   1 n j ijeff uieff d   ik  xk n j ijeff d   ik  xk , j ijeff d   ij   ijeff d    B  B  B B B B

  ij  ijeff

B

  ijeff

B

 ijeff

B

eff  cijkl  kleff

B

 ijeff

(2)

B

Let us assume this energy to be the same as the average deformation energy resulting from the microstructural solution in the RVE, i.e.

2U   ijms ijms  

1 B



B



1 1 1  ijms ijms d     ijms uims, j d    B B B B B

n j ijms uims d  

0

 ik

x B 

k, j

 ijms d  

B

 ij B

 



ms ms ij i ,j

u

d 

B

 1 timsuims d   ik  xk tims d    2B  B 



ms ij

d    ij  ijms

(3)

B

B

hence,

 ijms

B



1 B

x t

ms j i



d 

1 2B

x t

ms j i

ms  xi t ms j  d  : Tij .



Using the Voigt notation, the constitutive law in linear effective medium takes the form

(4)

  1eff  eff  2   3eff  eff  4   eff  5eff   6

   11eff   eff    22    33eff  :  eff    23    eff   13eff     12

  c11eff   eff   c21 eff   c31    eff   c41   c eff   51   c eff   61

c12eff eff c22 eff c32 eff c42 eff c52 eff c62

c13eff eff c23 eff c33 eff c43 eff c53 eff c63

c14eff eff c24 eff c34 eff c44 eff c54 eff c64

c15eff eff c25 eff c35 eff c45 eff c55 eff c65

c16eff eff c26 eff c36 eff c46 eff c56 eff c66

  1eff   eff  2    3eff    eff  4    eff   5eff     6

  1eff1   eff    22    33eff  :  eff   2 23   2 eff   13eff   2   12

    .     

(5)

Recall that in the plate bending theories  33eff  0 and  33eff does not enter the formulation. Thus, we could confine to the constitutive law in form   1eff  eff  2   4eff  eff  5   eff  6

   11eff   eff    22 eff  :   23   eff    13    eff   12

  c11eff   eff   c21 eff    c41   eff   c51   c eff   61

c12eff eff c22 eff c42 eff c52 eff c62

c14eff eff c24 eff c44 eff c54 eff c64

c15eff eff c25 eff c45 c5eff5 eff c65

c16eff eff c26 eff c46 eff c56 eff c66

  1eff   eff  2     4eff   eff   5    eff   6

  11eff   eff    22 eff  :  2 23   eff   213   2 eff   12

   .    

(6)

Now, we shall specify 5 particular microstructural BVP in the RVE, in order to calculate the effective material coefficients occurring in Eq. (6). BVP (1):

 ik(1)   i1 k111 , 11   (1)  const

eff eff [(3), (4)  2U (1)  11T11ms (1) ; (2), (1)  2U (1)  c1111  T11ms (1) / 11 1111 ]  c11eff  c1111

(6)  caeff1   aeff

(1) B

/ 11 for a  2, 4,5,6

Question is, how to determine  aeff

(1) B

? We can use two alternative approaches: (1)

 aeff

(i) ASE (average stress equivalence) approach:

B

  ams

(1) B

, hence and from (3), (4),

we have caeff1  Tams (1) / 11 (ii) BDE (boundary density equivalence) approach: uieff   ijeff

(1) B



1 2B

 

   kjeff (1) ki  d  

eff (1) ik kj

B

1 2B



 uims

 

eff (1) ik



 tieff



 tims



x j ,k   kjeff (1) xi ,k  d  

B



1  eff (1) 1  ik x j ,k   kjeff (1) xi ,k  d   nk  ikeff (1) x j   kjeff (1) xi  d      2 B B 2 B 



1 2B

 t

eff (1) i



(1) x j  t eff xi  d   j

1 2B

 t

ms (1) i



(1) x j  t ms xi  d   Tijms (1) j

(1)

Thus, caeff1   aeff

B

/ 11  Tams (1) / 11

 22   (2)  const

BVP (2):  ik(2)   i 2 k 2 22 ,

eff eff eff [(3), (4)  2U (2)   22T22ms (2) ; (2), (1)  2U (2)  c2222  22 22 ]  c22  c2222  T22ms (2) /  22

(6)  caeff2   aeff (2)

(i) ASE:  aeff (ii) BDE: uieff

B

(2)

/  22 for a  1, 4,5,6

B

  ams

 uims





(2) B

 Tams (2)  caeff2  Tams (2) /  22

 tieff



 tims

BVP (4):  ik(4)   i 2 k 3   i 3 k 2   23 ,



  aeff

(2) B

 Tams (2)  caeff2  Tams (2) /  22

2 23   (4)  const

eff [(3), (4)  2U (4)   23 T23ms (4)  T32ms (4)   2 23T23ms (4) ; (2), (1)  2U (4)  4c44  23 23 ] 

eff c44  T23ms (4) / 2 23 (4)

(6)  caeff4   aeff (4)

(i) ASE:  aeff (ii) BDE: uieff

B

B



/ 2 23 for a  1, 2,5,6 (4)

  ams

 uims



B

 Tams (4)  caeff4  Tams (4) / 2 23

 tieff



BVP (5):  ik(5)   i1 k 3   i 3 k1  13 ,

 tims



  aeff

(4) B

 Tams (4)  caeff4  Tams (4) / 2 23

213   (5)  const

eff [(3), (4)  2U (5)  13 T13ms (5)  T31ms (5)   213T13ms (5) ; (2), (1)  2U (5)  4c55 1313 ] 

eff c55  T13ms (5) / 213

(6)  caeff5   aeff (5)

(i) ASE:  aeff (ii) BDE: uieff

B



(5) B

/ 213 for a  1, 2, 4,6

  ams

 uims



(5) B

 Tams (5)  caeff5  Tams (5) / 213

 tieff



 tims



  aeff

(5) B

 Tams (5)  caeff5  Tams (5) / 213

BVP (6):  ik(6)   i1 k 2   i 2 k1  12 , 212   (6)  const

eff [(3), (4)  2U (6)  12 T12ms (6)  T21ms (6)   212T12ms (6) ; (2), (1)  2U (6)  4c66 1212 ] 

eff c66  T12ms (6) / 212 (6)

(6)  caeff6   aeff (6)

(i) ASE:  aeff (ii) BDE: uieff

B

B



/ 212 for a  1, 2, 4,5

  ams

 uims



(6) B

 Tams (6)  caeff6  Tams (6) / 212

 tieff



 tims



  aeff

(6) B

 Tams (6)  caeff6  Tams (6) / 212

One can see that both the ASE and BDE approaches yield the same expression for the effective material coefficients. Furthermore, the formula for effective material coefficients is formally the same for all considered BVP (b) and given as eff cab  Tams (b ) /  (b ) a, b 1, 2, 4, 5, 6

(7)

with Tams (b ) being given by Eq. (4), and for evaluation of Tams (b ) we need to solve particular microstructural BVP ( b 1, 2, 4, 5, 6 ) in the RVE. Let us consider periodically distributed empty pores in a macrostructure. Around each pore, we consider a cube cell. There is a question if the geometrical periodicity combined with above mentioned KUBC applied to the RVE with more pores leads to periodically symmetric solution in the RVE. For this purpose, let us consider the RVE shown in Fig.2.

Fig. 2 The RVE with 8 of periodically distributed cells

The boundary conditions applied on the surface of the RVE are given as

uims

 RVE

  ik xk

 RVE

, ti

 pore

  ikms nk

 pore

.

For the proof of periodical symmetry of the solution, let us consider  ik   i 2 k 2 22 . Then, the boundary conditions applied on the outer boundary of the RVE and the geometrical symmetry yield the following boundary conditions (b.c.) on two neighbour cells shown in Fig. 3

b.c. on sides of the L-cell: x1  0(l ) : uims   i 2 22 x2

b.c. on sides of the R-cell: x1  0(l ) : uims   i 2 22 x2

x3  0(l ) : uims   i 2 22 x2 x2  l : uims   i 2 22l x2  0 : uims  0

x3  0(l ) : uims   i 2 22 x2 x2  l : uims   i 2 22l x2  0 : uims  0

shift of c.s.

RBM

ui  ui  li 2 22

b.c. on sides of the L-cell: x1  0(l ) : uims   i 2 22 ( x2  l )

b.c. on sides of the L-cell: x1  0(l ) : uims   i 2 22 x2

x3  0(l ) : uims   i 2 22 ( x2  l ) x2  l : uims  0

x3  0(l ) : uims   i 2 22 x2 x2  l : uims   i 2 22l

x2  0 :

x2  0 :

uims   i 2 22l

uims  0

Fig. 3 Boundary conditions on two neighbour cells Finally, one can see that the boundary conditions on sides of the L-cell are RBM (rigid body motion) equivalent with the boundary conditions on sides of the R-cell. Similarly, it can be shown also for the b.c.  ik   i1 k111 . Thus, we know the b.c. on particular cells and it is sufficient to solve the BVP in one cell instead of the BVP in the RVE with more cells, as long as

 ik   i1 k111 and/or  ik   i 2 k 2 22 . Unfortunately, we cannot find the complete b.c. on the outer boundary of particular cells if  ik is not diagonal and we need to solve the BVP in the RVE with more cells, though it does not mean that the solution is not periodically symmetric even in such a case. Recall that the material coefficients of the linear elastic skeleton are assumed to be constant in order to justify the concept of the RVE in macrostructure with periodically distributed pores.

2.1 Numerical evaluation of effective material coefficients In order to evaluate the effective material coefficients, it is sufficient and necessary to solve numerically the five micro-structural boundary value problems specified in the previous passage. For this purpose, we have utilized the finite element method software provided by the commercial package COMSOL [19] with using the second order Lagrange polynomials as the shape functions [17], [18]. The effective material parameters were calculated from expression (7) for different porosities defined as p  Vp / Vt where V p is the volume of the pore or void space and Vt is the total or bulk volume of the material. To show the limit cases, we present the numerical values of the effective coefficients for two extremes of considered values of porosity: for the lowest value of porosity p  4.18 109

(i)

eff cab

2.979 -3.013 -3.4240   2.826e11 1.211e11   2.826e11 5.494 -1.792 -3.203    8.0769e10 -0.963 -0.819    8.0769e10 1.634    8.0769e10  

(ii)

eff cab

for the highest value of porosity p  0.40

1.976e11 6.113e10 19.663e5  1.976e11 7.917e5   5.678e10    

2.026e4 2.303e6 5.867e5 5.678e10

-1.511e6   1.643e6  2.775e5   -2.885e5  5.678e10 

It can be seen that the anisotropic nature of the homogenized elastic material is significantly more pronounced for higher levels of porosity.

Furthermore, in the limit p  0 , the

homogenized material becomes isotropic like the skeleton. 3. Unified formulation of plate bending problems in homogenized linear elastic media with empty micro-voids

Let us consider a straight plate structure occupying the 3D domain V

{ ( x1, x2 , x3 )

3

;x

( x1, x2 )

, x3

[ h / 2, h / 2]}

[ h / 2, h / 2] .

In order to unify the formulation for three plate bending theories (KLT- Kirchhoff-Love theory, FSDPT- 1st order shear deformation plate theory and TSDPT-3rd order shear deformation plate theory), the general expression for displacements is utilized in the form [7], [8].





vi (x, x3 )  i u (x)  c1 ( x3 )  x3  w, (x)  c1 ( x3 ) (x)  i3w(x) ,

(8)

In which u (x) ,  (x) , w(x) are in plane displacements, rotations, and deflection fields respectively. With assuming the expansion of in-plane displacements up to the third power of the transverse coordinate x3 , where  ( x3 ) : x3  c2 ( x3 ) ,  ( x3 ) :

4( x3 )3 3h 2

,

and  0, KLT , c1    1, SDPT

 0, c2    1,

FSDPT, KLT , TSDPT

(9)

are two factors introduced for switching among tree various theories incorporated in the unified formulation. Now, in view of (8), the strain tensor can be written as e (x, x3 )   (x)  c1 ( x3 )  x3  w, (x)  c1 ( x3 ) (x) ,

c e 3 (x, x3 )  1 ,3 ( x3 )  w, (x)   (x)  , 2

where

e33 (x, x3 )  0

(10)

 



1  ,   , 2



and  



1 u ,  u , 2



and ,3 ( x3 )  1  c2 ,3 ( x3 ) ,  ,3 ( x3 )  4  x3 / h  . 2

Recall that c1,3 (h / 2) vanishes and hence also e 3 (x, h / 2)  0   3 (x, h / 2) only if c1  0 (KLT) and/or c1  c2  1 (TSDPT). In the KLT, however, the shear stresses and strains are vanishing e 3 (x, x3 )  0   3 (x, x3 ) throughout the plate thickness x3 [h / 2, h / 2] , which is the shortcoming of the KLT. On the other hand, the FSDPT suffers from the shortcoming that the shear strains (stresses) are constant within the whole plate thickness, e 3 (x, x3 )  const for x3 [h / 2, h / 2] .

Recall that in the case of plate structures the trace of the strain tensor (divergence of the displacement field) is written as ekk (x, x3 )  e (x, x3 )  vk ,k (x, x3 )  v , (x, x3 )   u , (x)  c1 ( x3 )  x3  w, (x)  c1 ( x3 ) , (x) .

(11)

In order to derive the governing equations for the plate bending problem, we utilize the principle of virtual work  U   We  0 . Variations of the deformation energy and the work of external forces applied to the plate structure are given as 

h /2



   h /2



 U      ij (x, x3 ) eij (x, x3 )dx3  d  ,

 We   q(x) w(x)d  .

(12)



It can be seen

 ij eij     e  2  3 e 3

   e     u ,  c1 ( x3 )  x3  w,  c1 ( x3 ) ,  2  3 e 3    3c1,3 ( x3 )  w,    eff eff eff    c e  c 3 2e 3  c u ,   c1 ( x3 )  x3  w,  c1 ( x3 ) ,   eff  c1c 3,3 ( x3 )  w,   

(13)

  3  ceff3 e  ceff3 3 2e 3  ceff3 u ,  c1 ( x3 )  x3  w,  c1 ( x3 ) ,    c1ceff3 3,3 ( x3 )  w,    eff Where c is 4th rank tensor of effective elastic coefficients given in the matrix form by Eq. (6)

with utilizing the Voigt notation. Since the dependence on the x3 -coordinate is known a priori, the integration with respect to this transversal coordinate can be performed in advance and we obtain the formulation in 2D domain  . In view of (12) and (13), we have





( ) ( w) )    U   T (x) u , (x)  M (x) , (x)  M (x) w,  T( w 3 (x)  w, ( x)   ( x)  d  

(14)

in which the integral fields (average values over the plate thickness) are defined as h /2

eff eff T T T (x) :    (x, x3 )dx3  c   (x)  c 3  3 ( x)

(15)

 h /2

h /2

T   (x) : 

 h /2

u , (x)  c1 ( x3 )  x3  w, (x)  c1 ( x3 ) , (x) dx3

 T3 (x) :  w, (x)   (x)  c1  ,3 ( x3 )dx3 h /2

 h /2

h /2  M   x  c eff   M   x ( ) eff M (x) : c1   ( x3 )  (x, x3 )dx3  c      3  3  

(16)

 h /2

 M   x : h /2 c  ( x ) u (x)  c  ( x )  x w ( x)  c  ( x ) ( x) dx       1 3   ,  1 3 3  ,  3 1 3  ,  h /2 h /2

 c1 

 h /2

 ( x )u 3

 , ( x)  c2 ( x3 ) ( x3 ) w, ( x)   ( x3 ) 

2



 , (x) dx3

h /2 M   3   x  :  w, (x)   (x)  c1   ( x3 ),3 ( x3 )dx3  h /2

h /2

( w) M (x) : 

 h /2

 Mw x  ceff   Mw x eff      3  3    x3  c1 ( x3 ) (x, x3 )dx3  c

(17)

 Mw x : h /2 x  c  ( x ) u (x)  c  ( x )  x w (x)  c  ( x ) (x) dx        3 1 3   ,  1 3 3  ,  3 1 3  ,  h /2 h /2

 

 h /2

 x  c  ( x )u 3

1

3



 , (x)  x3  c1 ( x3 )

2 w, (x)  c1c2 ( x3 ) ( x3 ) , (x) dx3

h /2 Mw   3   x  :   x3  c1 ( x3 ) c1,3 ( x3 )  w, (x)   (x)  dx3   h /2





h /2

 w, (x)   (x) c1c2   ( x3 ),3 ( x3 )dx3 h /2



 h /2



) eff  (T ) (Tw ) T( w ( x)   3 ( x) :  c1   c2 ,3 ( x3 )   3 ( x, x3 ) dx3  c 3 c1  ( x)      h /2

(18)

)  ceff3 3 c1 (T3 ) (x)   (Tw ( x)  3   (Tw )   (x) : c1c2   ,3 ( x3 ) u , (x)   c1 ( x3 )  x3  w, ( x)  c1 ( x3 ) , ( x) dx3  h /2

 h /2 h /2

 c1c2 

 h /2

 ,3 ( x3 )u , (x)  c2 ( x3 ) ,3 ( x3 )w, (x)   ( x3 ) ,3 ( x3 ) , (x) dx3

)  (Tw (x) :   w, (x)   (x)  c1c2   ,3 ( x3 ),3 ( x3 )dx3 3 h /2

 h /2

The factor   (1  c2 )k  c2 was introduced by Reissner as modification of the shear stresses by the shear correction factor k in the FSDPT. The integrations defined in Eqs. (15)-(18) are accomplished in the Appendix. Eq. (14) can be rearranged as 

  w  ( w)   V  w d    n  

( )  U   n T  u  n M    M ( w)   



  T ,  u  



( ) M ,

)  T( w 3

  

( w) M  ,

)  T( w 3,

 w d 

(19)

where ( w) ( w) , T ( w) : t n M  , M ( w) : n n M 





( w) ( w ) V ( w) : n M    T ( w) (x c )  (x  x c )  ,   T 3 c

T ( w) . t

(20)

Now, in view of (12) and (19), the principle of virtual work  U   We  0 results into the governing equations T ,  0 ( w) ( w ) M ,  T 3,  q

(21)

( ) ( w ) M ,  T 3  0

and the boundary restrictions n T  u



 0  n T



 0 or [  u



 0 , i.e. u



is prescribed ]

 w   w   w  M ( w)    0  M ( w)  0 or [     0 , i.e.   is prescribed ]   n    n    n   ( ) n M  

V ( w) w





( )  0  n M 

 0  V ( w)





 0 or [ 



 0 , i.e. 



(22)

is prescribed ]

 0 or [  w   0 , i.e.  w  is prescribed ]

Substituting (A.11)-(A.14) into (21), we obtain the governing equations expressed in terms of the primary fields







eff eff c hd(0)0u , (x)  c 3c1h d(0)0  4c2 d(0)2 w, (x)   , (x)  0

(23)

  16 4 4   eff c h3  (c1  1)d(0)2  c1c2 d(0)6  w, (x)  c1c2  d(0)4  d(0)6   , (x)   9 3 3      ceff3 c1h  d(0)0  4c2 d(0)2  u , (x) 











 ceff3 3c1h  d (0)0  4c2 d(0)2  4c2 d(0)2  4d(0)4  w, (x)   , ( x)  q  

(24)   8 16 4  4   eff c c1h3  d (0)2  c2 d(0)4  c2 d(0)6   , (x)  c2  d(0)4  d(0)6  w, (x)   3 9 3  3      ceff3 c1h  d(0)0  4c2 d(0)2  u , (x) 







(25)





 ceff3 3c1h  d (0)0  4c2 d(0)2  4c2 d(0)2  4d(0)4  w, (x)   (x)  0  

In the governing equations (23)-(25), we have assumed the plate thickness to be constant. 4. Numerical solution and numerical examples In order to study a bending response of the plate which contains empty voids (characterized by porosity p ) the clamped square plates are considered with the side length L  1 [m] and two various values of the thickness h ( L / h  50 for thin plate and L / h  5 for thick plate). The plate is subjected to uniform transversal loading with q  107 [Pa] and the effective material coefficients are computed numerically by employing the proposed homogenization procedure for several values of porosity. Obviously, the standard C0 continuous finite elements are insufficient for discretization of the derived governing equations and approximation of filed variables such as deflections w and rotations  . In order to eliminate high order derivatives of field variables, the original system of the 4th order PDE (23)-(25) can be decomposed into the system of the 2nd order PDE by introducing new field variables defined by the 2nd order derivatives of the relevant primary fields m (x) : w, (x)

f (x) :  , (x)

,

.

(26) Substituting (26) into the governing equations (23)-(25), we obtain the system of the PDE for





field variables u , w,  , m , f and the original boundary value problems for plate bending are represented by the system of the 2nd order PDE, which can be solved via the finite element method provided by the commercial software COMSOL [19] with using the second order Lagrange polynomials as the shape functions [17], [18] and the Weak Form module.

Fig. 4 Deflection distributions for thin plates by three different theories (KLT, FSDT, TSDT) and various levels of porosity The influence of the level of porosity on deflections of thin and thick plates is shown in Fig. 4 and Fig. 5, respectively. It can be seen that the deflections are increasing with increasing the level of porosity (the bending stiffness is decreasing). The difference between the results by three various theories are negligible in the case of thin plates while they are observable in the case of thick plates. Apparently, the KLT fails in the case of thick plates.

Fig. 5 Deflection distributions for thick plates by three different theories (KLT, FSDT, TSDT) and various levels of porosity As the thickness of the plate increases the influence of shear strains on deflections is more expressive and the differences between the results by different plate theories become visible (Fig. 6). The KLT theory neglects the shear strains and is insensitive to the change of the plate thickness. Therefore the KLT is unreliable for bending analysis of thick plates. Remembering the shortcoming of the FSDT, one should be careful even with using the results by the FSDT. Thus, the results by the TSDT appear to be the most reliable. The anisotropy properties of the homogenized medium give rise to coupling between the in-plane deformation and bending modes, which is documented by finite in-plane displacements shown in Fig.7.

Fig. 6 Comparison of deflections computed by three different theories (KLT, FSDT, TSDT) for one value of porosity in: (i) thin plates, (ii) thick plates p=4.18e-9

p=0.11

p=0.26

p=0.40

Fig. 7. The contour maps of the in plane displacement field u1 (x) for thick plate (TSDT results).

5. Conclusion The main outcome of the paper is the study of bending of thin and/or moderately thick plates containing the empty micro-voids with using the multi-scale modelling. The plate bending is considered as the macro-structural elasticity problem in the homogenized medium with effective material coefficients. The complexity of the heterogeneous structure is treated via the concept of homogenization, when the effective material coefficients are calculated by solving properly selected boundary value problems with modelling microstructure in the representative volume elements. It appears that the homogenized medium can exhibit different material properties than the skeleton and it is not reduced to simple change of values of the non-vanishing coefficients in skeleton. Since in a plate bending theory the real three dimensional problem of elasticity is reduced to the two dimensional one, physically plausible is only the calculation of effective material coefficients from three dimensional representative volume element. Having known the material properties and effective coefficients of the homogenized medium, the unified formulation was developed for plate bending problems with assuming alternatively the assumptions of three plate bending theories, such as the Kirchhoff-Love theory for bending of thin plates as well as the 1st and 3rd order shear deformation plate theory. It has been shown that homogenized medium exhibits anisotropy, even if the skeleton was isotropic and the effect of anisotropy is becoming more pronounced with increasing the level of porosity. The coupling between in-plane deformation and bending modes arises due to the non-zero terms in full anisotropic matrix. In numerical simulations, it has been confirmed that the deflections of porous plates are increasing with increasing the level of porosity. In this paper, we assumed periodic distribution of voids. This assumption allows simplified evaluation of effective material coefficients in homogenized elastic continuum. Of course, the variability in the distribution of voids with certain density can affects the effective coefficients, but such a fine effect is not studied in this paper.

Appendix In this appendix, we perform the integrations defined in Eqs. (15)-(18). All the integrals can be expressed in terms of 1/2

d(0) a :



z a dz 

1/2

1  (1) a . (a  1)2a 1

(A.1)

Since d(0) a  0 for a  2n  1 , ( n  0,1, 2,... ), one can obtain

u , (x)  c1 ( x3 )  x3  w, (x)  c1 ( x3 ) , (x) dx3  hd(0)0u , (x)

(A.2)

 T3 (x)   w, (x)   (x)  c1  ,3 ( x3 )dx3  c1h  d(0)0  4c2 d(0)2   w, (x)   (x) 

(A.3)

h /2

T   ( x)  

 h /2

h /2

 h /2





 M   x  c h /2  ( x )u (x)  c  ( x ) ( x ) w (x)   ( x ) 2  ( x) dx      1   3   , 3  , 2 3 3 , 3  h /2

3 

 8 16 4  4    c1h  d (0)2  c2 d(0)4  c2 d(0)6   , (x)  c2  d(0)4  d(0)6  w, (x)  3 9 3  3     h /2 M   3   x  :  w, (x)   (x)  c1   ( x3 ),3 ( x3 )dx3  0

(A.4)

(A.5)

 h /2

 Mw x  h /2     

 h /2

 x  c  ( x )u 3

1

3



 , ( x)  x3  c1 ( x3 )

2 w, (x)  c1c2 ( x3 ) ( x3 ) , (x) dx3 

3 

 16 4 4    h  (1  c1 )d(0)2  c1c2 d(0)6  w, (x)  c1c2  d(0)4  d(0)6   , ( x)  9 3 3     h /2 Mw   3   x    w, (x)   (x)  c1c2   ( x3 ),3 ( x3 )dx3  0

(A.6)

(A.7)

 h /2

(Tw )   (x)  c1c2   ,3 ( x3 )u , (x)  c2 ( x3 ) ,3 ( x3 ) w, (x)   ( x3 ) ,3 ( x3 ) , (x) dx3  h /2

 h /2

 4c1c2 hd(0)2u , (x)

(A.8)

)  (Tw (x)    w, (x)   (x)  c1c2   ,3 ( x3 ),3 ( x3 )dx3  3 h /2

 h /2





 4c1c2 h d(0)2  4c2 d(0)4 w, (x)   (x)

(A.9)



Finally, in view of Eqs. (15)-(18) and (A.2)-(A.9), we obtain the explicit expressions for the integral fields in terms of the primary fields ( u , w,  ) and their derivatives as





eff eff T (x)  c hd(0)0u , (x)  c 3c1h d(0)0  4c2 d(0)2 w, (x)   (x)



(A.10)

  8 16 4  4   ( ) eff M (x)  c c1h3  d(0)2  c2 d(0)4  c2 d(0)6  , (x)  c2  d(0)4  d(0)6  w, (x)  3 9 3  3    

(A.11)

  16 4 4   ( w) eff M (x)  c h3  (c1  1)d(0)2  c1c2 d(0)6  w, (x)  c1c2  d(0)4  d(0)6   , (x)  9 3 3    

(A.12)

) eff   T( w 3 (x)  c 3 c1h  d (0)0  4c2 d(0)2  u , ( x) 









 ceff3 3c1h  d(0)0  4c2 d(0)2  4c2 d(0)2  4d(0)4  w, (x)   (x)  



(A.13)

Note that there is coupling between the in-plane deformation modes and the bending modes, eff eff provided that c 3  0 and/or c 3  0 .

Acknowledgment. The authors gratefully acknowledge the supports by the Slovak Science and Technology Assistance Agency APVV-14-0440. References [1] Reddy JN. Mechanics of laminated composite plates: theory and analysis. Boca Raton: CRC Press: 1997. [2] Reddy JN. A general non-linear third-order theory of plates with moderate thickness. Int J Nonlinear Mech 1990;25:677-686 [3] Reddy JN. A Simple higher-order theory for laminated plates [4] Barber JR, Elasticity, New York; Springer Science+Business Media B.V.; 2010. [5] Lurie AL. Theory of elasticity, Berlin: Springer-Verlag; 2005.

[6] Sladek V, Sladek J, Sator L. Physical decomposition of thin plate bending problems and their solution by meshfree methods. Eng Anal Boundary Elem 2013; 37:348-65 [7] Sator L, Sladek V, Sladek J. Coupling effects in elastic analysis of FGM composite plates by mesh-free methods. Compos Struct 2014; 115: 309–322. [8] Sator L, Sladek V, Sladek J. Elastodynamics of FGM plates by meshfree method. Compos Struct 2016; 140: 100–110. [9] Carrera E. On the use of transverse shear stress homogeneous and non-homogeneous conditions in third-order orthotropic plate theory. Compos Struct 2007;77:341-352 [10] Mantari JL, Oktem AS, Guedes Soares C. Bending response of functionally graded plates by using a new higher order shear deformation theory. Compos Struct 2012; 94(2):71423 [11] Ferreira AJM, Batra RC, Roque CMC, Qian LF, Martins PALS. Static analysis of functionally graded plates using third-order shear deformation theory and a meshless method. Compos Struct 2005;69(4):449-57 [12] Hill R. A self-consistent mechanics of composite materials. J Mech Phys Solids 1965; 13:213-222. [13] Sladek V, Musil B, Sladek J. Effective elasticity coefficients in dry porous materials. Numerical and semi-analytical approaches, In: Proc. of the VII European Congress on Computational Methods in Applied Sciences and Engineering - ECCOMAS Congress 2016 (M. Papadrakakis, V. Papadopoulos, G. Stefanou, V. Plevris, eds.), Crete Island, Greece 5-10 June 2016; https://doi.org/10.7712/100016.1941.7510 [14] Sladek V, Musil B, Sladek J, Kasala J. Microstructural Evaluation of Effective Elasticity Coefficients in Materials with Micro-voids. Procedia Engineering 2017;190:170-177 [15] Aboudi J, Arnold, B.A. Bednarcyk, Micromechanics of Composite Materials, Amsterdam, Elsevier;2013 [16] Dong B, Li C, Wang D, Wu CT. Consistent multiscale analysis of heterogeneous thin plates with smoothed quadratic Hermite triangular elements. Int J Mech Mater Des 2016; 12:539-562 [17] Johnson C, Numerical Solution of Partial Differential equations by the Finite Element Method, New York :Dover Publications; 2009

[18] Hughes TJR. The finite element method, New Jersey: Prentice-Hall Inc; 1987 [19] COMSOL Multiphysics Reference Manual, COMSOL, Inc, www.comsol.com