Modeling of Size Effects in Bending of Perforated Cosserat Plates

This paper presents the numerical study of Cosserat elastic plate deformation based on the parametric theory of Cosserat plates, recently developed by the authors.The numerical results are obtained using the Finite ElementMethod used to solve the parametric system of 9 kinematic equations.We discuss the existence and uniqueness of the weak solution and the convergence of the proposed FEM. The Finite Element analysis of clamped Cosserat plates of different shapes under different loads is provided. We present the numerical validation of the proposed FEM by estimating the order of convergence, when comparing the main kinematic variables with an analytical solution. We also consider the numerical analysis of plates with circular holes. We show that the stress concentration factor around the hole is less than the classical value, and smaller holes exhibit less stress concentration as would be expected on the basis of the classical elasticity.


Introduction
A complete theory of asymmetric elasticity introduced by the Cosserat brothers [1] gave rise to a variety of beam, shell, and plate theories.The first theories of plates that take into account the microstructure of the material were developed in the 1960s.Eringen proposed a complete theory of plates in the framework of Cosserat (micropolar) elasticity [2], while independently Green et al. specialized their general theory of Cosserat surface to obtain the linear Cosserat plate [3].Numerous plate theories were formulated afterwards; for the extensive review of the latest developments we recommend referring to [4].
The first theory of Cosserat elastic plates based on the Reissner plate theory was developed in [5] and its Finite Element modeling is provided in [6].The parametric theory of Cosserat plate, presented by the authors in [7], includes some additional assumptions leading to the introduction of the splitting parameter.The theory provides the equilibrium equations and constitutive relations and the optimal value of the minimization of the elastic energy of the Cosserat plate.The paper [7] also provides the analytical solutions of the presented plate theory and the three-dimensional Cosserat elasticity for simply supported rectangular plate.The comparison of these solutions showed that the precision of the developed Cosserat plate theory is compatible with the precision of the classical plate theory developed by Reissner [8,9].
The numerical modeling of bending of simply supported rectangular plates is given in [10].We develop the Cosserat plate field equations and a rigorous formula for the optimal value of the splitting parameter.The solution of the Cosserat plate was shown to converge to the Reissner plate as the elastic asymmetric parameters tend to zero.The Cosserat plate theory demonstrates the agreement with the size effect, confirming that the plates of smaller thickness are more rigid than expected from the Reissner model.The modeling of Cosserat plates with simply supported rectangular holes is also provided.
The extension of the static model of Cosserat elastic plates to the dynamic problems is presented in [11].The computations predict a new kind of natural frequencies associated with the material microstructure and were shown to be consistent with the size effect principle known from the Cosserat plate deformation reported in [10].
Since [10] was restricted only to the case of rectangular plates, the current article represents an extension of this work for the Finite Element modeling of the Cosserat plates for the different shapes, under different loads and different boundary conditions.We discuss the existence and uniqueness of the weak solution, convergence of the proposed FEM, and its numerical validation.
We also present the modeling of size effects in Cosserat plates with holes.We evaluate the stress concentration factor around the hole and show that the factor is smaller than the one based on the Reissner theory for simple elastic plates.The Finite Element comparison of the plates with holes confirms that smaller holes exhibit less stress concentration than the larger holes.

Parametric Theory of Cosserat Plates
In this section we provide the main equations of the theory presented in [7].
Throughout this article Greek indices are assumed to range from 1 to 2, while the Latin indices range from 1 to 3 if not specified otherwise.We will also employ the Einstein summation convention according to which summation is implied for any repeated index.
Let us consider a thin plate  of thickness ℎ and  3 = 0 representing its middle plane.The sets  and  are the top and bottom surfaces contained in the planes  3 = ℎ/2 and  3 = −ℎ/2, respectively, and the curve Γ is the boundary of the middle plane of the plate.The set of points  = (Γ × [−ℎ/2, ℎ/2]) ∪  ∪  forms the entire surface of the plate.Γ  × [−ℎ/2, ℎ/2] is the part of the boundary where displacements and microrotations are prescribed, while Γ  × [−ℎ/2, ℎ/2] is the part of the boundary where stress and couple stress are prescribed.
The equilibrium system of equations for Cosserat plate bending is given as where   are micropolar couple moments, all defined per unit length.The initial pressure  is represented here by the pressures p1 =  opt  and p2 = (2/3)(1 −  opt ), where  opt is the splitting parameter [7].
The constitutive formulas for Cosserat plate are given in the following form [7]: where  and  are symmetric constants and , , , and  are the asymmetric Cosserat elasticity constants [5].The vector of kinematic variables is given as where The system of equilibrium equations is accompanied by the zero variation of the work density with respect to the splitting parameter .This allows us to split the bending pressure on the plate ( 1 ,  2 ) into two parts corresponding to different orders of stress approximation. where and S() and E() are the Cosserat plate stress and strain sets The optimal value  op of the splitting parameter is given as in [10]  opt = 2W (00) − W (10) − W (01) 2 (W (11) + W (00) − W (10) where The approximation of the components of the threedimensional tensors   and   is where The approximation of components of the three-dimensional strain and torsion tensors   and   is The components of the three-dimensional displacements   and microrotations   [7] are where  = 2 3 /ℎ and ,  ∈ {1, 2} as before.

Modelling and Simulation in Engineering
We consider the vertical load and pure twisting momentum boundary conditions at the top and bottom of the plate The expressions for the three-dimensional stress, strain, and displacement components satisfy the equilibrium Cosserat elasticity equations the strain-displacement and torsion-rotation relations and the constitutive equations where   is the Levi-Civita symbol.
In order to obtain the micropolar plate bending field equations in terms of the kinematic variables, the constitutive formulas in the reverse form are substituted into the bending system of (1).The obtained Cosserat plate bending field equations can be represented as an elliptic system of nine partial differential equations in terms of the kinematic variables [10]: where  is a linear differential operator acting on the vector of kinematic variables V (unknowns): and () is the right-hand side vector that in general depends on : The operators   are given as follows: The coefficients   are given as

Finite Element Method for Cosserat Plates
The right-hand side of the system (17) depends on the splitting parameter  and so does the solution V that we will formally be denoted as V  .Therefore the solution of the Cosserat elastic plate bending problem requires not only solving the system (17), but also an additional technique for the calculation of the value of the splitting parameter that corresponds to the unique solution.The use of the splitting parameter adds an additional degree of freedom to the modeling.The optimal value of this parameter provides the highest level of approximation to the original 3dimensional problem.Considering that the elliptic systems of partial differential equations correspond to a state where the minimum of the energy is reached, the optimal value of the splitting parameter should minimize the elastic plate energy [12].The minimization corresponds to the zero variation of the plate stress energy (5).The Finite Element Method for Cosserat elastic plates is based on the algorithm for the optimal value of the splitting parameter.This algorithm requires solving the system of equations (17) for two different values of the splitting parameter , numerical calculation of stresses, strains, and the corresponding work densities.We will follow [10] in the description of our Finite Element Method algorithm.
(1) Use classical Galerkin FEM to solve two elliptic systems: for V 0 and V 1 , respectively.
(2) Calculate the optimal value of the splitting parameter  opt using (8).
(3) Calculate the optimal solution V  0 of the Cosserat plate bending problem as a linear combination of V 0 and V 1 : 3.1.Weak Formulation of the Clamped Cosserat Plate.Let us consider the following hard clamped boundary conditions similar to [13]: where n and ŝ are the normal and the tangent vectors to the boundary.These conditions represent homogeneous Dirichlet type boundary conditions for the kinematic variables: Let us denote by L 2 ( 0 ) the standard space of squareintegrable functions defined everywhere on  0 and by H 1 ( 0 ) the Hilbert space of functions that are squareintegrable together with their first partial derivatives: Let us denote the Hilbert space of functions from H 1 ( 0 ) that vanish on the boundary as in [14]: The space H 1 0 ( 0 ) is equipped with the inner product: Taking into account the fact that the boundary conditions for all variables are of the same homogeneous Dirichlet type, we look for the solution in the function space H( 0 ) defined as The space H is equipped with the inner product ⟨, V⟩ H : and relative to the metric induced by the norm ‖‖ = √⟨, ⟩ H , the space H is a complete metric space and therefore is a Hilbert space [15].Let us consider a dot product of both sides of the system of the field equations ( 17) and an arbitrary function V ∈ H: and then integrate both sides of the obtained scalar equation over the plate  0 : Let us introduce a bilinear form (, V) : The expression for (, V) Modelling and Simulation in Engineering 7 is a summation over the terms of the form where V  ∈ H  ,   ∈ H  , and L is a scalar differential operator.
There are 3 types of linear operators present in the field equations ( 17)-operators of order zero, one, and two, which are constant multiples of the following differential operators: where  is a 2 × 2 real matrix.For example, for the order two operator  1 ( 2 / 2 1 )+ 2 ( 2 / 2 2 ) the corresponding matrix  would be  = [  1 0 0  2 ].These operators act on the components of the vector  and are multiplied by the components of the vector V and the obtained expressions are then integrated over  0 : where V  ∈ H  and   ∈ H  .The weak form of the second-order operator is obtained by performing the corresponding integration by parts and taking into account the fact that the test functions V  vanish on the boundary  0 : The expression for  () (V) represents a summation over the terms of the form: Taking into account the fact that the optimal solution of the field equations (17) such that no vertex of the triangular element lies on the edge of another triangle (see Figure 1).Let us introduce the mesh parameter ℎ as the greatest diameter among the elements   : which for the triangular elements corresponds to the length of the longest side of the triangle.We now define the finite dimensional space Ĥℎ as a space of all continuous functions that are linear on each element   and vanish on the boundary: By definition H ℎ  ⊂ H  , and the Finite Element space H ℎ is then defined as The approximate weak solution  ℎ can be found from the Galerkin formulation of the clamped Cosserat plate bending problem [16].
Galerkin Formulation of the Clamped Cosserat Plate.Find all  ℎ ∈ H ℎ and  ∈ R that minimize the stress plate energy  S  ( ℎ , ) subject to The description of the function V ℎ  ∈ H ℎ  is provided by the values V ℎ  (  ) at the nodes   ( = 1, ).Let us define the set of basis functions { 1 ,  2 , . . .,   } of each space H ℎ  as excluding the points   on the boundary  0 .Therefore and the functions   are nonzero only at the node   and those that belong to the specified boundary and the support of   consists of all triangles   with the common node   (see Figure 2).Since the spaces H ℎ  are identical they will also have identical sets of basis functions   ( = 1, ).Sometimes we will need to distinguish between the basis functions of different spaces assigning the superscript of the functions space to the basis function; that is, the basis functions for the space H ℎ  are    .For computational purposes these superscripts will be dropped.

Calculation of the Stiffness Matrix and the Load Vector.
The bilinear form of the Galerkin formulation (48) is given as Since  ℎ  ∈ H ℎ  then there exist such constants Since ( 51) is satisfied for all V ℎ  ∈ H ℎ  then it is also satisfied for all basis functions  ()   ( = 1, ): where We define the block stiffness matrices   (,  = 1, 9): For computational purposes the superscripts of the basis functions can be dropped and the block stiffness matrices   can be calculated as Let us define the block load vectors   () ( = 1, 9) and the solution block vectors   corresponding to the variable  ℎ  ( = 1, 9): Equation ( 48) of the Galerkin formulation can be rewritten as The global stiffness matrix consists of 81 block stiffness matrices  i , the global load vector consists of 9 block load vectors   (), and the global displacement vector is represented by the 9 blocks of coefficients   .The entries of the block matrices   and the block vectors   () can be calculated as The block matrix form of equation ( 48) is given as ] .(61)

Existence and Convergence
Remarks.We will follow [17,18], where the analysis of the analytic regularity for linear elliptic systems and their general treatment are recently presented.
Employing integration by parts for the second-order operators   the bilinear form (⋅, ⋅) can be rewritten in the following form: where  and  are multi-indices.Since coefficients   are constant and therefore bounded on  0 , the bilinear form (⋅, ⋅) is continuous over H [17]; that is, there exists a constant  > 0 such that The strong ellipticity of the operator  was shown in [10].Since the operator  is strongly elliptic on  0 the bilinear form (⋅, ⋅) is -elliptic on H [12,17]; that is, there exists a constant  > 0 such that The existence of the solution of the weak form (43) and its uniqueness are the consequences of the Lax-Milgram Theorem [19,20].Note that the existence and uniqueness of the Galerkin weak problem (48) are also a consequence of the Lax-Milgram theorem, since the bilinear form (⋅, ⋅) restricted on H ℎ obviously remains bilinear, continuous, and -elliptic [12].Lax-Milgram theorem also states that the solution is bounded by the right-hand side which represents the stability condition for the Galerkin method.
The convergence of the Galerkin approximation follows from Céa's lemma and an additional convergence theorem [12,21].On the polygonal domains the sequence of subspaces of H = H 1 0 ( 0 ) 9 can be obtained by the successive uniform refinement of the initial mesh using the midpoints as new nodes thus subdividing every triangle into 4 congruent triangles.Therefore H  ⊂ H +1 for every  ∈ N and the sequence of spaces H  is dense in H [22], and thus and   converges to  as  → ∞ [12,23].It was shown that there exists a sequence of triangulations that ensures optimal rates of convergence in H 1 -norm for the FEM approximation of the second-order strongly elliptic system with zero Dirichlet boundary condition on polyhedron domain with continuous, piecewise polynomials of degree  [24].For the details on different norms and rates of convergence of the FEM we refer the reader to [20].

Validation of the FEM for Different Boundary Conditions
Let us consider the plate  0 to be a square plate of size [0, ] × [0, ] with the boundary Modelling and Simulation in Engineering the hard simply supported boundary conditions similar to [13], written in terms of the kinematic variables in the mixed Dirichlet-Neumann type: where The existence of a sequence of triangulations that ensures the optimal rates of convergence for the Finite Element approximation of the solution of a second-order strongly elliptic system with homogeneous Dirichlet boundary condition on polyhedron domain with continuous piecewise polynomials was shown in [24].For the case of piecewise linear polynomials the optimal rate of convergence in H 1norm is linear.
We propose using the uniform refinement to form the sequence of triangulations and estimate the order of the error of approximation of the proposed FEM in H 1 -norm and  2norm.
In our computations we consider plates made of polyurethane foam-a material reported in the literature to behave Cosserat-like and the values of the technical elastic parameters are presented in [25]: Taking into account the fact that the ratio / is equal to 1 for bending [25], these values of the technical constants correspond to the following values of symmetric and asymmetric parameters: which automatically satisfies homogeneous Dirichlet boundary conditions.Substituting the solution (70) into the system of field equations (17) we can find the corresponding righthand side function .The results of the error estimation of the FEM approximation in H 1 and  2 norms performed for the elastic parameters corresponding to the polyurethane foam are given in Tables 1 and 2, respectively.Let us consider mixed Neumann-Dirichlet boundary conditions.Simply supported boundary conditions represent this type of boundary conditions and therefore the FEM approximation can be compared with the analytical solution developed in [7] for some fixed value of the parameter .The results of the error estimation of the FEM approximation in H 1 and  2 norms performed for the elastic parameters corresponding to the polyurethane foam are given in Tables 3 and 4, respectively.
The boundary condition for the variable Ω 3 is a Neumann type boundary condition: and thus we will look for Ω 3 in the space H 1 (Δ,  0 ), where The boundary condition for the variables  and  * is a Dirichlet type boundary condition: and thus we will look for  and  * in the space H 1 0 ( 0 ) defined as [14] The boundary condition for the variables Ψ 1 , Ω 0 2 , and Ω0 2 is of mixed Dirichlet-Neumann type: and thus we will look for Ψ 1 , Ω 0 2 , and Ω0 2 in the following space [14]: The boundary condition for the variables Ψ 2 , Ω 0 1 , and Ω0 1 is of mixed Dirichlet-Neumann type: and thus we will look for Ψ 2 , Ω 0 1 , and Ω0 1 in the following space [14]: Therefore we will look for the solution of the Cosserat plate field equations (17) in the space H defined as where The space H is a Hilbert space equipped with the inner product ⟨, V⟩ H defined on H as follows: where ⟨, V⟩ H  is an inner product defined on the Hilbert space H  , respectively.Taking into account the essential boundary conditions we define the Finite Element spaces H ℎ  as follows: The finite dimensional space H ℎ is then defined as (87) We solve the field equations using described Finite Element Method and compare the obtained results with the analytical solution for the square plate made of polyurethane foam derived in [7].
The initial distribution of the pressure is assumed to be sinusoidal: The estimation of the error in H 1 norms shows that the order of the error is optimal (linear) in H 1 -norm for the piecewise linear elements for the simply supported plate.The results of the error estimation of the FEM approximation in H 1 and  2 norms performed for the elastic parameters corresponding to the polyurethane foam are given in Tables 5 and 6, respectively.
The comparison of the maximum of the displacements   and microrotations   calculated using Finite Element Method with 320 thousand elements and the analytical solution for the micropolar plate theory is provided in Table 7.The relative error of the approximation of the optimal value of the splitting parameter is 0.09%.

Numerical Results
In this section we study the Cosserat elastic plates with different shapes and boundary conditions.We investigate the size effect of circular clamped plate and of plate with holes.We also consider the plates with holes and calculate the stress concentration factor around the hole.The size effect of the Cosserat plate is illustrated in Figures 9 and 10.As in [10] we perform computations for different levels of the asymmetric microstructure by reducing the values of the elastic asymmetric parameters.For example, 1% of the microstructure means that 1% of the original values of the asymmetric parameters were used.The influence of the microstructure in thinner plates results in a more rigid plate than would be expected on the basis of the classical elasticity.

Stress Concentration.
In this section we will consider the plates with holes and calculate the stress concentration factor around the hole as a ratio of the highest stress  max to a reference stress  ref of the gross cross section: Let us consider a square plate of size [0, ] × [0, ] for  = 10.0 of thickness ℎ = 1.0 m and a hole of diameter  varying from 1.0 m to 5.0 m located in the center of the plate (note that in the notation of Figure 11, the diameter  defined as the largest distance in the interior of the hole is equal to 4).
The exterior boundary of the plate is  =  1 ∪  2 ∪  3 ∪  4 and the boundary of the hole is  =  1 ∪  2 ∪  3 ∪  4 .The plate is loaded under the uniform load and the following mixed hard simply supported and clamped boundary conditions written in terms of the kinematic variables:    (91) The stress distribution for the polyurethane foam plate and the stress concentration around the hole is provided in Figure 12.
The Finite Element calculations of the stress concentration factors around the hole show that the stress concentration factor for the Cosserat plate is always smaller than the classical value [26].The comparison of the stress concentration factors for Cosserat and Reissner plates is given in Figure 13.
The ratios of the Cosserat plate stress concentration factor to the Reissner plate concentration factor are plotted for different values of / and given in Figure 14.The presence of the length scale for Cosserat materials results in the predicted size effect.This graph confirms that the small holes exhibit less stress concentration than the larger holes as it is expected on the basis of the classical theory [26] and results in the enhanced toughness [27].Similar size effect results were  obtained earlier for the vertical deflection of Cosserat plates for the case of bending and for the main eigenfrequencies in the case of free vibration [10,11].

Conclusion
The article provides numerical results obtained using the Finite Element Method for solving the parametric system of 9 kinematic equations.The Finite Element modeling of the clamped Cosserat plates of different shapes under different loads is provided.The article also presents the numerical validation of the proposed FEM by estimating the order of convergence, when comparing the main kinematic variables with the analytical solution.Numerical simulation of size effects of Cosserat plates has been presented.It was shown that the stress concentration factor around the hole in the plates is less than the classical value, and smaller holes exhibit less stress concentration than would be expected on the basis of the classical elasticity.

Figure 1 :Figure 2 :
Figure 1: Example of the Finite Element triangulation of the domain  0 .

Figure 3
represents the Finite Element modeling of the bending of the simply supported square plate made of polyurethane foam.The comparison of the distribution of the vertical deflection of the clamped and simply supported plates is given in Figure4.The bending of a clamped circular plate under the uniform load is given in Figure5.Figures6-8represent the Finite Element modeling of the bending of the clamped polyurethane foam plates of different shapes under different loads including plates with holes.

Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :
Figure 3: Hard simply supported square plate 3.0 m × 3.0 m × 0.1 m made of polyurethane foam: the initial mesh and the isometric view of the resulting vertical deflection of the plate.

Figure 11 :
Figure 11: Square plate of size  ×  and thickness ℎ with a hole located in the center of the plate.

Figure 12 :
Figure 12: The distribution of the stress  11 for the square polyurethane foam plate 10.0 m × 10.0 m × 1.0 m demonstrates the stress concentration around the hole (160,000 elements used).

Figure 13 :Figure 14 :
Figure 13: Comparison of the stress concentration factors of the Cosserat plate with Reissner plate for different ratios /: Cosserat stress concentration factor is smaller than the classical value.

Table 1 :
Order of convergence in H 1 -norm for homogeneous Dirichlet BC.

Table 2 :
Order of convergence in  2 -norm for homogeneous Dirichlet BC.

Table 3 :
Order of convergence in H 1 -norm for mixed Neumann-Dirichlet BC.

Table 4 :
Order of convergence in  2 -norm for mixed Neumann-Dirichlet BC.

Table 5 :
Order of convergence in H 1 -norm for simply supported plate.

Table 6 :
Order of convergence in  2 -norm for simply supported plate.

Table 7 :
Relative error of the maximum values of the displacement and microrotations.
Figure 9: Comparison of the maximum vertical deflections of the circular Cosserat plate   with Reissner plate   for different values of width to thickness ratio /ℎ, which illustrates the influence of the thickness of the Cosserat plate.Figure 10: Comparison of the maximum vertical deflections of the circular plate with 4 holes,   (Cosserat) with   (Reissner) for different values of width to thickness ratio /ℎ, which illustrates the influence of the thickness of the Cosserat plate.