Numerical Simulations of Reinforced Concrete Behaviors Using Smoothed Particle Hydrodynamics
2024 Volume 5 Issue 2 Pages 40-44
Details
2024 Volume 5 Issue 2 Pages 40-44
In this study, the smoothed particle hydrodynamics (SPH) code for solid mechanics is extended to simulate the three-dimensional (3-D) behaviors of reinforced concrete (RC) structures. The reinforcing bar is discretized and modeled as a 3-D truss element. To save the total degrees of freedom in the system domain, an implicit representation of the reinforcement is introduced, where the displacements in the reinforcement are interpolated using the displacements of adjacent SPH particles representing the matrix material. An isotropic damage model is introduced to reproduce the concrete damage behaviors. Thereafter, these numerical methods of the SPH, reinforcement model, and the damage model are integrated and validated through simulations of tension stiffening behaviors of a RC member. It is demonstrated that the proposed approach effectively predicts the longitudinal force and strain curve, as well as the multiple cracking behaviors of the RC specimen.
The particle-based numerical approaches are attractive for engineering simulations involving the large deformations or failures. The system domain can be discretized according to the particle arrangements and each representative material domain is assigned for each particle domain. Generally, the particle-based numerical approaches (e.g., Smoothed Particle Hydrodynamics1),2),3), Discrete Element Method4),5), and Material Point Method6),7)) have the Lagrangian nature such that the severe deformations or fragmentations in solid domain are accurately captured with associated particle movements without requiring additional numerical treatments, such as the elemental technologies in nonlinear Finite Element Method8). Since the SPH is originally developed for the simulations of astrophysics1), it has been extensively applied for the various fluid2),3) and structural9),10) simulations. In particular, there are various coupling approaches of SPH with the FEM for the simulations of reinforced concrete (RC) structures11),12).
In this study, the previously developed in-house SPH code3),13) for the open-channel flows and fluid- structure interaction (FSI) is extended to accommodate the simulations of structural deformations and damage behaviors of composite structures (e.g., RC, fiber-reinforced composites). The elastic solid behaviors are simulated based on the algorithmic SPH formulations10), and concrete-like quasi-brittle materials are modeled according to the isotropic damage model and related tensile strain- softening law14),15). For reinforcement models, the coupling approach of SPH and beam element was developed and validated12), where the discretized momentum equations are derived from the total energy concept including the coupling energy in the overlapping zone of SPH particles and beam elements. Using the moving least square (MLS) interpolations, it was shown that the nodal displacements are accurately calculated12). In addition, the semi-discrete reinforcement model was developed and integrated with the Voronoi-cell lattice models16),17) to alleviate the numerical burden without increasing the number of computing nodes within the system domain.
In the present study, these numerical concepts are reflected for the reinforcement model. Therefore, without accounting the additional nodal degrees of freedom for the reinforcing element, the background adjacent SPH nodes are utilized to define the relationship between the displacements and forces in the reinforcement unit. In the following sections, the proposed combination of the SPH method, damage model, and reinforcement model is validated for simulating the RC behaviors. As a benchmark example, the direct tensile test of RC member is selected for numerical simulations. It is demonstrated that the proposed approach can accurately reproduce the tension stiffening and cracking behaviors of RC specimen.
Meanwhile, one of potential applications of particle-based approaches like SPH is evident in the development of physics-informed neural network (PINN) models18),19). The PINN model adopts both the data-based and physics-based training procedures at each discretized particle location within the system domain. The SPH method, a strong-form based approach, solves the momentum equations using the SPH formulations at each particle location. Therefore, the SPH methodologies to solve the partial differential equations and related numerical data can be primarily utilized for each of the physics-based and data-based learning procedures. In the near future, the developed SPH program code will be used for real-time prediction of the stress or strain fields in structures using thte PINN and related health monitoring methods.
(1) Total Lagrangian SPH formulation
Fig.1 shows the planar view of particle distribution and particle interaction between the particle i and surrounding particles j within the radius of kernel support domain kh, where each paricle interaction is defined according to the kernel function Wij as follows2),
where αd is dimensional parameter, set as 1/πh3 in 3-D space, and q(=|χi-χj|/h) is the normalized distance between the particle coordinates χi and χj with the smoothing length h.
The total Lagrangian formulations are introduced to describe the geometric nonlinearity in the structural deformations. The current position vector x is defined as the sum of the referencial (or original) position vector X and displacement vector υ. Then, the relationship between χ and X is described by the deformation gradient F as8),
where ι is the identity matrix. Considering the elastic material behaviors, the following stress-strain relation is defined as,
where λ and μ are the Lamé constants, E (=(FTF-ι)/2) is Green-Lagrange strain matrix, S is the second Piola-Kirchhoff stress matrix, J is the sum of the diagonal elements of E. To account for the material damage, the isotropic damage parameter d is introduced, which changes from the intact (d = 0) to the fully damaged states (d = 1) with respect to the historical maximum strain κ as14),
where E is the elastic modulus, ε0 is the initial strain limit, Gf is the fracture energy, and lp is the characteristic length. Then, using Eqs. (3) and (4), the stress state can be updated as,
In this study, the algorithmic formulations of SPH based on the total Lagrangian framework10) are utilized, where the momentum equations are discretized as follows,
where ρ and m are the particle density and mass, respectively, ν are the particle velocities, P is the the first Piola-Kirchhoff stress matrix, b is the body force vector, and Pν is the artificially viscous stress matrix,
where the coefficient β is set to 0.04, and cs is the sound speed. For the time marching analyses, the explicit time integration scheme is considered, where a Verlet-like time integration is consistently utilized as in the previous validations of SPH code3),13).
(2) Reinforcement model
In previous studies16), 17), as an implicit representation of the reinforcement, the modeling method, socalled semi-discrete approach, was utilized to model the fiber-reinforced cementitious composites or RC strategy is introduced, and the nodal sites of the reinforcing elements are not included in the total degrees of freedom in the system domain. Instead, the nodal displacements of each reinforcing element are intropolated from the surrounding particle displacements. Using the Shepard filter2), the nodal displacements υi can be computed as,
where V̅ is the particle volumn. The reinforcing bar is modeled and discretized as the truss element in the finite element method20). Considering the cross-sectional area As, unit length ls, and elastic modulus Es of the reinforcing bar, the stiffness coefficient κs(=Es As/ls) is calculated. Then, the relationship between the elemental displacements υe(=[υi υj] and the internal forces fe(=[fi fj]) is defined as,
where κs(=κs[1 -1; -1 1]) is the local stiffness matrix, and τe is the coordinate transformation matrix from the global and local coordinate axes20). Thereafter, the internal forces in the reinforcement are transferred and distributed to the adjacent SPH nodes with accounting each kernel weight contribution as in Eq. (8). The bilinear plastic behavior is assumed to model the reinforcement behavior. The slippage between the reinforcing bar and concrete is not accounted and the perfect bond behavior is considered.
For numerical validations, an example of the direct tensile simulation of RC beam21) is selected. Fig.2 shows the planal and cross-sectional views of the RC specimen and boundary conditions. The axial movement is fixed at the left end of the reinforcing bar, and a displacement-controlled tensile load is applied at the right end of the bar. The loading velocity remains constant, and the magnitude is sufficiently low to satisfy the quasi-static loading condition.
The material properties for concrete and reinforcement are determined from the data provided in the previous numerical study21). The elastic modulus Ec and Poisson’s ratio ν of concrete are defined as 20 GPa and 0.2, respectively. For tensile softening behavior, the tensile strength ft and fracture energy GF are set as 2.4 MPa and 100 N/m, respectively. In this study, the elastic-perfectly plastic behavior of the reinforcement is asummed such that the hardening modulus H is set to zero. The elastic modulus Es and yield strength fy of the steel reinforcement are defined as 200 GPa and 500 MPa, respectively.
Fig.3 shows the RC beam model, discretized as 8000 particles. As described in Section 2.1, the reinforcing bar is implicitly represented within the solid domain, where the reinforcement is not directly shown in the figure. As depicted in Fig.2, one reinforcing bar is located at the center of the beam, which is discretized as the truss elements. The nodal forces and displacements at each truss element are updated through the associated physical values of the adjacent SPH particles.
Basically, the tension stiffening phenomenon is due to the dual contribution of the stiffnesses of concrete and steel reinforcement. As the tensile loading increases, both the concrete and the reinforcing bar are stretched. After the tensile strain exceeds the elastic limit of the concrete, the concrete damage is progressed according to the softening relation in Eq. (4). At the same time, the contribution of the reinforcing bar to the internal resultant force increases linearly with the external loads until plastic yielding occurs. Such tension stiffening behaviors can be found in Fig.4. The figure shows the numerical results of the longitudinal load and average strain curve, including the analytical result of the pure reinforcement model serving as a baseline in the curve. In addition, the other numerical result21) of the finite element analysis is provided for comparison. It is demonstrated that the numerical result in the present study follows well the general trend in the other numerical result. Such curve trajectories can also be found in the conventional results in the tension stiffening experiments22).
Figs.5a through 5d show the simulated results of the damage distributions at each time moment of the average strains of 1.01×10-4, 2.52×10-4, 5.05×10-4, and 1.01×10-3 mm/mm. In the figures, it is shown that the concrete damages are symmetrically distributed about the vertical center line of the beam, which demonstrates the fulfillment of quasi-static loading condition during the numerical simulations. Because the external loads are transferred to the RC beam through each end of the reinforcement, the associated damage is initially developed near each end of the reinforcement regions (Fig.5a). As the tensile loading increases, a region in the middle of the beam experiences quasi-uniform strain. Once the tensile strain exceeds the elastic limit, the concrete damages are initiated, leading to the formation of two distinct vertical cracks (Figs.5b and 5c). Thereafter, the tensile stress fields are developed in the remaining concrete regions, which contributes to the two more vertical cracks (Fig.5d). The figure shows that the left main crack is located at the center of the region ωl between the left boundary surface and first cracked surface, which holds for the right side of the beam due to the symmetric condition. These fractured regions across the cross-sectional area of the specimen are clearly shown in Fig.6. The similar pattern of beam fractures can be observed in both experimental22) and numerical21) studies.
In this study, the combined approach of the SPH method, isotropic damage model, and reinforcement model is newly proposed and validated for simulating the RC behaviors. The solid domain is modeled as the SPH particle arrangements, and the stiffness of the reinforcing bar is reflected through the mechanical behaviors of the surrounding SPH particles. Therefore, the reinforcing bar can be modeled without dependence on the geometric locations of the reinforcements relative to the SPH particle locations. The simulated results agree well with other numerical and conventional experimental results in terms of the force-strain curve and associated multiple cracking behavior.
Future studies should include more validations, such as numerical applications to the impact behaviors of RC structures and convergence studies based on particle size. Additionally, the output data within the system domain can be easily extracted and reorganized from the physical information of each SPH particle. Therefore, further research will explore the application of the proposed approach to the development of a data-based prediction model using the PINN.
This research was supported by Korea Institute of Science and Technology Information (No. (KISTI) K24L2M1C5)
