Scale Effects in Orthotropic Composite Assemblies as Micropolar Continua: A Comparison between Weak- and Strong-Form Finite Element Solutions

The aim of the present work was to investigate the mechanical behavior of orthotropic composites, such as masonry assemblies, subjected to localized loads described as micropolar materials. Micropolar models are known to be effective in modeling the actual behavior of microstructured solids in the presence of localized loads or geometrical discontinuities. This is due to the introduction of an additional degree of freedom (the micro-rotation) in the kinematic model, if compared to the classical continuum and the related strain and stress measures. In particular, it was shown in the literature that brick/block masonry can be satisfactorily modeled as a micropolar continuum, and here it is assumed as a reference orthotropic composite material. The in-plane elastic response of panels made of orthotropic arrangements of bricks of different sizes is analyzed herein. Numerical simulations are provided by comparing weak and strong finite element formulations. The scale effect is investigated, as well as the significant role played by the relative rotation, which is a peculiar strain measure of micropolar continua related to the non-symmetry of strain and work-conjugated stress. In particular, the anisotropic effects accounting for the micropolar moduli, related to the variation of microstructure internal sizes, are highlighted.


Introduction
Complex composite materials are characterized by the presence of a heterogeneous and discontinuous internal structure, if observed at some length scales. A basic issue of mechanics of complex materials is the derivation of suitable constitutive models accounting for the presence of such internal structures (conventionally referred to as microstructures), as well as of its own evolution, due to plasticity, damage, fracture, contact with or without friction, and other nonlinear phenomena [1][2][3][4][5][6][7][8].
Discrete approaches, relying on direct modeling at smaller scales, were extensively used in the literature for simulating the mechanical behavior of composite materials, as well as other heterogeneous materials, including polycrystalline materials, jointed rock systems, and block masonry with periodic microstructures [9][10][11][12][13][14]. Such approaches, although numerically accurate in predicting the mechanical where ε ij and χ ij denote the (non-symmetric) strain and curvature tensors, respectively, and e ijk is the usual third-order permutation tensor. It follows that two stress measures work-conjugated to ε ij and χ ij must be considered in the balance equations, i.e., the non-symmetric stress and couple-stress tensors, denoted by σ ij and µ ij , respectively. Under the simplifying hypothesis of neglected body couples, the balance equations to be satisfied for each point of the micropolar body are σ ij,j + f i = 0, µ kj,j − e ijk σ ij = 0, (2) where f i denotes the body forces, while body couples are considered null. Moreover, from equilibrium considerations at the external boundary, the surface tractions t i and moment tractions m i are expressed in terms of σ ij and µ ij as t i = σ ij n j and m i = µ ij n j , respectively. The general linear anisotropic stress-strain relations of the micropolar continuum can be expressed as The total number of coefficients in Equation (3) is equal to 324; however, owing to the major symmetry requirements related to the existence of a well-defined strain energy function, the number of independent coefficients reduces to 171. Specific material symmetries imply further reduction of the number of elastic constants, and the representation of the above constitutive law can be simplified based on the symmetry properties of the considered microstructure. It can be shown for instance that, for centrosymmetric materials, the fourth-order tensors B ijkl and C ijkl traducing coupling between classical and curvature deformation effects vanish.
In the remainder of the paper, the attention is restricted to the case of a two-dimensional (2D) reduced model. Thus, each material particle of the continuum has only three degrees of freedom (DOFs), consisting in two in-plane displacement components, i.e., u 1 and u 2 , and one out-of-plane micro-rotation component, i.e., φ 3 . Coherently, the out-of-plane strain and stress components, as well as the in-plane curvature and couple-stress components, are not considered in the stress-strain relations, which assume the following matrix form: In the considered case of hyperelastic materials having major symmetries, the following relations hold: and the resulting number of independent elastic constitutive parameters turns to be 21.
For centrosymmetric materials B ijk = C kij = 0, the number of independent constitutive components is reduced to 13. Moreover, for orthotropic materials, which is the case of common running bond masonries, the additional conditions A 1112 = A 1121 = A 2212 = A 2221 = 0 and D 12 = 0 hold, and the number of material constants becomes 8.

Numerical Formulations for Anisotropic Micropolar Models
In this section, the theoretical formulations for the two numerical approaches employed for comparison purposes in the present work, i.e., the (classical) finite element and the strong finite element methods, are discussed. In particular, a non-standard anisotropic micropolar finite element was formulated and subsequently implemented within the finite element environment COMSOL Multiphysics ® , whereas the strong-form finite element method (SFEM) was implemented, for the same continuum model, via an in-house code written in MATLAB ® . In this section, the main implementation details of the above formulations are provided.

Finite Element Formulation
The governing equations of the micropolar linear elasticity problem in a 2D setting can be discretized via a standard displacement-based finite element approach, after introducing the following vectors collecting the relevant displacements, and strain components in the 2D framework: Similarly, the corresponding stress and couple stress vectors are introduced as follows: Via standard variational arguments, the boundary value problem defined in Section 2 can be reformulated in weak form, which is based on the virtual work principle for a 2D anisotropic micropolar continuum under the assumption of zero body couples and considering u and φ as a set of kinematically admissible displacement and rotation fields, respectively, such that where δ denotes the variational operator, f is the body force vector, and t and m are the traction and couple-traction vectors applied on the boundary Γ N . It is useful to note that the curvature vector χ contains the first-order partial derivatives of the micro-rotations. This means that the strong continuity requirement is not necessary and C 0 shape functions can be adopted for the dependent kinematic variables u and φ. It follows that the components of dependent variables u and φ in Equation (6) can be interpolated at any point in terms of their nodal values u and φ, referred to as primal unknowns of the present finite element formulation as where N u and N φ are the shape function matrices for u and φ. It is worth noting that different shape functions for displacements and micro-rotations are used, as shown in Reference [55]. Coherently, different Cosserat finite element formulations can be derived, by simply changing the order of these shape functions. In the present work, nine-node quadrangular elements were used, characterized by bi-quadratic and bi-linear interpolation functions for displacements and micro-rotations, respectively. It follows that all nine nodes possess displacement-type DOFs, whereas micro-rotation DOFs are referred only to the four corner nodes. The above shape function matrices can be collectively represented in the following matrix forms: where N u j (ξ 1 , ξ 2 ) and N φ j (ξ 1 , ξ 2 ) are the biquadratic and bilinear shape functions for the jth node, respectively. It is worth noting that these functions are expressed in terms of natural coordinates −1 ≤ ξ i ≤ 1.
Consequently, the micropolar strains defined in Equation (1) can be expressed in the following matrix form: where ∇ is the gradient operator in the 2D setting. If Equation (10) is substituted into Equation (12), the strain and curvature matrices can be expressed as where d = u T φ T T is the unknown vector of nodal (generalized) displacements. The matrices B ε and B χ collect the spatial derivatives of the shape functions N u and N φ . It is worth noting that the computation of the inverse of Jacobian matrix is required to evaluate these shape functions, since they are functions of natural coordinates. Therefore, the constitutive relations in Equation (4) can be expressed in the following compact form: where Note that, in the presence of centrosymmetric materials, the matrices D εχ and D χε , responsible for the coupling between classical and bending effects, turn to be zero as in the present work.
Finally, the discretized version of the weak statement governing the equilibrium problem of a micropolar body (neglecting also the body forces for the sake of simplicity) reads as δd T K d = δd T F for any virtual displacement δd, where denote the stiffness matrix and the nodal force vector of the adopted finite element for describing a 2D linearly elastic and anisotropic micropolar body. Moreover, a standard Gauss integration technique is adopted for the computation of the stiffness terms appearing in Equation (16). The above-described anisotropic Cosserat finite element was implemented within the COMSOL Multiphysics®environment, using the integrated Physics Builder functionality. Such a tool offers a user-friendly graphical interface for extending COMSOL's built-in finite element library without requiring any in-depth programming skill [56].

Strong-Formulation Finite Element Method
Strong-form finite element method is a numerical approach based on a domain-decomposition technique merged with an extremely accurate pseudo-spectral collocation approach. As with any other strong-form methodology, both governing and boundary equations have to be implemented and solved at the same time. This leads to continuous stress and displacement fields across the element boundaries. In particular, the aforementioned pseudo-spectral approach considered here is the well-known generalized differential quadrature (GDQ) method. The GDQ method is able to approximate functional derivatives as a weighted linear sum of the functional values at a given regular grid of points. This required regularity limits the application of the GDQ to rectangular (or regular) domains only. Therefore, in order to study complex domains of arbitrary shape, the domain-decomposition technique (with mapping) has to be employed. Initially, the geometrical domain is divided into macro-elements according to the given geometrical or material discontinuities. Inside each element, a GDQ grid is placed according to a mapping transformation which might follow Lagrangian-or NURBS-based (Non-Uniform Rational Basis Spline) mapping [50][51][52]. In the present paper, the irregularity is due to a concentrated load which is applied on a linear edge; therefore, Lagrangian mapping with four nodes is sufficient to have an exact geometrical approximation. It was observed that the most common source of errors in the geometrical mapping is due to the presence of curvilinear boundaries [53,54]. In most of the published references, the GDQ method is presented in one-dimensional (1D) form because the numerical structure does not change much when the problem is 2D. In fact, for any given 2D function f (ξ 1 , ξ 2 ), the following derivatives can be reported: where the repeated indices k and l indicate a sum up to the number of points N ξ for a fixed grid of N ξ × N ξ , and the indices i and j indicate the point on the grid along the two directions ξ 1 and ξ 2 , respectively. The indices n and m stands for the derivative order for the coordinate ξ 1 and ξ 2 , respectively. It is remarked that the weighting coefficients c ij along ξ 1 and ξ 2 are different if a different number of points and collocation in the two directions are considered; otherwise, they are the same. Thus, for the sake of generality, Equation (17) is reported in general form, but the same number of points and grid collocation are considered in the numerical applications below.
Equation (17) has to be implemented into a computer code matrix form of the equations to be solved. Thus, Equation (17) can take the form below.
where I is the identity matrix and c ξ 1 , c ξ 2 are the weighting coefficient matrices of Equation (17). Please note that the definitions in Equation (18) depend on how the collocation points are ordered in the computer code [52]. Equation (18) does not include the element mapping which modifies the algebraic equations according to selected procedure [52]; in fact, they are defined in the master (computational) element. Without losing generality, Equation (18) is mapped into Cartesian coordinates and can be presented as Now, Equation (19) can be used to carry out the discrete form of the governing equations as follows: A 1111 D (2) x 1 + 2A 1112 D (2) x 1 x 2 + A 1212 D (2) x 2 U 1 + A 1121 D (2) x 1 + (A 1122 + A 1221 )D (2) x 1 x 2 + A 2212 D (2) x 2 U 2 + B 111 D (2) x 1 + (B 112 + B 121 )D (2) x 1 x 2 + B 122 D (2) The boundary conditions in terms of displacements are where N 1 = n 1 I, N 2 = n 2 I. The boundary stresses in terms of displacements take the forms below.
The complete set of governing (Equations (20)-(22)) and boundary (Equations (24)-(26)) equations can be collected (after re-ordering) into the following algebraic system reported in matrix form KU + Q = 0 by separating boundary and domain grid points as where K bb , K bd indicate the discrete form of the boundary equations, K db , K dd represent the matrices of the fundamental equations, Q b , Q d are the vectors of boundary and domain forces, and U b , U d are the rearranged vectors of the displacement parameters (U 1 , U 2 , Ω) divided into boundary and domain degrees of freedom. Finally, to improve the numerical stability, the static condensation is performed so that where Clearly, Equation (28) can be solved using any numerical tool using, for instance, Cholesky decomposition as used by MATLAB. It is from Equation (23) that the total number of degrees of freedom (DOFs) in each problem for the SFEM can be computed as 3n e N ξ − 2 2 , where 3 is the number of DOFs per grid point, n e is the number of elements in the mesh, and N ξ − 2 2 is the total number of grid points per element.

Numerical Simulations
The present study aimed to compare the two different numerical approaches adopted in the modeling of orthotropic micropolar continua. The problem illustrated below considers a square domain/wall of width/side L = 4 m, fixed at the bottom edge and subjected to several top loads acting on lengths of different size a according to three ratios a 1 /L = 1, a 2 /L = 0.5, and a 3 /L = 0.25. A general sketch of the present geometry is depicted in Figure 1a,b, representing the three half-wall geometries termed Case 1, Case 2, and Case 3 used in the computations with evidence of the top load and bottom boundary condition used. For the sake of comparison, the resultant of the top load is kept constant for all three cases above as F = 10 MN. Thus, the intensities of the equivalent distributed force for the three geometrical configurations are q 1 = 2.5 MPa, q 2 = 5 MPa, and q 3 = 10 MPa, respectively. It is remarked that the physical problem is studied for the three different configurations (Case 1, Case 2, and Case 3) for different values of the load size footprint a 1 , a 2 , and a 3 , associated with decreasing value of the ratio a/L, which is responsible for the "structural size effects". As such a ratio tends to a value much smaller than 1, the considered distributed load behaves as a concentrated top vertical force F acting at the symmetry line of the panel. The load sketches considered in Figure 1b imply the mesh patterns selected as reported in Figure 2, not for numerical convergence reasons, but for capturing the load discontinuity given. three half-wall geometries termed Case 1, Case 2, and Case 3 used in the computations with evidence of the top load and bottom boundary condition used. For the sake of comparison, the resultant of the top load is kept constant for all three cases above as 10 MN F = . Thus, the intensities of the equivalent distributed force for the three geometrical configurations are 1 2.5 MPa q = , 2 5 MPa q = , and 3 10 MPa q = , respectively. It is remarked that the physical problem is studied for the three different configurations (Case 1, Case 2, and Case 3) for different values of the load size footprint 1 a , 2 a , and 3 a , associated with decreasing value of the ratio a L , which is responsible for the "structural size effects". As such a ratio tends to a value much smaller than 1, the considered distributed load behaves as a concentrated top vertical force F acting at the symmetry line of the panel. The load sketches considered in Figure 1b imply the mesh patterns selected as reported in Figure 2, not for numerical convergence reasons, but for capturing the load discontinuity given.   Finite element mesh for the FEM, macro-element discretization, and mesh-free points for the SFEM are shown.
In order to have a comparison, the medium was modeled as a classical and micropolar equivalent continuum. The adopted material constants are summarized in Table 1 and were evaluated using the coarse-graining approach presented in Reference [26]. With  as the significant microstructure internal length of the composite assembly, dependent on the brick size, a scale ratio / L  governs the so-called "material size effects" of the considered panel. Three material cases In order to have a comparison, the medium was modeled as a classical and micropolar equivalent continuum. The adopted material constants are summarized in Table 1 and were evaluated using the coarse-graining approach presented in Reference [26]. With as the significant microstructure internal length of the composite assembly, dependent on the brick size, a scale ratio /L governs the so-called "material size effects" of the considered panel. Three material cases corresponding to running bond sequences of bricks of increasing size were analyzed, named Material 1 ( 1 /L = 0.005), Material 2 ( 2 /L = 0.05), and Material 3 ( 3 /L = 0.5) [26]. Note that, as the Cauchy model does not depend on the size, the corresponding constitutive parameters do not vary for these three material cases. Moreover, the variation of the brick size does not affect the independent micropolar constants A1111, A1122, A1212, and A2121, while the bending moduli, D11 and D12, are strongly affected and play a fundamental role in retaining memory of the original discrete microstructure. In all the effected simulations, for the sake of simplicity, the out-of-diagonal terms are set equal to zero, and this corresponds to neglect dilatant effects in the joints between the bricks of the 2D composite solid considered.
The final aim was to show the capability of the micropolar model to retain memory of the original composite behavior under the action of a load applied on a limited area, as well as to numerically investigate the related mechanism of diffusion, using both FEM and SFEM approaches.
Due to the symmetry of the problem, only half of the domain was numerically studied, and the discretizations correspond to the aforementioned geometries ( Figure 2) for both FEM and SFEM. The finite element discretization is based on the macro-element decomposition. For the SFEM, two representations are provided: the one with macro-elements (domain decomposition) and the other with both elements and grid points. With this discretization, it was easier to determine a closer match in terms of degrees of freedom between FEM and SFEM. It can be noted that, in Case 1, the load is uniform (a 1 /L = 1); thus, only two elements are needed for SFEM. Actually, SFEM could handle it with just one element, but that would lead to a larger mesh distortion in the FEM counterpart. On the contrary, the other two problems, Case 2 (a 2 /L = 0.5) and Case 3 (a 3 /L = 0.25), consider four elements with the load applied on the top-right element of the given macro-element mesh pointing downward.
In order to show the numerical stability of the present SFEM numerical approach, convergence analysis of Case 3 is presented with material configuration Material 1 ( 1 /L = 0.005) as indicated in Table 1. A Chebyshev-Gauss-Lobatto grid [49] was considered by varying the number of grid points inside each element. The vertical displacement of the point on the symmetry axis of the geometry Case 3 under vertical pressure was considered as a reference value. The "error" between the displacement, computed by varying the number of points and the reference one with 15 × 15 points, was considered and is plotted in Figure 3. A negligible "error" was measured for such a 15 × 15 mesh, but it is clear that, for a lower number of points, larger differences emerge. Evidence of the numerical stability is clear from Figure 3, where it is clearly shown that the numerical technique becomes stable for 15 × 15 points. Therefore, such a selection was considered thereafter. Nevertheless, the presentation of the detailed numerical accuracy of both FEM and SFEM is out of the scope of the present study; thus, it is only marginally presented. In the authors' previous work [48], it was already shown that SFEM reaches numerical stability prior to FEM due to its strong-form nature. Moreover, stress plots in SFEM are continuous among elements; instead, FEMs have stress jumps because stresses are post-computed in the integration (interior) points of the element and not enforced a priori as in SFEM. For letting FEM contour maps be continuous as in SFEM, stress interpolation was applied among elements. It was observed in the literature that Chebyshev-Gauss-Lobatto grid [49] distribution provides the most accurate results and a uniform convergence by increasing the number of grid points. FEM models have uniformly distributed mesh elements, in particular, 21 × 21 for Case 1, and 15 × 15 for Case 2 and Case 3 for each macro-element. According to the elements used and by considering that FEM models are made of nine-node elements, the number of DOFs for Case 1 is 1014 for SFEM and 7503 for FEM, and the number of DOFs for Case 2 and Case 3 is 2028 for SFEM and 7339 for FEM.
Chebyshev-Gauss-Lobatto grid [49] distribution provides the most accurate results and a uniform convergence by increasing the number of grid points. FEM models have uniformly distributed mesh elements, in particular, 21 21 × for Case 1, and 15 15 × for Case 2 and Case 3 for each macro-element. According to the elements used and by considering that FEM models are made of nine-node elements, the number of DOFs for Case 1 is 1014 for SFEM and 7503 for FEM, and the number of DOFs for Case 2 and Case 3 is 2028 for SFEM and 7339 for FEM. The results are presented in terms of the two continuum models (classical or micropolar), for the three geometries ( Figure 1b) and material cases considered (Table 1). Therefore, a total of 12 configurations are provided. It can be noted that the bending moduli, responsible of the scale effects, also indirectly affect the relative rotations between the local rigid rotation, i.e., the macro-rotation , and the micro-rotation, φ , which is associated with non-symmetric angular strain components; therefore, they also have influence on the anisotropic features of the constitutive relations. Therefore, even if we did not consider a direct variation in the material symmetry (orthotropy) of the microstructure, the scale effect indirectly affects the results in terms of relative rotation. Conversely, the shear parameters, 1212 A and 2121 A , are directly associated with the resulting micropolar orthotropy of the panel, meaning that their values significantly influence the relative rotations. The direct variation in the material orthotropy is the object of a future work.  The results are presented in terms of the two continuum models (classical or micropolar), for the three geometries ( Figure 1b) and material cases considered (Table 1). Therefore, a total of 12 configurations are provided. It can be noted that the bending moduli, responsible of the scale effects, also indirectly affect the relative rotations between the local rigid rotation, i.e., the macro-rotation θ = 0.5(u 2,1 − u 1,2 ), and the micro-rotation, φ, which is associated with non-symmetric angular strain components; therefore, they also have influence on the anisotropic features of the constitutive relations. Therefore, even if we did not consider a direct variation in the material symmetry (orthotropy) of the microstructure, the scale effect indirectly affects the results in terms of relative rotation. Conversely, the shear parameters, A 1212 and A 2121 , are directly associated with the resulting micropolar orthotropy of the panel, meaning that their values significantly influence the relative rotations. The direct variation in the material orthotropy is the object of a future work.  Figure 4 shows the comparison between the two numerical approaches in terms of the vertical displacement, u 2 , along the symmetry axis of the problem. Figure 4a represents the classical material configuration, together with the three geometrical cases, whereas Figure 4b-d report each case separately by varying the micropolar material in each figure. From Figure 4a, it can be noted that a reduction of the load footprint size leads to a local stress concentration (as deduced from the nonlinear behavior of the vertical displacements for Cases 2 and 3), which is absent in Case 1 (characterized by a uniform vertical displacement state). FEM solutions are indicated with solid lines and SFEM ones are indicated with markers. It is evident that the two solutions accurately match; therefore. it can be said that both numerical tools are able to describe these phenomena. As expected, Case 1 (a 1 /L = 1) (Figure 4b) has a uniform top load and a uniform fixed boundary at the bottom and gives the same results for any (Cauchy or micropolar) material configuration, highlighting that the micropolar model shows different solutions from the classical elastic approach only in the presence of force concentrations, whereby it is able to activate rotational DOFs and also produce relative rotations between macro-and micro-rotations. In other words, Case 1, without any stress concentration, is not sensitive to the "scale effect" due to a change in the material constants D 11 and D 22 , and the results for Cosserat and Cauchy materials are coincident in terms of displacements and stresses. case separately by varying the micropolar material in each figure. From Figure 4a, it can be noted that a reduction of the load footprint size leads to a local stress concentration (as deduced from the nonlinear behavior of the vertical displacements for Cases 2 and 3), which is absent in Case 1 (characterized by a uniform vertical displacement state). FEM solutions are indicated with solid lines and SFEM ones are indicated with markers. It is evident that the two solutions accurately match; therefore. it can be said that both numerical tools are able to describe these phenomena. As expected, Case 1 ( 1 1 a L = ) (Figure 4b) has a uniform top load and a uniform fixed boundary at the bottom and gives the same results for any (Cauchy or micropolar) material configuration, highlighting that the micropolar model shows different solutions from the classical elastic approach only in the presence of force concentrations, whereby it is able to activate rotational DOFs and also produce relative rotations between macro-and micro-rotations. In other words, Case 1, without any stress concentration, is not sensitive to the "scale effect" due to a change in the material constants 11 D and   On the contrary, Case 2 and Case 3, reported in Figure 4c,d, respectively, differ among themselves as a function of the configuration (classical or micropolar) and material properties (Material 1 ( 1 /L = 0.005), Material 2 ( 2 /L = 0.05), and Material 3 ( 3 /L = 0.5)). The effect due to the applied load is obvious: concentrated loads (a 2 /L = 0.5 for Case 2 and a 3 /L = 0.25 for Case 3) give larger displacements, whereas the material configurations give different results if Case 2 and Case 3 are observed. In particular, the increase of the internal material length, through D 11 and D 22 , which is an intrinsic feature of the micropolar constitutive behavior, alleviates the stress concentration at the top, thus reducing the vertical displacements with respect to the classical case of Figure 4a. Figures 5-7 represent the vertical displacement u 2 , vertical stress σ 22 , and relative rotation,

Cauchy model Micropolar model
for Case 1 using both FEM and SFEM. This case has a uniform top load and a uniform clamped boundary condition at the bottom so that no load or geometrical discontinuity is present. Therefore, as expected and already observed above (Figure 4b), there is no difference among classical and micropolar solutions. The given results, for all cases, show a linear displacement field ( Figure 5) with maximum values at the top where the pressure load is applied and zero values at the bottom where the boundary condition is enforced. Due to the particular boundary and loading conditions for the present case, the vertical stress field σ 22 is constant in the whole domain ( Figure 6). Please note that, due to very small numerical oscillations, the color map in the representation is not of a single color but two colors, or else the same scattered contour lines might appear. This does not change the fact that the whole stress map is of a constant value σ 22 = −2.5 MPa as indicated. Finally, similar comments can be reported for Figure 7, where relative rotation, θ − φ, is shown. The relative rotation for Case 1 is zero, as expected, because there is no micropolar effect in structures without (geometric or material) discontinuities. It is also evident from these plots that not only are classical and micropolar cases equal but so are the solutions provided by FEM and SFEM, in terms of contour plots. The relative rotation is not shown for the classical model because, for such a case, the micropolar rotation φ is not present, as only the macro-rotation can be computed in the Cauchy model.    In Case 2 and Case 3, the effect of concentrated loads and difference materials is shown. Case 2 is depicted in Figures 8-10. For all cases, FEM and SFEM also agree very well. It should be pointed out that, in terms of displacements for the Cauchy case (Figure 8a), the map has contour lines which tend to lay horizontally on the whole medium except for the zone directly influenced by the load, whereas micropolar materials (Figures 8b-d and 8f-h) have contour lines which are wider-spread in the domain. It is remarked that Material 1 ( 1 /L = 0.005) has higher displacements close to the area where the load is applied and almost zero displacement near the free boundary (left boundary), whereas Material 3 ( 3 /L = 0.5) has the lower displacement value with strongly not-zero displacement on the free boundary, in agreement with the results in Figure 4c,d and Figure 4g,h. In fact, Material 3 is the one with the higher micropolar effect, due to higher values of the D 11 and D 22 parameters, resulting in the significant reduction of strain (and stress) gradients. The contour plot solutions in terms of stress σ 22 (Figure 9) reflect the aforementioned results in terms of the displacement components. In fact, the highest stress contour map is given by the classical configuration, and the micropolar Material 3 has the lowest stress field with the small stress distribution within the medium. In other words, Material 3 is associated with the most significant reduction of stress (and strain) concentrations. The micropolar material models are able to re-distribute the stresses within elastic media better than classical elasticity [48]. This effect is tailored by the micropolar mechanical properties and is stronger when D 11 and D 22 have higher values, which are associated with larger ratios, where /L is the material intrinsic length. The micropolar effect can be easily observed from the relative rotation contour plots (Figure 10). In fact, all micropolar materials show a strong effect between the area where the load is applied and the free surface at the top boundary; in order words, the relative rotation has high values close to the discontinuity. Such an effect is much stronger for the micropolar solid with the highest bending constants (Material 3) with respect to the others. However, Material 1 and Material 2 show a relative rotation that goes down toward the boundary edge, whereas Material 3 has a relative rotation closer to the top free surface. In Case 2 and Case 3, the effect of concentrated loads and difference materials is shown. Case 2 is depicted in Figures 8-10. For all cases, FEM and SFEM also agree very well. It should be pointed Analogous comments can be reported for Case 3 (a 3 /L = 0.25), where the concentrated load is much stronger (applied on a smaller area) in Figures 11-13. Obviously, the corresponding values of vertical displacements u 2 , vertical stresses σ 22 and relative rotations,θ − φ, are higher than the previous cases, even though the resultant force has equal magnitude. This is due to the fact that the geometrical discontinuity is stronger in Case 3 ( Figure 4). Moreover, it can be observed that the gradients represented in Figures 11-13 are stronger with respect to Figures 8-10, as confirmed by the reported contour lines. Figure 13 in particular shows that Material 1 (associated with the smallest scale ratio 1 /L = 0.005) presents the most remarkable anisotropic features, providing the greatest peak for the relative rotations (related to the antisymmetric part of the strain tensor), as well as for the strain (curvature) and stress components. In this case, the load diffusion is less evident, and the vertical deflection is larger than in the other considered cases, thus highlighting smaller overall flexural stiffness properties. Another observation can be made on Figure 13, related to the relative rotation, whereby contour lines are slightly rotated with respect to Figure 10. This is due to the fact that the geometrical discontinuity in Case 3 is not aligned with the center of the half-medium presented here as in Case 2. In fact, in both cases, the relative rotation for Material 1 and Material 2 points toward the center of the medium; however, since in Case 2 the discontinuity is already aligned with the center, the contour plot shape is almost symmetric, whereas, for Case 3, the latter is slightly deformed.
The results here obtained are in agreement with the results already presented in previous works [20,26] and the extended simulations contained herein. The comparison with new experimental results will be the object of a future work. displacement on the free boundary, in agreement with the results in Figures 4c,d and 4g,h. In fact, Material 3 is the one with the higher micropolar effect, due to higher values of the 11 D and 22 D parameters, resulting in the significant reduction of strain (and stress) gradients. The contour plot solutions in terms of stress 22 σ (Figure 9) reflect the aforementioned results in terms of the displacement components. In fact, the highest stress contour map is given by the classical configuration, and the micropolar Material 3 has the lowest stress field with the small stress distribution within the medium. In other words, Material 3 is associated with the most significant reduction of stress (and strain) concentrations. The micropolar material models are able to re-distribute the stresses within elastic media better than classical elasticity [48]. This effect is tailored by the micropolar mechanical properties and is stronger when 11 D and 22 D have higher values, which are associated with larger ratios, where L  is the material intrinsic length. The micropolar effect can be easily observed from the relative rotation contour plots ( Figure 10). In fact, all micropolar materials show a strong effect between the area where the load is applied and the free surface at the top boundary; in order words, the relative rotation has high values close to the discontinuity. Such an effect is much stronger for the micropolar solid with the highest bending  σ and relative rotations,θ φ − , are higher than the previous cases, even though the resultant force has equal magnitude. This is due to the fact that the geometrical discontinuity is stronger in Case 3 ( Figure 4). Moreover, it can be observed that the gradients represented in Figures 11-13 are stronger with respect to Figures 8-10, as confirmed by the reported contour lines. Figure 13 in particular shows that Material 1 (associated with the smallest scale ratio 1 0.005 L =  ) presents the most remarkable anisotropic features, providing the greatest peak for the relative rotations (related to the antisymmetric part of the strain tensor), as well as for the strain (curvature) and stress components. In this case, the load diffusion is less evident, and the vertical deflection is larger than in the other considered cases, thus highlighting smaller overall flexural stiffness properties. Another observation can be made on Figure 13, related to the relative rotation, whereby contour lines are slightly rotated with respect to Figure 10. This is due to the fact that the geometrical discontinuity in Case 3 is not aligned with the center of the half-medium presented here as in Case 2. In fact, in both cases, the relative rotation for Material 1 and Material 2 points toward the center of the medium; however, since in Case 2 the discontinuity is already aligned with the center, the contour plot shape is almost symmetric, whereas, for Case 3, the latter is slightly deformed.
The results here obtained are in agreement with the results already presented in previous works [20,26] and the extended simulations contained herein. The comparison with new experimental results will be the object of a future work.  σ and relative rotations,θ φ − , are higher than the previous cases, even though the resultant force has equal magnitude. This is due to the fact that the geometrical discontinuity is stronger in Case 3 ( Figure 4). Moreover, it can be observed that the gradients represented in Figures 11-13 are stronger with respect to Figures 8-10, as confirmed by the reported contour lines. Figure 13 in particular shows that Material 1 (associated with the smallest scale ratio 1 0.005 L =  ) presents the most remarkable anisotropic features, providing the greatest peak for the relative rotations (related to the antisymmetric part of the strain tensor), as well as for the strain (curvature) and stress components. In this case, the load diffusion is less evident, and the vertical deflection is larger than in the other considered cases, thus highlighting smaller overall flexural stiffness properties. Another observation can be made on Figure 13, related to the relative rotation, whereby contour lines are slightly rotated with respect to Figure 10. This is due to the fact that the geometrical discontinuity in Case 3 is not aligned with the center of the half-medium presented here as in Case 2. In fact, in both cases, the relative rotation for Material 1 and Material 2 points toward the center of the medium; however, since in Case 2 the discontinuity is already aligned with the center, the contour plot shape is almost symmetric, whereas, for Case 3, the latter is slightly deformed.

Conclusions
This work proposed a detailed numerical investigation of the scale effects in orthotropic composites, such as brick/block assemblies, modeled as micropolar continua, under the action of localized loads. The numerical solution of the underlying boundary value problems was performed using two different numerical methods, i.e., the finite element method and the strong-formulation finite element method (SFEM). The results show that FEM and SFEM approaches provide comparable results both for classical and micropolar materials. The effect of orthotropic micropolar mechanical properties was shown as a function of geometrical discontinuities and the material properties, derived from the description of the composite microstructure. In particular, two length scales were considered in the present study, i.e., the scale ratio a L between the characteristic size a of the loaded area and the overall structural size L , chosen as a suitable measure of the load concentration, and the scale ratio L  , with  being the internal material length that directly affects the scale-dependent micropolar bending moduli. It was observed that the relative rotation represents a significant measure of the micropolar effect of an elastic medium subjected to concentrated loads. For the trivial solution with constant top pressure (i.e., fixing 1 1 a L = ), no micropolar effect was also carried out and shown. In fact, for such a configuration, no difference could be observed between classical and micropolar models. By

Conclusions
This work proposed a detailed numerical investigation of the scale effects in orthotropic composites, such as brick/block assemblies, modeled as micropolar continua, under the action of localized loads. The numerical solution of the underlying boundary value problems was performed using two different numerical methods, i.e., the finite element method and the strong-formulation finite element method (SFEM). The results show that FEM and SFEM approaches provide comparable results both for classical and micropolar materials. The effect of orthotropic micropolar mechanical properties was shown as a function of geometrical discontinuities and the material properties, derived from the description of the composite microstructure. In particular, two length scales were considered in the present study, i.e., the scale ratio a/L between the characteristic size a of the loaded area and the overall structural size, chosen as a suitable measure of the load concentration, and the scale ratio /L, with being the internal material length that directly affects the scale-dependent micropolar bending moduli.
It was observed that the relative rotation represents a significant measure of the micropolar effect of an elastic medium subjected to concentrated loads. For the trivial solution with constant top pressure (i.e., fixing a 1 /L = 1), no micropolar effect was also carried out and shown. In fact, for such a configuration, no difference could be observed between classical and micropolar models. By tailoring the micropolar mechanical properties. the effect on orthotropic media changes both in terms of displacements and stresses. With reference to the cases of nonhomogeneous stress states (i.e., Case 2 a 2 /L = 0.5, and Case 3 a 3 /L = 0.25) associated with localized loads (similarly to the case of other geometry or boundary discontinuities), it was shown that the contribution of the relative rotation between the macro-and the micro-rotation can be increased by using higher values of the micropolar bending moduli corresponding to increasing values of the scale ratios, /L. Its beneficial role allows the diffusion of concentrated loads, thus alleviating the associated severe stress gradients, especially for Case 3 (i.e., that associated with the smallest ratio a 3 /L = 0.25), thus highlighting the capability of the micropolar model to retain memory of the underlying microstructure response under the action of localized loads.
Moreover, the micropolar bending moduli were proven to be influenced by the anisotropic features of the constitutive relations, which directly affect the micropolar shear parameters, kept fixed in the present study. In the cases of nonhomogeneous stress states under investigation, Material 1 (associated with the smallest scale ratio 1 /L = 0.005) showed the most remarkable anisotropic features. In this case, the load diffusion was less evident, the vertical deflection was larger than in the other considered cases, and the relative rotations (related to the antisymmetric part of the strain tensor) showed the greatest peak, thus highlighting smaller overall flexural stiffness properties.
From a numerical point of view, it was also shown that, even if in the cases analyzed the two solutions provide the same results, SFEM is a more reliable and simple technique to apply in terms of discretization procedure, because the mesh can be constructed starting from the discontinuities present in the physical problem. Many finite elements have to be used in all computations, and mesh refinements have to be applied for refining the numerical solutions. On the contrary, once the mesh is defined in terms of macro-elements in SFEM, a simple re-population in terms of grid collocation points can be easily performed. Moreover, SFEM needs only to employ derivative discretization, whereas FEM implements both integrals and derivatives, which leads to slower calculations in terms of computational cost.