Comparative Analysis of the Stiffness Characteristics of Masonry Wall Based on Experimental Data and an Orthotropic Model

1. Introduction

A critical challenge in designing composite materials lies in predicting their elastic properties. Various methods exist for predicting the elastic properties of composites [1,2,3,4,5,6], which can be broadly categorized into analytical and numerical approaches. Analytical methods, primarily represented by micromechanical models [7,8,9], provide explicit or semi-explicit solutions and thus require relatively low computational effort. However, these methods often rely on specific assumptions and simplifications. On the other hand, numerical methods, predominantly based on finite element analysis of representative volume elements [10,11,12,13], do not require special assumptions to obtain solutions but typically demand significant computational resources. A reliable analytical method is highly desirable in composite material design, as the design process usually involves iterative evaluations of numerous intermediate configurations before achieving the desired properties [14,15].

Consequently, predicting the mechanical behavior of macro-heterogeneous materials—accounting for their constituent properties and spatial arrangement—is of paramount importance. When designing structural elements using such materials, simple and robust computational methods are needed to determine their physico-mechanical properties, stress-strain conditions, and strength estimates. The problem of assessing the stress-strain state of structural elements made of heterogeneous materials [16,17] is addressed through solid mechanics, which often requires complex mathematical formulations. While extensive literature exists on composite material strength [18,19], a unified theory for formulating constitutive laws of structurally heterogeneous media remains lacking.

Research on the mechanical properties of brick masonry is crucial for reliable structural design and simulation. The work of Galyalyuk and Demchuk [20] presents both experimental and numerical data for three series of masonry specimens subjected to compressive loading. In this study, we reproduce and utilize these data to construct the stiffness matrix D_ij based on an orthotropic model, identifying patterns and differences between the specimen series.

Recent studies on masonry mechanics, including works presented in [21,22,24,25], highlight the importance of accounting for anisotropy and structural heterogeneity. As demonstrated in [22], proper consideration of boundary conditions in masonry modeling is essential for accurately describing its anisotropic properties: weak mortars require an increased number of courses in the model to mitigate edge effects. For reinforced masonry [24], it has been established that mesh reinforcement compensates for anisotropy, reducing the required number of modeled courses. The study [21] proposes a homogenization method based on numerical simulations of 9-course masonry fragments in Abaqus [22], enabling the identification of a representative volume element (RVE) and revealing the stiffness matrix asymmetry (D₁₂≠D₂₁).

Research [25] demonstrates that the stiffness ratio between brick and mortar (E_b/E_m) significantly influences strain distribution in masonry (Figure 1a). This is corroborated by calculations of the stiffness matrix D_ij (Table 2) and aligns with experimental data for the KPO-1–KPO-3 series. The results illustrate how E_b/E_m affects strain distribution: for KRO-1 (E_b/E_m=1.12), the difference Δ₁−Δ₂ is minimal, whereas for KRO-3 (E_b/E_m=2.58), maximum strain divergence is observed. This is quantitatively confirmed by the increase in the anisotropy coefficient D₁₁/D₂₂ from 1.13 to 1.25 (Table 2).

This study aims to develop the stiffness matrix D_ij for three distinct series of brick masonry specimens and validate the computational results against experimental data. A complementary objective involves evaluating the applicability of the Abaqus finite element analysis system for simulating the structural performance of masonry, which is critical for predicting its mechanical behavior under various loading conditions.

2. Masonry Materials

2.1. Brick

The study utilized solid ceramic bricks with dimensions of 250×65×120 mm. The brick material exhibited the following mechanical properties:

• Elastic modulus, E_b = 11850 MPa
• Shear modulus, G_b = 4740 MPa
• Poisson’s ratio, μ_b=0.113

The compressive strength of the bricks, determined according to EN 772-1 standard, was f_b = 44.1 MPa.

2.2. Mortar

Three types of masonry mortar with different strength characteristics were employed in the study:

• Series KRO-1:
  • Mortar strength, f_m=10.9MPa
  • Elastic modulus, E_m=10580 MPa
  • Shear modulus, G_m= 4232 MPa
  • Poisson’s ratio, μ_m=0.17.
• Series KRO-2:
  • Mortar strength, f_m=7.9 MPa
  • Elastic modulus, E_m=9210 MPa
  • Shear modulus, G_m= 3684 MPa
  • Poisson’s ratio, μ_m=0.19.
• Series KRO-3:
  • Mortar strength, f_m=3.1 MPa
  • Elastic modulus, E_m=4600 MPa
  • Shear modulus, G_m= 1840 MPa
  • Poisson’s ratio, μ_m=0.23.

The mortar was prepared using dry mix No. 111/11 M100, with adjusted proportions to achieve the required strength. The thickness of mortar joints was maintained at 10 mm for all series.

3. Research Methodology

Masonry Specimens

Figure 1b shows the experimental setup for compression testing of brick masonry specimens. The tests were conducted using a P-125 hydraulic press in compliance with the EN 1052-1 standard [28]. The specimens were subjected to stepwise loading applied perpendicular to the horizontal mortar joints. Displacements and strains were measured using mechanical dial gauges and strain gauges.

where:

Figure 1b (a) $-$ diagram of the arrangement of devices (I1-I4 – clock-type movement indicators IGP-1 with a division value of 0.001 mm; h_s = 515 mm; l_s = 510 mm; t_s = 120 mm)

Figure 1b (b) −diagram of the arrangement of devices (I1-I4 – clock-type movement indicators IGP-1 with a division value of 0.001 mm; h_s = 780 mm; l_s = 305 mm; t_s = 120 mm)

The current study employs the following methods:

a.

Experimental Investigations

Compression tests were conducted on masonry specimens of series KRO-1, KRO-2 and KRO-3 (Figure 1b), determining the following mechanical properties: elastic modulus E_i, Poisson’s ratios (μ_xy, μ_yx) and shear modulus G_xy, calculated using the formulas:

$$E_i={\displaystyle \frac{F_{i,max}}{3\cdot A_i{\epsilon }_i}},{\mu }_{xy}={\displaystyle \frac{{\epsilon }_x}{{\epsilon }_y}},{\mu }_{yx}={\displaystyle \frac{{\epsilon }_y}{{\epsilon }_x}},G_{xy}={\displaystyle \frac{{\tau }_{xy}}{{\gamma }_{xy}}}$$

(1)

where: F_i,max = ultimate failure load, A_i = cross-sectional area of the specimen, ε_i = relative longitudinal strain, “3” = load equal to one-third of the ultimate failure load, ε_x = transverse strain, ε_y = longitudinal strain.

b.

Analytical method for determining the stiffness matrix of anisotropic brick masonry.

This study employs an orthotropic model to construct the stiffness matrix D_ij based on the known elastic moduli E_x, E_y and Poisson’s ratios μ_xy, μ_yx obtained from references [20,23]. The masonry behavior is assumed to follow a linear orthotropic model under plane stress conditions, where the stress-strain relationship is expressed through the compliance matrix:

$$\left\{\begin{array}{c}{\epsilon }_x={\displaystyle \frac{{\sigma }_x}{E_x}}-{\displaystyle \frac{{\mu }_{yx}}{E_y}}{\sigma }_y\\ {\epsilon }_y={\displaystyle \frac{-{\mu }_{xy}}{E_x}}{\sigma }_x+{\displaystyle \frac{{\sigma }_y}{E_y}}\\ {\epsilon }_{xy}={\displaystyle \frac{{\sigma }_{xy}}{G_{xy}}}\end{array}\right.$$

(2)

To determine the components of the stiffness matrix DD, we solve the system of equations for stresses:

$$\left\{\begin{array}{c}{\sigma }_x=D_{11}{\epsilon }_x+D_{12}{\epsilon }_y\\ {\sigma }_y=D_{21}{\epsilon }_x+D_{22}{\epsilon }_y\\ {\sigma }_{xy}=D_{33}{\epsilon }_{xy}\end{array}\right.$$

(3)

Let us examine two independent loading cases:

Case 1: Axial Loading in x-direction

$${\sigma }_x\ne 0,{\sigma }_y=0=\left\{\begin{array}{c}{\epsilon }_x={\displaystyle \frac{{\sigma }_x}{E_x}}\\ {\epsilon }_y={\displaystyle \frac{-{\mu }_{xy}}{E_x}}{\sigma }_x\end{array}\right.$$

(4)

Case 1: Axial Loading in Y-direction

$${\sigma }_y\ne 0,{\sigma }_x=0=\left\{\begin{array}{c}{\epsilon }_x={\displaystyle \frac{-{\mu }_{yx}}{E_y}}{\sigma }_y\\ {\epsilon }_y={\displaystyle \frac{{\sigma }_y}{E_y}}\end{array}\right.$$

(5)

Formulation of the Equation System

By applying the boundary conditions to the governing equations, we obtain:

Case 2: Axial Loading Along the y-Direction (σₓ ≠ 0, σ_y = 0, τ_xy = 0)

$$\begin{array}{c}{\sigma }_x=D_{11}{\displaystyle \frac{{\sigma }_x}{E_x}}+D_{12}\left({\displaystyle \frac{-{\mu }_{xy}}{E_x}}{\sigma }_x\right)\\ 0=D_{21}{\displaystyle \frac{{\sigma }_x}{E_x}}+D_{22}\left({\displaystyle \frac{-{\mu }_{xy}}{E_x}}{\sigma }_x\right)\end{array}$$

(6)

Case 2: Axial Loading Along the y-Direction (σₓ = 0, σ_y ≠ 0, τ_xy = 0)

$$\begin{array}{c}0=D_{11}\left({\displaystyle \frac{-{\mu }_{yx}}{E_x}}{\sigma }_x\right)+D_{12}{\displaystyle \frac{{\sigma }_x}{E_y}}\\ {\sigma }_y=D_{21}\left({\displaystyle \frac{-{\mu }_{yx}}{E_y}}{\sigma }_y\right)+D_{22}{\displaystyle \frac{{\sigma }_y}{E_y}}\end{array}$$

(7)

Solution for the Stiffness Matrix Components D_ij

$$D_{11}={\displaystyle \frac{E_x}{1-{\mu }_{xy}{\mu }_{yx}}},D_{12}={\displaystyle \frac{{\mu }_{yx}{E}_x}{1-{\mu }_{xy}{\mu }_{yx}}},D_{22}={\displaystyle \frac{E_y}{1-{\mu }_{xy}{\mu }_{yx}}},D_{21}={\displaystyle \frac{{\mu }_{xy}{E}_y}{1-{\mu }_{xy}{\mu }_{yx}}}$$

(8)

We obtain the stiffness matrix D in the following form:

$$Dij=\left[\begin{array}{ccc}{\displaystyle \frac{E_x}{1-{\mu }_{xy}{\mu }_{yx}}}& {\displaystyle \frac{{\mu }_{yx}{E}_x}{1-{\mu }_{xy}{\mu }_{yx}}}& 0\\ {\displaystyle \frac{{\mu }_{xy}{E}_y}{1-{\mu }_{xy}{\mu }_{yx}}}& {\displaystyle \frac{E_y}{1-{\mu }_{xy}{\mu }_{yx}}}& 0\\ 0& 0& G_{xy}\end{array}\right]$$

(9)

c.

Numerical Modeling

Finite element modeling of a 7-course unreinforced masonry fragment was performed in Abaqus (micromodeling) using C3D20R elements, with a mesh size of 10 mm for bricks and 5 mm for mortar joints, including local refinement in stress concentration zones. The contact interaction was defined using surface-based cohesive behavior, while boundary conditions consisted of full fixation at the bottom rigid beam and uniformly distributed compressive loading on the top plate, applied incrementally with energy convergence control (ETOL < 1×10⁻⁶) perpendicular to 10 mm thick horizontal mortar joints. For macromodeling, an orthotropic material model with Tie constraints between rigid beams and masonry was implemented using a 20 mm mesh, followed by experimental validation (Figure 2).

Two elastic-domain computational models were sequentially implemented: (1) an Engineering Constants model (Formula 1) with parameters derived from experimental elastic characteristics, and (2) an orthotropic model with stiffness matrix D_ij (Formula 9). However, elastic-domain simulations failed to reproduce the experimentally observed stiffness matrix asymmetry (D₁₂ ≠ D₂₁). Accurate elastic-domain representation of this effect would require a fully anisotropic model with 36 independent stiffness components, whereas the experimental study only determined five key parameters (E_x, E_y, μ_xy, μ_yx, G_xy), making complete anisotropic characterization unachievable within the existing experimental dataset limited to elastic behavior.

The presented orthotropic model describes only the linear-elastic behavior of masonry, which is its main limitation; as experiments [20,23] show, in masonry with low-strength mortar (KRO-3), microcracks develop under loads exceeding 30–50% of the ultimate, leading to progressive stiffness degradation, and future studies are proposed to address this by introducing scalar damage parameters or employing multistage homogenization that accounts for the crack formation stage.

Finite element modeling in Abaqus was conducted using two approaches: (1) micromodeling of a 7-course masonry fragment (510×515×120 mm for loading perpendicular to joints and 305×780×120 mm for loading parallel to joints) with explicit brick/mortar representation, and (2) macromodeling of equivalent dimensions using an orthotropic material model. This resulted in four distinct models: two micromodels (one per loading direction) and two macromodels (one per loading direction), maintaining geometric consistency with experimental specimens.

Two types of macromodels were used in this study:

• Macromodel EXP – A homogenized model based on experimentally measured elastic properties (Young’s moduli E_x, E_y, Poisson’s ratios μ_xy, μ_yx and shear modulus G_xy) from reference [20].
• Macromodel D_ij – A refined model incorporating the asymmetric stiffness matrix (calculated via Eqs. 8–9), where D₁₂≠D₂₁.

In Abaqus, the contact is modeled using Cohesive Behavior with separation into normal (K_nn) and shear (K_ss) stiffnesses. The elastic stage is defined by the interface stiffness matrix (10), which relates normal and shear stresses to relative displacements [29]:

$$\overrightarrow{t}=\left[\begin{array}{c}{t}_n\\ t_s\end{array}\right]=\left[\begin{array}{cc}{k}_{nn}& 0\\ 0& k_{ss}\end{array}\right]\left[\begin{array}{c}{\delta }_n\\ {\delta }_s\end{array}\right]$$

(10)

where: t_n, t_s — normal and shear stresses at the interface, δ_n, δ_s — relative displacements (normal and shear), K_nn, K_ss — stiffness coefficients for tearing/compression and shear, respectively. To describe the behavior of “brick-mortar” contact surfaces, a micro-model approach is used that takes into account the equivalent stiffness of the interfaces. The stiffness characteristics of the contact surfaces are determined through the mechanical properties of the materials and the geometric parameters of the weld according to formula (11):

$$k_{nn}={\displaystyle \frac{E_b{E}_m}{h_m\left(E_b-E_m\right)}},k_{tt}=k_{ss}={\displaystyle \frac{G_b{G}_m}{h_m\left(G_b-G_m\right)}}$$

(11)

where h_m is the mortar thickness of 10 mm.

Table 1. Initial mechanical characteristics of materials for calculating contact stiffness in Abaqus.

Series	contact stiffness in Abaqus
	Eb (MPa)	Em (MPa)	Gb (MPa)	Gm (MPa)	${\mathit{\mu }}_{\mathit{b}}$	${\mathit{\mu }}_{\mathit{m}}$	${\mathit{k}}_{\mathit{n}\mathit{n}}$ (MPa)/мм	${\mathit{k}}_{\mathit{s}\mathit{s}}$ (MPa)/мм
KRO -1	11850	10580	4740	4232	0.113	0.17	987.2	394.9
КRО-2		9210		3684		0.19	413.4	165.4
KRO -2		4600		1840		0.23	75.2	30.1

Table 1. Initial mechanical characteristics of materials for calculating contact stiffness in Abaqus.

Series	contact stiffness in Abaqus
	Eb (MPa)	Em (MPa)	Gb (MPa)	Gm (MPa)	${\mathit{\mu }}_{\mathit{b}}$	${\mathit{\mu }}_{\mathit{m}}$	${\mathit{k}}_{\mathit{n}\mathit{n}}$ (MPa)/мм	${\mathit{k}}_{\mathit{s}\mathit{s}}$ (MPa)/мм
KRO -1	11850	10580	4740	4232	0.113	0.17	987.2	394.9
КRО-2		9210		3684		0.19	413.4	165.4
KRO -2		4600		1840		0.23	75.2	30.1

For the studied series of masonry (KRO-1, KRO-2, KRO-3), the rigidity parameters were calculated on the basis of experimental data:

4. Results and Discussion

4.1. Results

The results of calculating the stiffness matrix D_ij, presented in Table 2, were obtained on the basis of experimental data and shear moduli using the orthotropic model.

The shear modulus G_xy in the xy plane is calculated according to the data of [20], E_x, E_y: Elastic moduli in principal material directions, D₁₁, D₂₂: Diagonal components of the stiffness matrix.

Key Observations:

• Model Validation:

The ratios E_x/E_y and D₁₁/D₂₂ are nearly identical for each series (difference < 0.001), confirming the validity of the orthotropic model.

2.

Anisotropy Trends:

• Maximum anisotropy occurs in KPO-2 (E_x/E_y=1.264)
• Minimum anisotropy is observed in KPO-1 (E_x/E_y=1.134)

3.

Unexpected Stiffness Behavior:

For KPO-3, both E_x/E_y and D₁₁/D₂₂ exhibit a slight decrease compared to KPO-2, despite a more pronounced difference between E_x and E_y.

The discrepancy between experimental stiffness values and the orthotropic model (Eqs. 8–9) arises due to directional loading effects in experiments (measuring deformation modulus along load axes) versus the model’s comprehensive strain-stress coupling (including Poisson’s ratios μ_xy, μ_yx). Additionally, micromechanical heterogeneity at brick-mortar interfaces induces stiffness matrix asymmetry (D₁₂≠D₂₁), unaccounted for in classical orthotropy, necessitating modified homogenization approaches with asymmetry correction factors (e.g., α=D₁₂/D₂₁) and representative volume element (RVE) upscaling to mitigate edge effects.

Analysis and Comparison of Data in Table 2

The analysis revealed that KRO-1 specimens exhibited the highest stiffness characteristics in both X and Y directions, attributable to the use of high-strength mortar (grade M109). KRO-2 demonstrated intermediate values, while KRO-3 showed the lowest performance, consistent with experimental findings.

The D₃₃ shear stiffness component progressively decreased from KRO-1 to KRO-3, indicating degradation of brick-mortar interfacial mechanical behavior with reduced mortar strength. Comparative analysis confirmed the correlation between stiffness matrix components (D₁₁, D₂₂) and elastic moduli (Eₓ, Eᵧ) from reference [23], demonstrating a consistent trend: mortar strength reduction directly diminishes masonry stiffness. This effect was particularly pronounced for D₃₃, highlighting the critical role of material bonding in structural performance.

Using Table 2 data, anisotropy coefficients (k, n), parameter m, and averaged symmetric model Dsym were calculated through relationships (12, 13) from references [26,27]:

$$k=\sqrt{{\displaystyle \frac{E_x}{E_y}}},m={\displaystyle \frac{E_x}{G}}-2{\mu }_{xy},n=\sqrt{2k+m}$$

(12)

The parameters k and n, related to m through the equation m = n² - 2k, represent independent characteristics of an orthotropic plate. Their introduction simplifies stress calculations in formula (10) and generalizes results for plane strain problems through reduced coefficients, where:

• k quantifies the degree of stiffness anisotropy between orthogonal directions (x - along mortar joints, y - across joints), with k=1 indicating isotropic material (Eₓ=E_y) and k≠1 confirming anisotropy;
• m reflects the combined influence of shear stiffness (Gₓᵧ) and transverse deformations on anisotropy, where m=0 for isotropic materials and m>0 for masonry due to low shear stiffness of mortar joints;
• n serves as a comprehensive parameter combining longitudinal and shear anisotropy effects, with n=2 for isotropic cases and n≠2 demonstrating anisotropy.

Formula (13) introduces the averaged symmetric stiffness matrix D_sym that approximates the asymmetric components D₁₂ and D₂₁ of the classical orthotropic model through their arithmetic mean:

$$D_{sym}={\displaystyle \frac{D_{12}+D_{21}}{2}}$$

(13)

For an asymmetric stiffness matrix (D₁₂ ≠ D₂₁), a correction factor α is introduced, reflecting the degree of asymmetry:

$$a={\displaystyle \frac{D_{12}}{D_{21}}}$$

(14)

$a$ $-$ characterizes the degree of deviation from the classical orthotropic model.

Modified parameters (15) are introduced to account for the stiffness matrix asymmetry D₁₂≠D₂₁:

$$k_{mod}=\sqrt{{\displaystyle \frac{D_{11}}{D_{22}}}}a,m_{mod}={\displaystyle \frac{D_{11}}{G}}-({\mu }_{xy}+{\mu }_{\mathrm{y}\mathrm{x}}),n_{mod}=\sqrt{2k_{mod}+m_{mod}}$$

(15)

Analysis of the data presented in Table 3 leads to the following conclusions:

• Effect of stiffness matrix asymmetry:

The consideration of asymmetric stiffness matrix components (D₁₂ ≠ D₂₁) results in a significant 30-65% increase in anisotropy coefficients compared to the symmetric model.

2.

Critical case for KRO-2 series:

The most substantial deviations were observed for KRO-2 specimens, where the asymmetry coefficient α reaches 1.65 and the modified parameter nₘₒ𝒹 attains a value of 2.38.

3.

Theoretical implications:

The obtained results demonstrate the necessity to modify the classical Lekhnitsky approach [26] through the introduction of a correction factor α, which is particularly relevant for materials with pronounced structural heterogeneity (Figure 3a).

The error is calculated as:

$$error={\displaystyle \frac{D_{ij}-D_{sym}}{D_{sym}}}\times 100\%$$

(16)

For α, the error is not specified, since the symmetric model by definition assumes α=1.

Table 4 clearly demonstrates that the asymmetric model (D₁₂≠D₂₁) more accurately captures masonry behavior, particularly for parameter k where symmetric model errors reach 65%, confirming the necessity of accounting for anisotropy in stiffness calculations; while differences are less critical for parameter n, simplifying its estimation in engineering applications (Figure 3b), with series KRO-2 showing critical error levels (65%) that mandate asymmetric model implementation, as derived from stiffness matrix D_ij analysis.

The numerical simulation results presented in Figures 4–6; demonstrate the distribution of vertical compressive stresses and horizontal strains in perpendicular and horizontal mortar joints for three experimental series (KRO-1, KRO-2, KRO-3), obtained through: (1) micro-modeling, (2) EXP macro-modeling of a homogeneous anisotropic medium (based on experimental values from [20]), and (3) D macro-modeling of a homogeneous anisotropic medium (using the calculated stiffness matrix D_ij), conclusively demonstrating the significant influence of mortar joint strength on the formation of structural stress states.

In Figure 4a (loading perpendicular to the joints), a uniform distribution of vertical stresses and horizontal deformations is observed in the KRO-1 series masonry, which is attributed to the high strength of the mortar (E_b>/E_m=1.12), whereas in Figure 4b (loading parallel to the joints), the stress and strain distribution is less uniform but remains stable, confirming the anisotropy of the masonry. At the same time, as demonstrated in [22], to fully eliminate edge effects, modeling of masonry fragments with at least 9 vertical rows and 4 bricks in length is required, highlighting the importance of accounting for edge effects and the micromechanics of “brick-mortar” interactions for accurate prediction of masonry behavior, especially under complex loading conditions.

Table 5, Table 6 and Table 7 present calculated errors in vertical stress and horizontal strain determinations for KRO-1, KRO-2, and KRO-3 masonry under different load orientations relative to mortar joints, with the strain error (ε) calculation formula for each macromodel versus the micromodel given by:

$$error(\mathit{\%})={\displaystyle \frac{{\epsilon }_{macro}-{\epsilon }_{micro}}{{\epsilon }_{micro}}}$$

(17)

Stres error (σ):

$$error(\mathit{\%})={\displaystyle \frac{{\sigma }_{macro}-{\sigma }_{micro}}{{\sigma }_{micro}}}$$

(18)

The error quantification methodology employs: ${\epsilon }_{macro}$ (macromodel-calculated strain), ${\epsilon }_{micro}$ (micromodel reference strain), ${\sigma }_{macro}$ (macromodel stress), and ${\sigma }_{micro}$ (reference micromodel stress), with percentage deviation calculated analogously for both stress and strain components using the unified formula (17, 18). For the analysis, control isopoles 2, 6, and 10 were selected, uniformly distributed along the specimen height, to: (1) eliminate local edge effects in contact zones with support plates, (2) cover all characteristic stress state regions (bottom, middle, and top), and (3) ensure data representativeness by analyzing zones with established stress fields.

Analysis of Figure 5a and Figure 5b:

• Figure 5a (loading perpendicular to joints, KRO-2): A uniform distribution of vertical stresses (σ_y) is observed, with localized concentrations (up to 15%) in edge zones, while horizontal strains (ε_x) exhibit elevated values near boundaries (EXP errors: up to 29%, D-model: up to 18.3%), confirming the influence of edge effects, partially mitigated by the orthotropic model.
• Figure 5b (loading parallel to joints, KRO-2): Stress distribution is less uniform, with pronounced concentrations near vertical joints (σ_y errors up to 15.1%), and strains (ε_x) in edge zones show significant deviations (EXP errors: up to 35.7%, D-model: up to 28.6%), highlighting the critical role of brick-mortar contact micromechanics and the need for larger RVEs to minimize inaccuracies.

As seen in Figure 6a and Figure 6b, masonry with low-strength mortar (KRO-3) exhibits significant stress heterogeneity, which correlates with high macro-modeling errors (Table 7), confirming the need to account for the micromechanics of contacts in computational models; moreover, as shown in [22], accurate description of edge effects requires modeling masonry fragments with a minimum of nine rows, which is particularly important for specimens with low mortar strength.

Stress-strain visualization (Figures 4–6) reveals pronounced masonry anisotropy, evidenced by fundamentally different mechanical responses under loading perpendicular (Figure 4a, Figure 5a and Figure 6a) versus parallel (Figure 4b, Figure 5b and Figure 6b) to mortar joints, with KRO-1 (high-strength mortar) exhibiting uniform stress distribution (≤15% variation) while KRO-3 (low-strength mortar) shows acute stress concentrations (local peaks ≤40%) and marked strain field heterogeneity.

Comparative experimental-calculated data analysis (Table 5, Table 6 and Table 7) demonstrates satisfactory agreement (≤7% error) for KRO-1 using orthotropic modeling, whereas KRO-3 displays significant discrepancies (≤65%) due to unaccounted microstrains in weak joints, substantiating the need for modified orthotropic models incorporating stiffness matrix asymmetry and edge effect correction factors.

4.2. Discussion

Experimental investigations revealed significant stiffness matrix asymmetry (D₁₂ ≠ D₂₁) [21,24], demonstrating the necessity to account for micromechanical effects, including localized brick-mortar interfacial deformation processes and structurally-induced anisotropy. The classical orthotropic model’s assumption of elastic symmetry (D₁₂ = D₂₁) produced statistically significant deviations from experimental data, particularly pronounced in specimens with weak mortar joints (KRO-3), where joints acted as stress-concentrating weak links. While macromodeling (orthotropic approach) simplified computations, it failed to capture microstructural influences on joint stress distribution.

To enhance predictive accuracy, we propose a modified anisotropic model incorporating asymmetry corrections through an anisotropy coefficient α = D₁₂/D₂₁ (Formula 14) for masonry series (KRO-1, KRO-2, KRO-3), with optimal E_b/E_m ratios identified for edge effect mitigation.

The proposed orthotropic model with asymmetry correction provides accurate predictions for uniaxial compression of unreinforced brick masonry with medium-to-high mortar strength (series KRO-1–KRO-3). However, its applicability is limited by:

(a)

Material and geometric constraints: The model requires validation for reinforced masonry, thin-bed joints, or alternative brickwork patterns (e.g., English bond).

(b)

Loading scenarios: Nonlinear effects (cracking, creep) and multidirectional stresses were not considered.

(c)

Scale effects: The RVE size (4 courses for vertical loading) may not capture behavior in full-scale structures under complex boundary conditions.

5. Comparison of Voigt and Reuss Methods for Masonry Stiffness Assessment Considering RVE

In this section, a comparative analysis is presented of methods for evaluating masonry stiffness based on the classical Voigt and Reuss approaches, taking into account the determination of the representative volume element (RVE). These methods are widely used for predicting the effective mechanical properties of heterogeneous materials; however, their applicability to anisotropic structures such as brick masonry requires further investigation [30,31,32,33].

The conducted analysis demonstrates that classical methods based on the assumption of isotropy or symmetry of properties can lead to significant errors in masonry stiffness evaluation, particularly for specimens with low-strength mortar (KRO-3). To improve calculation accuracy, a modification of the orthotropic model has been proposed through the introduction of an asymmetry coefficient α = D₁₂/D₂₁, which enables accounting for micromechanical effects at the “brick-mortar” interfaces.

5.1. RVE Determination via Indicator Zones

The Representative Volume Element (RVE) was established using experimental indicator zone data (Figure 1b). For brickwork with brick dimensions 250×65×120 mm and 10 mm joint thickness, RVE selection followed:

(a)

Y-axis loading (vertical, perpendicular to joints): RVE encompasses 4 brick courses (indicator zones)

(b)

X-axis loading (horizontal, parallel to joints): RVE comprises 2 brick rows (indicator zones)

5.2. Stiffness Evaluation Methods

Two classical approaches were applied for effective stiffness calculation:

• Voigt Method: Assumes uniform strain across components, computing stiffness as volume-weighted average:

$$D^v=A_m{D}_m+A_b{D}_b$$

(19)

where $A_m$ and $A_b$ represent the volume fractions of mortar and brick, and $D_m$, $D_b$ denote the stiffness tensors of the isotropic mortar and brick materials, respectively. The area fractions within the RVE were determined for two loading directions:

• Y-axis (perpendicular to joints): $A_b$ = 86.7%, $A_m$ = 13.3%
• X-axis (parallel to joints): $A_b$ = 86.7%, $A_m$ = 13.3%

with calculations based on cross-sectional area ratios in the RVE according to the formula.

$$R_i={\displaystyle \frac{A_i}{A_{total}}},i=\mathrm{b}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{k},\mathrm{m}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{a}\mathrm{r}$$

(20)

where A_i − cross-sectional area of component in RVE, A_total = A_b+ A_m− total cross-sectional area сечения.

2.

Method Reuss:

Assumes uniform stress across components, calculating compliance as:

$$S^R=R_m{S}^m+R_b{S}^b,D^R={\displaystyle \frac{1}{s^R}}$$

(21)

where ${\mathrm{S}}^{\mathrm{m}},{\mathrm{S}}^{\mathrm{b}}$ − compliance tensors for mortar and brick

Stiffness tensors for isotropic materials are defined via Young’s modulus (E) and Poisson’s ratio (μ):

$$D_m={\displaystyle \frac{E_m}{\left(1+{\mu }_m\right)\left(1-2{\mu }_m\right)}}\left[\begin{array}{ccc}1-{\mu }_m& {\mu }_m& 0\\ {\mu }_m& 1-{\mu }_m& 0\\ 0& 0& {\displaystyle \frac{1-2{\mu }_m}{2}}\end{array}\right]$$

(22)

$$D_b={\displaystyle \frac{E_b}{\left(1+{\mu }_b\right)\left(1-2{\mu }_b\right)}}\left[\begin{array}{ccc}1-{\mu }_b& {\mu }_b& 0\\ {\mu }_b& 1-{\mu }_b& 0\\ 0& 0& {\displaystyle \frac{1-2{\mu }_b}{2}}\end{array}\right]$$

(23)

Table 8 compares two classical methods for evaluating masonry stiffness—the Voigt method and the Reuss method—against the calculated stiffness matrix components D_ij derived from experimental data and the orthotropic model. The goal is to assess the accuracy of these methods for different sample series (KRO-1, KRO-2, KRO-3) and identify their limitations.

The error for each method was calculated as the percentage deviation from the experimental D_ij values:

$$error={\displaystyle \frac{D_{method}-D_{ij}}{D_{ij}}}\times 100\%$$

(24)

Analysis of Table 8 reveals that the Voigt method demonstrates lower error for diagonal stiffness matrix components D₁₁, D₃₃, particularly for KRO-1 and KRO-2 series (error ≤ 5.6%), while the Reuss method shows better accuracy for off-diagonal terms (D₁₂, D₂₁) in KRO-3 series (error ≤44.9%). Both methods, however, significantly underestimate the matrix asymmetry (D₁₂, D₂₁), confirming the necessity to modify classical approaches to account for masonry anisotropy and micromechanical effects.

Figure 7 shows a graphical comparison of the errors obtained using two classical methods for estimating masonry stiffness — the Voigt method and the Reuss method — relative to the orthotropic model, with the errors calculated for the components of the stiffness matrix D_ij (Table 8) and expressed as percentage deviations.

Although the Voigt and Reuss estimates provide reasonable results for isotropic composites, their underlying assumptions of (a) uniform stress/strain fields and (b) material symmetry make them fundamentally unsuitable for modeling the orthotropic behavior of masonry, particularly the experimentally observed stiffness matrix asymmetry (D₁₂/D₂₁>1.3) arising from the micromechanics of interfacial boundaries.

The Figure 8 illustrates the dependence of the effective elastic moduli on the stiffness ratio E_b/E_m for the three series of specimens (KRO-1, KRO-2, KRO-3).

Both E_x and E_y decrease significantly as (E_b/E_m) increases, reflecting the progressive weakening of mortar relative to brick. Specifically, E_x decreases from approximately 11.2 to 6.5 GPa, while E_y decreases from about 9.8 to 5.1 GPa.

The stiffness matrix components obtained from Eqs. (8–9) (D₁₁, D₂₂) closely follow the experimental trends and remain consistently bounded by the classical Voigt and Reuss estimates. The Voigt method provides an upper bound that slightly overestimates stiffness, while the Reuss method gives a lower bound that underestimates it. The widening gap between the bounds with increasing (E_b/E_m) indicates a stronger influence of structural heterogeneity when weaker mortar is used, accompanied by a marked increase in anisotropy between the x- and y-directions.

Thus, the comparison confirms that the orthotropic model derived from Eqs. (8–9) accurately captures masonry behavior. Moreover, mortar stiffness is identified as a governing factor controlling both the magnitude of effective moduli and the degree of anisotropy. These findings support the applicability of the modified orthotropic model for numerical calibration (e.g., finite element simulations) and for predicting the in-plane mechanical response of masonry walls under compressive loading.

6. Conclusions

This study established fundamental relationships between mortar joint mechanical properties and anisotropic behavior of unreinforced brick masonry. Experimental and numerical results (including micro- and macro-modeling using orthotropic formulation) demonstrated significant stress distribution dependence on the brick-to-mortar elastic modulus ratio (E_b/E_m). The identified stiffness matrix asymmetry (D₁₂ ≠ D₂₁) necessitates accounting for structural heterogeneity and micromechanical processes at brick-mortar interfaces in computational models.

The highest anisotropy occurred in low-strength mortar specimens (KRO-3 series) with anisotropy coefficient D₁₁/D₂₂ reaching 1.253. KRO-1 series (high-strength mortar) exhibited more uniform stress distribution, indicating reduced anisotropy effects. Nevertheless, even in this case, classical orthotropic models assuming matrix symmetry showed statistically significant deviations from experimental data.

The research results highlight the critical need to account for the following factors to enhance masonry behavior prediction accuracy:

a)

Material anisotropy driven by brick-mortar stiffness disparity

b)

Edge effects, particularly pronounced in low-strength mortar specimens

c)

Elastic property asymmetry, necessitating classical orthotropic model modifications

d)

Real service conditions, including loading direction and masonry component interactions

The stiffness matrix asymmetry represents a fundamental physical characteristic of masonry rather than an artifact, demanding model refinements. While the correction factor α can currently be employed for engineering calculations, comprehensive understanding requires additional experimental and numerical investigations.

Future work should address:

• Extension to reinforced masonry and non-standard geometries.
• Incorporation of nonlinear material models for damage analysis.
• Experimental validation under shear and combined loading.