Coupled DEM/CDF analysis of impact of free water on the static and dynamic response of concrete in tension regime

In this paper, the quasi-static and dynamic behavior of partially fluid-saturated concrete under two-dimensional (2D) uniaxial tension at the mesoscale is investigated numerically. It was calculated what effect the free pore fluid content (gas and water) had on the fracture process and strength of concrete in a tension regime. An improved pore-scale hydro-mechanical model based on DEM/CFD was used to simulate the behavior of fully and partially fluid-saturated concrete in quasi-static and dynamic conditions. The foundation of the fluid flow idea was a network of channels in a continuous area between discrete elements. In very porous, partially saturated concrete, a two-phase laminar fluid flow was proposed. To track the liquid/gas content, the position and volumes of the pores and cracks were considered. Numerical simulations of bonded granular specimens of a simplified spherical mesostructure that resembled concrete were carried out under dry and wet circumstances for two different strain rates. Investigations were conducted on the effects of fluid pore pressures, fluid saturations, and fluid viscosities on the strength and fracture process of concrete. The quasi-static tensile strength dropped nonlinearly with increasing fluid saturation and fluid viscosity during fluid migration through pores and cracks which contributed to a promoted fracture process. However, during rapid dynamic tensile deformation, a fracture process attenuated due to a fluid migration constraint from insufficient time to exit the pores. It caused the dynamic tensile strength to rise nonlinearly as fluid saturation and fluid viscosity increased.


Introduction
Concrete constructions are often in contact with water.Because concrete is porous and has many macropores, capillary pores, fissures, and other openings, free water may seep through it.Concrete's static and dynamic mechanical characteristics may be significantly impacted by this (Rossi and Boulay, 1990;Rossi, 1991;Li, 2004;Wittmann et al., 2007).Therefore, it is essential to research how free water affects concrete's mechanical response while building concrete structures like hydraulics, marine engineering, bridges, and tunnels.However, no thorough explanation consensus has yet been established regarding the physical mechanism of strength changes under varying water saturation.Concrete contains chemically and physically linked water in addition to free water.The effect of the moisture variation in concrete during loading is usually disregarded in practice.
It is commonly recognized that under quasi-static stress conditions in laboratory testing, wet concrete materials behave much differently than under dynamic loading situations (Reinhardt et al., 1990;Ross et al., 1996;Yan et al., 2005;Zhang et al., 2010;Deng et al., 2014).The impact of water saturation on the compressive strength of concrete increases with initial porosity (Shen and Xu, 2019;Wang et al., 2022).In comparison to concrete with a low water-cement ratio, the mechanical characteristics of concrete with a high water-cement ratio are also more susceptible to water saturation.Generally speaking, wet concrete responds to loads similarly to wet mortars and rocks.
The physical mechanism leading to the decreased strength of wet concrete in quasi-static tests can be interpreted in a variety of ways, such as the loosening of a molecular system in ITZs (Zhang et al., 2019) and the growth of an internal humidity gradient (Xu et al., 2017;Wang et al., 2022).According to some scholars, the primary causes are water seeping into fractures (Wang et al., 2022) and the increase in pore pressure brought on by pore closure (Oshita and Tanabe, 2000).It was also discovered that the weaker van der Waals forces and cohesiveness between microscopic particles (Zhang et al., 2019a;Wang et al., 2008) were the causes of the strength decline.Ultimately, the primary causes of the strength decrease were identified as capillary suction (Rossi, 1991;Chen et al., 2012), viscosity (Rossi, 1991;Fu et al., 2021), viscous friction (Fu et al., 2021), and temperature-induced expansion (Weiss et al., 2004;Fu et al., 2021).
One of the most crucial elements influencing the impact of material strain rate is water content (Cadoni et al., 2001;Wang and Li, 2007;Wang et al., 2009;Zhou et al., 2011;Wang et al., 2016;Fu et al., 2021;Chiu et al., 2021;Wang et al., 2022;Wang et al., 2023;Zhang et al., 2021).The dynamic strength of cement-based materials is enhanced by free water when the strain rate reaches a particular threshold.As the loading rate rises, so does this increase.The confinement effect of pore water to the cement matrix due to water migration constraint, the meniscus cracking effect, and the Stefan effect (viscosity effect) − induced delay of crack initiation and propagation are some of the existing explanations for how moisture affects porous materials in a dynamic range, leading to increased strength (Wang et al., 2022;Fu et al., 2021;Chiu et al., 2021).Because wet concrete is prepared differently, it should be noted that there can be some discrepancies in the data that are currently available (Sun et al., 2020).
Our current article aims to model wet concrete behavior under 2D quasi-static and dynamic conditions in uniaxial tension.In quasi-static and dynamic tensile mesoscale circumstances, the paper seeks to: a) explain the physical mechanism for strength dependence of fluid (water and gas) saturation in concrete, and b) quantify the effect of free fluid migration on the fracture process and concrete strength.The simulations were carried out in quasi-static (the strain rate was 0.01 1/s) and dynamic conditions (the strain rate was 1.0 1/s).The solid-fluid interaction was studied using an improved fully coupled hydro-mechanical approach based on 2D DEM/CFD (Krzaczek et al., 2020(Krzaczek et al., , 2021)).In the approach, two domains coexisted in one physical system: a discrete solid domain and a continuous fluid domain.A direct numerical simulation (DNS) technique was used to solve fluid flow equations in a fluid domain between discrete elements.Since the 3D DEM/CFD model is still undergoing testing, a 2D model was utilized instead.The suggested model considers immiscible and compressible fluids flowing in a two-phase laminar flow.The fluid can be moved only by changes in pore volumes, which can cause changes in local pressure without the requirement for external pressure.This allows the suggested model to simulate the way that fluids and solids interact during compression/tension tests.
The proposed DEM/CFD model employs a direct numerical simulation method.DEM was utilized to simulate the mechanical behavior of concrete specimens with discrete spherical elements interacting locally through normal and tangential contacts that caused cracks.The literature has produced several discrete models for concrete mechanical behavior, including discrete element approaches, rigid-body-spring models, lattice discrete particle models, and classical lattice models (Bolander et al., 2021).Since DEM realistically depicts a fracture process in concrete (Skarżyński et al., 2015;Nitka andTejchman, 2015, 2018;Suchorzewski et al., 2018aSuchorzewski et al., , 2018b)), we used it.Furthermore, it depicts mesostructure with the use of simple dynamic equations to represent contact forces.
The laminar viscous two-phase fluid flow (water and gas) through pores and cracks in a continuous domain between the discrete elements was represented by CFD using a flow network composed of specified channels.In this research stage, calculations were carried out under isothermal conditions on 2D small-size bonded granular specimens of a simple mesostructure that resembled concrete under uniaxial tension.Without making a distinction between aggregate, mortar, ITZs, and macropores (Nitka and Tejchman, 2015;Suchorzewski, 2018aSuchorzewski, , 2018b)), concrete was seen as a single-phase material consisting of bonded spheres with varying sizes.In addition, the range of particle diameter sizes was narrow and the number of spheres was quite low.As a result, the mechanical behavior of concrete (stress-strain curve and fracture process) in laboratory experiments was roughly mimicked.However, those limitations did not qualitatively affect the nature of a coupled hydro-mechanical problem investigated (Krzaczek et al., 2023a).The paper primarily examined the impact of pore fluid migration and pressure on the tensile strength of the specimen.The various fluid saturation and fluid viscosity were the factors investigated in a series of coupled DEM/CFD simulations.It should be noted that 2D modeling is an intermediate step to gain some understanding of the physical model and basic phenomena but it does not provide quantitatively accurate modeling of real 3D behavior.
For the sake of simplicity, the capillary pressure, which is the difference between the partial pressures of the wetting phase (water) and the non-wetting phase (water vapor or moist air), was disregarded in the current DEM analysis.This definition states no capillary movement when the concrete is completely saturated.Because of the viscous fluid flow, capillary flow can happen in partially saturated concrete even in the absence of an external fluid source as during uniaxial tension.Certain capillaries may be filled with water by the liquid phase flow to the point where the snap-off mechanism (Singh et al., 2017;Shams et al., 2021) is triggered (the only capillary flow mechanism that could occur inside the specimen).This mechanism instantly fills capillaries with water.This mechanism can solely occur locally under particular conditions.The proposed DEM-CFD model (Krzaczek et al., 2023a) included capillary pressure and capillary flow mechanisms.
The literature has proposed several coupled hydro-mechanical DEM-CFD models (Cundall, 2000;Hazzard et al., 2002;Al-Busaidi et al., 2005;Lathama et al., 2008;Shimizu et al., 2011;Chareyre et al., 2012;Yoon et al., 2014;Catalano et al., 2014;Xiao-Dong et al., 2015;Zeng et al., 2016;Ma et al., 2017;Zhang et al., 2017;Papachristos et al., 2017;Liu et al., 2019;Zhang et al., 2019b;Caulk et al., 2020).A review of those coupled models based on DEMs was provided in (Krzaczek et al., 2020).The DEM/CFD-based mesoscopic approach for fluid flow in partially fluid-saturated porous cohesive-frictional materials with very low porosity, which is reported in this paper, offers various advantages over the existing coupled DEM-CFD models.These advantages include.1.The direct numerical simulation approach is used to solve fluid flow equations in a two-dimensional continuous fluid domain between discrete elements.2. The fluid and discrete domains of solids reside in one physical system.Both domains are discretized into one triangular mesh.3. Due to the triangular mesh in the fluid domain, the variable geometry, size, position, and volume of pores/cracks are considered to correctly track the liquid/gas content.Hence, the effect of material skeleton deformation on the fluid pressure distribution can be studied.4.An effective method is developed to automatically mesh and re-mesh these domains to account for changes in the geometry and topology of particle and fluid domains.5. Coarse meshes of solid and liquid domains are utilized to build a virtual fluid flow network (VPN).6.The two-phase immiscible fluid contains both a liquid (water) and a gas (water vapor or moist air).
In summary, the fully coupled DEM-CFD model used in the paper is significantly more accurate than all existing ones regarding fluid pressures and fluid and solid velocities.
The existing DEM-CFD models using co-simulation techniques, based on unresolved and resolved methods (Kuruneru et al., 2019;Zhou and Zhao, 2022) have several shortcomings such as the empirical drag force, the dependence on the relative velocity and volume fraction, the condition on the number of particles compared to the cell size, the disregard for the flow behavior surrounding the particle, the simple particle shape, etc.Thus, those models are not able to realistically capture fluid flow between discrete elements, in particular when the initial porosity of Fig. 1.Two domains coexisting in one physical system: a) co-existing domains before projection and discretization and b) solid and fluid domains after discrete elements projection and discretization (fluid domain is in red and solid domain is in grey) (Krzaczek et al., 2020(Krzaczek et al., , 2021)).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)M. Krzaczek et al.

Downloaded from mostwiedzy.pl specimens is very low.
There also exist hydro-mechanical models, based on lattice and lattice discrete particle models (Bolander and Berton, 2004;Grassl and Bolander, 2016;Luković et al., 2016;Athanasiadis et al., 2018;Eliáš and Cusatis, 2022;Mašek et al., 2023).The models are simplified since the pore fluid pressure does not depend on volume changes which are ignored.The fluid flow is solely caused by external fluid pressure or relative humidity difference.However, during loading, the fluid flow is

Downloaded from mostwiedzy.pl
mainly driven by changes in pore pressures, which are affected by changes in pore volumes.For the problem of the effect of free water in concrete, there is also a DEM solution in the compression situation (Benniou et al., 2021); however, it does not account for the true interaction between solid (DEM) and fluid mechanics (CFD).Moreover, several numerical coupled hydromechanical continuum models were also used to study the effect of free water content on the global mechanical behavior of concrete under different stress levels in static (Forquin et al., 2015;Bian et al., 2018) and dynamic loading scenarios (Zhao and Wen, 2018;Huang et al., 2020).In those simulations, a fracture process was ignored (Forquin et al., 2015;Bian et al., 2018) or simplified empirical equations were taken into account (Zhao and Wen, 2018;Huang et al., 2020).
The authors conducted a numerical analysis of hydraulic fracturing in rocks (Krzaczek et al., 2020(Krzaczek et al., , 2021) ) and hydraulic/capillary flow in unsaturated mortars and concrete (Krzaczek et al., 2023a) using the improved fully coupled DEM-CFD model.The coupled DEM/CFD technique was confirmed with detailed 3D CFD simulations employing the Reynolds-averaged Navier-Stokes equations in the continuous domain between particles (Abdi et al., 2022(Abdi et al., , 2023)).It was shown that in the simplified fluid flow model, the turbulent kinetic energy and turbulent dissipation energy are negligible and can be ignored.The authors' coupled DEM/CFD model now includes heat transmission (Krzaczek et al., 2023b;Krzaczek and Tejchman, 2023).
The current paper is organized as follows.Following the introduction in Section 1, a mathematical model of the coupled hydro-mechanical approach based on DEM/CFD is presented in Section 2. The data used for DEM-CFD simulations is described in Section 3. Section 4 compiles the findings of the pure DEM simulation.The behavior of fully and partially saturated concrete under uniaxial tension is shown in Section 5 through numerical static and dynamic simulation findings using DEM-CFD.Some closing thoughts are included in Section 6.The authors implemented the improved fully coupled DEM-CFD technique for partially saturated bonded granular materials in the open-source YADE software application (Kozicki and Donze, 2008;Šmilauer et al., 2021).

Two-dimensional DEM/CFD-based model
The improved coupled DEM-CFD model was explained in detail by Krzaczek et al., 2020Krzaczek et al., , 2021. .Only the most important features of the model are included in Section 2 for clarity.

DEM for cohesive-frictional materials
The discrete element open-source software YADE's 3D explicit solver is used to do DEM simulations (Kozicki and Donze, 2008;Šmilauer et al., 2021).Particles in a DEM interact with one another during translational and rotational motions through the application of Newton's second law of motion and an explicit time-stepping technique (Cundall and Strack1979).The model predicts a cohesive bond at the grain contact and brittle failure below the critical normal tensile force.Shear cohesion failure under typical compression results in sliding, which is governed by the Coulomb friction law.Appendix A contains the fundamental DEM controlling equations (Eqs.A1-A7) (Skarżyński et al., 2015;Nitka andTejchman, 2015, 2018;Suchorzewski et al., 2018aSuchorzewski et al., , 2018b)).In simulations, non-viscous damping is selected (Cundall and Hart, 1992) (Eq.A7) to accelerate convergence.The possibility of particle overlap in DEM allows for the possibility of achieving an arbitrary micro-porosity.For DEM simulations, the ensuing material constants are necessary: E c , υ c , μ c , C, and T (Appendix A).The parameters R s , ρ (mass density), and α d (damping factor) are also required.It is standard practice to set the damping factor at α d = 0.08 (Nitka and Tejchman, 2015).To accurately replicate the distribution of shear and tensile cracks, the ratio between the uniaxial compressive and tensile strength, and the failure mode of specimens (brittle or quasi-brittle), the particle contact ratio C/T must be carefully considered (Tomporowski et al., 2023).
To find the material constants, a series of DEM simulations are usually performed and the outcomes are compared with experimental data from basic tests including uniaxial compression, triaxial compression, and simple shear.Material softening is not taken into account by the DEM model.The authors have effectively utilized the model for simulating engineering particle materials, namely granular materials (Widulinski et al., 2011;Kozicki et al., 2013Kozicki et al., , 2014;;Kozicki and Tejchman, 2018), concretes (Skar żyński et al., 2015, Nitka and Tejchman, 2015, 2018, 2020;Suchorzewski et al., 2018aSuchorzewski et al., , 2018bSuchorzewski et al., , 2019) ) and rocks (Krzaczek et al., 2020(Krzaczek et al., , 2021;;Tomporowski et al., 2023).The current DEM models did not account for grain damage.While their time will be significantly increased, this process may be easily included in DEM simulations (by employing collections of particles to represent the grains allowing for intra-granular fracturing).One of DEM's drawbacks is how long industrial-scale problems take to calculate.

Fluid flow model
A fundamental idea for streamlining the fluid flow model was put out by Krzaczek et al., 2020Krzaczek et al., , 2021. .In the authors' developed 2D DEM-CFD model: a discrete 3D domain (solid) and a continuous 2D domain (fluid) coexisted in one physical system (Fig. 1).This is not the case with the basic concept of conventional fluid flow networks.Since using a 3D DEM solver was necessary for the YADE software, it was assumed that the discrete domain was initially 3D.The solid domain was composed of a single layer of discrete three-dimensional objects (spheres), whereas the initial fluid domain was two-dimensional.The centers of gravity of discrete elements were found on a plane (2D surface).To create circles, discretized 3D elements were first projected onto a plane.These circles, along with voids, were then discretized into a 2D triangular mesh (Fig. 1).
In both domains, the meshing/remeshing processes automatically produced triangular meshes (Krzaczek et al., 2020).If one or more discrete elements were moved beyond the specified maximum distance, the operation was initiated.Initially, a collection of points in two dimensions was produced.The points were generated along discrete elements and contact edges and in the voids between discrete elements.Next, the first mesh of convex objects was created using Delaunay triangulation methods.Afterward, the mesh was refined using the Alpha Shapes theory to replicate more realistically the geometry (e.g.concave voids).Although coarse, the produced mesh was of acceptable quality.The identification process, which classified each mesh element into a solid or fluid domain was next completed.In the last stage, under the assumption that the mass was a topological invariant, the remeshing process transferred the computation results from the previous mesh to the current one.The remeshing procedure increased the calculation time by factor 3.
After projection and discretization, there were two 2D domains left: the fluid domain and the solid domain.Note that the movement of the spheres was fixed in the direction of OZ, that was, perpendicular to the plane.As a result, the mechanical problem was effectively 2D, much like the fluid flow model, even when using YADE's 3D DEM solver.All definitions relevant to the geometry of the two domains (fluid and solid) were handled as two-dimensional in the CFD model, which had a unit size in its third dimension.
The fluid flow network in the 2D fluid domain is formed by the channels connecting the centroids of the triangles.A re-meshing technique discretizes the overlapping circles, fixes the contact lines, and removes the overlapping areas to provide a more accurate distribution of pressures, fluid-phase fractions, and densities (Krzaczek et al., 2020(Krzaczek et al., , 2021)).Based on the pressure and shear stress for a specimen thickness equal to the mean particle diameter, the contact forces are computed.Appendix B contains the fundamental equations for the fluid flow model (Eqs.B1-B28).

Input data for 2D DEM/CFD simulations
DEM and CFD are usually calibrated independently using mechanical (DEM) and permeability and sorptivity testing (CFD) (Krzaczek et al., 2020(Krzaczek et al., , 2021(Krzaczek et al., , 2023a)).Due to 2D simplified coupled simulations, the calibration procedure was also suggested for the 2D conditions.For the numerical calculations in the first calculation stage, a simple 2D onephase model for concrete made out of spheres with varying diameters was employed to roughly depict concrete's performance in uniaxial tension as compared to laboratory tests.The diameters of spheres d were arbitrarily set to range from 1.0 mm to 1.8 mm, with a mean diameter of d 50 = 1.4 mm as during uniaxial compression (Krzaczek et al., 2024a(Krzaczek et al., , 2024b)).For uniaxial tension tests with dry and wet specimens (containing about 1600 spheres), a quadratic specimen of 50 × 50 mm 2 was assumed with two opposite initial weak elements at the mid-height to induce a tensile crack (see Figs. 3 and 4 in Section 4).The spheres had no cohesion in weak elements.The specimen thickness was always equal to the mean grain diameter.Horizontal boundaries at the top and bottom were frictionless and smooth.Vertical boundaries were free to move.
One grain at the bottom midpoint was fixed.There were no pre-existing macro-pores or cracks in the specimen.A continuous rise of the top specimen boundary caused the tensile deformation in the specimen.The strain rate was chosen as ∊ = 0.01 1/s in static and ∊ = 1.0 1/s in dynamic simulations.The initial micro-porosity was the same as during uniaxial compression p = 5 % (being equivalent to typical concrete (Suchorzewski et al., 2018a;Nitka and Tejchman, 2020)).The predicted tensile strengths, the corresponding vertical normal strains, the elastic moduli, and the crack patterns were in agreement with the experimental ones for usual concrete (Tejchman and Bobinski, 2013).
For permeability 2D simulations with pure CFD, a quadratic nondeformable specimen of 10 × 10 mm 2 was chosen (Krzaczek et al., 2023).The macroscopic permeability coefficient κ of the concrete specimen, calculated with Darcy's law, was 4.0 × 10 − 16 m 2 for the Fig. 9. DEM-CFD simulation results for one-phase concrete specimen during uniaxial tension with different initial fluid volume fraction α q : relationship between vertical normal stress σ y and vertical normal strain ε y : a) pure DEM, b) DEM-CFD with α q = 0.3, c) DEM-CFD with α q = 0.5, d) DEM-CFD with α q = 0.7 and e) DEM-CFD with α q = 1.0 (A) quasi-static analyses and B) dynamic analyses).
The basic material constants assumed for concrete specimens in the 2D coupled DEM-CFD calculations are given in Table 1.The adaptive time step in DEM and CFD was always defined (Krzaczek et al., 2020).The maximum time step was limited to 1⋅10 − 8 s.The computation time of one simulation was about 5 h on a computer with two Intel Xeon Platinum processors 8280 (2.70 GHz).The computational cost of the simulation was relatively high because the existing DEM-CFD model was parallelized on threads only but not in a distributed mode (on cluster computer nodes).

Pure DEM simulation results
For quasi-static and dynamic mechanical simulations Section 3, a concrete specimen of 50 by 50 mm 2 was assumed, which is comparable to the specimen used in the uniaxial compression experiments that were discussed by Krzaczek et al., 2024aKrzaczek et al., , 2024b.The numerical results for concrete, which DEM defines as a one-phase material consisting of spheres with sizes ranging from 1.0 to 1.8 mm (without differentiating between aggregate, mortar, ITZs, and macropores (Skarżyński et al., 2015;Nitka andTejchman, 2015, 2018)) are shown in uniaxial tension in Figs.2-8.The stress-strain graphs in both static (strain rate ∊ = 0.01 1/s) and dynamic ( ∊ = 1.0 1/s) conditions are shown in Fig. 2. The evolution of the fracture patterns and broken particle contacts for the different vertical normal strains are displayed in Figs.3-6.The yellow lines in Figs. 5 and 6 connect the particle centers and the red lines correspond to the broken contacts.Fig. 7 shows the evolution of the relative broken contact number.The distribution of the tensile and compressive contact forces at the maximum vertical normal stress is displayed in Fig. 8.The contacts can be broken either by tensile or shear forces.During uniaxial tension, tensile contact damage took place only.
The behavior of the static and dynamic specimens was comparable.The predicted tensile strengths were approximately f t = 2.70 MPa and the elastic moduli were E = 16 GPa for a vertical normal strain of ε y = 0.018-0.022% (Fig. 2).The range of elastic region was ε y = 0.014-0.015%.The stress-strain curve was only roughly reproduced, especially in a quasi-static regime, as a result of various simplifications assumed in the calculations to speed up the simulations (e.g.one-phase material composed of spheres, narrow particle diameter size range, a small number of spheres, and 2D conditions).The computed specimen had a too-linear, zigzag pre-peak response under quasi-static conditions, and an excessively brittle post-peak response (Fig. 2).The pre-peak dynamic response of the specimen was also excessively torn (Fig. 2).In simulations of uniaxial compression, splitting tension, and bending, the concrete behavior (stiffness, strength, brittleness, and fracture) can be realistically simulated without those three initial simplifications (Skarżyński et al., 2015;Nitka andTejchman, 2015, 2018).Crack formation was the reason behind jumps on the stress-strain curves in Fig. 2.
The specimen failed due to a single macrocrack of a sinusoidal shape between two opposing material flaws (Figs.3-6).In quasi-static and dynamic situations, the specimen experienced the onset of micro-cracks far in advance of the peak stress, at approximately ε y = 0.013 % (Fig. 7).They began to change at ε y = 0.015 % and continued to do so until the damage rate increased to the point of failure.Regardless of the strain rate, at the maximum vertical stress (ε y ≅ 0.02 %), about 0.6 % of contacts were damaged.The external vertical load was transferred by a network of vertical normal contact forces (Fig. 8).Although some compressive normal contact forces were also present beyond the macrocrack, tensile normal contact forces predominated (Fig. 8b).Specimen failure was caused by damage of vertical tensile contact forces (Fig. 8b).The empty spaces seen in Fig. 8 are due to the scale assumed.

DEM/CFD simulation results
For hydro-mechanical calculations in tension, the same 50 × 50 mm 2 concrete-like specimen (Fig. 3a and 4a) was used.The specimen's initial fluid pressure (water and gas) was P o = 0.1 MPa.Fig. 9 illustrates the effect of different initial fluid volume fractions α q (α q = 0.3, 0.5, 0.7, and 1.0) on the stress-strain curve during uniaxial tension.Fig. 10 illustrates the relationship between the tensile strength f t and α q .
Regardless of the parameter α q , the numerical curves σ y = f(ε y ) Fig. 10.DEM-CFD simulation results for one-phase concrete specimen during uniaxial tension with different initial fluid volume fraction α q : relationship between tensile strength f t and α q (a) quasi-static tests and b) dynamic tests).exhibited a comparable shape for quasi-static and dynamic situations (Fig. 9).The quasi-static tensile strength f t and the corresponding vertical normal strain ε y decreased as the initial fluid volume fraction α q rose (Fig. 10).This result agrees with laboratory tests where the splitting tensile/tensile strength continuously decreased with the growth in water saturation (Zhao et al., 2023;Chen et al., 2012).Because the specimens' original porosity was low (5 %), there were little variations in their tensile strength.
As shown in Fig. 10, the expected dynamic strength f t and the associated vertical normal strain ε y increased steadily.In the laboratory studies, when the water saturation rose, the dynamic tensile strength increased either non-linearly (Cadoni et al., 2001) or linearly (Wang et al., 2023).Because of the lack of cracks and low specimen porosity (Sun et al., 2020), low pore fluid pressures (see Fig. 16) maintained the calculated initial material stiffness up to ε y = 0.015.This numerical result is consistent with the compressive experiment by Zhang et al., 2019a but it differs from several lab studies in which the elastic modulus increased (Zhu et al., 2012;Wang et al., 2022;Zhao et al., 2023) or decreased (Shoukry et al., 2009;Boxu et al., 2018) with a greater saturation level.Compared to a dry specimen (α q = 0), a fully saturated specimen (α q = 1.0) had a 20 % lower vertical normal strain and an approximately 8 % lower tensile strength for the quasi-static deformation (Fig. 9A).In dynamic simulations, a fully saturated specimen (α q = 1.0) had 5 % higher vertical normal strain and around 5 % higher tensile strength (Fig. 9B).The dynamic brittleness of a partially saturated specimen was marginally higher.The variations in f t with α q were marginally less than in uniaxial compression test outcomes (Krzaczek et al., 2024a(Krzaczek et al., , 2024b)).

Mechanical results for wet specimens
Similar to the pure DEM results (Figs. 3-7), the DEM-CFD findings for the fluid-saturated specimen (α q = 1.0) are shown in  show how the fracture patterns and broken contacts developed.The damaged contact number changed with specimen deformation, as seen in Fig. 15.
The failure of the wet specimen was similar to that of the dry specimen (Figs.3-6).On the other hand, compared to the dry specimens, the macro-crack in the wet specimens (Figs. 13 and 14) was flatter and more irregular (had more branches).The wet specimen experienced a faster breakage of sphere contacts during quasi-static simulations (ε y = 0.014 %) compared to the dry specimen (ε y = 0.015 %) (Fig. 15A).The number of broken contacts in the dynamic simulations was similar to ε y = 0.018 %; however, in the dry specimen, the broken contacts occurred more quickly in the range of ε y = 0.018-0.021% (Fig. 15B).

Fluid flow results for fully saturated specimen
Fig. 16a shows the evolutions of the average pore fluid pressures during uniaxial tension, while Fig. 16b shows the evolutions of the average pore fluid pressures during dynamic simulations.In Fig. 17a (quasi-static simulations) and Fig. 17b (dynamic simulations), the average pore fluid velocities are evolving.The distribution of pore fluid velocities in the specimen in quasi-static (Fig. 18a) and dynamic (Fig. 18b) analyses at the peak stress is depicted in Fig. 18.
During quasi-static deformation, the average fluid pore pressure rose to roughly 0.104 MPa (Fig. 16a).It exceeded the initial value of 0.10 MPa.The dynamic average fluid pore pressure likewise increased up to ε y = 0.012 % (long before ε y peak = 0.02 %) before declining (Fig. 16b).Its minimum value was 0.098 MPa, and its maximum value was 0.103 MPa, which was comparable to quasi-static analyses.For ε y > 0.012 %, the pore fluid pressure was less than the initial one of 0.1 MPa.As a result, the specimen experienced suction.
Other DEM-CFD results using a partially saturated specimen for the two different initial fluid volume fractions α q (α q = 0.3 and 0.7) (not shown here) showed a similar material behavior.For α q = 0.7 and α q = 0.3, the maximum average fluid pressures were around 0.1002 MPa and 0.10003 MPa, respectively.The quasi-static average fluid pressure (0.104 MPa) was higher than the dynamic one (0.101 MPa) at the peak stress.
The dynamic fluid migration was more limited, as shown by the maximum quasi-static average fluid pore velocity (0.0018 m/s) being much higher than the dynamic one (0.0003 m/s) (Fig. 17).In quasistatic and dynamic simulations, the maximum pore fluid velocity at the peak stress was determined to be 0.16 m/s and 0.003 m/s, respectively (Fig. 18).Fig. 13.Coupled DEM-CFD simulation results for one-phase wet (α q = 1.0) concrete specimen during quasi-static uniaxial tension: distribution of broken contacts in non-deformed specimen for different vertical normal strain ε y (a) ε y = 0.013 %, b) ε y = 0.015 %, c) ε y = 0.0175 % (at stress peak stress), d) ε y = 0.019 % (after peak stress)) (yellow lines connect particle centers, red lines denote broken contacts between spheres due to tension).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 18 shows the distribution of pore fluid velocity vectors in the specimen for the peak stress during uniaxial tension for both quasi-static and dynamic analysis.For the peak stress in the quasi-static analysis, the distribution of the pore fluid velocity vectors showed large differences.The fluid velocity distribution (as well as the maximum velocity) in a concrete specimen might vary significantly between time steps.Sudden changes in fluid velocity were caused by rapid and random changes in fluid volume in pores and cracks.The fluid continued to flow through the entire specimen volume but its velocities at the specimen bottom were smaller compared to the maximum fluid velocity at this time step and hence invisible in Fig. 18a.

Effect of viscosity
The influence of viscosity μ q on the stress-strain curve during tension was examined in the case of totally fluid-saturated concrete (α q = 1).The specimen's tensile strength increased by 4 % in the quasi-static simulation and decreased by 1 % in the dynamic simulations with a factor 2 decrease in viscosity.The effect of viscosity resembled this during compression (Krzaczek et al., 2024a(Krzaczek et al., , 2024b)).
In summary, the physical mechanism for strength dependence of fluid saturation was affected by changes in the fluid pore pressure during loading which variously impacted a cracking process.If loading was quasi-static (slow), the fluid pore pressure sped up cracking and contributed to the material strength decrease.Flowing out of pores the fluid promoted higher contact forces and faster material degradation.If loading was dynamic (fast), the fluid pore pressure retarded cracking.The fluid had no time to freely flow out of pores and cracks, causing a delay in its migration.Thus, the fluid confinement appeared in pores and cracks causing the material strength to increase.DEM is more suitable than continuum models for describing the physical mechanism for strength dependence since it takes the development of cracks into account from the beginning of deformation through strongly diverse local strength criteria (Chodkowski et al., 2022).

Summary and conclusions
This work uses an improved fully coupled DEM/CFD-based hydromechanical technique to predict the behavior of partially fluidsaturation concrete in uniaxial tension under simplified 2D isothermal circumstances.The method discretized the geometry of pores and cracks while accounting for the two-phase laminar fluid flow in the network of designated channels.The numerical model was validated based on 2D uniaxial tension and permeability test results.The technique proved to

Downloaded from mostwiedzy.pl
be effective in assessing the tensile strength, fracture process, and pore fluid pressures of wet concrete specimens.Significant research was done on the effects of fluid saturation and viscosity on pore fluid pressures, concrete strength, and the fracture process.In terms of strength change with fluid saturation and viscosity, the numerical results demonstrated qualitative agreement with laboratory studies reported in the literature.The following findings are drawn from the 2D simulations for concrete specimens with an initial porosity of 5 %: -Compared to dry specimens, the quasi-static/dynamic tensile strength of concrete for totally fluid-saturated specimens was approximately 8 %/5 % lower/greater (Fig. 10).Under quasi-static or quasi-dynamic situations, it decreased/rose with higher fluid saturation (Fig. 10).The decrease exhibited nonlinearity as the fluid saturation increased (Fig. 10).When the fluid viscosity was two times lower, the tensile strength of concrete for fully fluid-saturated specimens increased/decreased by 4 %/1 % in the quasi-static/ dynamic regime.This was qualitatively similar to the response under uniaxial compression.-Changes in tensile strength were caused solely by deformations of the material skeleton forced by an external load in the absence of external pressure sources in the fluid area.Those deformations led to changes in the volume of pores and cracks along with the fluid contained in them.Changes in fluid volume caused a change in fluid pressure.That contributed to fluid flow.The fluid had more time to reach a higher velocity in the static tensile test and force the fluid pressure to change over a larger area.For this reason, the impact of fluid flow in the pores and cracks affected the tensile strength differently in static and dynamic tests.-The tensile strength was reduced by the quasi-static deformation, which promoted a fracture process during fluid migration through pores and cracks.The loss in tensile strength was raised with increasing pore fluid pressure.-Due to rapid dynamic tensile deformation, the fluid migration constraint resulted from insufficient time for the fluid to exit the pores.The raised dynamic tensile strength was facilitated by an attenuated fracture process due to the fluid confinement in pores.

Appendix A
Figure A1 illustrates the mechanical response of DEM (Kozicki and Donze, 2008;Šmilauer et al., 2021).The DEM equations are listed below: where

Appendix B
The fluid flow network in the 2D fluid domain is formed by the channels connecting the centroids of the triangles.A remeshing procedure discretizes the overlapping circles, establishes the contact lines, and removes the overlapping areas to provide a more accurate distribution of pressures, fluid-phase fractions, and densities.The mass change in triangular cells is correlated with the density change in a fluid phase, which results in pressure changes.Triangles defy the equation for the conservation of momentum as a result, but their mass is conserved across the whole volume.This process is performed by employing a Virtual Pore Network (VPN), which is an explicit formulation for each VP in the fluid flow network.The gas and liquid may initially be found in the material and any pre-existing discontinuities.VPN introduces two different channel types (Krzaczek et al., 2020): (A) artificial 'S2S' channels are those that connect discrete concrete elements that are in contact with one another and (B) actual 'T2T' channels are those that link grid triangles in pores that are in contact with one another along a common edge.
The distance between the gravity centers of neighboring grid triangles is assumed to be the channel length.Solving the continuity and momentum equations for the laminar flow of an incompressible fluid yields the mass flow rate in channels.Only the fluid flow via the network's face (edge) is used to estimate the mass flow rate of the fluid.
A single gas phase flow with a gas phase fraction (a), a single liquid phase flow with a liquid phase fraction (b), and a two-phase flow (liquid and gas) (c) are the three flow regimes found in the channels (Krzaczek et al., 2020).The mass flow rate in flow regimes a) and b) (single phase fluid flow) is calculated using the Poiseuille equations (Reynolds 1883, Batchelor, 2000).The two-phase fluid flow caused by a pressure gradient in neighboring VPs behaves similarly to a two-phase fluid flow of two immiscible and incompressible fluids in channels (Fig. B1).
The gas phase fraction and the liquid phase fraction are denoted by the symbols α p and α q , respectively, and they sum equals 1.These are the mesh's single cell's local parameters.The criteria for setting the initial conditions can be the same for every grid cell.The liquid-gas interface is parallel to the channel plates and is calculated as the average of the liquid phase fraction in adjacent cells.However, the fraction of liquid and gas phases changes over time, as well as the position of the phase interface in the channel.Gravity's effects are disregarded.While the volumetric flow rates of the fluid phases are unknown, the interface position is connected to fractions of the fluid phases in adjacent VPs.Equations of continuity and momentum describe the flow in each phase.(1) mass flow rate for each phase of fluid flowing through the cell faces (in channels surrounding VP) can be calculated by solving momentum and continuity equations, (2) phase fractions and their densities in VPs can be calculated using equations of state and continuity, (3) pressure in VPs can be calculated using the equation of state, and (4) material properties (such as liquid and gas densities in each cell grid) can be updated.
A modified empirical method created by Hökmark et al. (2011) is used to calculate the hydraulic aperture h of the artificial channels 'S2S': ) where h inf − the hydraulic aperture for the infinite normal stress, h 0 − the hydraulic aperture for the zero normal stress, σ n − the effective normal stress at the particle contact and β − the aperture coefficient.The geometry of the nearby triangles has a direct bearing on the hydraulic aperture of the actual channels 'T2T' (Fig. B2) as h = γecos(90 where e is the edge length between two adjacent triangles, ω denotes the angle between the edge with the length e and the center line of the channel that connects two adjacent triangles, and γ is the reduction factor established to maintain the maximum Reynolds number Re along the main flow route at a value below the critical value for laminar flow (Re = 2100).The parameter γ was determined in parametric tests at fluid pressure up to 140 MPa.Following Barmak et al. (2016), the continuity and momentum equations are developed that are rendered dimensionless where u j = ( u j , v j ) and p j are the velocity and pressure of the fluid phase j, ρ j and μ j are the corresponding density and dynamic viscosity.The Reynolds number is Re p = ρ p u i h p /μ p and the density and viscosity ratios are r = ρ q /ρ p and m = μ q /μ p .The lower and upper phases in the dimensionless formulation, are located in the regions − n d ≤ y ≤ 0 and 0 ≤ y ≤ 1, where n d = h q /h p .The velocities satisfy the no-slip boundary condition at the channel walls.The continuity of velocities and tangential stresses is required by the boundary conditions at the interface (Barmak et al., 2016) u q (y = 0) = u p (y = 0) (B6) and The mass flow rates M q,x and M p,x for both fluid phases are derived by solving Eqs.B3 and B4 with the boundary conditions (Eqs.B5-B7), as well as the shear stress τ f0 at the channel surfaces.
Unlike the model of fluid flow in the channels, the fluid in VPs (triangular cells) is assumed to be compressible.The mass conservation equation for the liquid phase in discretized form is where V n+1 i and V n i are the volume of VP i at a time increment n + 1 and n (the third dimension is the unit dimension), respectively, f is the face (edge) number, U n f denotes the volume flux through the face [m 3 /s], based on the average velocity in the channel, α n q,f is the face value of the liquid phase volume fraction [-], t is the time step [s], n denotes the time increment and i is the VP number [-].The same equations are defined for the gas phase.
For both fluid phases in VPs, the Peng-Robinson equation of state is used to describe fluid behavior above the critical point where P is the pressure [Pa], R denotes the gas constant (R = 8.314 J/(mol K)), V q/p is the molar volume of liquid (q) and gas (p) fraction [m 3 /kmol] and T denotes the temperature [K].The parameters in Eq.B10 are: ) n q/p = 0.37464 + 1.54226ω q/p − 0.26992ω 2 q/p , (B12) a q/p,0 = a c,q/p β(T), (B13) a c,q/p = 0.457247R 2 T 2 q/p,c P q/p,c , (B14) b q/p = 0.07780RT q/p,c P q/p,c , (B15) ) where T q/p,c is the critical temperature of phase [K], P q/p,c denotes the critical pressure of phase [Pa], ω q/p is the acentric factor of phase [ − ], and T r denotes the reduced temperature T Tc .When c 1 = c 2 = c 3 = 0, the original model is obtained.The extra factors help connect vapor pressure data from highly polar liquids like water and methanol.Equations B11-B16 provide a good fit for the vapor pressure, however predicting molar volumes for liquid phase can be very inaccurate.The forecast of saturated liquid molar quantities might deviate by l0-40 % (Mathias et al., 1989).Peneloux et al. (1982) proposed an effective correction term where s is the small molar volume correction term that is component dependent; V q is the molar volume predicted by Eq.B9 and V corr q refers to the corrected molar volume.The value of s is negative for higher molecular weight non-polar and essentially for all polar substances.The molar volume correction term is considered to be 0.0 m 3 /kmol and − 0.0034 m 3 /kmol for the gas phase and liquid phase (water), respectively.The Peng-Robinson equation of state has the advantage of being able to describe the behavior of supercritical fluids at extremely high fluid pressures and temperatures.For each phase, the mass conservation equation is used.The mass transfer between phases and the grid velocity is ignored when there is no internal mass source.The discretized form of the mass conservation equation for the liquid phase is with where f is the face (edge) number, U n f denotes the volume flux through the face [m 3 /s], based on the average velocity in the channel, α n q,f is the face value of the fluid phase volume fraction [-], t is the time step [s], n denotes the time increment and i is the VP number [-].The explicit formulation is used instead of an iterative solution of the transport equation during each time step since the volume fraction at the current time step is directly computed from known quantities at the previous time step.Similarly, the mass conservation equation for a gas phase is introduced.The product ρ q U n f α n q,f in Eq.B18 is the mass flow rate M q,f of the liquid phase flowing through the face f (edge of a triangle) of VP i .The density of the liquid phase can be calculated by solving the mass conservation equation for both phases The density of the gas phase can also be computed in the same way.It should be noted that the molar volume V (q/p) is related to the gas density and to the liquid density where w p and w q are molar weights of gas and liquid phases, respectively.Due to the fact that the fluid phases share the same pressure the fluid phase fractions are computed.Inserting Eq.B21 for the gas phase and Eq.B22 for the liquid phase into Eq.B23, a polynomial equation is obtained with respect to the liquid fraction α n+1 q,i .The gas-phase fraction is computed as α n+1 p,i = 1 − α n+1 q,i .Equation B10 is used to calculate the new pressure P n+1 i in VP i .The explicit formulation is utilized instead of an iterative solution of the transport equation during each time step.The edge of the channel experiences shear stress due to the passage of a viscous fluid.In the case of immovable parallel plates with no-slip boundary conditions (zero velocity), the fluid's shear stress profile is triangular.The shear stress at the liquid-gas interface is computed as The forces F → P,j and F → S,j acting on spheres are derived from the fluid pressures in VPs and 'S2S' channels.The fluid pressure forces acting on the sphere are computed, for simplicity, by multiplying the fluid-solid contact area b the pressure in the cell ( F → P,j ) or channel ( F → S,j ).Instead of calculating the contact area as a portion of the sphere's surface, it uses a section of the cylinder's surface with a height equal to the sphere's diameter.In doing so, the computing procedure is made simpler and the forces are only somewhat overestimated.
where n → − the unit vector normal to the discretized sphere's edge, P i − the pressure in VP, i − the VP number, j − the sphere number, and A k − the contact area between the fluid in VP i and sphere.
A k = 2r j e k , (B27) with r j − the sphere radius and e k − the sphere edge length.Ultimately, the forces acting on spheres are derived from the shear stresses as where I → − the unit vector parallel to the channel wall and oriented in the fluid flow direction, τ f0,i − the shear stress in the channel, i − the channel number, j − the sphere number and A k − the contact area between the channel and sphere, and L k − the channel length.

Fig. 2 .
Fig. 2. Pure DEM simulation results for one-phase dry concrete specimen (α q = 0) during uniaxial tension: relationship between vertical normal stress σ y and vertical normal strain ε y (a) quasi-static deformation and b) dynamic deformation).

Fig. 3 .
Fig. 3. Pure DEM simulation results for one-phase dry concrete specimen (α q = 0) during quasi-static uniaxial tension: evolution of fracture pattern (a) ε y = 0 %, b) ε y = 0.018 %, c) at peak stress (ε y = 0.022 %), and d) after peak stress (ε y = 0.023 %)) (weak particles with 0-cohesion are in pink, blue particles have constant vertical velocity and green particles are blocked in vertical direction, displacements are magnified by factor 100). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .Fig. 6 .
Fig. 4. Pure DEM simulation results for one-phase dry concrete specimen (α q = 0) during dynamic uniaxial tension: evolution of fracture pattern (a) ε y = 0 %, b) ε y = 0.016 %, c) at peak stress (ε y = 0.018 %), and d) after peak stress (ε y = 0.023 %)) (weak particles with 0-cohesion are in pink, blue particles have constant vertical velocity and green particles are blocked in vertical direction, displacements are magnified by factor 100). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Pure DEM simulation results for one-phase dry concrete specimen (α q = 0) during uniaxial tension: evolution of relative number N of broken contacts in [%] (a) static simulations and b) dynamic simulations).

Fig. 8 .
Fig. 8. Pure DEM simulation results for one-phase concrete specimen (α q = 0) during uniaxial tension: distribution of contact forces at stress peak: A) static simulation, B) dynamic simulation (a) compressive forces and b) tensile forces, force thickness denotes force magnitude).

Fig. 11 .
Fig. 11.Coupled DEM-CFD simulation results for one-phase wet (α q = 1.0) concrete specimen during quasi-static uniaxial tension: evolution of fracture pattern (a) ε y = 0 %, b) ε y = 0.015 %, c) at peak stress (ε y = 0.0175 %), and d) after peak stress (ε y = 0.019 %)) (weak particles with 0-cohesion are in pink, blue particles have constant vertical velocity and green particles are blocked in vertical direction, displacements are magnified by factor 100). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 .
Fig. 12. Coupled DEM-CFD simulation results for one-phase wet (α q = 1.0) concrete specimen during dynamic uniaxial tension: evolution of fracture pattern (a) ε y = 0 %, b) ε y = 0.016 %, c) at peak stress (ε y = 0.020 %), and d) after peak stress (ε y = 0.022 %)) (weak particles with 0-cohesion are in pink, blue particles have constant vertical velocity and green particles are blocked in vertical direction, displacements are magnified by factor 100). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 14 .
Fig. 14.Coupled DEM-CFD simulation results for one-phase wet (α q = 1.0) concrete specimen during dynamic uniaxial tension: distribution of broken contacts in non-deformed specimen for different vertical normal strain ε y (a) ε y = 0.014 %, b) ε y = 0.016 %, c) ε y = 0.020 % (at stress peak stress), d) ε y = 0.022 % (after peak stress)) (yellow lines connect particle centers, red lines denote broken contacts between spheres due to tension).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 15 .
Fig. 15.Coupled DEM-CFD simulation results for one-phase concrete specimen during uniaxial tension: evolution of relative number N of broken contacts in [%] (a) α q = 0.0 and b) α q = 1.0,A) static simulations and B) dynamic simulations).

Fig. 18 .
Fig. 18.Coupled DEM-CFD results for fully fluid-saturated specimen (α q = 1.0): distribution of pore fluid velocity vectors in specimen for peak stress during uniaxial tension (a) quasi-static analysis and b) dynamic analysis).
F → n − the normal contact force, U − the overlap between discrete elements, N → − the unit normal vector at the contact point, F → s − the tangential contact force, F → s,prev − the tangential contact force in the previous iteration, X → s − the relative tangential displacement increment, K n − the normal contact stiffness, K s − the tangential contact stiffness, E c − the elastic modulus of the particle contact, v c − the Poisson's ratio of particle contact, R s − the particle radius, R A and R B contacting particle radii, μ c − the Coulomb inter-particle friction angle, F s max − the critical cohesive contact force, F n min − the minimum tensile force, C − the cohesion at the contact (maximum shear stress at zero pressure), and T − the tensile strength of the contact, F → k damp − the dampened contact force, F → k and v → k p − k th − the components of the residual force and translational particle velocity v p and α d − the positive damping coefficient smaller than 1 (sgn(•) that returns the sign of the k th component of velocity).

Table 1
Basic material constants assumed for concrete, fluid and gas in DEM/CFD calculations.