Orthotropic Hyperelastic Model of Concrete Deformation Under Combined Multiaxial Loading

Full Text

1. Introduction The main causes of nonlinear behavior of concrete are microcracking under tension and failure (microcracking and crushing) under compression, as well as concrete creep. The need to take nonlinear properties of concrete into account [1] arises when solving various problems related to the evaluation of strength of reinforced concrete structures. Therefore, in 2003, instead of SNiP 2.03.01-84[2] building code introduced in the USSR, new design standards for concrete and reinforced concrete structures - SNiP 52-01-2003[3]- were approved and put into effect. For the first time in the Russian Federation, these standards declare the analysis based on a nonlinear deformation model using stress-strain diagrams for concrete and reinforcement as the primary analysis method. This approach, based taking nonlinear deformation into account, was preserved in the updated version of SP 63.13330.2018[4], accepted in 2019. When describing the behavior of concrete in triaxial stress state, nonlinearity can be taken into account either in terms of nonlinear elasticity theory [1; 2-9] or plasticity theory [1; 9-25], as well as using material models that account for damage [9; 26-50]. Detailed reviews of models of elastic and inelastic deformation of concrete are presented in [1; 51]. A number of critical elements in reinforced concrete structures (near-support zones of box-section members of bridge superstructures, anchorage zones of prestressed reinforcement in bridge superstructure beams, joints between columns and slabs and beams, diaphragms of bridge superstructures, pylons of cable-stayed bridges, etc.) are subjected to multiaxial stress. The development of a three-dimensional, nonlinear, anisotropic model of concrete deformation remains an unresolved problem. Existing models of nonlinear deformation of concrete [10-13; 19; 31; 40; 52, etc.] are generally limited to the consideration of isotropic material. Orthotropy, which reflects the nature of microcracking in concrete, is taken into account, for example, in the model proposed by N.I. Karpenko [3]. A shortcoming of this model is the absence of elastic potential. In this study, the nonlinearity, multiaxiality, and anisotropy of concrete deformation are taken into account by using a model of locally orthotropic hyperelastic material with orthotropy axes coinciding with the directions of the principal stress axes. Due to the presence of potential, the model under consideration offers a number of advantages in the numerical solution of boundary value problems, such as symmetry of stiffness matrices and improved convergence of iterative procedures for finding nonlinear solutions. It should be noted that the hyperelastic model also possesses the property of thermodynamic consistency. Under uniaxial tension, the initiation and propagation of cracks during microcracking occur in the plane orthogonal to the direction of the applied stress. In a combined multiaxial stress state, the initiation and propagation of cracks during microcracking occur predominantly in the plane orthogonal to the direction of the maximum principal tensile stress. When there are two (or three) positive principal stresses, cracks form in two (or three) mutually orthogonal planes, the orientation of which is determined by the principal planes of the stress tensor [1]. The formation of a system of identically oriented microcracks leads to a change in the mechanical properties of concrete. Its stiffness in certain directions decreases. It can no longer be considered an isotropic body, and an appropriate approximation of its anisotropy is, in general, an orthotropic material. Local orthotropy is attained as the load increases due to a non-uniform stress state resulting from the directed, depending on the orientation of the principal stress planes, development of internal microcracks. In accordance with the principle of symmetry formulated by P. Curie - “the symmetries of the causes are to be found in the effects” [53] - the constitutive equations of the material must account for the elements of symmetry in the microstructure of concrete. Orthotropy has been taken into account using an elastic material model in [2; 4; 6-8, etc.], and when using plasticity and damage models in [41; 47-49, etc.], but it should be noted that elastic potential was not introduced in any of the cited works. The second factor contributing to the accuracy of the nonlinear solution, along with the consideration of anisotropy, is the choice of a function that approximates the uniaxial experimental curves under compression and tension. Various empirically motivated analytical approximations of stress-strain diagrams of concrete have been proposed in [54-63, etc.]. A detailed review of research in this area is presented in [51; 60; 62]. Examples of such relationships for describing stress-strain curves for short-term uniaxial compression include: ¡ Eurocode recommendation (EN 1992-1-1[5]) , (1) where is the peak stress (strength) of concrete in compression; is the strain corresponding to the peak stress ; ; is the initial (secant) modulus of elasticity of concrete, determined at the stress level of . ¡ N.I. Karpenko’s approximation [1] , (2) where is the initial modulus of elasticity of concrete in compression; is the value of the rate of change coefficient of the secant modulus at the top of the diagram; for the ascending branch , , and for the descending branch - , . In expression (2), the plus sign corresponds to the ascending branch of the stress-strain diagram of concrete, and the minus sign corresponds to the descending branch. ¡ Exponential-power approximation of G.V. Murashkin [60] , (3) where coefficients A, b, c are determined from conditions , , . ¡ A.V. Benin’s approximation [59] , (4) where, is the residual strength, and are the stress and strain at the beginning of crack formation. Along with uniform analytical relationships for the entire diagram or a portion of it (1)-(4), etc., piecewise-linear approximations have become widely used in numerical analysis. The simplest variants are bilinear and trilinear approximations. With a large number of data points, multilinear approximations provide highly accurate description of the analytical relationship or experimental data (which are typically obtained in tabular form). A piecewise-linear approximation of the stress-strain diagram allows a mathematical representation , (5) where is the Heaviside function; Ai are the constants determined based on the slopes of stress-strain diagrams for tension and/or compression. The aim of this study is to develop a general orthotropic hyperelastic model of nonlinear deformation of concrete under uniaxial and multiaxial loading, with the ability to simplify the identification of material parameters based on piecewise-linear approximations of the stress-strain diagram. 2. Methods 2.1. Constitutive Equations of Orthotropic Hyperelastic Model The development of a general three-dimensional nonlinear anisotropic model of concrete deformation under combined multiaxial loading has not yet been fully resolved. Existing models of inelastic concrete deformation are limited to the case of isotropic material. Orthotropy, which reflects the nature of directional microcracking in concrete, is taken into account in the model by N.I. Karpenko [1-3]. One of the shortcomings of this model is the absence of elastic potential. This study attempts to account for the nonlinearity, multiaxiality, and anisotropy of concrete deformation by using a model of locally orthotropic hyperelastic material with orthotropic axes coinciding with the directions of the principal stresses. In general case, when analyzing the combined multiaxial stress state of reinforced concrete structural members, it is necessary to use the general three-dimensional form of the constitutive equations, which represents the relationship between the six components of the stress tensor and the six components of the strain tensor, which becomes nonlinear at high stress levels. One of the rational material models for describing physically nonlinear deformation of concrete is the orthotropic hyperelastic material. Hyperelasticity implies the possibility of introducing scalar potential (free energy, Helmholtz thermodynamic potential) for the strain tensor, which allows the stress tensor to be defined as its derivative , or in the form of direct tensor calculus used hereinafter . (6) Introduction of a convex potential ensures the existence, uniqueness of the solution, and the possibility of using efficient finite-element methods to solve (elliptic) boundary value problems. For isotropic material, the potential can be considered as a function of the invariants of the strain tensor. For linearly elastic material, the potential is a quadratic function of the first and second invariants. The symmetry group of orthotropic material includes reflections with respect to three orthogonal planes, characterized by the normals to the microcrack growth planes . For cracked concrete, these directions correspond to the directions of the principal stresses and may vary depending on the combined non-proportional loading program. It is assumed that elastic potential ψ for an orthotropic material is a function of the invariants of the strain tensor and the second-rank tensors characterizing the anisotropic properties of the material [64]: ; (7a) ; (7b) , (7c) where symbol represents dyadic (tensor) product of vectors. When solving specific problems, it is sufficient to consider only two anisotropy tensors and , since the following relationship holds: , where is a unit tensor. Given this, 7 invariants can generally be the arguments of the elastic potential: , , ; (8) , ; (9) , . (10) In the orthotropic hyperelastic material model, the definition of the elastic potential (11) fully determines the form of the constitutive equations and the elastic moduli tensors. In the case of small strains, the simplest representation of elastic potential ψ is the quadratic form of the components of the strain tensor: , (12) where λ and μ are the Lamé parameters; are the material constants. Isotropic material corresponds to the case when . Elastic potential (12) is a generalization of the variant proposed in [65], which takes into account invariants and . To take into account the physical nonlinearity of concrete behavior, it is proposed to generalize (12) by introducing nonlinear functions of the invariants: , (13) where are twice differentiable piecewise-continuous functions. The constitutive equation, which expresses the relationship between stresses and strains, is obtained by substituting (13) into (6): (14) Differentiating this expression yields a quasi-linear relationship between the stress rates and strain rates: , (15) where the tensor of tangent elastic moduli of rank 4 (16) can be expressed as: (17) The symbols for direct and cross dyadic (tensor) products are used here: , . (18) In the case of piecewise-linear approximation of stress-strain diagrams (5), for i = 1,2; for j = 3, 4,…,7 and, as a consequence, , which leads to substantial simplification of expression (17): (19) In this case, the tensor of tangent moduli will be piecewise-constant. As follows from (19), the nine independent components of the tangent tensor of elastic moduli are associated with (along intrinsic anisotropy axes ) the elastic potential constants introduced in (12) by the relations (20) The elastic potential constants can be expressed in terms of the components of the tangent tensor of the elastic moduli as the solution to the linear system of algebraic equations (20): (21) Taking (20) into account, expression (15) can be written in matrix form: . (22) The structure of the matrix in (22) allows to clarify the physical significance of the introduced material constants. Constants and characterize the influence of microcracking processes on the relationship between normal stresses (acting along and ) and coaxial longitudinal strains. Constants and characterize the influence of microcracking processes on the relationship between shear stresses (acting in planes with normals and ) and the corresponding shear strains. Constants and characterize the influence of microcracking processes on the relationship between the normal stresses and transverse strains. The reciprocal form of constitutive equations (22), which allows the strain rates to be expressed in terms of the stress rates, is as follows: , (23) where the introduced coefficients in the tangent compliance matrix (23) (which represent the tangent Young's moduli in the three directions along orthotropy axes , three Poisson’s ratios and three shear moduli ) are expressed in terms of the elastic potential coefficients using the relations obtained by inverting the upper-left 3×3 block of the matrix of elastic moduli from (22): (24) In this respect, the symmetry conditions for the compliance matrix (23) are assumed to be satisfied: . (25) 2.2. Identification of Material Parameters and Functions In general case, to determine the 9 material parameters that appear in the elastic potential expression (13), it is necessary to conduct 6 experiments for hard or soft loading: 3 for uniaxial tension/ compression along 3 different axes of orthotropy and 3 for shear in different orthogonal planes. Such methods for determining material parameters in full-scale or computational experiments are typical for composites and materials with technologically induced orthotropy and are unlikely to be applicable to concrete and similar brittle materials. At the same time, it is necessary to take into account the specific nature of the difference in resistance under tension and compression (sensitivity to the type of stress state). The simplest case of parameter identification, where tension and compression curves are only available for uniaxial loading, along with data on transverse strain, is considered further. For uniaxial tension along the x-axis, the principal stresses are determined as follows: , . Normal to the plane with the maximum principal stresses is oriented along the x-axis. It is clear that under uniaxial tension, the transverse strains will be negative and equal to each other: . Function is assumed to be known based on experimental data. The expressions for the invariants for the case of tension are the following: , , , , , . (26) Considering the assumptions made, the system of nontrivial equations (22) in the case where the normal to the plane of maximum principal stresses is oriented along the x-axis, and lies along the z-axis, is transformed into: (27) where is the piecewise-constant function of the tangent moduli, determined from the stress-strain curve for uniaxial tension. Due to symmetry (indistinguishability of axes y and z under uniaxial tension from the natural state along the x-axis, corresponding to ) in addition to (27), the following conditions should be added:,, which lead to equalities (28) Parameters and are determined based on the initial values of the Young’s modulus and Poisson’s ratio of an initially isotropic body at zero strain. Due to the limited availability of experimental data, the lack of equations is compensated by making assumptions about approaching the isotropic case in the form of . Then, from equations (27) and (28), all the basic material constants for the tension case are found: (29) For uniaxial compression along the x-axis, the principal stresses are determined by the following equations: . Normal to the plane with the minimum principal stresses is oriented along the x-axis. Under uniaxial tension, the transverse strains will be positive and equal to each other: . Function is assumed to be known based on experimental data. The expressions for the invariants under compression are the following: , , , , , . (30) Considering the assumptions made, the system of nontrivial equations (22) in the case when the normal to the plane of minimum principal stresses is oriented along the x-axis, and - along the z-axis, transforms into: (31) where is the piecewise-constant function of the tangent moduli, determined from the stress-strain diagram for uniaxial compression. Due to symmetry (indistinguishability of axes y and z under uniaxial compression along the x-axis, corresponding to ) in addition to (31), the following conditions should be added: , which lead to equalities (32) Due to the limited availability of experimental data, the lack of equations is compensated by making assumptions about approaching the isotropic case in the form of . Then, from equations (31) and (32), all the basic material constants for compression are found: (33) The nonzero coefficients obtained in (29) and (33) are piecewise-constant functions. Invariants serve as the parameters that define the interval within which the coefficients remain constant: (34) where is the Heaviside function, are the Macaulay brackets. The proposed model of orthotropic hyperelastic material was implemented in CES version 5.26 (Constitutive Equation Studio) [66], using which the results presented below were obtained. 3. Results and Discussion The integration of the constitutive equations in the case of rigid loading is performed directly based on (22). In the case of soft loading, when the components of the stress tensor are specified and the components of the strain tensor need to be determined, it is necessary to invert the matrix of tangent moduli. The structure of the tangent modulus matrix (22), which contains a large number of zero components, corresponds to a representation in the coordinates of the axes of intrinsic anisotropy . They are governed by microcracking processes and are not directly related to the directions of the principal stress axes. Therefore, in the general case of arbitrary multiaxial loading, solving these problems requires preliminary determination of the principal axes of the stress tensor. Along the axes of intrinsic anisotropy, the tangential components of the stress tensor are absent, and for a simplified set of constants (see Section 2.2, (29), (33) and (34)) when and , the following representation is obtained: . (35) In the considered case of piecewise-constant coefficients (34), equation (35) can be easily integrated. This results in piecewise-linear approximations of the stress-strain diagram (5). 3.1. Uniaxial Loading Figure 1 shows a comparison of the results calculated using the orthotropic hyperelastic material model (6)-(35) with experimental data [59] for uniaxial compression and tension of B25 and B40 concrete. The results were obtained for the case of hard loading, for which it is possible to get both ascending and descending (post-peak) branches of the stress-strain diagram. A total of 32 points were specified as input data for generating a complete stress-strain diagram, including both the compression curve and the tension curve, taking into account the ascending and descending branches. In this case, it is possible to obtain a design curve with a good degree of accuracy. The characteristic values of the parameters for B25 and B40 concrete [59], used in the calculations, are presented in Table. Надпись: σ, MPa Надпись: Figure 1. Comparison of calculation results for the orthotropic hyperelastic material model with experimental data under uniaxial compression and tension for B25 and B40 concrete S o u r c e: made by A.S. Semenov. Experimental data obtained by A.V. Benin [59]. Надпись: σ, MPa Надпись: a Надпись: b Figure 2. Stress-strain diagrams for B25-concrete under: a - tension; b - compression, obtained on the basis of the orthotropic hyperelastic material model with a different number of segments N of a piecewise linear approximation S o u r c e: made by A.S. Semenov. Values of parameters of B25 and B40 concrete used in calculations Concrete parameters В25 В40 , MPa 18.5 29.1 , MPa 1.55 2.1 , % 0.199 0.201 , % 0.0123 0.0124 , MPa 30000 36000 , MPa 9343 1455 , MPa 12602 16934 S o u r c e: made by A.V. Benin, A.S. Semenov. Надпись: σ, MPa Depending on the number of segments N, the stress-strain diagram yields approximations with varying degrees of accuracy. Figure 2 shows multi-segment approximations of the stress-strain diagram for B25 concrete in tension and compression for various values of doubling number of segments. A comparison of the calculation results for B25 concrete in compression using the hyperelastic orthotropic model (6)-(35) with N.I. Karpenko’s orthotropic model (2) and experimental data [59] is shown in Figure 3. On the ascending branch of the diagram, both models demonstrate high accuracy in describing the curve; however, on the descending branch, a systematic deviation is observed for the N.I. Karpenko’s model. The use of a multi-segment piecewise-linear approximation of the stress-strain curve ensures high accuracy in approximating the hyperelastic orthotropic model. A 16-point approximation was used to construct the compression curve. 3.2. Multiaxial Loading Here, monotonic triaxial compression of concrete with varying degrees of transverse confinement is examined. The ratios of the three axial components of the stress tensor took the following values: · 10%-confinement with · 6.9%-confinement with Experimental data corresponding to the specified loading conditions are presented in the paper by Yu.N. Malashkin [67]. The stress-strain diagrams for the concrete under consideration are taken from [1]. A comparison of the results calculated using the hyperelastic orthotropic model (6)-(35) with N.I. Karpenko’s orthotropic model (2) and experimental data [67] for various degrees of confinement is shown in Figures 4 and 5. Both orthotropic models demonstrate satisfactory predictive accuracy (deviation from the experiment of less than 15%) compared to the experiment. 4. Conclusion As an alternative to the currently dominant empirical approach to formulating the constitutive equations for concrete, a thermodynamically consistent, locally orthotropic, hyperelastic model is proposed to describe the nonlinear deformation of concrete under monotonic multiaxial load, taking into account the tension-compression anisotropy. Local orthotropy develops as the load increases due to the growth of internal microcracks, the direction of which is determined by the orientation of the principal stress planes. Hyperelastic potential was defined based on a piecewise quadratic approximation using 7 joint invariants of the strain and anisotropy tensors. Надпись: σ, MPa Надпись: Figure 3. Comparison of results calculated using the hyperelastic orthotropic model with the orthotropic model of N.I. Karpenko and experimental data for uniaxial compression S o u r c e: made by A.S. Semenov. Experimental data obtained by A.V. Benin [59]. Надпись: Надпись: σ, MPa Figure 4. Comparison of results calculated using the hyperelastic orthotropic model with the orthotropic model of N.I. Karpenko and experimental data for multiaxial loading S o u r c e: made by A.S. Semenov, A.V. Benin. Experimental data obtained by Yu.N. Malashkin [67]. Надпись: Надпись: σ, MPa Figure 5. Comparison of results calculated using the hyperelastic orthotropic model with the orthotropic model of N.I. Karpenko and experimental data for multiaxial loading S o u r c e: S o u r c e: made by A.S. Semenov, A.V. Benin. Experimental data obtained by Yu.N. Malashkin [67]. In general case, the model contains 9 material parameters. A method is proposed for simplified identification of material parameters based on experimental stress-strain diagrams for compression and tension, with measurements of longitudinal and transverse strains. The obtained results allow to draw the following conclusions: 1. A comparison of the results obtained using the proposed orthotropic hyperelastic deformation model with experimental data for B25 and B40 concrete demonstrated good prediction accuracy under uniaxial loading (deviation from the experiments of less than 1%) and multiaxial loading (deviation from the experiments of less than 15%). 2. The proposed model allows to describe both the ascending and descending branches of the concrete stress-strain curve. 3. When describing the behavior of concrete on the descending branch of the strain diagram, the proposed model demonstrated higher accuracy compared to N.I. Karpenko’s model of orthotropic elastic material. Further development of the model involves explicit consideration of damage accumulation during loading and development of methods for identifying material parameters based on experimental data obtained under multiaxial loading.

About the authors

Andrey V. Benin

Emperor Alexander I St. Petersburg State Transport University

Email: nich@pgups.ru
ORCID iD: 0000-0001-5646-0354
SPIN-code: 8251-4345

Candidate of Technical Sciences, Head of the testing laboratory “N.A. Belolyubsky Mechanical Laboratory,” Associate Professor of the Department of Mechanics and Strength of Materials and Structures

9 Moskovsky Pr., Saint Petersburg, 190031, Russian Federation

Artem S. Semenov

Peter the Great St. Petersburg Polytechnic University

Author for correspondence.
Email: semenov.artem@googlemail.com
ORCID iD: 0000-0002-8225-3487
SPIN-code: 8879-2244

Doctor of Physical and Mathematical Sciences, Professor of the Department of Mechanics and Control Processes

29 Polytechnicheskaya St, Saint Petersburg, 195251, Russian Federation

References

Supplementary files