Double scalar elastic damage constitutive model considering shear damage

The concrete open-web sandwich slab is composed of the vierendeel girders intersecting each other. The vierendeel girder are composed of the upper rib, lower rib and shear key. The upper and lower rib mainly bear the axial force, and the shear keys mainly bear the shear force. A double scalar elastic damage constitutive model considering shear damage was derived to accurately simulate the mechanical property of the shear key in which one scalar reflected the volume damage, and the other one reflected the shear damage. The damage constitutive model in the study, the uniaxial tension and compression damage constitutive model (GB50010 (2015)), and the double scalar elastic damage constitutive model of tension and compression (Liu (2021)) were compared by simulating the shear key. The results showed that the elastic damage constitutive model proposed in the study can be used for the mechanical property analysis of the shear key and can better evolve the damage development of the shear key. The damage constitutive calculation model can be provided for the analysis of the vierendeel girder and the open-web sandwich slab.


INTRODUCTION
In recent years, due to the high integration of industrial production and the requirement of saving land resources, the demand for multi-story large-span building structures has increased rapidly.In the 1980s, the reinforced concrete open-web sandwich slab structure was proposed by Ma et al. (2006)(Figure 1), which solved this problem well.The reinforced concrete open-web sandwich slab is composed of vierendeel girders intersecting each other, which has the characteristics of lightweight, low steel consumption, and small structure height (Heilongjiang Provincial Department of Housing and Urban-Rural Development, 2014).The basic components of the vierendeel girder are the upper rib, the lower rib, and the shear key (Figure 2(a), 2(b)).The shear key and the lower rib are cast-in-situ concrete as a whole, the upper section of the shear keys are reserved for the construction joints.After the erection of the upper rib formwork is completed, the upper rib concrete is poured.The effective floor height of the building is increased for the electrical and plumbing conduit can pass through the gap between the upper and lower ribs.The upper and lower ribs of the vierendeel girder mainly bear the axial force, and the shear keys mainly bear the shear force (Figure 2(c)).The stiffness and force condition of each component of the vierendeel girder had different effects on its overall force state and deformation, which would further affect the mechanical properties of the open-web sandwich slab (Ma et al. 2000 andWang et al. 2019).As an important component connecting the upper and lower rib, the shear key performance played an important role in the bearing capacity of the vierendeel girder and the open-web sandwich slab (Fan et al. 2020).Concrete as an important engineering building material, the research on the constitutive model of which has never stopped.Some scholars had done a lot of research on isotropic damage models by defining one or two scalar damage variables (Kim and Abu Al-Rub. 2011, Zhou et al. 2013, Zhao et al. 2018).In addition, the anisotropic damage model was studied by incorporating second order, fourth order or eighth order damage tensors (Abu Al-Rub and Kim. 2010, Badreddine et al. 2015, Badreddine and Saanouni. 2017).Due to damage tensor is difficult to implement in nonlinear analysis, many authors used scalar isotropic damage variable because of its simplicity.
In the 1980s and 1990s, the elastic damage constitutive model of concrete was studied by many scholars.The Ladevèze-Mazars model (Mazars. 1986) was the most representative, which gave the whole damage formula comprised of the tensile and compressive damage.Later, the model was revised and improved by some scholars, among which the representative elastic damage constitutive models were established by Papa and Taliercio (1996), Comi and Perego (2001).Babu et al. (2010) proposed a phenomenological constitutive model for elastic damage of concrete and selected scalar damage parameters to quantify the damage.In recent years, there were few studies on the elastic damage constitutive model of concrete, and most of the studies focused on the elastoplastic damage constitutive model (Li et al. 2017, Huang et al. 2017, Feng et al. 2018, and Kim et al. 2021).In addition, the research object was not only concrete but also included rock (Wang et al. 2017, Wang et al. 2018).After introducing strain softening, the mechanical behavior of quasi-brittle materials such as concrete and rock can also be better described.Although the elastic damage constitutive model does not consider plastic deformation, the calculation formula of this model is simple, the damage evolution parameters are few, and the calculation efficiency is high, and which makes the calculation accuracy have a certain guarantee and veracity.Therefore, the elastic damage constitutive model is still having great significance for the response calculation of large concrete structures under complex load conditions.
With the deepening of research, people realized that the essential cause of material damage is the energy dissipation during the generation and the development process of micro-cracks.In general, under the action of external Latin American Journal of Solids and Structures, 2021, 18(6), e387 3/19 load, the deterioration of material properties is an irreversible thermodynamic process.Therefore, irreversible thermodynamics related theories have a pivotal position in the development of damage theory and the establishment of damage constitutive models.Leng and Wu (2013) established a concrete damage constitutive model based on thermodynamics, which can better describe the tensile damage and compression damage of concrete, and it follows the second law of thermodynamics.A novel four-parameters damage constitutive model for concrete in strain space was proposed in the framework of the thermodynamics of irreversible process by Qi et al. (2020), which can be used to determine the damage of materials under uniaxial and multiaxial loading conditions, and describe the non-linear behavior of materials after damage.The study of Wu et al. (2006) showed that during tensile loading the tensile damage and shear damage occurred and during compressive loading, the compressive consolidation and shear damage occurred.Ahmed et al. (2020) proposed a new elastic-plastic-damage constitutive model for concrete which considered six damage criteria to account for each type of damage, and the corresponding damage threshold for shear and pure tension and compression were derived from the uniaxial tension and compression damage threshold.
Liu (2012) derived a double scalar elastic damage constitutive model considering tensile and compressive damage through damage mechanics and continuum thermodynamics.Liu and Shen (2021) optimized the formula based on Liu (2012).The elastic damage constitutive model was established by the UMAT interface in Abaqus and the bearing capacity of the shear key of the vierendeel girder was analyzed.The shear keys of the vierendeel girder were mainly shear deformation (Jiang et al. 2019 andShen et al. 2020).It was difficult to reflect the mechanical performance of the shear key truly with the aforementioned double scalar elastic damage constitutive model.Therefore, to analyze the shear capacity of shear keys more accurately (Shang et al. 2019), an elastic damage constitutive model dominated by shear deformation was required.
The objective of the present paper is to establish the elastic damage constitutive model of concrete considering shear damage, and apply it to the force analysis of shear keys.Based on the work of Liu (2012), Liu and Shen (2021), a double scalar elastic damage constitutive model considering shear damage is derived in Section (2), in which one scalar is used to reflect the volume deformation damage, and the other is used to reflect the shear deformation damage.In Section (3) the FORTRAN77 is used to compile the UAMT subroutine in the interface of Abaqus, and the damage constitutive subroutine flowchart is given.In Section (4) through the normal stress-line strain curve of C45 concrete in uniaxial tension and uniaxial compression in the Chinese Code for Design of Concrete Structures (China Building Industry Press, 2015), and the shear stress-shear strain curve of Guo (2013), the damage evolution parameters of the double

Basic formula
Base on the small deformation hypothesis, the Helmholtz free energy function 0  is defined as follows when there is no damage: where: 0 D is the initial elastic stiffness tensor,  is the mass density, ε is the second-order strain tensor.
When damage occurs, Eq. ( 1) is defined as a function of strain tensor ε and damage scalar  , as follows: Calculating derivative of Eq. ( 2) and bring it into the Clausius-Duhem inequality, the outcome is: where: σ is the second-order stress tensor.
It can be known from Eq. (3): The damage energy release rate is defined as: The non-damage free energy is defined as Eq.(1), which is divided into two parts: volume effect and shear effect: where: 0 K is the initial bulk modulus, G are the initial elastic modulus and the initial shear modulus respectively;  is the Poisson's ratio; v  is the body strain ( ); e is the strain deviator.
Latin American Journal of Solids and Structures, 2021, 18(6), e387  5/19 To characterize the difference in damage between volume deformation and shear deformation, the free energy of damage is defined as follows: where: 1  reflects volume deformation damage,

2
 reflects shear deformation damage.Eq. ( 7) is reorganized to obtain: Eq. ( 4) and ( 5) are combined on the definition of stress tensor and damage energy release rate, the stress tensor can be obtained as: where: I is ij


. The damage energy release rate is: The total damage scalar  is defined as follows: Due to 0 D is constant tensor; 0 K is constant; the damage scalars are functions of strain tensor; and combining the expression of v  , calculating derivative of Eq. ( 9), the rate form of the stress tensor was written: where the continuous consistency tangent modulus D is: According to the chain derivation rule, we have: The evolution expressions of volume damage variable 1  and shear damage variable 2  are defined as respectively: Latin American Journal of Solids and Structures, 2021, 18(6), e387 6/19 where: Taking the partial derivative of Eq. ( 15), we obtain: The partial derivative of Eq. ( 10): Substituting Eqs. ( 16) and ( 17) into Eq.( 14):

Derivation of the relationship between damage variables
1 2 ,   Writing Eq. ( 9) in the form of tensor expansion: Expanding the elastic Jacobian matrix 0ijkl D , which can be written form in the following: where: Substituting Eq. ( 20) into Eq.( 19), we obtain: To simplify the expression, letting: Eq. ( 21) can be organized as: Expanding Eq. ( 10) damage energy release rate 2 Y , organizing it to get:

Damage evolution
The damage constitutive model considers volume deformation and shear deformation because the two types of damage evolution are not the same, so their damage criteria are also different.The hardening factor is considered in the damage loading surface, then the volume damage loading surface and the shear damage loading surface are defined as: where: 1 R is volume damage hardening function; 2 R is shear damage hardening function, defined as: Latin American Journal of Solids and Structures, 2021, 18(6), e387 11/19

Damage flow law and consistency conditions
The damage flow law defines the evolution criterion of damage variables and hardening parameters, and the damage development direction is consistent with the normal direction outside the damage loading surface, can be written as: where: 1 2 ,     are the volume and shear damage multipliers respectively.Through the damage loading surface to determine the damage: 1 0 H  : volume damage develops, the volume damage variable is updated; 1 0 H  : volume damage does not develop, the volume damage variable is not updated; 2 0 H  : shear damage develops, the shear damage variable is updated; 2 0 H  : shear damage does not develop, the shear damage variable is not updated.Substituting the expression of  into Eq.( 54), the relationship between 1 2 ,   is as follows after sorting out:

Numerical realization when
To facilitate iteration, letting: A F I J into Eq.( 57), sorting out, and letting: Latin American Journal of Solids and Structures, 2021, 18(6), e387 12/19 Calculating derivative of the above Eq.( 58): Using the Newton iteration method to find 2  , the formula is: Substituting Eqs. ( 58) and ( 59) into Eq.( 60), the iterative formula of 2  can be obtained.

Numerical realization when
In double scalar damage considering shear effect under shear force loading, shear damage 2  is the main damage, and volume damage 1  is secondary damage, there is a connection between the two damages.It can be obtained 1  from 2  according to the Newton iteration method.To facilitate iteration, letting: , , , , , A F I J into Eq.( 57), sorting out, and letting: Calculating derivation of the above Eq.( 61): Using the Newton iteration method to find 1  , the formula is: Substituting Eqs. ( 61) and ( 62) into Eq.( 63), the iterative formula of 1  can be obtained.

Program process
Generally, once the iteration step starts, the program will interact with data through the UMAT interface.The initial value of the state variables of the current element integration point can be transferred to the corresponding variables via UMAT interface.At the end of UMAT, the updated value of the state variables will be returned to Abaqus through the UMAT interface.
The flowchart of the UMAT subroutine is shown in Figure 4. Volume damage or shear damage should be selected as the dominant damage in damage constitutive subroutine after calculating and judging the numerical magnitude relationship between energy release rate 1 Y , 2 Y .uniaxial compression loading.To ensure the accuracy of the analysis, it was finally determined that the mesh size used was 10mm after the mesh sensitivity analysis.To avoid the problem of hourglass control, the C3D8R unit was used and set to enhance the hourglass control.After a lot of trial calculations and adjustments, the material parameters and evolution parameters of uniaxial tension and compression were given, as shown in Table 1.
Table 1 Material properties and evolution parameters of uniaxial tension and compression considering shear damage.Figure 6 shows the uniaxial tensile stress-strain curve calculated by the UMAT subroutine, the uniaxial tensile curve of China Building Industry Press (2015), and the uniaxial tensile damage evolution calculated by UMAT.It can be seen from Figure 6 that in the ascending phase it is in good agreement with China Building Industry Press (2015), but there is almost no descending phase after the peak load.The peak stress calculated by the UMAT subroutine is basically the same as the value given by China Building Industry Press (2015), and at this time, the corresponding strain is about 1.0×10 -4 , and the corresponding tensile damage variable develops to 0.75.In the double scalar damage constitutive model considering shear damage, the damage development is relatively complete and the damage development speed is relatively fast in uniaxial tension loading, so there is no decline section basically.To consider the overall stress field of the model when the uniaxial tensile stress reaches the peak stress, a comparison of stress S11 calculated by constitutive theory in China Building Industry Press (2015) and UMAT subroutine was given, as shown in Figure 7. Since there is a certain conversion relationship between tensile stress and tensile damage, the tensile damage distribution can be indirectly obtained from the tensile stress distribution, to predict the area and location of the first damage.It can be seen from Figure 7(a) that using the constitutive calculation model in China Building Industry Press (2015), the stress S11 of the model is mostly around 2.31MPa; while in Figure 7 Figure 8 shows the uniaxial compression stress-strain curve calculated by the UMAT subroutine and the constitutive theory in China Building Industry Press (2015), as well as the damage evolution of uniaxial compression.At the stage when the strain is small, the stress-strain curve calculated by UMAT is in good agreement with China Building Industry Press (2015), and the calculated strength value of UMAT is slightly lower than the calculated value of China Building Industry Press ( 2015), but after the peak, the UMAT calculation has only one short phase.In the UMAT calculation, when the strain corresponding to the peak stress is 1.81×10 -3 , the damage variable is 0.75; when the strain reaches 2.14×10 -3 , the damage variable is 0.9.The damage evolution changes rapidly before the stress peak point and has a slowing down after the stress peak point, but overall, the damage development is relatively complete.Figure 9 shows the comparison of the stress S11 distribution calculated by the UMAT subroutine and constitutive theory in China Building Industry Press (2015) when compressive stress reaches the maximum value.It is generally believed that the damage starts at the location of the maximum stress, according to the compressive stress distribution, the compression damage distribution can be obtained indirectly, and then predict the damage location and damage zone.It can be seen from Figure 9(a) that when the peak stress is reached, the stress S11 of the model almost reaches the standard value of 29.60MPa for C45 concrete axial compression, but in Figure 9 A pure shear test model of 20mm×20mm×20mm was established, as shown in Figure 10.The displacement of the left section of the model in the three directions x, y, and z were constrained; the coupled point is set at the right section, the mesh size was 5mm, and a z-direction displacement of 0.5mm was applied.After many trial calculations and adjustments, the material properties and evolution parameters of the pure shear state were given, as shown in Table 2.  Guo ( 2013) obtained a shear stress-shear strain curve based on a large number of pure shear tests and regression analyses.Figure 11 shows the shear stress-shear strain curves calculated by Guo (2013) and UMAT subroutine, as well as the shear damage calculated by UMAT.It can be seen from Figure 11 that when the shear strain is less than 1×10 -4 , the two constitutive curves are almost the same.When it is greater than 1×10 -4 , the shear strain growth rate of Guo ( 2013) is accelerated, and the shear strain growth rate calculated by UMAT is smaller than that of Guo (2013), but the maximum shear stress calculated by UMAT is the same as Guo (2013) basically.The shear damage calculated by UMAT finally reached 0.21, which shows that the shear damage has not fully developed.

Uniaxial loading of shear key
After determining the material property parameters and damage evolution parameters of uniaxial tension, compression, and pure shear loading, according to Loland (1980) and Construction Department of Guizhou Province (2005), the unit cell of vierendeel girder was established.The force diagram of the model is shown in Figure 12 and the finite element model is shown in Figure 13.The displacement in the y, z directions of the left and right sections in the upper rib and the left section in the lower rib were constrained; the displacement in the x, y, and z directions of the right section in the lower rib was constrained.The right section of the upper rib was coupled and displacement loading system was applied at the coupling point with 10mm of the displacement in the negative x-direction.The C3D8R element was used and the element size was 100mm.The comparison of shear key shear forces calculated by the uniaxial tension-compression constitutive model of China Building Industry Press (2015), Liu and Shen (2021), and the UMAT subroutine in the study, are shown in Figure 14.It can be seen that at the primary peak point, the shear capacity of China Building Industry Press ( 2015) is 370kN, the peak bearing capacity calculated by Liu and Shen (2021) is 347kN, which is about 6.2% lower than China Building Industry Press (2015).While the shear capacity calculated by UMAT in the study is 413kN, which exceeds the calculation result of China Building Industry Press (2015) by about 11.6%.After the primary peak point, the three curves have declined to a certain extent.Among them, China Building Industry Press (2015) has the most decline, Liu and Shen (2021) has the least decline, and UMAT in the study is the second.When the displacement is less than 4.5mm, the curves calculated by UMAT and China Building Industry Press (2015) basically the same changes, and there are secondary peak points in the two curves, and the shear force at the secondary peak points differs by about 26kN, but there is no peak point in the latter stage of the curve of Liu and Shen (2021).In general, the elastic damage constitutive model by Liu and Shen (2021) cannot accurately analyze the shear key, such kinds of the main shear members, which will cause the shear capacity of the member to be low and overestimate the shear force after the peak point.The elastic constitutive model considering shear damage in the study improves the shear capacity of the shear key, and it has appropriate shear capacity after the primary peak point.Therefore, the damage constitutive model in this paper can better evolve the damage development of the shear key and has good applicability.

CONCLUSION
(1) The double scalar elastic damage constitutive model of concrete considering shear damage is derived through the basic theory of thermodynamics, of which, one scalar is used to reflect the volume deformation damage, and the other is used to reflect the shear deformation damage.And the data interaction with finite element software Abaqus was realized by the UMAT subroutine in UMAT interface.
(2) Based on the existing uniaxial tension, uniaxial compression and pure shear constitutive model, the material properties and evolution parameters of the double scalar elastic damage constitutive model of concrete considering shear damage are determined.The damage development of the damage constitutive model in the study is Latin American Journal of Solids and Structures, 2021, 18(6), e387 18/19 concentrated in the ascending phase before the peak point for uniaxial tension loading, uniaxial compression loading, and pure shear loading.
(3) The double scalar elastic damage constitutive model considering shear damage can be used for the stress analysis of concrete subjected to tension, compression, and shear deformation, which can evolve the damage development of the shear key and has good applicability.It provides a constitutive calculation model for the overall analysis of vierendeel girder and open-web sandwich slab.

Figure 1
Figure 1 Engineering examples: (a) Xingyi City Office Canteen Roof; (b) Natatorium of Guizhou Aluminum Factory Company.

YY
are the volume energy release rate, shear energy release rate; are the volume energy release rate threshold and shear energy release rate threshold respectively; a , b are the volume damage evolution parameters of material; m ,n are the shear damage evolution parameters of material.

Figure 3
Figure 3 Pure shear stress state.

1
must be different from the development of shear damage variable 2  , because volume damage and shear damage are essentially different.According to the different values of the energy release rate 1 Y and 2 Y , the compilation of the damage constitutive program was divided into two situations to reflect the different development of damage under different stress states.

1Y is greater than 2 Y 1 
In double scalar damage considering shear effect under uniaxial loading, volume damage 1  is the main damage, and shear damage 2  is the secondary damage, there is a connection between the two damages.It can be obtained 2  from according to the Newton iteration method.

Figure 4 Figure 5
Figure 4 Double scalar elastic damage constitutive subroutine considering shear damage program flowchart.

Figure 7
Figure 7 Stress S11 distribution comparison of uniaxial tensile peak point: (a) Stress S11 distribution calculated by China Building Industry Press, 2015; (b) Stress S11 distribution calculated by UMAT.
(b), most of the stress S11 of the model is 2.39MPa.The comparison between the two shows that given the rapid development of damage, the maximum stress of 2.46 MPa obtained by the constitutive calculation model considering shear damage is slightly less than the maximum stress of 2.55 MPa calculated by constitutive theory in China Building Industry Press (2015), but the overall calculation results of the two models are relatively close.Latin American Journal of Solids and Structures, 2021, 18(6), e387 15/19

Figure 9
Figure 9 Stress S11 distribution comparison of uniaxial compression peak point: (a) Stress S11 distribution calculated by (China Building Industry Press (2015); (b) Stress S11 distribution calculated by UMAT.

Figure 10
Figure9shows the comparison of the stress S11 distribution calculated by the UMAT subroutine and constitutive theory in China Building Industry Press (2015) when compressive stress reaches the maximum value.It is generally believed that the damage starts at the location of the maximum stress, according to the compressive stress distribution, the compression damage distribution can be obtained indirectly, and then predict the damage location and damage zone.It can be seen from Figure9(a) that when the peak stress is reached, the stress S11 of the model almost reaches the standard value of 29.60MPa for C45 concrete axial compression, but in Figure9(b), most of the stress S11 of the model is less than the standard value of C45 concrete axial compression, and its value is 26.20MPa.In general, UMAT calculation results are relatively close to China Building Industry Press (2015) calculation results.4.2 Pure shear loading

Figure 12 19 Figure 13
Figure 12 Force diagram of the unit cell of vierendeel girder.

Figure 14
Figure 14 Comparison of shear force of shear key.

2 DERIVATION OF DOUBLE SCALAR ELASTIC DAMAGE CONSTITUTIVE MODEL CONSIDERING SHEAR DAMAGE
constitutive model are determined.Further, the elastic damage constitutive model is used to analyze the shear resistance of the shear key.The results show that the elastic damage constitutive model in the study can better evolve the damage development of the shear key, and can be used for simulation of concrete under tension, compression, and shear deformation.The research in the study can provide a constitutive calculation model for the overall analysis of the vierendeel girder and the open-web sandwich slab.

Table 2
Material properties and evolution parameters under pure shear state.