There are two main types of titanium alloys used for marine structures:TC4 and TA31 (see Table 1 for the main alloy composition). TC4 titanium alloy is a dual-phase alloy, which is widely used in the field of aerospace, and a lot of related research has been done (Luo et al., 2008; Allahverdizadeh et al., 2015). Compared with TC4 titanium alloy, TA31 titanium alloy has the advantages of high plasticity, high fracture toughness and good weldability. However, there are few studies on the constitutive relation model of TA31 titanium alloy. In this paper, Experimental and numerical investigation of TA31 titanium alloy will be committed. Figure 1 shows the microstructure of the tested TA31 titanium alloy, which consisted of equiaxed primary α phase in a transformed β matrix. The α+β phases are clearly visible and the average grain size is approximately 50 μm. The characteristic of the bi-modal structure is that the equiaxed primary α phase with a content of less than 50% is distributed on the matrix of the β transformation structure. The bi-modal structure is a combination of the basket-weave structure and the equiaxed structure, which has the characteristics of both types of structure and has better comprehensive mechanical properties.

The aim of this experiment is to study the plastic hardening behavior and failure mechanism of TA31 alloy. Four types of specimens are designed to represent different stress states, as shown in Figure 2, including smooth round bar (SRB) tensile test, notched round bar tensile test with different notch radius R (NRB) and shear-induced specimens test (Shear).

In the test, an electronic extensometer with a gauge distance of 50 mm is used to obtain the load-displacement curve until the fracture, and the test load is controlled by displacement. In order to capture the fracture position within the gauge length, a slow loading rate (2 mm/min) is adopted.

The measured yield strength and ultimate strength of TA31 alloy are 813MPa and 881 MPa respectively, and Young’s Modulus is 118 GPa There is no obvious yield platform for TA31 alloy, and notch size has great influence on bearing capacity and elongation of samples. Figure 3 shows the different macro-fracture phenomena, such as SRB, NRB with different notch radius and shear-induced specimens. The size of oblique shear lip decreases with the decrease of notch radius, which indicates that the influence of the shear fracture mechanism is gradually weakened. On the contrary, very small and almost zero dimples appear on the smooth fracture surface of the shear-induced specimens, and the necking phenomenon is not observed.

The plasticity model and damage model will be formulated in terms of invariants of stress tensor, such as the equivalent stress σ ¯ , the stress triaxiality η and the Lode parameter L. σ m = σ 1 + σ 2 + σ 3 / 3 (1) σ ¯ = σ 1 − σ 2 2 + σ 2 − σ 3 2 + σ 3 − σ 1 2 / 2 (2) η = σ m / σ ¯ (3) L = 2 σ 2 − σ 1 − σ 3 / σ 1 − σ 3 (4)

Where σ ₁, σ ₂ and σ ₃ denote principal stresses. The plane stress state is dominant in the most of pressure-hulls for deep-sea equipment. It was shown by Xue (Xue and Wierzbicki, 2008) and Wierzbicki that the condition σ ₃ = 0 relates the parameters η and L as shown in Eq. 5. η = 3 + L / 3 L 2 + 3 (5)

Eq 5 calculates the corresponding value for the stress triaxiality under plane stress with a given Lode parameter between [-1, 1]. In this paper, the fracture loci under plane stress will be depicted in the space of (η, σ ¯ ).

The Johnson-Cook plasticity model (Johnson and Cook, 1985) is an empirical model for metals subjected to large strains, high strain rates and high temperatures. σ = A + B ⋅ ε ¯ p n ⋅ 1 + C ⁡ ln ε ˙ / ε ˙ 0 ⋅ 1 + T m (6)

Where ε ¯ p is the equivalent plastic strain, σ is the Von Mises tensile flow stress, ε ˙ is the effective current stain rate, ε 0 ˙ is reference stain rate, and T is temperature. These five material parameters are A, B, C, n and m. This paper mainly studies the failure behavior at room temperature under quasi-static condition. Therefore, the J-C model could be simplified as Eq. 7. σ = A + B ⋅ ε ¯ p n (7)

Based on the J-C model, Zhang (Zhang et al., 2022c) developed a new plasticity model of TA31 titanium, which was short for Modified Johnson-Cook (MJC) model. In order to consider the pressure sensitivity effect, a parameter accounting for stress triaxiality was introduced into the strain hardening modulus, B, and hardening exponent, n, as described in Eq. 8. σ = A + B ⋅ ε p n B = b 1 ⁡ exp − η ¯ a v g / b 2 + b 3 n = n 1 ⁡ exp − η ¯ a v g / n 2 + n 3 (8) Where b ₁∼b ₃, n ₁∼n ₃ are parameters to be determined. Table 2 shows the obtained parameters for MJC model of TA31 titanium alloy after the calibration.

Damage is considered as a thermodynamic state variable, which characterizes the deterioration of the material. Assuming that the distribution of damage is isotropic, the value of damage in all directions can be shown by the scalar factor D. The progressive damage of a ductile material is based on the change of the shape and the number of voids within the material, which can also be represented by the change of the Young’s modulus of material, as shown in Eq. 9. According to the strain equivalence hypothesis (Lemaitre, 1985), the constitutive equations of a damaged material are the same of the virgin material with no damage where the stress is simply replaced by the effective stress, and the value of the strain in the damaged material is equal to the value of the strain in the undamaged material with a stress value, as shown in Eq. 10. D = 1 − E eff E 0 (9) σ eff = σ 1 − D (10)

Where σ _eff is the effective stress, σ is the stress value in the damaged material, and E ₀ and E _eff are the Young’s modulus of the undamaged and damaged material, respectively. In the absence of any coupling between the damage and plasticity dissipation potentials, the total dissipation potential is given by linear superposition of the dissipation potential with plastic deformation, f _p, and the damage dissipation potential, F _D, for the Bonora damage model (Bonora, 1997). The damage dissipation is assumed to depend on the deformation history, where the damage, D, variable is coupled with plastic strain, as formulated in Eq. 12. F = f p + F D (11) F D = 1 2 − Y S 0 2 S 0 1 − D ⋅ D cr − D α − 1 α p n + 2 n (12)

Where S ₀ is a material constant, α is damage exponent that determines the shape of the damage evolution curve, n is the hardening exponent of plasticity model in a power-law form, D _cr is the value of damage at failure, and p is equivalent plastic strain. The damage associated variable Y, which is called damage strain energy release rate, is defined as Eq. 13. Y = 1 2 E σ 1 − D 2 ⋅ 2 3 1 + μ + 3 1 − 2 μ η 2 (13) Where μ is Poisson’s ratio, and η is stress triaxiality. By normality rule, the following expression for the kinetic law of damage evolution can be obtained in Eq. 14. D ˙ = − p ˙ 1 − D ∂ F D ∂ Y (14)

Under the assumption of proportional loading condition, the damage evolution is given according to the following expression. D = D cr 1 − 1 − ln p / p th ln ε ¯ cr / ε ¯ th ⋅ 2 3 1 + μ + 3 1 − 2 μ η 2 α (15)

Where pth is the threshold strain under multiaxial stress at which damage process are initiated, and ε ¯ t h is the threshold value of uniaxial equivalent plastic strain for damage process. Thomson and Hancock experimentally observed that the damage threshold strain is scarcely sensitive to stress triaxiality, so the effect of stress triaxiality on the threshold strain is neglected in this paper. ε ¯ c r is the critical value of uniaxial equivalent plastic strain to failure, and this value does not coincide with the failure strain of a standard tensile test, since ε ¯ c r is the maximum value of the specimen, on the contrary the failure strain of the test is the mean value. D _cr is the critical damage value at which complete failure occurs, which is directly related to the critical value of equivalent plastic strain that suddenly reduces the load carrying capability of the effective resisting section. In order to simplify the calculation, D _cr is set as one in this paper.

In this paper, The FEA models are established in the ABAQUS/explicit platform to simulate the foregoing quasi-static tests. Damage behavior and failure stain could be predicted by FEA simulation using an appropriate developed subroutine VUMAT. Instead of separating the elements, an element is removed if the critical damage value of the material point is reached. Figure 4 shows the FEA models used in this paper. Three dimensional hexahedral elements with reduced integration (C3D8R) are employed, and the minimum element size is kept to be 0.1 mm in the minimum section. Adaptive meshing is employed to avoid element distortion in the large plastic deformation. Taking the advantage of symmetry to reduce computing time, one-eighth of the model is adopted for round bar specimens, and a half of the model is used for shear specimens making a half of the thickness.

Damage threshold strain indicates the strain level at which damage process starts to take place, and voids nucleation are large enough to affect the material stiffness at this stage. In order to simplify, the damage threshold strain is equal to the uniaxial plastic threshold strain ε ¯ t h , and the influence of stress triaxiality on the damage threshold strain is ignored. According to the definition of damage, D _cr is set to 1, indicating that it becomes suddenly incapable of sustaining the external loading consequent failure of the material at the microscopic scale. In this paper, there are two parameters ( ε ¯ c r and α) in the Bonora model needed to be calibrated. A new approach for numerical determination of the Bonora’s damage parameter is used in this paper, which is proposed by Shamshiri et al. (2021) to determine the Lemaitre’s damage parameter. Based on the initial guessed value of the parameters, The foregoing tests is numerically simulated through the MJC model and the Bonora damage model. After the finite element simulations, the force-displacement curve and true stress-strain diagrams are numerically obtained and compared with the tests data. If the numerical and experimental diagrams are in good accordance and show adequate correlation, the trial values are accepted as the true for the Bonora’s damage parameters ( ε ¯ c r and α). Else, a new guess value is tried until the comparison of numerical and practical results exhibits sufficient adaptation. After reaching a proper agreement between the numerical and tests data, the last trial value is finalized for the Bonora’s damage parameter.

The mechanical responses and failure mechanisms of specimens in different stress states are quite different as shown in Figure 5, which shows the load-displacement curves predicted by FEA simulation and those obtained in the tests. A similar tendency is observed between the experimental and numerical results for the tensile tests with the round bar specimens, and the critical points are located in the center region of the minimum section, and the maximum load occurs at the beginning of necking. Whereas the maximum load of shear-induced specimens occurs until complete fracture, and the critical point is located at the outer surface, as shown in Figure 5D.

Owing to the absence of the softening behavior of material, the MJC model can only accurately predict the load-displacement relationship before necking from this figure. It can be seen that the numerical results based on the Bonora model are in good agreement with the results of tests. However, the damage parameters ( ε ¯ c r and α) are sensitive to the stress triaxiality, and the calibrated parameters of Bonora damage model for different specimens are reported in Table 3.

Failure strain is the strain at the end of plastic deformation, or right before final crack at the point where fracture initiates. Measure of failure strain is a big challenge in the tests. Because fracture initiates inside of metals, it is impossible to measure the moment of fracture. The failure strain obtained from the foregoing tests is the mean value of the specimen, which is lower than the strain in the interior of necking section as shown in Figure 6.

Currently, a hybrid experimental-numerical approach has been widely employed to obtain failure strain. In this method, tests are committed at first and load-displacement curves are recorded. Then these tests are simulated based on FEA model and the failure strains are extracted from simulations at the fracture onset. In this paper, this hybrid experimental-numerical method is adopted to measure failure strains in smooth round bar specimens, notched round bar specimens and shear-induced specimens for TA31 alloy. As shown in Figure 5 above, it is observed that the predicted fracture onset based on the Bonora damage model agrees well with tests for all of the specimens. Simultaneously, the evolution of stress triaxiality with plastic strain of fracture onset can be obtained from numerical analysis of the forgoing tests as shown in Figure 7. It can be found that as the stress triaxiality increases, the failure strain ε ¯ c r decreases for the round bar specimens. On the contrary, for the shear-induced specimen the stress triaxiality of the element at the critical point remain minor, and the failure strain is lower than the strain at the critical point of smooth round bar specimen. The major reason for this difference is known to be the diverse stress states in the tensile tests with different mechanism.

The stress state plays a significant role in determining the parameter of Bonora damage model. As shown in Figure 8, both ε ¯ c r and α decrease as η increases for the round bar specimens, but for the shear-induced specimens, there is a different trend.

The changes in stress triaxiality and equivalent plastic strain are non-proportional with increasing load, and the distributions are predicted by FEA simulation with Bonora damage model in Figure 9 and Figure 10. Figure 9A shows that the ε _p-value and the η-value remain 0.0508 and 0.333 respectively across the minimum section before damage initiates for the smooth round bar specimens, whereas the value of the notched round bar specimens exhibit a larger variation across the minimum section. Figure 9B shows that the ε _p-value and the η-value decrease from the center point to the outer surface along the minimum section, which indicates that the crack is predicted to start from the center of the section and to propagate towards the outer surface for the round bar specimens.

For the shear-induced specimens, the η-value varied greatly around zero, which indicates that pure shear deformation does not exist in practice. However, As shown in Figure 10, the η-value remain minor across the outer surface, which is approximate to the pure shear state. The critical point is located at the shear slip band which is at the outer surface close to the edge of the shear gauge.

Scanning electron microscope (SEM) images of the fracture surface have been acquired to analyze the fracture mechanism. In the round bar specimens, many dimples of varying size are present on the fracture surface as shown in Figures 11(A–C). It can be found that the size of the dimple varies with the notch radius. The SRB specimens with the smallest triaxiality produces the largest dimple, whereas the NR5 specimen with the highest triaxiality produces the smallest dimple. However, in the shear-induced specimens, as shown in Figure 11D, very minor microvoids appear on the smooth fracture surface. Qi (Qi et al., 2017) pointed out that the growth of cavity is the comprehensive effect of shear stress and tension stress. Hence, higher stress triaxiality, where tension stress plays a dominant role, will promote the isotropic enlargement of microvoids, which will lead to the constriction of the ligament between two voids before shear aggregation. In addition, it seems that the high hydrostatic stress, on behalf of high stress triaxiality, will enhance the expansion of spherical cavity. In the shear-induced specimens, the stress triaxiality is close to zero, and the shear stress dominates, resulting in the elongation of the voids in the fracture along the shear fracture direction. The above research shows that the deformation mechanism of micro voids of TA31 alloy is greatly influenced by the stress state, which is also the result of competition between shear mechanism and tension mechanism in the plastic deformation process.

The evolution of the internal microvoids is varying with stress triaxiality, which is also illustrated briefly in Figure 12. The SRB specimens are under uniaxial load with η ¯ a v g ≈1/3. Firstly, the preexisting voids inside evolve and the new ones nucleate in one direction. Subsequently, the localized plastic flow between voids will lead to the formation of cracks through internal necking, as shown in Figure 12A. The NRB specimens exhibit three-dimensional stress state with higher stress triaxiality, which will accelerate voids nucleation and growth, leading to the increase of the number of voids and the interaction between the mechanical fields around primary voids, as shown in Figure 12B. In this case, the macroscopic displacement along the uniaxial loading direction of the NRB specimens are smaller than the SR specimen, indicating that a higher η-value will restrain the size of dimples. However, Figure 12C shows that the deformation characteristics of shear specimen, where shear mechanism dominates, is different from the tensile tests. Here, coalescence occurs by extensive localized shear deformation in the ligaments between voids. It is shown that a higher η representing triaxial stress state will accelerate void initiation and inhibit coalescence between voids at the beginning of damage and that the material becomes softer and suddenly cracks along the weakest direction, which will lead to a lower ε ¯ c r and α as shown in Figure 8. On the contrary, the shear mechanism that η-value is nearly zero will accelerate coalescence between voids, resulting in a different trend for the parameters of ε ¯ c r and α shown in Figure 8. When η is in the range of [0, 1/3], both shear mechanism and tension mechanism coexist and compete with each other.