A multi-phase mechanical model of biochar–cement composites at the mesoscale

This study presents a five-phase mesoscale modeling framework specifically developed to investigate crack propagation and mechanical properties of biochar–cement composites. The multi-phase model includes porous biochar particles with precise geometric construction, sand aggregates, cement matrix, and interfacial transition zone adjunct to both the biochar particles and sand aggregates. The 3D porous biochar library was first proposed and established in this study, which could provide an external interface for describing different pore shapes, wall thicknesses, and pore areas. All the simulation results were experimentally validated using a digital image correlation. Through precise geometric modeling, the unique failure modes and timing of biochar particles within the mortar were identified. This is analogous to the “strong column– weak beam” concept, accounting for the enhanced ductility observed in the biochar–cement composites under compression test. This work can advance the geometric modeling of porous aggregates broadly and elucidate their mesoscopic failure mechanisms in cementitious materials, thus providing new insights for developing high-ductility

are essential (N.Wang & Adeli, 2014;Zavadskas et al., 2018).Within this context, biochar, a carbon-negative material recognized by the Intergovernmental Panel on Climate Change, has attracted considerable attention for partially replacing cement and/or aggregates in cementitious composites and redefining low-carbon construction materials.The existing research on biochar-cement composites has predominantly been confined to experimental analysis, focusing on their compressive strength (L.Chen et al., 2022;Gupta et al., 2018) and flexural strength (Ahmad et al., 2015).These studies have demonstrated the capacity of biochar to mitigate concrete cracking, making it a promising material for low-carbon construction.However, these results fall short of providing a detailed analysis of the full stress-strain curve under compression, which is essential for a holistic understanding of the mechanical responses of cement composites.An in-depth exploration of the inherent properties of biochar, such as crack formation and propagation, is still lacking.The absence of numerical models for biochar-cement composites further hinders a comprehensive understanding of their failure mechanisms.
For large-scale concrete structures, the neural dynamics model has been increasingly adopted to facilitate the optimal design of building components (Ahmadkhanlou & Adeli, 2005;Aldwaik & Adeli, 2016).Furthermore, the macroscopic mechanical properties of concrete can now be predicted using machine learning (Rafiei et al., 2016), significantly aiding the concrete mix design and reducing manpower and financial costs (Rafiei et al., 2017).In mesoscale modeling, the Lattice Discrete Particle Model is widely utilized to investigate the failure behavior of concrete (Cusatis et al., 2011;Shen et al., 2020).Concurrently, models such as the Discrete Lattice and Rigid-Body-Spring Model have been developed for applications in durability (Kuntal et al., 2021;M. Liang et al., 2022) and creep (M.Liang et al., 2024) evaluation, as well as 3D concrete printing (Chang et al., 2021).The development of mesoscale models has been facilitated by the evolution of geometric representations, starting from 2D polygon geometry (Z.M. Wang et al., 1999) and progressing to 3D ellipsoid geometry (Häfner et al., 2006) and spherical geometry (Wriggers & Moftah, 2006).With advancements in modeling technology and algorithms, there has been a shift toward irregular geometry that closely resembles the actual shape of sand aggregates (Naderi & Zhang, 2021).However, there is currently a lack of a methodology for establishing the intricate geometry of porous aggregates, which poses a significant challenge and remains the primary barrier to the mesoscale modeling of cement composites with lightweight aggregates such as biochar.Currently, available simulations conducted for porous aggregates, such as asphalt aggregates (Y.Wang et al., 2021) and coral aggregates (Guo et al., 2023), are commonly relied on homogeneous theories.They simplified the porous structure as solid geometry and incorporated a reduction factor in material strength to account for possible effects of porosity.Although this method, with adjustment of material properties, could align the stress-strain curves with numerical simulation data, it falls short in accurately capturing the strain state and failure mode of porous aggregates and the overall interactions between the porous aggregates and the surrounding matrix.Hence, generating a precise 3D geometry for porous aggregates is essential to accu-rately simulate the specific stress distribution and failure mechanisms and predict their overall impact on structural integrity.Such advanced modeling is a solid foundation for achieving a realistic and mechanistic understanding of the mechanical behavior of biochar-cement composites.
To address the aforementioned shortcomings and knowledge gaps, this paper introduces a novel approach to constructing porous aggregate entities, marking the first instance of complex porous structures being modeled and simulated within a mesoscale framework.By establishing a customizable porous aggregate library, the model allows for precise geometrical representation of diverse materials such as biochar aggregates, coral aggregates, and various types of lightweight aggregates through the adjustable parameters of pore shapes, pore walls, and pore areas.As a porous material, biochar has a non-negligible influence on the mechanical properties of biochar-cement composites (Maljaee et al., 2021).Hence, the geometric information established for biochar is essential to accurately simulate the mechanical behavior of mortar.This research incorporates a five-phase mesoscale modeling framework that includes porous biochar particles with precise geometric construction, sand aggregates, cement matrix, and interfacial transition zone (ITZ) adjunct to both the biochar particles and sand aggregates.The precise geometric parameters of the biochar particles were determined via field emission scanning electron microscopy (FE-SEM) and ImageJ, and the multi-phase model was analyzed using finite element method.The cohesive interface element (CIE) was used for the biochar particles (Wimmer et al., 1997;Zhu et al., 2023), where the parameters were obtained from the nanoindentation tests (Zickler & Schöberl, 2006).The concrete damage plasticity (CDP) was adopted in the matrix and ITZ region in this model.In addition, significant emphasis in this study was placed on establishing a congruent mesoscale material structure between real and virtual mortar specimens.The failure modes and stress-strain relationships in concrete are influenced by the specimen size (Bažant, 2000).To counteract the size effect in the numerical simulations and experiments, the volume and dimensions of the specimens are identical across both setups after careful evaluation (see Appendix A).For experimental verification, digital image correlation (DIC) technology was utilized to obtain the strain-stress responses and detect the crack initiation and propagation during compression tests.Finally, a comprehensive investigation was conducted to examine the influence of biochar placement and orientation within the mortar matrix on the mechanisms and extent of biochar damage, thereby providing new insights into the underlying mechanisms of failure modes.

Individual aggregate and ITZ generation
Due to the significant complexity and diversity of irregular aggregate shapes, 3-D polygon geometry is usually generated based on relatively regular shapes.In this model, the regular polyhedron (with 28 vertices and 48 faces) is adopted as the base geometry, and three random fluctuations (  ,   ,   ) are introduced to the vertices in polar coordinates.Then, the procedures below are followed to randomly generate aggregate shapes.
Step 1: Generate three random numbers, presenting the aggregate diameter r, azimuth angles θ, and φ.These three values are determined by the random fluctuations and follow a uniform distribution, which can be expressed as where  stands for a uniformly distributed random number between 0 to 1;   ,   ,   stand for the minimum radius of the aggregate grading segment and angle degree at each layer, respectively.
Step 2: Randomly generate 28 points based on regular tetra octahedron under polar coordinates (, , ), and transfer the polar coordinates to Cartesian coordinates (  ,   ,   ), which can be formulated as where   ,   ,   represent the radius and azimuth angles of th point, respectively.
Step 3: Connect each point layer by layer and generate a surface accordingly, followed by creating the solid geometry based on the shell model in ABAQUS using Python coding.
The random fluctuation and aggregate size significantly influence the convergency.If the random fluctuation is too large, the surface generated by the lines may overlap and cause geometric errors.The suggested random fluctuation value of radius and angles lies in (0-0.3) and (0-30), respectively.Figure 1, as an example, shows the polygonal aggregate geometry at different random fluctuation parameters.
The ITZ is formed between aggregates and cement, and it has been recognized as the weakest region in mortar (He et al., 2022).The thickness of ITZ is usually 10-50 μm (Barnes et al., 1979), considerably smaller than the aggregates used in this model, which range in size from 2.36 to 4.75 mm.If the thickness of ITZ is set as small as the reality, the number of elements will be enormous and increase computation difficulties.Based on the previous models, the thickness of ITZ is usually set between 0.1 and 0.5 mm (Kim & Abu, 2011), which is approximately 1/10 of the diameter of the aggregate (H.Chen et al., 2018).However, the selected ITZ thickness significantly exceeds the typical ranges observed in empirical studies.Existing observations indicate that the influence on computational outcomes is considerably diminished when the ITZ thickness is reduced to less than 0.3 mm (Maleki et al., 2020).Therefore, in this study, the ITZ is set as 0.07 of the aggregate diameters (i.e., less than 0.3 mm) to satisfy both modeling accuracy and computational efficiency.The generation of the ITZ region is based on the offset point method (Naderi & Zhang, 2021) and coded in Python.The modeling process mainly follows three steps: (a) find the centroid of the generated aggregate in ABAQUS, (b) offset a total of 28 vertices by 0.07 of aggregate diameter, and (c) sequentially generate lines, surface, and solid.

Biochar library generation
The granular-shaped biochar pyrolyzed from wood waste at 400  (Figure 2a-c).To better define the geometric properties, including wall thickness over cross-section diameter ratio, roundness, and pore area, approximately 200 pores were analyzed using ImageJ (Collins, 2007).The wall thickness of biochar lies mainly at 2.5−4 μm regardless of the pore area, and the pore area is related to the ratio of wall thickness to edge (Figure 2d).Another important geometric parameter, the roundness of pores, shows the exponential relationship to the ratio of long-to-short edge (Figure 2e).
It is important to note that these analyses were conducted on 2D images, which can introduce systematic effects such as angular deviations and potential measurement biases due to the projection of 3D structures onto 2D planes.
The SEM images are carefully selected with pores directly facing the camera to minimize these effects.
Voronoi cells are widely used for geometric construction and mesh generation (Park et al., 2024;Q. Zhang et al., 2024), and the generated shapes are very similar to the porous structure of biochar.Therefore, based on the experimental observations, the pore structure of biochar can be generated based on Voronoi cells as demonstrated in Figure 3. Four control parameters, namely, roundness (ℛ), porosity (№), circle radius (Ω), and thickness factor (ℱ), are defined to govern the random distribution of the biochar structure, including the wall thickness, roundness, and pore distribution as discussed above.First, a solid cube is generated as the reference geometry, and various Voronoi cells are meshed based on different values of roundness (ℛ) and porosity (№) in this domain as shown in Figure 3a.Next, two circles with random radii and positions are used to determine the pore area and distribution (Figure 3b).Meanwhile, the wall thickness factor (ℱ) is introduced to control the ratio of wall thickness to the long edge of the pores (Figure 3b).The generation of these Voronoi cell-based shapes is based on the geometric features of biochar visually confirmed under FE-SEM as illustrated in Figure 2a-c.In this study, the fine aggregate is replaced with granular biochar with particle sizes ranging from 2.36 to 4.75 mm, which was sieved before preparing the mortar.To ensure the simulations are realistic and aligned with the experimental setup, an innovative and first-ever biochar library is generated based on Voronoi cells, containing 12 different particle sizes and pore structures based on SEM images (Figure 3c).The particle size of each model aligns with that of the fine aggregate employed in the experiments, falling within the specified range for fine aggregate.It is essential to emphasize that, given the results of FE-SEM, the pores distribution of biochar is continuous in the length direction and does not show significant variations.Therefore, the 3D geometry is established from the cross-section of the pores.
The model-building process of the porous biochar can then be summarized into three steps: (a) exporting the porous assemblies from Rhino and Grasshopper into IGES files that are compatible with ABAQUS for further processing, (b) segmenting the imported geometries (square Voronoi cells) into circular shapes of varying radii to align with the sieved particle size distributions as illustrated in Figure 3c, and (c) using Python coding to locate the centroid of each assembly and move it to the origin position (0, 0, 0) for subsequent particle placement.The orientation of each biochar particle within the composite is controlled through a randomized algorithm, which assigns a random rotation angle during placement.
Subsequently, four different biochar volume contents (1%, 3%, 5%, and 10%) are chosen to assess the feasibility of complex geometry placement.Examples in Figure 4 demonstrate the placement and interaction details of the biochar particles, which are generated with a Fuller distribution (Wriggers & Moftah, 2006) ranging from 2.36 to 4.75 mm, measured by their cross-section diameter.After deposition, the particles are randomly rotated by a certain angle to better match and simulate real-world scenarios.These biochar granules and sand aggregates with different volume content are integrated into the cement matrix and simulated in ABAQUS.

Inserting cohesive elements
When introducing biochar as fine aggregate into the mortar, the ITZ region surrounding the biochar particles is  formed due to the high water retention capacity of biochar (Zhu et al., 2023).Due to the irregular geometric shape of biochar particles, it is challenging to establish an ITZ in the same manner as for sand.The cohesive element in ABAQUS, which is designed to simulate the crack propagation and delamination behavior at material interfaces, is particularly well-suited for modeling biochar particles (Z.Zhang et al., 2007).Therefore, a program specifically designed to handle the intersection of cohesive elements for complex geometric models is developed.This process is implemented through Python coding by modifying the INP file, and the specific operational procedure is documented as follows.
Step 1: Mesh the whole composites using the 4-node linear tetrahedron (C3D4) element type.Read the keywords "*Node" and "*Element" of the mesh information in the INP file, Step 2: Generate the cohesive element geometry set by iterating all the node and element information.
Store the node number, Cartesian coordinates (x, y, z), and element connectivity information in their respective Python dictionaries.
Step 3: Number the nodes through the variable "Node appearance" in the INP file, which defines the occurrence count of common nodes for elements.This is important as cohesive elements are generated between adjacent element faces and require multiple splits of common nodes (Figure 5a).The following formula is used to record the node numbers: where l represents the total number of digits for node numbering, δ represents the current node number, and γ represents the number of node splits.
Step 4: Write the nodes and update the element connectivity to generate the CIE, which is a zerothickness six-node triangular prism.Notably, when updating the element connection, the node order of cohesive elements should follow the thickness direction as demonstrated in Figure 5a.
Figure 5b shows the cohesive elements for the biochar (CIE_BC) and the interface of biochar-ITZ (CIE_BCITZ), which are simulated using the COH3D6 element in ABAQUS.The material law of the cohesive element will be further explained in Section 2.4.2.

Mesoscale model generation
There are two different mesoscale models considered in the numerical simulations: one is the reference group that only contains fine aggregates, and the other is the BC-10 group with 10 vol% biochar replacement of fine aggregates.In the previous study, the aggregate distribution size is often characterized by the Fuller curve (Wriggers & Moftah, 2006) as shown below: where () is the cumulative percentage passing as a sieve with an aperture diameter ,  max is the maximum size of aggregates, and  was usually taken as 0.5.
The above method is commonly used for concrete mesoscale modeling, where the mortar is assumed to be a homogeneous material, and only coarse aggregates are considered (Thilakarathna et al., 2020;Wu et al., 2022;Xu et al., 2022).The aggregate volume is usually 35%, which is not large enough when only counting coarse aggregate (over 5 mm).However, when simulating the 50-mm cubic representative volume element (RVE) under mortar scale, the number of aggregates will be over 10,000 and make calculation impossible if considering the aggregate size below 2.56 mm.An effective solution to achieve consistency between simulated assumptions and experimental designs on a mortar scale is to employ a gap-graded aggre-gate design.Similar to the approach by H. Zhang et al. (2020) in their study on the tensile strength of mortar, they prepared a 10-mm cube and only considered aggregate particle sizes between 1 and 2 mm for both experiments and modeling.However, the 10-mm size is too small to adequately investigate the compression behavior (strainstress relationship) and cracking propagation using DIC.Therefore, the suggested sample size should be a minimum of 50-mm cube.To bridge the numerical simulations and experiments more efficiently, only an aggregate size of 2.36−4.75mm is considered both in tests and modeling based on volume fraction.
Various placing algorithms are adopted from the previous study, including the take-and-place method (Meng et al., 2020; Z. M. Wang et al., 1999), graphic-based method (Qiu & Dai, 2021;Zheng et al., 2021), and size scaling method (H.Zhang, et al., 2020;Quey et al., 2011).The takeand-placed method is adopted for randomly placing aggregate due to its high efficiency and adaptability, as it allows for the simultaneous generation and placement of realistic aggregates without relying on computed tomography.The flow chart for the method is shown in Figure 6.Initially, biochar particles were placed before sand aggregate due to their irregular geometric shapes and lower relative volume fraction.From an established biochar library (Figure 3), a particle is randomly selected and positioned within the domain with a random position (x, y, z coordinates).This step involves an iterative process of 600 cycles to ensure no intersections occur; if intersection criteria are still not met after these iterations, a new particle is chosen from the library.Once positioned, each biochar particle is rotated at a random angle.The process continues until the total volume of biochar particles meets the preset requirements.Subsequently, the placement of sand aggregates begins, following a similar procedure.It is important to note that the generation of sand aggregate and ITZ is implemented via the algorithm described in Section 2.1.1,while the biocharcement ITZ is created through the insertion of cohesive elements.
Once all particles are generated and placed in the domain, all aggregates and ITZ are converted as one entity region through a merging operation.Figure 7a shows the procedures of three-phase mesoscale modeling, with 30 vol% aggregates and its ITZ according to Section 2.1.1.To ensure that the different phases have common element nodes after meshing, two steps of the Boolean operation are executed: The first step involves combining the ITZ and aggregates to form the aggregate-ITZ entity, and the second step involves combining the matrix with aggregate-ITZ entity to form the three phases (aggregate, ITZ, and matrix) entity model.All three phases are meshed using solid elements, and the unstructured mesh element tetrahedral (C3D4) is generated for meshing due to the  Two additional phases should be considered in the BC-10 group, that is, the biochar phase and the ITZ near biochar particles, together with previous phases to form the five-phase mesoscale model.From the observation of SEM and backscattered electron images (Maljaee et al., 2021;Zhu et al., 2023), the cement matrix is assumed to fill the void of biochar particles in the five-phase model.The modeling process of the five-phase model is similar to the three-phase model, but the generation algorithm is more complex due to the incorporation of biochar particles as shown in Figure 7b.The biochar particles are randomly selected from the biochar library (as shown in Figure 3), then placed and rotated to varying degrees at random before the aggregates are positioned.The cohesive element is adopted to model the cracking behavior for biochar particles and their ITZ.All the components are meshed using tetrahedral elements, with the characterized element size   = 1 mm to balance the accuracy and computational cost.For cement matrix, sand aggregate, and sand ITZ, a 4-node linear tetrahedron (C3D4) element type is employed, utilizing free meshing techniques.For biochar and biochar-cement ITZ, a 6-node three-dimensional cohesive element (COH3D6) is chosen.The sensitive analysis can be found in Section 2.3.All code and datasets accompanying this manuscript are available on GitHub at https://github.com/LMduo.

Loading and boundary conditions
As shown in Figure 8a, two rigid plates are used for the uniaxial compression simulations.The specimen modeled contains 1242 aggregates and 145 biochar particles in total.Two reference points (RFs) are established and linked to the upper and bottom surface through kinematic coupling interaction to control the movement of the loading and supporting surface.The displacement control loading is applied on the upper RF to simulate the quasi-static situation.The supporting condition is governed by the RF on the lower surface, where the movement along all the directions is fixed.The mechanical behavior (cracking propagation and compressive strength) is significantly influenced by the friction coefficient between the specimen's surface and plate during the experiments (Bandeira et al., 2022), which is discussed in Appendix B, guiding this study to select a coefficient of 0.2.It should be noted that a longer prism would be more effective in reducing the frictional effect (Fascetti et al., 2022).However, using the prism would significantly increase the mesh size and computational time due to the large amount of sand aggerate in the mortal scale, which led us to choose a 50-mm cube for this study.To cover a larger range of strain-stress responses and compare the results with DIC, the maximum loading displacement applied is set as 0.25 mm (0.5% strain).The loading time is set to be 0.02 s, which is long enough to minimize the potential dynamic effect (see details in Section 2.3).

Sensitivity analysis
It is widely acknowledged that the numerical results are significantly influenced by the finite element mesh size.
A very small element size results in unnecessary computation expenditure, whereas a large mesh element size produces imprecise results.Moreover, the loading time is one critical parameter in the explicit quasi-static analysis, which addresses the convergence problem related to crack propagation.Theoretically, the loading time should be as long as possible to minimize the dynamic effect, but the computational cost drastically increases.To find a suitable mesh size and loading time that satisfy both accuracy and efficiency, the reference model with 30% volume fraction aggregate is meshed with four different sizes (0.8, 1, 1.2, and 1.4 mm) and loading under four different periods (0.05, 0.01, 0.02, and 0.03 s).In the current simulations, conducted on a system equipped with a CPU 13900K, 64GB of RAM, and an RTX 4090 GPU (NVIDIA), the computational time varied significantly with mesh size.Specifically, when the finite element mesh size was set to 1.4 mm, the simulation was completed in approximately 4 h, while a finer mesh size of 0.8 mm extended the computational time to about 15 h.The stress-strain curve for every loading time and finite element size under uniaxial compression is illustrated in Figure 8b.The stress is obtained from the total nodal reaction force of the RF divided by the surface area (50 × 50 mm 2 ).The engineering strain is calculated as the average nodal displacement of the RF divided by the length of the RVE.It can be found that the mesh size presents a negligible effect on the stress-strain curves within 1.2 mm.Substantial violation of the stress-strain responses only appears (Figure 8b) when the time step is 0.005 s.Hence, to balance the accuracy and computational cost, the mesh size and time step are set as 1 mm and 0.02 s, respectively.

Concrete damage plasticity law
The CDP model is widely used to simulate the damage for both cement matrix and ITZ (Huang et al., 2015), where the total strain rate .

𝜀 𝑝𝑙
(5) The compression and tension damage parameters used in this study are proposed by Jankowiak and Lodygowski (2005), where the stress-strain relationship is defined below: where   and   present compressive and tensile stress;   and   are the damage parameters for compression and tension;    ,    ,    , and    are the total strain and plastic strain for both compression and tension; and  0 is the initial elastic modulus for undamaged concrete.
The states of failure and damage are governed by the yield function  in Equation ( 7) by representing a surface in the effective stress space.Lubliner et al. (1989) first proposed the yield function in the CDP model, which was further modified by Lee and Fenves (1998): where The σmax represents the maximum principle effective stress, and this value of ITZ is 0.5-0.9 for mortar (Jin et al., 2017).σ (ε   ) and σ (ε   ) are the effective compressive and tension cohesion stresses, respectively.
The degree of damage is indicated through DAM-AGEC for compression, DAMAGET for tension, and SDEG presents the stiffness degradation, which are used to evaluate the crack propagation, where "0" and "1" stand for undamaged and completely damaged states, respectively.Although the CDP model cannot define the cracking depth accurately because its theory is based on continuum mechanics, it can effectively simulate and predict the crack propagation.In this study, the model assumes cracking failure when the SDEG value is over 0.9 (Naderi & Zhang, 2021).

Cohesive interface element law
The biochar particles and the ITZ near biochar are modeled using cohesive elements, which are widely used for material cracking (Naderi & Zhang, 2022).The bilinear traction-separation law is adopted in biochar modeling, which consists of three main definitions, namely, linear elastic (Equation 11), damage initiation (Equation 13), and damage evolution regions (Equation 16).The linear elastic behavior is represented as follows: where  and  are the nominal traction stress vector and separation vector, respectively, and  is the elastic matrix consisting of three components that should be defined as, namely, the normal stiffness   and two shear stiffness   and   .
Based on our previous study (Zhu et al., 2023), the elastic modulus of the ITZ surrounding the biochar is slightly higher than that of the biochar itself.To relate this to the data obtained from nanoindentation (see Section 3.2), the cohesive stiffness matrix can be approximated using Turon's approach below (Turon et al., 2007): where  is the coefficient affecting the compliance of the composite, which should be set to be large enough (Le et al., 2017;Neville, 1995); the   is the characteristic element size, and  is the passion ratio of the material, which is taken as 0.37 for biochar (Gibson & Ashby, 1997;Wimmer et al., 1997).
The damage initiation is defined based on the quadratic nominal stress criterion in Equation ( 13): where    denotes the peak stress in the normal direction;    and    are the peak stress in two shear directions.The symbol ⟨ ⟩ represents the Macaulay bracket: The damage evolution is defined based on the fracture energy (), which is calculated under the bilinear tractionseparation law, as calculated below: where   is the effective traction at damage initiation, and   is the effective displacement at failure.A damage evolution law defines the progressive stiffness degradation of the cohesive elements.In this study, a linear damage evolution is adopted with the effective displacement (  ) (Camanho et al., 2007) and scalar damage variable , governed by the equations below: where   is the normal relative displacement,   and   represent the two shear relative displacements; the  max  is the maximum effective displacement that can be obtained during the loading; the    and    represent the effective displacement corresponding to damage initiation and complete failure.The progressive stiffness degradation can be estimated from variable  from intact (0) to complete failure (1) in SDEG.

Input parameters
The input parameters of the model include geometric definition (Table 1) and material law (Table 2).Table 1 summarizes the input conditions for geometric control parameters, including the size of the RVE, aggregate particle size distribution, ITZ thickness, biochar particle size distribution, and the control parameters for porous biochar based on Voronoi tessellation.The sand aggregate used in the mix design is river sand, which is resistant to damage during compression and will not break in the current experimental conditions.In this model, the material law of sand is assumed to be elastic, with an elastic modulus of 70 GPa.(Naderi & Zhang, 2021).
For the matrix, a water-to-cement (w/c) ratio of 0.4 results in a 28-day compressive strength of approximately 42 MPa (Abd elaty, 2014).The descending branch of the stress-strain curve and the elastic modulus are defined according to GB50010.Typically, the maximum principal effective stress for ITZ ranges between 0.5 and 0.9 for mortar (Jin et al., 2017).Other parameters include the plastic potential eccentricity, ratio of initial equiaxial compressive yield stress to initial uniaxial compressive yield stress, ratio of the second stress invariant on the tensile meridian   , and the viscosity is taken as 35 • , 0.1, 1.16, 0.667, and 0.0005 for ITZ and mortar (Sümer & Aktaş, 2015).For the biochar particles, the cohesive stiffness is calculated through Equation ( 12), where the elastic modulus is obtained by nanoindentation (see Section 3.2), and the characteristic element size   is taken as 1 mm.The elastic modulus for the biochar-cement ITZ, obtained from our previous study (Zhu et al., 2023), indicates that the ITZ around biochar exhibits slightly

Material and sample preparation
The mortar mixture was prepared using CEM I Portland cement, manufactured by Hong Kong Green Island cement company (see Table 3 for its chemical composition), deionized water, and river sand with a particle size ranging from 2.36 to 4.75 mm.A total of 12 samples were cast, with three samples allocated to each group for compressive testing.The sample was prepared using a w/c ratio of 0.4, and extra water (50 wt% of biochar) was added to the mortar for pre-soaking of biochar.As only 30 vol% sand was used in the mortar, the sand/cement mass ratio was calculated to be 0.75 after multiple iterations (see Appendix A for estimation).Biochar, characterized by a specific gravity of 0.58 kg/m 3 , was introduced as a partial replacement for sand in quantities ranging from 10 to 30 vol%, as detailed in Table 4.The mixing procedures entailed several steps: Initially, dry cement was blended with sand for 30 s while simultaneously pre-absorbing biochar with water at 50 wt%.Subsequently, the biochar was mixed with sand and cement for another 30 s, followed by the addition of deionized water and 90 s of mixing (30 s at a low speed and 60 s at a high speed).An additional 90-s mixing step at a low speed was included, resulting in a total mixing time of approximately 3 min.After mixing, the fresh mixture was poured into three 50-mm cube molds and vibrated 20 times to minimize air entrapment.The cubes were then sealed and cured for 28 days under controlled laboratory conditions, maintaining a temperature of 20 ± 2 • C.

Nanoindentation test of biochar
Due to limited studies on the mechanical properties of biochar, the nanoindentation test (Hysitron TI Premier, Bruker) was applied to provide a robust and comprehensive baseline for numerical simulation.One piece of biochar sample (∼ 10 × 10 × 5 mm 3 ) was selected and impregnated with resin epoxy (Buehler EpoxiCure 2).
To remove the open voids and air bubbles, the vacuum chamber (Cast N'Vac, Buehler) was used for five cycles (10 min/cycle).Then, the sample was polished with SiC paper (#300, #600, and #1200 grit) to remove the hardened epoxy resin on the surfaces of the samples, followed by polishing with suspensions containing diamond (MicroPolish II Alumina, Buehler) with the particle size of 9, 3, and 1 µm in sequence.Each polishing step was carried out for 30 min, applying an 8.8 N force.Finally, the ultrasonic bath filled with alcohol was adopted to remove the residual alumina on the surfaces of the samples.
To investigate the mechanical properties of biochar under nanoindentation more accurately, different indentation loads, indentation rates, and holding times were designed (Figure 9) in this study.The indentation load varied from 0.2 to 10 mN; the indentation rate ranged from 10 to 300 µN•s; and the holding time was explored from 0 to 100 s to investigate their effects on the reduced elastic modulus (Er) of biochar.The pores in a single biochar sample can be clearly visualized from a Surface Topographical Image (shown in Figure C1).To prevent potential probe damage and erroneous results associated with a mesh grid approach, the point-by-point indentation approach was adopted in this study.For enhanced reliability of the results, over 30 different locations were selected for indentation, and more than 10 experiments were conducted for each loading scheme.This approach involved collecting and analyzing over 300 points across various loading schemes as detailed in Section 4.1.

Digital image correlation
For the uniaxial compression test, four cubic specimens (BC-10, BC-20, and BC-30) with dimensions of 50 × 50 × 50 mm 3 were fabricated for reference.The test was conducted using a 1000-kN universal hydraulic testing machine (HCT-A Series, WANCE).The compression test was executed in the "displacement control mode," maintaining a constant displacement rate of 0.2 mm/min to obtain stable and comprehensive stress-strain curves.Throughout the test, the 3D-DIC system (Figure 10) was used to monitor both the crack propagation process and the surface strain field.The speckles were applied using a hand painting in advance: first by evenly spraying a matte white paint onto one side of the sample, serving as the background; subsequently, matte black spray paint was applied from a 30-50 cm distance at a 45-degree angle relative to the sample, resulting in the uniform distribution of speckles across the sample surface.Two industrial cameras (Q-400 GenTL/GenICam) combined with a 25-mm fixed focal length lens were operated to capture speckle images at a rate of 5 Hz.Three LED lights were employed to provide ample illumination, and an 80-mm calibration target (SS-HP-AL) was used to calibrate two cameras before the test.The software ISTRA 4D provided by DantecTM was applied for correlating two cameras and analyzing the speckle images.

Micro-mechanical properties of biochar
One technical difficulty in accurately probing the mechanical properties is finding a suitable experimental approach because the bulk compression test does not work specifically on the nature of a porous material.The nanoindentation technique enables the tip to indent several micrometers and to record the load-displacement behavior (p-h) simultaneously based on the pre-defined loading stage (Pharr et al, 1992).From the p-h curves (see Figure 11a), the indentation modulus, M, and the hardness, H, can be obtained based on the Oliver-Pharr justification (Oliver & Pharr, 1992;Pharr et al., 1992).
where  is the unloading stiffness,  max is the maximum indentation load, and  is the projected area estimated from the maximum indentation depth (ℎ max ) (Oliver & Pharr, 1992).Indentation modulus (i.e., reduced elastic modulus), , depicts the combination of the elastic characteristic of the material and the indenter tip.For the homogeneous material, the elastic modulus, E, can be calculated through Equation ( 19) (Oliver & Pharr, 1992): where   and   are the elastic modulus and Poisson's ratio of the indenter tip, respectively, and ν is the Poisson's ratio of the material.
Biochar pyrolyzed from wood waste exhibited exceptional elasticity (Figure 11a) without any residual depth observed within an indentation load ranging from 0.2 to 10 mN.As illustrated in Figure 11c, the elastic modulus of biochar decreased gradually with the increase in indentation load, which can be proven in the relationship shown in Figure 11b.A possible explanation for this phenomenon is that the pores in the porous structure of biochar resulted in stress concentration around the pores.As the indentation depth increased, the stress concentration effect became more pronounced, resulting in a lower elastic modulus of the porous material (Gong et al., 2011;Kasyap & Senetakis, 2021).Moreover, the elastic modulus of biochar was independent of the loading rate but increased slightly due to insufficient holding time, which may be attributed to the viscoelastic behavior of biochar, requiring a sufficient holding time to ensure energy dissipation (S.Liang & Wei, 2020).A detailed geometric representation of the porous structure of biochar was established in computational simulations.In addition to recognizing the influence of porosity on the reduction of the elastic modulus, the precise modeling of the geometric structure in this study can adopt the elastic modulus values obtained under low indentation loads (specifically, within the range of 0.2-0.3mN).Under such minimal loads, the behavior of materials closely resembles a non-porous solid state.

Strain-stress curve and cracking detection
The complete uniaxial strain-stress curves of different specimens from both experiments (plain mortar, BC-10, BC-20, and BC-30, as solid curves) and numerical simulations (plain mortar and BC-10, as scattered plots) are presented in Figure 12a.The stress values were calculated by dividing the reaction forces at the RF by the initial area (50 × 50 mm 2 ).The strain values were obtained by dividing the applied displacement load by the size of the RVE (50 mm).A small discrepancy was observed for the REF group at the end of the curve, where the simulated values slightly exceeded the experimental data, which might stem from the CDP model limitation in accurately replicating the nonlinear and stochastic process of complete material failure (Lee & Fenves, 1998;Lubliner et al., 1989;Naderi & Zhang, 2021).In the ductility analysis, the computation of strain ductility typically employs a strain value at 85% of the compressive strength (Y.Zhang et al., 2019).This corresponded to strain values of 0.0021 and 0.0026 for the REF group and the BC-10 group, respectively.Within this specific range, the simulation results could exhibit a high degree of concordance with the experimental data, both in terms of curve patterns and stress magnitudes.To gain deeper insights into the influence of biochar as a partial substitute for fine aggregates on the material performance, the mechanical responses were examined across four distinct stages: elastic stage, crack initiation stage, crack development stage, and failure stage.Figure 12b showcases the corresponding crack propagation at each stage, obtained through DIC and ISTRA4D software.Certain regions of the crack propagation exhibited blank areas due to the occurrence of speckle fracture, rendering the correlation process infeasible.During the elastic stage (A), the specimens demonstrated negligible damage, and the stress exhibited the most rapid growth, characterized by a slope equivalent to Young's modulus.Subsequently, in the crack initiation stage (B), stress growth gradually decelerated, accompanied by the emergence of yield stress and the initiation of small-scale crack propagation as visualized in Figure 12b-ii.As the materials progressed to the crack development stage (C), the ultimate compression strength was observed, coinciding with a pronounced increase in the number of cracks and a rapid decline in stress.Finally, at the failure stage (D), the stress reduction became relatively stable, while the number of cracks escalated swiftly until complete material failure as depicted in Figure 12b-iv.For the BC-20 and BC-30 groups, the crack patterns are almost symmetric, and the final failure mode exhibits a pyramid shape.However, the symmetry is less apparent in the REF group and BC-10 group, where some cracks appear parallel to the loading directions.This could be attributed to the limitation of DIC in capturing the full failures and the difficulties in controlling boundary restrictions (friction coefficient) during the experiments.
The energy-dissipation coefficient, that is Equation (20), and strain ductility coefficient, that is, Equation ( 21), are introduced in this study to assess the ductility and energy dissipation characteristics of biochar-cement composites (Z.Chen et al., 2014): where   is the energy-dissipation coefficient,   denotes the area enclosed by the stress-strain curve at a peak stress of 0.85, and   denotes the area enclosed by the rectangle formed by the strain value corresponding to 0.85 of peak stress at the descending branch times the peak stress.
The strain ductility coefficient is obtained by the following equations (Liu et al., 2020;Y. Zhang et al., 2019) where β is the strain ductility coefficient,  0.85 and  0.5 represent the strain corresponding to 85% and 50% of compressive strength after failure, respectively, and   represents the strain at the yield point.Table 5 lists the strain ductility and energy dissipation coefficients of the biochar-cement composites at different biochar replacement ratios.With the progressive substitution of fine aggregates by biochar particles, there was a reduction in the compressive strength and elastic modulus of the materials.However, this decrease was accompanied by an enhancement in both strain ductility and energy dissipation coefficient of the composites, particularly evident in the BC-30 group.The strain ductility coefficient at 0.85 and 0.5 of the peak stress in the BC-30 group exhibited an increase of 12% and 42.5%, respectively, compared to the reference group.The energy dissipation coefficient at BC-30 and BC-20 groups displayed a 16.3% increase, compared with the reference group.

Crack propagation
The cracking process in biochar-cement composites can be analyzed by studying the distribution of damage fields.In this study, failure of the mortar is assumed when its SDEG value exceeds 0.9.BC-10 was chosen for simulations to balance computational efficiency to observe the biochar cracking mechanisms, as the differences between 10, 20, and 30 were insignificant for this purpose.shaped damage forms.Notably, boundary restraints were idealized and controlled with a 0.2 friction coefficient in simulations.However, in the experiments (Figure 12b), boundary restraints are challenging to control and can significantly impact crack distribution, potentially explaining the discrepancies between simulation and experimental results.When the external strain load reached 0.0036, the REF model showed a higher matrix crack density (14.1%), compared to the BC-10 model (11.2%).
The sand ITZ and biochar-cement ITZ exhibit different modes of cracking.In comparison to the cement matrix, the sand aggregate ITZ undergoes initial damage, resulting in the formation of minor cracks around the sand and the extension toward the cement matrix (Figure 13).These minor cracks subsequently intersect and merge, leading to the formation of the primary cracks.While the biochar-cement ITZ experiences simultaneous damage with the biochar matrix, the severity and location of damage are closely related to its orientation (Figure 14).The discrepancy in the crack density of the two groups can be attributed to several factors.First, replacing 10% sand aggregates with biochar led to a reduction in the volume of the CDP element, which is calculated together with the matrix volume.As a result, the overall damaged element number could be reduced.However, this effect is minimal due to the low ITZ thickness, being only 0.5% of the total meshed elements.A more significant factor was the reduced stiffness and strength of biochar, compared to the aggregates and matrix materials.Moreover, the inherent elasticity and porous structure of biochar, as visually demonstrated in Figure 11, caused biochar to undergo failure at an earlier stage under external loads.This phenomenon can be further investigated in Figure 14a, where certain biochar particles displayed damage at a strain of 0.01 (Stage A), although the composite remained within the elastic stage.This failure pattern, where the weaker region fails before the others, resembles the concept of "strong-column weak-beam" in structural design (Mahini & Ronagh, 2010).This involves sacrificing compressive strength to mitigate crack formation and enhance ductility, which is proved by the strain-stress curve (Figure 12a).Furthermore, this behavior became more pronounced when there was a higher proportion of biochar replacement.As the ductility results are shown in Table 5, the peak strain is delayed, and the strain ductility is significantly increased as more biochar particles are added to the mix design.These experimental and simulation results confirmed that replacing aggregates with porous, low-stiffness materials (e.g., biochar) can increase the damping coefficient of cement-based materials (Chi et al, 2019;Kharitonov et al., 2015).Therefore, replacing aggregates with an appropriate volume of biochar can enhance the ductility of cementitious materials while maintaining the compression strength and slightly inhibiting the potential crack growth.

Anisotropic damage patterns of biochar particles
For an in-depth observation of the biochar failure process, the damage characteristics of the BC-10 group were specifically analyzed.The progression of biochar failure is illustrated by the fracture process of biochar as the strain load increases at an interval of 0.0005 (Figure 14a-i).The simulated results revealed the failure of biochar in a regional pattern (Figure 14a-ii).The upper portion of biochar experienced failure first, and afterward, the damage gradually propagated from the top and bottom toward the center.The major fracture of biochar happened when the strain load reached 0.003, while a small portion of biochar particles exhibited negligible damage throughout the entire loading process (Figure 14a-i).At the strain load of 0.005, most of the undamaged biochar particles were aligned parallel to the loading direction, whereas the biochar particles oriented perpendicular to the loading direction experienced some degree of damage.This behavior was closely related to the anisotropy caused by the porous structure of biochar particles.
Figure 14b illustrates the failure process of biochar particles with specific pore structures under strain load at various orientations.The failure mode of biochar was closely related to its position and angle within the mortar.When the biochar particles were horizontally placed (i.e., perpendicular to the direction of force, see Figure 14c), the cracks initiated inside the pores and rapidly propagated to both sides, ultimately resulting in extensive damage to the central region of the pore structure.When the biochar particles were vertically embedded in the mortar (i.e., parallel to the direction of force, see Figure 14c), the extent of damage was minimal.Cracks originated from the sidewalls of the pores and gradually propagated inward.Most cracks concentrated on the sidewalls, while the central pore structure was preserved.Biochar particles exhibited a mixed failure pattern when they were inclined at the 45 • angle to the direction of force (see Figure 14c).This configuration demonstrated the most severe level of damage, as cracks simultaneously propagated from the interior and exterior surfaces.Most biochar particles were completely destroyed when the strain load reached 0.003.The key reason underlying this phenomenon probably stems from the preservation of the inherent anisotropy of wood waste (Song et al., 2018) within the biochar after pyrolysis.As a result, substantial discrepancies appeared in the mechanical responses and failure characteristics of the biochar particles across different orientations.This phenomenon was consistent with the simulation results, as hardly any inclined biochar particles could be preserved when the strain load reached 0.005 (Figure 14a-i).
To conduct a more comprehensive investigation into the failure process of biochar particles as fine aggregates, the damage in different orientations of various pore structures was simulated (Figure 14c), and the extent of damage was quantified (Figure 15a), where the damage percentage was calculated by dividing the number of failed elements by the total number of elements.Notably, the damage progression exhibited similar trends across different biochar pore structures, with crack propagation showing no significant variations based on the differences in pore structures.
However, the orientation of biochar placement was confirmed as a critical factor influencing the extent of damage.Vertical orientation resulted in the least damage regardless of the biochar structure, as all three pore structures exhibited damage levels below 30% at the maximum loading.Conversely, when biochar was inclined at a 45-degree angle to the applied force, the damage exceeded or approached 75% for all three pore structures (Figure 15a).Furthermore, the occurrence of cracks varied depending on the angle at which biochar was subjected to the external load.Some pore structures only experienced damage when the strain load reached approximately 0.0018 under vertical loading, while the cracks appeared in the Degree45 group at a low strain load of 0.0006.Nevertheless, the biochar particles experienced earlier damage, compared to the cementitious matrix, as evident by the fact that most biochar particles inclined at a 45-degree angle reached 24% damage at a strain load of 0.0012, while the cementitious matrix remains largely undamaged (Figure 15b).By the time the strain load reached 0.016, varying degrees of damage were observed in all three orientations of biochar particles (38% horizontal), whereas the damage to the cementitious matrix only constituted approximately 3% of the total volume.

CONCLUSIONS
This study introduces a novel five-phase mesoscale modeling methodology, presenting a pioneering paradigm in the field of numerical simulation for biochar-cement composites with accurate geometric information of porous biochar in joint with homogenized cement matrix, irregular shape of sand aggregates, and their ITZ.The nanoindentation and FE-SEM were performed to characterize the relevant statistical measures of the pore structure and mechanical properties of biochar particles.Then, a uniaxial compression test with DIC was conducted to rigorously validate the numerical model.This research substantially advances the understanding of the mechanical and failure behavior of biochar-cement composites from both theoretical and experimental perspectives, offering valuable guidelines for the design and development of more resilient and durable construction materials.The key findings can be summarized as follows: 1.A novel approach is proposed for precise geometric modeling of porous aggregates, accompanied by a dedicated geometric library.By adjusting the parameters related to the pore size and shape, this library is suitable for various materials, including biochar, coral aggregate, and lightweight aggregate.Meanwhile, for the first time, the numerical model of biochar-cement compression simulations under varying fraction coefficients.Damage (cracks) exhibited significant variation across different fraction coefficients, with a notable concentration at the corners under frictionless conditions, gradually evolving into an "X-shaped" pattern as the coefficient increased (Figure B1a).The ultimate compression strength at maximum coefficient was approximately 30% higher than frictionless conditions, which was consistent with Mattews's experiments (Bandeira et al., 2022;Figure B1b).Experimental verification demonstrated that under various compression conditions, the compression strength of the reference group (rigid steel plate) surpassed the frictionless group by an increment of up to 20%.This corresponded to the simulated friction coefficients ranging from 0.15 to 0.25.Therefore, a friction coefficient of 0.2 was chosen for this simulation.

APPENDIX C SURFACE TOPOGRAPHY OF BIOCHAR
F I G U R E C 1 Surface topography contour plot of biochar.

F
Meshed sand aggregate at different random fluctuation values (with 1-mm mesh element size), varying according to (a) angle fluctuation and (b) radius fluctuation.

F
Generation process of the porous geometry using Voronoi cells: (a) transferring the basic cube geometry to Voronoi cells based on roundness and pores number control, (b) introducing radius and thickness factor, and (c) inputting the porous geometry to finite element model.F I G U R E 4 Numerical models with different biochar particles in a volume unit of 2500 mm 3 (1%, 3%, 5%, and 10% volume noted by P bc ).

F
Cohesive element insertion and node numbering for biochar particles and biochar interfacial transition zone (ITZ) region: (a) 2D scenario and 3D scenario and (b) cohesive element display for single biochar particles.CIE_BC, cohesive elements for the biochar; CIE_BCITZ, interface of biochar-ITZ.

F
I G U R E 6 Flow chart for generating biochar and aggregate particles in the model.F I G U R E 7 Procedure of establishing: (a) three-phase mesoscale finite element model (30% volume fraction of aggregate) and (b) five-phase mesoscale finite element model (30% volume fraction of biochar replacement of sand).
14678667, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mice.13307,Wiley Online Library on [25/08/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License complex and irregular shape of aggregates.The matrix and ITZ region are modeled using the CDP law.

F
I G U R E 8 Prerequisites of the simulation: (a) boundary conditions of the computational uniaxial compression test of 10 vol% biochar-cement composites model and (b) sensitive analysis of simulation parameters, including element number and time step on the stress-strain response using the reference model.

14678667, 0 ,
Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mice.13307,Wiley Online Library on [25/08/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License F I G U R E 9 Loading schemes of nanoindentation tests: (a) varied indentation load, (b) varied loading rate, and (c) varied holding time.

14678667, 0 ,
Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mice.13307,Wiley Online Library on [25/08/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License F I G U R E 1 0 Uniaxial compression test setup with DIC: (a) DIC setup and (b) specimen testing.F I G U R E 1 1 Nanoindentation results of biochar particles: (a) displacement-loading curves, (b) elastic modulus changing with indentation depth, and (c) comparison of the measured elastic modulus amongst different loading stages.

14678667, 0 ,
Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mice.13307,Wiley Online Library on [25/08/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License F I G U R E 1 2 (a) Typical compressive stress-strain curves for experiments and numerical simulations among different biochar replacement percentages and (b) DIC results for crack detection at different stages with varying biochar particle substitution percentages.
Figure 13a,b illustrates the propagation of exterior and interior cracks within 3D mesoscale models under uniaxial compression, respectively, with a biochar volume replacement of 0% (REF) and 10% (BC-10).The applied load was directed in the z-axis, and crack propagation was distinctly observed during the fracture stages defined by stress-strain (stage C & D: ε = 0.0036 & 0.005).Both models exhibited similar crack initiation patterns (Stage B), where cracks started from the sides and gradually extended toward the center as shown in Figure 13a-iii,b-iii, with the expected pyramid-F I G U R E 1 3 Crack propagation of 3D mesoscale model under uniaxial compression in Z direction among four stages (ε∈(0.002,0.005)) in (a) the reference group and (b) BC-10 group; (i) & (ii) 3D view and (iii) cross-section view (elements with SDEG values over 0.9 are deleted).

F
Crack propagation in biochar particles in (a) BC-10 group under compression in Z-direction; (b) specific porous structures with various orientations (i.e., horizontal, vertical, and 45-degree inclination); and (c) various porous structures.

14678667, 0 ,
Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mice.13307,Wiley Online Library on [25/08/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License F I G U R E 1 5 Damage quantification: (a) cumulative damage percentage of biochar particles with various pores structure at different placing directions and (b) cumulative damage percentage of the matrix.

14678667, 0 ,
Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mice.13307,Wiley Online Library on [25/08/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License F I G U R E B 1 (a) Failure modes of concrete under uniaxial compression at various friction coefficients: (i) full view and (ii) cross-sectional view, and (b) strain-stress curve of concrete under uniaxial compression at various friction coefficients.
TA B L E 1
TA B L E 4
TA B L E 5