A new homogenization scheme for beam and plate structures without a priori requirements on boundary conditions

This contribution picks up on a novel approach for a first order homogenization procedure based on the Irving-Kirkwood theory and provides a finite element implementation as well as applications to beam and plate structures. It does not have the fundamental problems of dependency from representative volume element (RVE) size in determining the shear and torsional stiffness for beams and plates, that is present in classic Hill-Mandel methods. Due to the possibility of using minimal boundary conditions whilst simultaneously reusing existing homogenization algorithms, creation of models and numerical implementation are much more straight forward. The presented theory and FE formulation are limited to materially and geometrically linear problems. The approach to determining shear stiffness is based on the assumption of a quadratic shear stress distribution over the height (and width in case of the beam), which causes warping of the cross-section under transverse shear loading. Results for the homogenization scheme are shown for various beam and plate configurations and compared to values from well known analytical solutions or computed full scale models.


Introduction
When modeling large structures, there are two main concepts to reduce model complexity and with it computational costs. One is the use of beam, plate and shell elements, if the geometric prerequisites are met by the structure of interest. The other is a multiscale approach, that uses two layers of models -a macro level with a coarse mesh, in which the structure is assumed to be at least partwise homogeneous and a micro/meso level that accounts for inhomogeneities. The micro/meso scale serves the purpose of generating a material law for every point of the macro scale by arguments of separation of scales and homogenization. We use the term micro/meso scale as for beams and plates the lengths of inhomogeneties can be of the same order of magnitude as the B Maximilian Müller mmueller.paper@gmail.com 1 Fachgebiet Festkörpermechanik, Technische Universität Darmstadt, Franziska-Braun-Str. 7, 64287 Darmstadt, Germany 2 Lehrstuhl für Baustatik und Baudynamik, Rheinisch-Westfälische Technische Hochschule Aachen, Mies-van-der-Rohe-Straße 1, 52074 Aachen, Germany lengths of the cross-sections. Therefore, the cross-section is modeled in full length in the micro/meso scale which makes it more of a meso than a micro scale for beams and plates. For more details, see Sect. 4

.2.2.
For a maximum reduction of computational costs, it is desirable to use both techniques at the same time. Using the well established Hill-Mandel FE2 methods, see e.g. [1,2], for multiscale modeling in combination with structural elements such as beams, plates and shells, leads to significant challenges, see e.g. [3][4][5][6][7]. Solutions for these difficulties have only recently been introduced for shells in [8,9] and for beams in [7]. In [8] a higher order homogenization scheme was introduced for shells with shear soft kinematics and successfully tested on complex structures. However, the authors did not provide any benchmark results for well known properties of homogeneous structures. In [9] a classic Hill-Mandel FE2 approach was considered for shells, but lead to transverse shear stiffness parameters, that were dependent on the size of the RVE, as was shown in [10]. In [7] the classic Hill-Mandel FE2 approach was applied to shear soft beams and additional inner constraints were introduced via interface elements on the micro/meso level. This method leads to well established cross-sectional values for linear elastic bench-mark tests and shows results in accordance with full scale models for complex structures with nonlinear geometry and materials. However, the method relies on interface elements for additional constraints and periodic boundary conditions on the micro/meso scale. Both are not easy to use in a practical sense as they have special requirements regarding meshing of the structure. They need to be taken into consideration by the user and are not easy to meet, depending on the nature of the modeled structure.
The aim of this contribution is to show a homogenization scheme for shear soft structures such as beams and plates, that is not based on the Hill-Mandel condition. Instead, it uses the more recently developed homogenization scheme based on the Irving-Kirkwood theory, that was presented in [11,12]. This approach does not need any interface elements when applied to shear soft beams or plates and has no apriori requirements concerning boundary conditions, as the micro/meso scale is no longer loaded via boundary conditions but rather via an additional, global constraint that links strains between macro and micro/meso scale.
Main properties of the presented multiscale approach for beams and plates are: • The macro structure is modeled with a shear soft kinematic as Timoshenko beam and Reissner-Mindlin plate. Therefore, basic elements with six degrees of freedom for spatial beams and five degrees of freedom for plates, including membrane deformation, are used. • On macro and micro/meso level all strains are assumed to be small. A linear geometry and linear elastic materials are considered. • RVEs on the micro/meso scale are modeled in three dimensions and allow the modeling of arbitrary inhomogenieties. • Results for homogenized properties are independent of the RVE's length without any further adjustments. • Benchmark tests show perfect accordance with wellknown values, except for shear correction factors in more complex geometry. The latter still show good estimates and reasons for their deviation are identified and discussed. • Inhomogenieties in beam axis direction and reference surface directions are considered and reactions of the macro structure show good accordance with full scale reference solutions.

Beam and plate theory
We limit the theory to small strains and thus linear geometry. Furthermore, we make the following assumptions and simplifications for beams and plates: beams: • prismatic with straight, untwisted reference axis, • reference axis through centroids of cross-sections, • cross-sections are doubly symmetric, • no warping of cross-sections, • stresses S x x , S xy , S xz are small and can be neglected. plates: • reference surface is plane and the mid surface, • inextensibility in thickness direction • no warping of cross-sections, • in-plane strains of reference surface are included, • thickness normal stresses S zz are small and can be neglected.
We denote the linear strain tensor byε. To avoid confusion with beam and plate strains, denoted by ε B and ε P , we rename its Voigt's notation with respect to a cartesian base system For a beam, the relevant parts are denoted by Analogous to the strains, the stress tensor is denoted byσ and its Voigt's notation renamed to S V = σ V = [S x x , S yy , S zz , S xy , S xz , S yz ] T = [σ x x ,σ yy ,σ zz ,σ xy ,σ xz , σ yz ] T to avoid confusion with beam and plate stress resultants, denoted by σ B and σ P . Relevant parts of the stresses for beams are denoted by S B = [S x x , S xy , S xz ] T and for plates by S P = [S x x , S yy , S xy , S xz , S yz ] T .

Kinematics
Let B 0 be the body of interest (beam or plate) and {e x , e y , e z } an orthonormal base system. For a beam, the base system is chosen in a way that e x aligns with the undeformed beam axis. For a plate, e x and e y span the plane of the undeformed reference surface.
In accordance with the Timoshenko beam theory, displacementsū = [ū x ,ū y ,ū z ] T of any point in the beam with coordinates x = [x, y, z] T are given with respect to the displacements and rotations u B = [u x , u y , u z , β x , β y , β z ] T of the reference axis as (1) The six independent beam strains ε B = [ε x , γ xy , γ xz , κ x , κ y , κ z ] T B are defined and connected to the relevant strains: For plates, a Reissner-Mindlin kinematic is used. The displacements and rotations of the reference surface are given by [ū x ,ū y ,ū z ] T of any point within the plate with coordinates x = [x, y, z] T are then given bȳ The relevant strains are The eight plate strains ε P = [ε x , ε y , γ xy , κ x , κ y , κ xy , γ xz , γ yz ] T P are defined and connected to the relevant strains: We solely include the in-plane strains of the plate to showcase the potential of the presented theory for extension to shell structures. The in-plane parts are not an area of focus in this contribution, but it can be shown that the corresponding

Constitutive equations for the stress resultants
The beam stress resultants are normal and transverse forces as well as torsional and bending moments, which are arranged in a vector Here, A denotes the area of the beam cross-section. A linear elastic isotropic material model is used to compute the stresses Where E and G = E/[2(1 + ν)] denote Young's and shear modulus, respectively, and ν the Poisson's ratio. Introducing (10) into (9) yields The material matrix D B for the beam is modified to with shear correction factors y and z . They are introduced to correct the overly stiff response to transverse shear loading. The area moments of inertia are denoted by I y and I z while I T denotes the torsional moment of inertia.
In an analogous way, stress resultants for a plate with thickness h are defined as Again, a linear elastic isotropic material model is used to connect stresses and strains and then introduced into (13): The material matrix D P contains the entries where Analogous to the beam, the shear part D s is modified to incorporate a shear correction factor to correct the overly stiff response to transverse shear loading.

Basics
This section briefly describes the relevant parts of the findings presented in [12]. There, a homogenization framework for solid bodies was presented, that is based on the Irving-Kirkwood theory (see [14]). It does not depend on any a priori assumptions on boundary conditions or loading conditions of the micro/meso scale, but only assumes volumetric averaging between micro/meso and macro scale solely for basic extensive mechanical quantities mass, linear momentum and energy. A body B is described in two scales. On a first, macro scale it is assumed to be homogeneous. On a second, micro/meso scale, a lot more details are incorporated, including inhomogeneities or a micro/meso structure with characteristic lengths, that are small compared to those of the macro structure of the body. For details regarding size requirements between macro and micro/meso scale and inhomogenizations, see for example [15][16][17].
The body B takes up the region with boundary ∂ . Its macro scale representation is described in an ortho-normal base system {e y 1 , e y 2 , e y 3 }, its micro/meso scale representation accordingly in {e x 1 , e x 2 , e x 3 }. A point within the body can then be identified with its coordinates y and x in macro and micro/meso scale, respectively.
With ρ the mass density, v the velocity of a material point, t the boundary force vector on ∂ and b the body force vector, the continuity equation for mass and the balance of linear momentum can be formulated independently from each other for the whole body in the macro scale and for a small subsection P with boundary ∂P of the micro/meso scale. Entities of the macro scale are denoted by superscript (·) M , those of the micro/meso scale with (·) m .
The continuity of mass then reads as and the balance of linear momentum as d dt To connect the entities of both scales, an averaging theorem for mass density and linear momentum is introduced For this averaging theorem, it is essential, that the subsection P is a small surrounding area of the point y. It is in fact the RVE corresponding to y. The authors in [11] and [12] showed that this theorem connects both scales in a consistent manner and that they bring about averaging equations for the body force vector and stress tensor, that can be described for the quasi-static case as For entities of the reference configuration, analogous homogenization equations can be derived.

Macro to micro/meso transition
On macro level, the strong form of the boundary value problem is described by Where with boundary ∂ is the region occupied by the body. σ M is the stress tensor, b M denotes the body force vector. The boundary ∂ is divided into a part ∂ u with given boundary displacementsū and a part ∂ t with given boundary forcest. The union of both parts makes up the full boundary: On micro/meso level, the equilibrium equation for stresses and body forces holds, but as mentioned before, there are no a priori boundary conditions so far (apart from restraining rigid body motions): In order to define a load on the micro/meso level consistent with the macro level, a homogenization of the strain tensor analogous to the homogenization Eqs. (22)-(25) is assumed [12]: Since the micro/meso region P corresponds to a single point in the macro scale, the macro strainsε M are constant throughout P. Therefore, holds. This is essentially an additional constraint for the micro/meso scale boundary value problem, that can be introduced into the weak form of (29) in Voigt's notation via Lagrange parameters: Here, g(u, δu) denotes the weak form of (29) and δ(·) denotes the variation of the quantity (·). A vector of Lagrange parameters λ V ∈ R 6 is used to enforce (31). They clearly show to be stresses.
We assume E M V , E m V and λ V as unknown. Assuming E M V as an unknown quantity leads to a symmetric system matrix in the finite element formulation and also makes it possible to use already existing algorithms to generate the homogenized entities later on. The given values of E M V will then be provided by means of boundary conditions. Equation (32) shows the need for appropriate ansatz functions for λ V in order to fulfill the side conditions in a weighted average instead of a point wise manner. A point wise fulfillment of the side condition would mean a uniform strain across all of P and simply lead to Voigt's bound for homogenized entitites, which does not account for the internal structure and material distribution within P. In essence, the micro/meso scale needs to allow for more complex deformation than the macro scale to make the concept of homogenization worthwhile.
The complete weak form for the micro/meso scale without boundary and body forces then reads

Required element formulations
A standard isoparametric, displacement-based finite element formulation is invoked. We use the same basic scheme for brick, plate and beam elements. The region , that is occupied by the body of interest B, is approximated with a finite number of elements e : e . The superscript h refers to the finite element approximation and numel is the total number of elements. Geometry, displacements and rotations, as well as virtual quantities are interpolated with functions N K . We use linear and quadratic functions N K (ξ ) for beam elements, bilinear/biquadratic N K (ξ, η) for plates and trilinear/triquadratic Here, nel denotes the number of nodes per element and L describes a differential operator that connects the displacement vector to the strain vector. It is different for each application (3D body/ plate/ beam). For a 3D body (brick element) a displacement vector u V = [u x , u y , u z ] T with three translational degrees of freedom is used. The strains are ε V = E V and the corresponding differential operator L V reads For a beam, a displacement vector u B = [u x , u y , u z , β x , β y , β z ] T with three displacements and three rotations of the reference axis is used. The corresponding differential operator reads For a plate, a displacement vector u P = [u x , u y , u z , ϕ x , ϕ y ] T with three translations and two rotations of the reference surface is used. The two rotations ϕ x = −β y and ϕ y = β x are used for easier application of boundary conditions and readability of results. The corresponding differential operator L P reads The approximations (34)-(37) are inserted into the weak form in vector notation on the macro scale. Omitting the subscripts (·) V , (·) P and (·) B for a general representation, standard procedure leads to The quantities K , v and F are generated by assembling all element quantities k e , v e , f e into one corresponding system quantity each.
For beam and plate elements, a selectively reduced integration scheme is used to mitigate the effects of transverse shear locking. One fewer Gauss point is used per element dimension to integrate the shear parts of the element stiffness matrix k e than is used for the other parts.

Elements on micro/meso level
On micro/meso scale, 3D brick elements are used. For approximation of displacement quantities and micro/meso strains in the weak form (33), we use the method described in Sect. 4.1. The approximation of Lagrange parameters is dependent on the type of elements used on macro level as they need to account for the relations between beam or plate on macro level and 3D body on micro/meso level. We first develop an implementation for use with 3D bodies on macro level to show differences between the approaches presented here and in [12]. We do not showcase any numerical examples for 3D bodies as this was already done there.

3D bodies
When brick elements are used on macro level to describe a three dimensional body, macro and micro/meso scale both use the same strains within their respective elements. Therefore, the macro strain vector ε M V = E M V is directly available for use in the weak form on micro/meso level (33) and none of its entries are zero by definition. This means, an approximation for E M V and λ V with constant functions as in is feasible. The assumption of constant Lagrange parameters for the whole micro/meso model leads to a fulfillment of the side condition (31) only in an averaged manner as desired (see Sect. 3.2). Inserting those approximations together with (34)-(37) into (33) leads tō and Assembly of element quantities into system quantities leads tō As mentioned before, the system matrix shows to be symmetric. The additional degrees of freedom for ε M V and λ V are associated with additional global nodes, that are shared by all elements on micro/meso level, see Fig. 3. Degrees of freedom, that are used for storage of ε M V i , are fixed at their respective values.

Beams
When using beam elements on macro level, the macro strains E M V needed for the weak form on micro/meso level are not directly accessible as those elements internally use the beam strains ε B . However, they can easily be converted by using Eq. (4): As there are only three strains in E M V , that are non-zero, a meaningful side condition can only be invoked for those three. The following is introduced for a consistent notation: This leads to an expression of the weak form in the micro/ meso scale as In this representation, the Lagrange parameters are stresses. This leads to some problems, that need to be addressed. When using a beam on macro level, the argument of separation of scales only holds true for one dimension -the direction of the beam axis. The cross-section of the beam needs to be represented in full within the micro/meso scale because there is no separation of length scales for the remaining two dimensions (see Figure 4). This also means, that there cannot be constant Lagrange parameters λ i j associated with averaged stress entities, that accurately describe the beam's reactions to the imposed macro strains without major information loss. For example, a beam experiencing bending would react with a linear normal stress S x x across its width or height. An over the cross-section averaged normal stress, and subsequently the corresponding Lagrange parameter, would be zero and all information on the beam's reaction would be lost. Therefore, the Lagrange parameters cannot be assumed constant, but need appropriate global ansatz functions with associated parameters to accurately describe the beam's reactions. The ansatz functions need to be global in order to fulfill the side condition (31) only in an averaged manner.
In the stress resultants, there is a quantity that is constant for each cross-section and can easily be connected to the stresses. Therefore, they are used as new parameters and ansatz functions are built around them. Thus the vector λ B = [λ N , λ Q y , λ Q z , λ M x , λ M y , λ M z ] T is introduced and connected to the Lagrange parameters bȳ wherē This leads to constant normal stresses for a normal force and linear normal stresses for bending moments. A torsional moment leads to shear stresses, that are proportional to the distance from the beam axis. Torsional warping is not accounted for. This limits the applicability of the proposed theory to either warp-free cross-sections or problems, where torsion is negligible. Shear stresses caused by transverse forces are assumed to be quadratic in an attempt to account for shear warping without giving up the Timoshenko theory on the macro level. They are assumed quadratic as this is the easiest enhancement from constant shear stresses and fulfills the stress boundary condition for free boundaries. The effects of this assumption are studied in Sect. 5.
The presented approximations are now introduced into the weak form and analogous equations and quantities to Sect. 4.2.1 are derived: where . . .
Assembly of element quantities into system matrices leads to the expression

Plates
Using plate elements on macro level leads to very similar problems as using beam elements. Subsequently, analogous strategies can be applied. First, the plate strains are converted: E M P = A P ε M P . Second, the number of Lagrange parameters is set according to the number of non-zero strains and an analogous notation is invoked: In a third step, an ansatz for the Lagrange parameters is chosen, that is based on the stress resultants and quadratic shear stresses, such that λ P ∈ R 8 with wherē Finally, the micro/meso problem can be described with same equations as the beam by substituting quantities with subscript (·) B with respective quantities with subscript (·) P , see Eqs. (51) -(53).

Consistent boundary deformations
Using the described theory and approximations and FE interpolations without any boundary conditions other than to restrict rigid body motions may lead to inconsistent boundary displacements. This is due to the fact, that the additional constraints and therefore the kinematic relations between macro and micro/meso scale are only fulfilled in an averaged manner, when implementing the Lagrange parameters as described. As this average fulfillment includes the boundary of the micro/meso scale, boundary displacements are not subjected to any hard, point wise requirements, so they deform as if they are free boundaries. But the micro/meso model is essentially a subsection of the whole body and therefore its boundaries are not free. They are interfaces to the remainder of the body's continuum, except for the outer boundaries of the cross-section for the beam and the top and bottom surfaces for the plate, which are real outer boundaries and therefore free in their deformation. These interfaces to the body's continuum cannot undergo arbitrary deformations if basic physics are to be respected. Especially they cannot deform in a way, that would lead to holes or interpenetration within the whole continuum. Rather, they need to be periodic. Boundary deformations of opposing boundaries (if they are continuum interfaces) need to fit together without any holes or interpenetration.

Displacement boundary conditions
The most obvious way to make opposing boundaries deform periodically is to apply periodic or linear displacement boundary conditions. The former are rather complicated to implement and need special consideration when discretizing the micro/meso problem, which contradicts the aim of an easy to use, one for all homogenization scheme. This leads to the consideration of linear displacement boundary conditions (DBC). They are easily implemented and impose no requirements on the used discretization.
For 3D bodies, DBC are applied onto the continuum interface boundaries of the micro/meso scale using the following expressions for the boundary displacements u R For beams, shear deformations are applied to u R x only and the expression arises. However, if DBC are applied like this to u R y and u R z , they lead to wrong results for the bending and extensional stiffness, similar to the findings in [9]. As those boundary conditions only apply torsional deformations and there are no incompatible boundary deformations caused by torsional loads because torsional warping is neglected, they can be omitted. Accordingly, DBC are only applied to u R x : For plates, analogous arguments hold for applying DBC and therefore only u R x and u R y are subject to prescribed boundary displacements using

Scaling of Lagrange parameters
A less obvious way to achieve consistent boundary deformations, but one without the need for boundary conditions, is a slight modification of the ansatz for the Lagrange parameters. However, it can only be applied if beam or plate elements are used on the macro scale and special requirements for the internal structure of the material are met.
If the body consists of a layered structure with each consisting of a homogeneous material, the previously proposed ansatz of overall constant λ x x for a given ε M x leads to nonuniform strains across the layers, but constant normal stress  Figure 5). If the ansatz is modified in a way, that the Lagrange parameter is not overall constant, but only constant within each layer and the value dependent on the stiffness of the layer, the strains would become uniform across all layers and thus consistent boundary deformations would be achieved (see bottom row of Figure 5).
This method can obviously only be applied if there is no change of material along the direction of the beam axis or within the directions of the reference surface of the plate. Would there be such an inhomogeneity, there would then be a jump in the normal stresses at the material interface, that would violate the equilibrium conditions. For the same reason it cannot be applied to 3D bodies.
In detail, the necessary modifications to the Lagrange parameters affect only the ones related to stresses in beam axis direction and the directions of the reference surface of the plate. Those Lagrange parameters will be multiplied with the associated entries of the material matrix C ∈ R 6×6 in Voigt notation for a 3D continuum:

Micro/meso to macro transition
In order to generate the homogenized quantities from the solution of the micro/meso problem, we adapt the algorithm presented in [9]. The overall discretized, coupled problem is described by its weak form as in Here, i denotes the averaging quantity dependent on the macro problem (RVE volume for 3D body, RVE mid surface area for plate, RVE length for beam) and N G P denotes the number of Gauss points per macro element. For every macro Gauss point, there is one micro/meso problem described by its weak formḡ i . Using the FE discretization described previously and combining all degrees of freedom, system matrices and load vectors of each micro/meso problem together with i into new generalized entitiesv m ,K m andF m , the coupled problem is represented by The element stiffness matrix k eM and element load vector f eM of the macro system are dependent on the solution of the micro/meso problems. But the entities describing the micro/meso problem are independent from each other, so they can be solved independently in parallel computations. The overall structure of the equation describing the micro/meso problem is the same for 3D bodies, plates and beams or with element representation: For this representation, it is assumed, that all rigid body motions are already restrained by boundary conditions and all corresponding degrees of freedom are accounted for by static condensation. This leads to K being regular.
In cases, where there are any additional displacement boundary conditions, displacements v e can be split into a part v e R subjected to boundary conditions and a free part v e F . The latter can be connected to the overall displacement vector v using the assembly matrix a e while v e R can be connected to the macro strains using the matrix R I according to Sect. 4.3: and where Combining the two ε M and δε M rows and columns and introducing Lastly, unknown degrees of freedom for displacements and Lagrange parameters are combined and correspondinglȳ are introduced, which then leads to Using [δv T , δε M T ] i = 0, static condensation is possible. The unknown displacements and Lagrange parameters can be obtained from the first row with This can be inserted into the second row of (75) and the homogenized entities can be obtained by comparing coefficients The homogenized material matrix D M i and the vector of homogenized stresses or stress resultants σ M i at Gauss point i in element e of the macro problem are given by

Examples
The developed models are programmed in FEAP [18]. All computations are done with FEAP, if not explicitly stated otherwise.
In this section, we use abbreviations according to Table 1. Normalized values W norm and relative errors / relative deviations W rel are computed by Here W hom denotes a value received from the presented homogenization scheme and W ref a reference value obtained from literature results or a full scale model. For beams, only plane macro problems are investigated as torsional warping is neglected and all other deformations can be linear superposed for the different directions, due to the linear elastic material and linear geometry.

Influence of RVE length, boundary conditions and order of ansatz functions
In contrast to classic Hill-Mandel homogenization scheme, where torsional and shear stiffness of beams (and plates) are dependent on boundary conditions and RVE length, if no further adjustments are made (see [7,10]), the homogenization scheme presented here does not show any such effects. Solely the use of linear displacement boundary conditions leads to expected boundary effects, because they stand in opposition to the assumption of quadratic shear stresses by using a quadratic ansatz funktion for the corresponding Lagrange parameters. The influence of those boundary effects decreases with increasing RVE length (see the following examples in this section).
In Sect. 5.2, results show that elements with quadratic ansatz functions have much better convergence rates over linear elements (see Fig. 6 and 7 ). Therefore, only quadratic elements are used in all following examples.

Homogeneous, isotropic beam with square cross-section
A homogeneous isotropic beam with a square cross-section is used as a benchmark to demonstrate basic functionality of the proposed homogenization method. We consider a square cross-section with a width of b = h = 1 cm and a linear elastic, isotropic material with E = 21 000 kN/cm 2 and ν = 0.3. This leads to a reference For the shear correction factor , a lot of literature can be found, e.g. [19][20][21][22][23][24], in which some depend on the value of ν, but their deviation from ν = 0 is minimal for a square cross- section and ν = 0.3. Therefore, we use the universal = 5/6 as reference value. The proposed beam kinematics do not account for torsional warping, which is why I P = I y + I z is used for reference instead of I T . Table 2 summarizes results obtained with the newly proposed homogenization scheme.
The correct diagonal structure for D B is obtained. There are no erroneous coupling terms. Extensional and torsional stiffness are always exact. Bending stiffness is always exact for quadratic elements and converges correctly with increasing discretization if linear elements are used (Fig. 6). The shear stiffness converges with increasing discretization against the correct value (Fig. 7). Although there is a boundary effect caused by application of DBC, its influence decreases with increasing RVE length (Fig. 8). This is due to the DBC hindering shear warping displacements near the boundary. With larger RVE lengths this effect becomes negligible.
Comparison with results from our previously proposed homogenization scheme [7], which is based on a classical Hill-Mandel approach, shows, that only the torsional stiffness differs. The previous approach yielded the Saint-Venant torsional stiffness G I T = 1 130 kN cm 2 instead of G I P = 1 346 kN cm 2 .

Beam with layered cross-section
As a more complex example, a layered cross-section with geometry according to Fig. 9 is studied. The fraction ρ C = h C /h describes the ratio of core height h C to total height h = h C + 2h L and is varied as well as the stiffness ratio α = E C /E L of Young's modulus E C of the core and Young's modulus E L of the face layers. The Young's modulus of the face layers is fixed at E L = 1 000 kN/cm 2 and Poisson's ratio for face layers and core is chosen as ν = 0.3.

Homogenized material matrix
Considering a typical use of such structures, focus lies in studying results for extensional stiffness D 11 = E A, bending stiffness D 55 = E I y and shear stiffness D 33 = z G A.
Analytical values for the considered stiffness parameters D ii are given by Here, s L = (h C + h L )/2 and z can be computed according to [25] with and 3 15 3ρ 2 C + 9ρ C + 8 , These analytical values are also achieved with our previous approach based on classical Hill-Mandel homogenization presented in [7] and henceforth referred to as reference values. Table 3 summarizes results of the computations.  Again, D B is computed correctly as diagonal matrix and the values for E A and E I y are exactly the reference values ( Fig. 10 and 11 ). The boundary effect of the DBC configuration is stronger than for the homogeneous cross-section and much higher RVE length and discretization are needed for the deviation in z between DBC and Scal configuration to become negligible (Fig. 13). The DBC configuration is therefore deemed to be not feasible for this example and all further investigations are made within the Scal configuration.
The computed shear correction factor z shows to be near the reference values for all considered combinations of ρ C and α but does not meet it exactly (Fig. 12). The relative deviation between computed and reference value increases with increasing core height and increasing stiffness ratio to up to 16.6 %. For a detailed study into the reason for this deviation, see Sect. 5.3.2.

Multiscale model
In order to better understand the influence of the error in shear stiffness, a macro problem of a clamped beam is considered (see Fig. 14). It is loaded with a concentrated load in   its center. We study the maximum displacement w and the stresses S x x and S xz for two exemplary configurations with α = 0.5, ρ C = 0.5 (I) and α = 0.01, ρ C = 0.8 (II). Considering symmetry of the system, only half of the beam is simulated. The macro problem is discretized with 20 linear Timoshenko beam elements according to Sect. 4.1. For the micro/meso problem 24 quadratic brick elements are used per direction with eight elements over the height for each layer. For reference, we use a full scale model with quadratic brick elements and same discretization density as for the micro/meso problem (15 360 brick elements overall). Table 4 shows results for the displacement w of both models. The error in the displacement is slightly lower than, but clearly increasing with, the error in shear stiffness.
Taking a look at stress distributions S x x (z) and S xz (z) at x = 7.75 cm in Fig. 15 and 16 , the normal stresses show perfect accordance with the full scale model. The shear stresses however show a significant deviation. While the RVE can only show a quadratic shear stress distribution due to the ansatz in Eqs. (50) and (61), the full scale model shows a more complex distribution. For configuration (I), deviation from a quadratic parabola is small, but for configuration (II) it is much more significant and hence the error in shear stiffness is much larger. Figure 18 shows the geometry of a beam with soft inclusions, that shall be studied in this section. Main material and inclusions are again linear elastic, isotropic with Young's modulus E 1 = 21 000 kN/cm 2 and Poisson's ratio ν 1 = 0.3 for the main material and E 2 = 10kN/cm 2 , ν 2 = 0 for the inclusions.  As the beam is inhomogeneous in direction of the beam axis, the RVE can only be modeled by means of the DBC configuration instead of the Scal configuration. Considering results from previous examples, a convergence study is conducted considering RVE length and RVE discretization. This shows a RVE consisting of three inclusions (L RVE = 60 mm) discretized with quadratic brick elements with edge length l e ≈ 0.9 mm is sufficient. One inclusion is then discretized with eight elements per direction and the main material accordingly. The used RVE can be seen in Fig. 19.

Beam with soft inclusions
As macro problem, a clamped beam is chosen (see Fig. 17). The midpoint of the beam is displaced by w = 18 mm and reaction forces are studied. For reference, a full scale model with 27-node quadratic brick elements is used. The full scale model is discretized with the same density as the RVE (41 472 elements overall). The macro model of the FE2 study is discretized with 16 linear beam elements according to Sect. 4.1.  Table 5 shows results for reaction force and reaction moment. Both quantities show good accordance between full scale and FE2 model with only 1.10 -1.85 % lower quantities in the FE2 model. As the inclusions lower reactions by 7.5 % (F z ) and 6.6 % (M y ) in comparison to a beam consisting only of the main material, the FE2 model only slightly overestimates the weakening effect caused by the inclusions.

Plate with inclusions
As examples for a homogeneous and a layered plate lead to very similar results as the homogeneous and layered beam, they are not discussed in detail. Instead, a more complex example of a plate with periodically distributed cube shaped inclusions is considered. The geometry of one such inclusion and its surroundings are displayed in Fig. 20. The matrix material of the plate is chosen as E 1 = 21 000 kN/cm 2 and ν 1 = 0.3. The inclusion material is chosen as E 2 = 100 kN/cm 2 and ν 2 = 0.45. The plate has a thickness of h = 50 cm, the distance between two inclusions is 30 cm for both x-and y-direction.
Again, convergence studies considering RVE length and discretization are carried out. They show a RVE length of 100 cm in both directions and an element edge length of l e ≈ 5.6 cm to be sufficient when using quadratic elements. This means, one inclusion is discretized with three quadratic elements per direction. The discretized RVE is displayed in Fig. 22.  For the macro problem, a clamped plate under constant loading is studied (Fig. 21). Using symmetry arguments, only a quarter of the plate is simulated. For reference, a full scale model with quadratic brick elements is used, where the load is applied on the upper surface of the plate. It is computed in ABAQUS CAE 6.14 [26] with elements of type 3D20 (20-node brick element with quadratic ansatz functions), because of memory issues in FEAP due to the size of the model. Convergence studies lead to a discretization of the full scale model with 5.19 million degrees of freedom. The macro model of the FE2 model is discretized with 50 linear elements per direction (elements according to Sect. 4.1).
Results in Table 6 show good accordance between FE2 and full scale model and only minimal deviation of 0.27 % in the vertical displacement at the midpoint of the plate. Comparison with a homogeneous plate shows that the inclusions cause a 1.0 % increase of the vertical displacement in the reference models. The weakening effect is not as strong as for the beam with inclusions, but still measurable and the FE2

Conclusions
A first order homogenization scheme based on the Irving-Kirkwood theory is presented and applied to shear soft beam and plate structures. The micro/meso scale is loaded via a side condition, that enforces the macro scale strains to be the volumetric average of the micro/meso scale strains. No boundary conditions are needed to enforce the macro deformations on the micro/meso scale. Nontheless, minimal boundary conditions must be applied to prevent rigid body motions and additional ones can be applied for other purposes if necessary, e.g. to incorporate boundary conditions in an RVE located at the support of the macro structure. The side condition and 3D modeling of the RVE allow for broader capabilities than just simple cross-sectional modeling within the Timoshenko or Reissner kinematic of the macro scale. A simple extension to quadratic shear stresses over the cross-section within the RVE leads to decent results for shear correction factors, even for non-homogeneous cross-sections. Macro quantities like reaction forces, reaction moments, in-plane stresses and displacements show good accordance between the proposed multiscale modeling and full scale models.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.