Selected Concrete Models Studied Using Willam’s Test

Willam’s test is a quick numerical benchmark in tension–shear regime, which can be used to verify inelastic (quasi-brittle) material models at the point level. Its sequence consists of two separate steps: uniaxial tension accompanied with contraction—until the tensile strength is attained; and next for softening (cracking) of the material—tension in two directions together with shear. A rotation of axes of principal strains and principal stresses is provoked in the second stage. That kind of process occurs during the analysis of real concrete structures, so a correct response of the material model at the point level is needed. Some familiar concrete models are selected to perform Willam’s test in the paper: concrete damaged plasticity and concrete smeared cracking—distributed in the commercial ABAQUS software, scalar damage with coupling to plasticity and isotropic damage—both implemented in the FEAP package. After a brief review of the theory, computations for each model are discussed. Passing or failing Willam’s test by the above models is concluded based on their results, indicating restrictions of their use for finite element computations of concrete structures with predominant mixed-mode fracture.


Introduction
Sophisticated and complicated simulations for concrete structures are presented in numerous papers. It happens sometimes that it is not convincing whether an employed material model, which can be commonly used, satisfies all basic requirements for proper nonlinear analysis. Appropriate results should actually be guaranteed both for primary stress states and consequently for their complex combinations. First of all, a considered phenomenological model for quasi-brittle materials like concrete should be examined by simple benchmarks illustrating softening/cracking response e.g., for uniaxial tension, also at the point level.
One such benchmark is the so-called Willam's test which was originally devised in order to study differences between fixed and rotating crack models [1]. This test is still useful for identifying differences between various concrete models in the case when shear induced cracking is important. It has been remarked in [1,2] that the results of particular models can differ even if these models exhibit a similar behavior in uniaxial tension. This tension-shear test with two loading steps is helpful to prevent undesirable effects in more advanced computations because a rotation of principal directions of the strain and stress tensors is very often observed. If a model fails the Willam's test, then there is no guarantee that it can render properly the structural behavior of concrete elements with predominant mixed-mode fracture. The test itself is quite demanding and even well-known models can fail it, see, e.g., [3].

Phase I.
Uniaxial horizontal tension with vertical contraction due to the Poisson's effect, according to the relation between the strain increments: where 11 and 22 are horizontal and vertical strain components, respectively, γ 12 is shear strain and ν is Poisson's ratio. In Figure 1a, the scheme of prescribed displacements corresponds to the uniaxial tensile strain state. Such relation is valid until the tensile strength is attained. Phase II.
Immediately after the tensile strength is reached, the change of configuration is enforced, Figure 1b. Now, the proportions for the strain increments are arranged in the following way: This relation induces tension in two directions with additional shear effect. As a result of such combination, a rotation of principal strain axes occurs; however, tension regime is preserved. At the beginning, the rate of rotation is fast, but, during the evolution, it goes down gradually.  It is suggested by [2] that this two-phased loading process is observed during the analysis of real reinforced concrete (RC) structures. Nevertheless, in the literature, there seem to be no experiments for this or similar tests with rotating principal directions [7]. It turns out that Willam's test is difficult to be reproduced in a laboratory. However, different authors verified their own proposals in modeling of quasi-brittle materials by means of this numerical test; for example, a comparison of models with respect to multi-surface plasticity is given by [3].
It should be emphasized that Willam's test is passed if two conditions are satisfied:

Condition 1:
The maximum principal stress is lower than or at most equal to the given uniaxial tensile strength.

Condition 2:
All stress components should converge to zero at the final stage.
Both conditions are physically motivated. Condition 1 is obvious. Concerning Condition 2, the principal strains 1 and 2 have positive values during Phase II. Their directions rotate and in the limit angle Θ between the principal strain 1 and the horizontal strain 11 takes value 52.018 • . The test is performed in the plane stress state without confinement, hence material can freely shorten in the out-of-plane direction. Therefore, Willam's test in Phase II is basically a biaxial tension test with changing directions and no confinement. In such a situation, dilatancy does not appear and in the limit the principal stresses σ 1 and σ 2 should converge to zero, which justifies Condition 2.
In this study, four models (mentioned above) are taken into account. The essential theory for these models is concisely described in Section 2. The paper does not attempt to describe a current development of material models for concrete as well as their more general historical context, so the reader can consult e.g., the following sources on this topic [8][9][10][11][12]. The nonlinear analysis is confined to statics and the assumption of small strains. The main Section 3 shows the corresponding data and the results obtained for each used concrete model. Furthermore, a comparison of diagrams of principal stresses and their directions for these models is shown in Section 4. The last Section 5 contains final conclusions on passing or failing Willam's test by all the models.

Concrete Damaged Plasticity (CDP)
The most popular model to simulate concrete behavior, which is delivered by the ABAQUS software [4], is the concrete damaged plasticity (CDP) model. This model was originally proposed by [13] and then modified by [14].
The theory of non-associated plasticity is linked with isotropic damage, where distinction between stress states for tension and compression is admissible. This nonlinear behavior of the actual body (the material point) is represented by a relationship between stress tensor σ and strain tensor . Moreover, effective stress tensorσ is assigned here to the fictitious undamaged counterpart of this body. The above concept is proposed together with the postulate of strain equivalence =ˆ [15,16] in the real and effective configurations.
The yield function is defined in the effective stress space: whereq is the effective Mises equivalent stress,p is the effective hydrostatic pressure and σ max is the positive part of the maximum principal effective stress. Parameters A, B and C decide about the shape of the yield surface; however, ABAQUS users define the following strengths for concrete: the uniaxial tensile strength denoted as f t and f c , f c , f bc are maximum uniaxial, initial uniaxial, and initial biaxial strengths for compression, respectively. These coefficients affect parameters A, B and C as well as the yield criterion, so that the shape of F p CDP can be determined. Its initial form in 2D for the data applied in Willam's test (see Section 3) is shown in Figure 2. This yield surface has first been introduced in [13]. It matches the experimental data well for concrete in stress states with prevailing tension and dominant hydrostatic compression. There are a few other yield surfaces of similar quality widely used for quasi-brittle materials like concrete, cf. [17][18][19]. The current stress state is described via stress-strain relationships separately for uniaxial compression and uniaxial tension: The primary argument of the above yield function is equivalent plastic strain˜ p , split into˜ p c for compression and˜ p t for tension. Of course, the standard additive decomposition of strain rate into elastic and plastic parts is assumed in the CDP model: More in-depth description concerning the determination of nonlinear functions for compressive crushing and tensile cracking can be found e.g., in [4,20,21]. The plastic flow potential function is determined in the following way: where ψ is the dilatancy angle and the so-called eccentricity E decides about the shape of the tip of surface G p CDP . Figure 3 illustrates the influence of dilatancy angle ψ on function G p CDP described in thep-q plane. Different aspects of the dilatancy definition connected with the plastic potential are addressed e.g., by [22]. It should be noticed that both functions F p CDP and G p CDP are formulated as modifications of the classical Burzyński-Drucker-Prager (BDP) yield criterion. If damage phenomenon is additionally activated in this plasticity-based model to consider a stiffness degradation of concrete, then the stress σ for the real material structure corresponds to the effective stressσ in the undamaged skeleton of the body via the damage parameter ω. It should also be recalled that ω evolves from 0 for the intact material to 1 for a full loss of stiffness. The constitutive relationship can be written in the form: where D is the fourth order tensor for the initial elastic stiffness. In a similar fashion to the distinction of uniaxial stress-strain relations in compressive and tensile regimes, the stiffness degradation as well as its recovery factors can be determined independently. Damage ω is split in such way: Variable ω t reduces the stiffness in tension and consequently ω c -in compression, while functions s c (σ) and s t (σ) are responsible for stiffness recovery. A wider discussion of the damage definition in the CDP model and of the crack closing phenomenon can be found in [23].
This model in the ABAQUS/Standard version is equipped with the regularization by the viscous term for the scalar plastic strain rate and also for the rate of degradation, according to the generalized approach of [24]. The sensitivity analysis of the viscosity parameter is performed e.g., by [25]. However, the regularization is not regarded at the material point level, so the CDP model has no viscosity in Willam's test and relaxation time is equal to zero.

Concrete Smeared Cracking (CSC)
Using the ABAQUS/Standard version for simulations of concrete structures another nonlinear model can also be employed. It is called the concrete smeared cracking (CSC) model [4] and allows one to compute the problem under monotonic loading. Hence, it works in the nonlinear analysis of quasi-brittle materials but without the crack closing effect.
The concept of smeared cracking in concrete comes from the late sixties, cf. [26]. In the CSC model, crack orientation is normal to maximum principal stress σ I reaching tensile strength f t and then held on fixed. Consequently, the elastic stiffness becomes reduced. In the region of predominant tensile stress, a yield surface F p CSC,t is introduced, called here the crack detection surface: whereq andp are introduced similar toq andp, but without stress components of σ associated with open cracks-this way, secondary cracks can form only perpendicular to already existing ones (no more than two cracks perpendicular to each other in 2D case). Relationσ t (λ t ) for uniaxial tension governs softening of the yield surface and parameter T influences the shape of the yield surface. In order to control softening of function F p CSC,t , the classical formalism of the associated elastoplasticity is introduced with strain rate decomposition:˙ =˙ e +˙ p t (10) and the associated flow rule: In the above equation, p t is the plastic strain tensor for the crack detection and λ t is the tensile plastic multiplier.
It has to be stressed that the elastoplastic formulation described above is only used for a proper description of softening of the yield surface F p CSC,t . Immediately after crack detection, a damage elasticity approach is introduced linking total values of stress and strain tensors according to the formula: σ = D cr : (12) where D cr is a secant elastic operator formulated in the spatially fixed coordinate frame n,t,s aligned with the crack plane (with axis n normal to this plane). The normal component of the stiffness operator D cr in direction n is defined as: where the value of the stress σ nn for a given value of total strain nn is computed using the stress-strain relation for uniaxial tension, a so-called tension stiffening curve [4], defined by a user. Additionally, Poisson's effect is neglected for an active crack leading to the formula: The shear stiffness decreases for cracking, so Kirchhoff (shear) modulus G is reduced according to the shear retention concept [1,27]: Function ρ( nn ) is illustrated in Figure 4 for cases computed in Section 3.2 and defined as follows:  In the region of predominant compressive stress in the CSC model, the standard associated elastoplastic approach is used with the compressive yield surface in the classical BDP form: where S is a material parameter specified by the ratio of biaxial to uniaxial compressive strengths and hardening/softening of the yield surface is governed by the cohesion function τ c (λ c ). The standard strain rate elastoplastic additive decomposition is adopted: but now the associated plastic flow is: where p c is the plastic strain connected with compression, λ c is the compressive plastic multiplier, parameter W depends on the ratio of stress components for biaxial and uniaxial compression as well as the ratio of respective strains, and σ c ( c ) is the hardening/softening curve for uniaxial compression, cf. Equation (4) 1 .
The diagram σ c ( c ) serves as well to define the cohesion function τ c (λ c )-the details are given in [4]. The compressive yield surface F p CSC,c as well the crack detection surface F p CSC,t in their initial state are presented in Figure 5 based on the data for Willam's test (see Section 3). Hence, the yield criterion is specified in the stress space and consists of two different surfaces similarly to [28]. As can be seen from the above description that the CSC model belongs to the family of total smeared fixed crack models [1,27,29], enhanced by the standard elastoplastic formulation for compressive stress. More detailed explanations for the CSC model, e.g., how to recalculate values of parameters T , W and S from the data entered by users are included in [4]. It should be mentioned that, for this model, mesh sensitive results can occur. Therefore, if the CSC model is applied in the computations of concrete structures, then the crack band theory [30] should be employed to set a proper size of the finite element corresponding to the expected cracking zone, see also [31,32]. Wide discussion concerning smeared cracking models can be found e.g., in [27,29]. Moreover, the concept of rotating crack model has also been included in the Modified Compression Field Theory (MCFT) developed for a combined description of cracked concrete and reinforcement at the material point level, see [33]. However, in that case, the descending branch of the stress-strain curve for concrete represents the tension stiffening phenomenon for reinforced concrete (RC) rather than the tension softening for plain concrete.

Damage-Plasticity (DAP)
A combination of the damage theory defined in strain space with hardening plasticity given in stress space leads to the damage-plasticity (DAP) model. The model presented here is implemented in the FEAP package [5] and is based on works of [16,34,35].
Starting from the pure scalar damage model, one damage measure ω is introduced in the classical way. The effective stressσ is distinguished as in the CDP model: The assumption of strain equivalence also remains valid. Damage is characterized by the following loading function: where˜ ( ) is an equivalent strain measure and κ d is a damage history parameter. This damage activation function F d DAP controls the behavior of the material after the damage threshold κ o is attained during the loading history. The equivalent strain measure˜ , which should demonstrate different behavior in tension and compression, can be defined for instance according to the idea by [36]: where I is the positive part of I-th principal strain I . The second proposal, employed in the DAP model, is the modified von Mises definition [37]: where I 1 is the first strain tensor invariant, J 2 is the second deviatoric strain invariant, ν is the Poisson's ratio and k = f c / f t is the ratio between compressive and tensile strength. Both functions are depicted in Figure 6 according to the data used in Willam's test (see Section 3). The damage growth function directly depends on the damage history parameter κ d and can be determined e.g., as the exponential softening relation. It properly reproduces tensile fracture phenomenon in concrete by the asymptotic function [38]: Parameter α corresponds to the residual stress (1 − α) E κ o , so, if α = 1, the total loss of the stiffness is attained. The ductility parameter η is connected with the rate of softening and the concrete fracture energy G f . The elastic constant E is Young's modulus.
During unloading, the secant stiffness (1 − ω)D results in a return to the origin, i.e., no residual strains are observed and damage does not grow. The damage-based model can be coupled to plasticity in order to include physically observed irreversible strains.
The yield function for the plastic component in the DAP model is formulated in the effective stress space: whereσ(σ) is an equivalent measure of effective stress and σ y (κ p ) is an isotropic hardening law with the yield strength limit. Function F p DAP can be described, for example, by the Huber-Mises-Hencky (HMH) or BDP criteria. Proportionality between the plastic multiplierλ and the plastic strain measure κ p is assumed. The plastic multiplier λ determines the magnitude of plastic strains p according to the non-associated flow rule: where m is the plastic flow direction and G p DAP is a plastic potential function. Based on the standard additive decomposition (5) and considering the plastic consistency condition, the elastic strain rate is written as: where h is the hardening (or softening) modulus and the gradient tensor n = ∂F p DAP /∂σ is associated with the yield function F p DAP . Finally, using the Sherman-Morrison formula, the tangential relation is derived:˙σ = D ep :˙ with the elastoplastic tensor: It is seen that, in the DAP model, the coupling of damage and plasticity is given in Equation (20). The stress rate during the evolution of damage and plasticity can be computed as: The rate of damage during loading (κ d =˜ ) is calculated in the way: and during unloadingω = 0. The stress-strain relation for the coupled model is: where the following definitions are used: There are two possibilities of coupling in the DAP model. The equivalent strain˜ can depend on the total strain tensor or its elastic part e . If the full coupling via˜ ( ) is employed, then the plastic strains also stimulate the damage growth. When the second option˜ ( e ) is selected, the coupling effect is weaker, and it seems to be more relevant in the modeling of quasi-brittle materials.
Beyond the point level, the DAP model can be enhanced to a nonlocal version by means of a gradient-type or an integral-type approach in order to ensure mesh-objective results, see, e.g., [35,[39][40][41][42].

Isotropic Damage (IDA)
Material (elastic stiffness) degradation can be introduced into the model using a fourth-order damage tensor or a second-order damage tensor for anisotropic description, see, e.g., [43,44]. The decomposition of damage effect into different tensile and compressive actions is described by the CDP model presented above or for instance by [38,[45][46][47]. Remaining within the isotropic description a volumetric-deviatoric split with two damage variables and one equivalent strain measure can also be applied by [48].
Even a simpler upgrade of the scalar damage theory is originated in [49], where one damage parameter influences a different decrease of stiffness for bulk modulus K and shear (Kirchhoff) modulus G. Moreover, the magnitudes of degradation of K and G depend on the sign of dilatancy. This approach will be called the isotropic damage (IDA) model here. Firstly, the elasticity operator D is written as: where I is the second order identity tensor, Q = I − I ⊗ I / 3 is a fourth order tensor and I is the fourth order identity tensor. Combining Equations (20) and (34), the following constitutive equation is obtained: Now, the deviatoric strain is dev = Q : and the dilatancy is θ = I : , so Equation (35) is rewritten: Referring to [49], where the anisotropic damage description is shown, positive and negative parts of hydrostatic pressure and stress components are distinguished. It is indicated that damage is connected with a micro-defect (micro-crack) pattern which evolves in a different way under tension and under compression. The action of damage can be reduced for the negative (compressive) effective stress and the volumetric part of stiffness when the material is subjected to compression. The influence of damage in the IDA model is governed by dilatancy θ, so two cases are considered.
If θ > 0, then Equation (36) is expressed as: where the power P regulates damage of a volumetric-deviatoric split. In fact, P > 1.0 accelerates the damage progress for the bulk modulus. If P = 1.0, then pure scalar damage is retrieved. It is possible P < 1.0, but, in that case, the development of volumetric degradation is slowed down. The rate of stress is: If θ < 0, Equation (36) is redefined: where R ∈ [0.0, 1.0] is a damage reduction factor. The closer to 0 the factor is, the less degradation of the bulk stiffness is involved. It should be noted that, in [50], a similar relation is introduced with P = 1.0 and R = 0.0, i.e., the bulk stiffness remains elastic for the negative dilatancy. Accordingly, for linearization, the rate of stress should be written as: The IDA model can also be implemented as nonlocal to prevent spuriously sensitive mesh discretization, cf. [49,51].

Testing of Considered Models
The overall data for each concrete model are as in [3]. The elastic constants are as follows: Young's modulus E = 32000 MPa and Poisson's ratio ν = 0.20. The next parameters of the particular model are tuned to initial uniaxial tensile strength f t = 3 MPa, maximum uniaxial compressive strength f c = 38.3 MPa, and tensile fracture energy G f = 0.11 N/mm. Small strains and plane stress conditions are assumed. The finite element size is equal to 100 mm. As mentioned in the Introduction, it is expected that the response of each model can vary despite the same set of global data. Furthermore, crucial features of the models are exposed in the presentation of results.

Concrete Damaged Plasticity (CDP)
In the CDP model, stress-strain relations have to be defined separately for tensile cracking and compressive crushing. For softening in the tension regime, starting from the point where a so-called cracking strain cr t = 0.0 (i.e., difference of total and elastic strains) corresponds to the tensile strength f t = 3.0 MPa, the curvilinear relationship between the uniaxial tensile stress σ t and the cracking strain cr t is determined as shown in Figure 7a. When uniaxial compression is considered, the stress-inelastic strain function similar to a parabola, cf. [52,53], is employed. Maximum compressive strength f c = 38.3 MPa is adopted for inelastic strain in c = 0.00172. This diagram is depicted in Figure 7b.   Figure 2 in 2D effective principal stress space. The ratio of equivalent stressq on the tensile and compressive meridians is K c and the default value 2 3 is used in the computations, see also [4]. For the CDP model without activated damage, the influence of dilatancy angle ψ is studied; hence, it can range from 5 • to 55 • . The plastic potential functions G p CDP for different angles are compared in Figure 3 with the same eccentricity E = 0.1. This default value is assumed in the numerical analysis.
In the case of activation of damage growth functions for tension and compression as shown in Figure 8, the test is run only for dilatancy ψ = 25 • . Moreover, damage ω c is insignificant in Willam's test, where compression is not present.   First of all, the influence of different dilatancy angles ψ in the plasticity model is verified. In the computations, the same values of ψ = 5, 25, 35 and 55 • as for the simulation of punching shear presented by [22] are introduced. It is clearly visible that, for each option, the values of all stress components are at most equal to the uniaxial tensile strength f t = 3 MPa, so the first condition of Willam's test is satisfied. It is also noticed that the stresses decrease in the second phase of the test. For ψ = 5 • and ψ = 25 • (see Figures 9a,b and 10), they approach zero at the final stage; hence, the second condition of Willam's test is also fulfilled. However, when the dilatancy angle ψ in the CDP model becomes larger than 35 • , negative values of horizontal stress σ 11 and minimum principal stress σ 2 are observed. It is evidently demonstrated in Figure 9d for ψ = 55 • . Therefore, for the CDP model, a similar observation as in [22] can be made, i.e., the dilatancy angle ψ should not be larger than about 35 • . Larger values produce a dilatancy effect (manifested by negative stresses), which could be justifiable in confinement conditions as e.g., in [54], but for Willam's test is simply non-physical and violates Condition 2. Additionally, the results illustrated in Figure 10 are prepared for the case where also damage functions are incorporated in the CDP model. The diagrams are similar to the one presented in Figure 9b. The difference is noticed only for σ 22 and σ 2 , where activated damage reduces their maximum values. Indeed, the presence of damage in this model for Willam's test is hardly relevant. It should be recalled that damage in the CDP model is basically introduced to simulate a stiffness degradation during unloading. Moreover, it occurs (results not included in the paper) that only damage evolution for tension, cf. Figure 8, is able to modify the solution.

Concrete Smeared Cracking (CSC)
The material relations in the CSC model are determined identically as in the CDP model. Figure 7a shows the softening branch of tensile stress σ t as the function of cracking strain cr t , while Figure 7b depicts the relationship of compressive stress σ t and inelastic strain in c . The associated strengths are also the same, but now two initial yield surfaces F p CSC,t for tension and F p CSC,c for compression are combined to obtain the yield criterion, see Figure 5.
ABAQUS users can decide about the shape of both surfaces as well as the values of parameters T , W, and S related to them via so-called failure ratios [4]: • the ratio of the biaxial compressive stress to the uniaxial compressive stress equals 1. 16  It should be emphasized that the associated flow is assumed in the CSC model. If no more options are stated, then full shear retention is given as the default. In the case of detailed specification of the shear retention option, the following parameters are defined: fraction ρ cl of shear modulus G for closed cracks in concrete (the default value is 1.0) and the maximum value of the strain nn normal to the crack plane is max -cf. Figure 4 (the default value is a large number).  The shear retention effect is analyzed in the computations for the CSC model. Analogically to the previous subsection, all stress components are depicted depending on strain 11 . Figure 11a shows results for ρ cl = 1.0 and max = 10.0, which corresponds to full shear retention. This value of max for 11 = 0.002 gives ρ = 0.9998 ≈ 1.0. It is visible that not only tensile strength f t = 3 MPa is exceeded many times by the maximum principal stress σ 1 and shear stress σ 12 , but they do not converge to zero stress. The minimum principal stress σ 2 from about 11 = 0.0002 becomes negative and further makes a mirror image of σ 1 relative to the horizontal axis for zero stress. For the case with ρ cl = 1.0 and much smaller max = 0.001, the effects of shear retention are partly active, see Figure 11b. Of course, formation of the first crack is correlated with the first peak when σ 1 = f t = 3 MPa, but, after that, the maximum value of σ 1 is larger than 4 MPa at the moment of formation of the second crack. Next, a sudden drop is noticeable in the diagrams of stress components. Both the increase of values of σ 1 and σ 12 as well as the decrease of σ 2 (negative values occur again) are connected with shear retention. For 11 ≈ 0.0005, this effect disappears, but the principal stresses rather deviate than go to zero stress. The results for the opposite case with small ρ cl = 0.1 and large max = 10.0 are presented in Figure 11c. Again, the formation of two cracks is observed, but further shear retention makes that σ 1 and σ 12 go to plus infinity, while σ 2 to minus infinity. The diagram of minimum principal stress σ 2 looks like a mirror image (with the negative sign) of the diagram of maximum principal stress σ 1 . Figure 11d illustrates the diagrams of stress components for the option with small ρ cl = 0.1 and also small max = 0.001. As previously, the crack formation is noticed and moreover the value of f t = 3 MPa is not exceeded, the presence of shear retention is almost imperceptible, but finally principal stresses σ 1 and σ 2 depart from zero value. Figures 11a and 12a indicate that the parameter ρ cl plays an important role in the tension regime, and its small value can lead to a substantial reduction of stress σ 12 . This observation is at odds with the CSC model description given in [4], according to which the parameter ρ cl should have no influence on the behavior of an open crack as it defines only the shear retention factor for a closed crack (in compression). It can be concluded that the actual implementation of the CSC model in the ABAQUS package departs at this point from its description [4]. The shear retention effect is suppressed for small ρ cl = 0.1. Summarizing, the CSC model fails Willam's test, although the results for the last considered case pass the first condition of the test. In fact, the behavior of the CSC model is very similar to the behavior of the standard fixed smeared crack model as observed by [1,27,29]. For both these models, the source of non-physical behavior seen in Figures 11a and 12a is quite obvious and has already been identified in [1]-retaining the full stiffness for shear in Phase II leads to the unbounded increase of shear stress for increasing shear strain.

Damage-Plasticity (DAP)
In damage theory, the threshold is calculated as quotient of the uniaxial tensile strength and Young's modulus, so κ o = f t /E = 0.00009375. In the test, exponential softening given in Equation (24) is taken into account. Parameters α and η are determined based on fracture energy G f = 0.11 N/mm. The residual stresses should asymptotically go to zero, so the first parameter α = 1.0 and a complete loss of stiffness is accepted. The ductility parameter η is estimated as equal to 4000, which seems to be unrealistically huge for the DAP model, but it is connected with G f and related to the element size, so this value is truly correct. If modified von Mises definition in Equation (23) is employed, the ratio between compressive and tensile strength is equal to k = f c / f t = 12.7667.
Plasticity with the HMH criterion for F p DAP is selected if the coupled model is turned on. The yield stress σ y is the uniaxial tensile strength, and f t = 3 MPa and isotropic linear hardening are applied. The hardening modulus h = 0.5 E is adopted. This value seems to be large, but it is known from [35] that, in the DAP model, hardening effects are connected with the fictitious (effective) configuration, i.e., with the material skeleton. The plastic part of the DAP model coupled with damage can also influence the development of microcracks. As shown by [55], the solution approaches the response as for pure damage when the value of h → ∞. Two ways of coupling can be considered: total or elastic. For h = 0.5 E, the differences in results for manners of the coupling are clear enough.  As previously, the stress components versus strain 11 are analyzed in the diagrams. The results for pure scalar damage model and the two different definitions of the equivalent strain measure are firstly compared. Figure 12a depicts the diagrams of stress components for the Mazars definition given in Equation (22), while Figure 12b illustrates the solution obtained for the modified von Mises measure defined in Equation (23). A more rapid decrease of stresses is noted for the second option. It is also observed that the maximum values of stresses σ 22 and σ 12 are about 50% smaller than for the Mazars definition. For both options, the uniaxial tensile strength f t = 3 MPa is kept and all stress components tend to zero, hence it can be concluded that, for pure scalar damage, Willam's test is passed.
The diagrams for two manners of coupling in the DAP model are presented in Figure 13. Now, only the modified von Mises definition is employed. The results have the same tendency as for pure scalar damage. The tensile strength is not exceeded and stresses tend to zero in the second phase of the test. However, after the peak, the descending paths run in such way that, for option ( ) shown in Figure 13a, they are below those resulting from pure damage, cf. Figure 12b. This is the case of total coupling in the DAP model. The reverse is the case when the weak coupling˜ ( e ) is considered, see Figure 13b. Now, all stress paths are above those presented for pure damage. It can generally be indicated that the DAP model passes Willam's test in each case.

Isotropic Damage (IDA)
The data for the IDA model are the same as for pure damage in the previous section, but additional parameters P and R have to be determined. They decide how the damage parameter ω degrades bulk modulus K and shear modulus G in different ways. In Willam's test, dilatancy only increases and θ ≥ 0, so the constitutive relation given in Equation (37) is employed or, in other words, the damage reduction factor R = 1.0, cf. also Equation (39). Other values of R are impossible in this test. The results for different R are shown by [51], where θ < 0 is admitted e.g., for a splitting test. In this analysis, the value of the power P is introduced as smaller or larger than 1.0 to demonstrate deviations of the response of the IDA model. The power P = 1.0 corresponds to the DAP model with pure damage.
Again, in this subsection, axial strain-stress component's diagrams are shown as previously. Figure 14a presents the results for P = 0.1. The volumetric degradation is significantly reduced, so, after the first peak for σ 11 and σ 1 = f t = 3.0 MPa and their quick and slight decrease, a second increase of stresses occurs up to value 7.69 MPa for 11 = 0.0008. After passing this point, all diagrams descend, probably to zero. All components, apart from shear stress σ 12 which is zeroed, run together. An analogical behavior is observed for the next case depicted in Figure 14b. However, it is found that, for P = 0.225, a second hump is reached when σ 1 equals 3.0 MPa, since the uniaxial tensile strength is not exceeded and the first condition of Willam's test is passed. The second condition is also fulfilled, because the stresses approach zero. The results for cases P = 0.5 and P = 4.0 illustrated in Figure 14c,d satisfy Willam's test as well. Increasing power P accelerates the process of stiffness degradation due to larger and larger reduction of bulk modulus K. In the case P = 4.0, the steepest slope of the softening path for σ 11 and σ 1 is noticed. Negative values of σ 22 and σ 2 are manifested, but finally they return to zero. It means that the second condition is always passed, even if quite large values of P are introduced. It can therefore be concluded that the IDA model passes Willam's test only if P ≥ 0.225.

Discussion-Comparison of Models for Principal Stresses and Their Directions
In this section, selected eight cases are compared in the diagrams prepared to demonstrate a change of maximum and minimum principal stresses as well as of principal stress directions. These cases have previously been presented; they are characteristic options for all the models discussed in the paper. The comparison is done in order to show directly the different behavior of the selected concrete models for the same general data. The list of selected cases together with their most important features is given in Table 1. In the first column, corresponding acronyms of all cases investigated below are written. The second column includes full names of the models and their details. For additional help, in the last column, the number of figure related to the considered option is noted.

Acronym
Model Crucial details Figure   CDP25 concrete damaged plasticity dilatancy angle ψ = 25 • Figure 9b CDP25dam concrete damaged plasticity dilatancy angle ψ = 25 • , tensile damage-Figure 8a       Table 1. Figure 15b illustrates analogical results for the minimum principal stress σ 2 . The diagrams for both options connected with the CDP model (CDP25-dark magenta dashed-dotted line and CDP25dam-dark green solid line) almost overlap for σ 1 and σ 2 , but, for σ 2 , the maximum value of stress is smaller for CDP25dam than for CDP25. It means that, for Willam's test, where an active process is only considered without any type of unloading, the damage component in the CDP model does not matter. The most deviating diagrams are obtained for the CSC model, see the green dotted line for CSCfull and brown dotted line for CSC in Figure 15. Moreover, when σ 2 is taken into account, then negative values can appear. The results for the DAP model prove that coupling of the damage model with plasticity influences the response. The blue dashed curve (DAPtot) for full coupling by˜ ( ) in the DAP model is below the black solid curve (DAP) for pure damage, while the red dashed one (DAPela) for weak coupling bỹ ( e ) is above the black one. A more ductile response is visible when only the elastic part of strains e stimulates the damage growth and, conversely, if the total strain tensor influences in the damage process, then a more brittle response is noticed. However, apart from the results for the CSC model which are unacceptable, the diagram for the IDA model (gray solid line) gives the most ductile solution. It is observed for both principal stresses σ 1 and σ 2 , cf. Figure 15a Figure 16 illustrates the evolution of principal stresses in first and fourth quadrants of the principal stress plane, cf. Figures 2 and 5. It is seen that, in phase I, which corresponds to the uniaxial tension, σ 1 grows from zero to 3.0 MPa for all cases, while σ 2 is equal to zero in that stage. In phase II, when softening (cracking) occurs, the values of σ 1 decrease. At the same time, the values of σ 2 initially increase, but finally tend to zero. This is observed for all cases except CSCfull and CSC. The values of the principal stress σ 1 for the CSC model exceed the uniaxial strength f t . After that, it seems that the stress paths tend to zero, but for about 2.0 MPa, this process is broken and they shoot up to infinity. For the case CSC, a second return is noticed, but near the origin this curve turns again and finally goes to infinity. The calculations for this case are interrupted. Hence, it is verified once more that the CSC model fails Willam's test. The maximum values of σ 2 are obtained for IDA and next for CDP25. It can also be noticed that the differences between CDP25 and CDP25dam as well as DAP, DAPtot, and DAPela are small. The CDP, DAP, and IDA models behave in the principal stress space in a similar way. The directions of principal stresses are not constant during the loading process, see also e.g., [2]. After the peak, at the beginning of phase II in Willam's test, the rotation of principal directions increases fast, then this change slows down and finally tends to the angle with value 52.018 • . The rotation of principal directions can be expressed by evolution of angle Θ for strains and Θ σ for stresses during the loading process. In Figure 17 for two options DAP and IDA, where the scalar or isotropic damage is employed, the angle of principal directions Θ σ for stresses and Θ for strains evolves in the same manner, since it is seen that, for pure damage (without any coupling), principal strains and stresses are coaxial. For the coupled version of the DAP model, no matter whether by total strains -case DAPtot or by elastic strains e -case DAPela, the change of angle Θ σ is faster than the change of Θ . Hence, the coaxiality of principal directions between strain and stress fields can be lost. It should also be noticed for options DAPtot and DAPela that both diagrams overlap and approach the final value 52.018 • for the angle of principal directions. The existence of plasticity accelerates the effect of rotation of the principal stresses. It is confirmed for the CDP model as well. For option CDP25, the principal directions for stresses grow very fast and reach the angle with limit value 52.018 • for 11 = 0.0002. When the damage component is added in the CDP model, then, for option CDP25dam, this value of angle is achieved for 11 = 0.0003. It seems that the value 52.018 • is attained immediately for the CSC model, but next, for both cases CSCfull and CSC, the angle drops to about 35.0 • , which corresponds to manifestation of the presence of the primary crack and then is slowly reduced till 11 ≈ 0.00022. After that, the change of the angle depends on whether stronger or weaker shear retention is assumed in the CSC model. For option CSCfull, the value of the angle slowly increases up to about 45.0 • . When option CSC is considered, the angle decreases almost to zero and afterwards increases to 40.0 • for 11 = 0.002, see the internal subfigure in Figure 17.

Conclusions
In the paper, popular concrete models are selected to verify if they pass or fail Willam's test [1]. This test is a one finite element benchmark containing two phases: uniaxial tension till the tensile strength is achieved and softening in the biaxial tension-shear regime. It is proved that the results of the test for each concrete model can be different even if the parameters are calibrated in such a manner that the models exhibit almost identical behavior in uniaxial tension.
The following models are tested: concrete damaged plasticity (CDP) and concrete smeared cracking (CSC) models delivered in the ABAQUS software [4], the da-mage-plasticity (DAP) model without or with coupling of both theories by elastic or total strain tensor and finally an isotropic upgrade (IDA) of the damage model where the volumetric-deviatoric split is applied. The DAP and IDA models are implemented in the FEAP package [5]. Table 2 summarizes the content of the paper. The CDP model passes Willam's test, but only when the dilatancy angle ψ is smaller than 35 • ; otherwise, exaggerated dilatancy is observed. Please note that a large dilatancy angle, i.e., ψ ≥ 49 • , is considered in several works, cf. [21,54,56,57]. The recommendation to use the model for ψ ≤∼ 35 • is similar to that given in [22], where punching shear in slabs is simulated. The CSC model fails Willam's test, even if the effect of shear retention is substantially suppressed. On the other hand, the DAP model passes the test independently of the presence and kind of coupling with plasticity. When the scalar damage is upgraded to the isotropic version as in [48] or in the fashion of the IDA model, then the parameters can decide about passing or failing Willam's test. Here, for the IDA model, the power P governs the degradation of the volumetric part of the stiffness, but it should be equal to or larger than 0.225 to satisfy the first condition and thus pass Willam's test. As stated in the Introduction, failing Willam's test by a given concrete model raises serious doubts concerning its ability to describe properly the structural behavior with predominant mixed-mode fracture, e.g., RC beams failing in shear as shown in [12]. Therefore, among the investigated models, the CSC model cannot be recommended for such structural analyses. The other models, i.e., CDP, DAP, and IDA, seem to be well suited for nonlinear FE computations of concrete structures with the predominant mixed-mode fracture; however, the restrictions mentioned above for the CDP and IDA models should be taken into account. As a general suggestion-in the authors' opinion, a verification using Willam's test should be mandatory when any new material model for concrete is proposed.