Bingnan Qian

Jinyong Zhang

Yangyang Fu

Fan Sun
email: fan.sun@chimieparistech.psl.eu

Yuan Wu

Jun Cheng

Philippe Vermaut

Frédéric Prima

In-situ microstructural investigations of the TRIP-to-TWIP evolution in Ti-Mo-Zr alloys as a function of Zr concentration

Keywords: Metastable beta (β) Ti-alloys, TRIP/TWIP effects, Deformation mechanisms, in-situ traction-EBSD, TEM observation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Titanium and its alloys have attracted much attention due to their high strength [1,2] , low density [3] , excellent properties of biocompatibility, and good corrosion resistance [1,4,5] . However, compared to traditional structural alloys, such as stainless steel and Co-Cr alloys, their low ductility and the lack of strain-hardening limited their applications where combination of strength and ductility are required. Recently, combination of high strength, large ductility, and excellent strain-hardening rate has been achieved in a new family of metastable β titanium alloys, such as Ti-12Mo [6] , Ti-9Mo-6W [7] , Ti-6Cr-4Mo-2Al-2Sn-1Zr [8] , Ti-12Mo-5Zr [9] , Ti-11Mo-5Sn-5Nb [10] and Ti-18Zr-13Mo [11] , with TWIP and/or TRIP effects, leading to greater potential applications.

Recent studies have reported complex deformation mechanisms mainly involving deformation twinning ({332}<113> mode and {112}<111> mode) and stress-induced martensite (SIM α″), in addition to dislocation slip. As a result of these simultaneously activated mechanisms, a high strain-hardening rate may be reached by a microstructural evolution during plastic flow. This microstructural evolution has been shown to give birth to: (1) the so-called dynamic Hall-Petch effect [12,13] ; (2) a strong mechanical contrast effect between the twins and their matrix due to high elastic anisotropy of BCC β phase in titanium alloys [14] ; (3) a dense multi-modal twinning network [11] ; (4) a hierarchical twin/SIM α″ microstructure [6,7,13,15,16] .

As well-known, the TRIP and TWIP effects and their influence on mechanical properties are highly dependent on the stability of β phase in Ti-alloys. Martensitic transformation (β → α″) can be favorably activated by applying external stress or strain in the β-metastable Ti-alloys when the chemical stability of β phase (usually estimated by Mo-equivalent content, Moeq wt.%) is relatively low at a rough range of 10.0% ≤ alloy Moeq ≤ 13% [17,18] . Mechanical twinning can be mainly activated when the β stability becomes higher (Moeq roughly between 12% and 15%) [6,19,20] . For alloys displaying intermediate β stability (overlapping zone between the two ranges), the simultaneous occurrence of TRIP/TWIP effects has been observed [6,7] . Therefore, a continuous transition of deformation mechanisms, from near-TRIP to near-TWIP effects, can be expected by increasing progressively the chemical stability of a TRIP/TWIP alloy.

The Zr element has been found to be a neutral element (in Ti-Zr system) or weak β-stabilizing element (in Ti-Nb-Zr system) [21,22] . Previous studies have reported that the nature of Zr in a Ti-M-Zr system (M = β-stabilizer alloying element) strongly depends on the β-stability and concentration of M [9,21,22] . A theoretical equivalence factor of Zr in the Moeq estimation has been proposed to be +0.47 by Wang et al. [23] , showing that the β stabilizing effect of Zr is about half of that of Mo. Therefore, using Zr could be two-fold advantageous since a wide composition range (up to tens of weight percent) can be added into the master alloy: (1) to achieve better tolerance of alloys composition for TRIP/TWIP fine-tuning; (2) to obtain the solution strengthening effect. Furthermore, the exact effect of different Zr additions in the TRIP/TWIP Ti-Mo system has not yet been reported with systematical experimental results. In order to study the effect of Zr addition on the TRIP to TWIP transition in Ti-Mo system, a comparative investigation is performed between the TRIP/TWIP Ti-12Mo-xZr (x= 3, 6, 10) (wt. %) alloys. This work focuses on in-situ investigations of the deformation mechanisms in β phase with different stability levels. The investigation is performed at various length scales using the Ti-12Mo-xZr system. The deformation sequence, orientation relationships and the volume fraction of SIM α″ at different strain levels are traced by in-situ EBSD analysis during tensile loading and after unloading. The in-situ observations reveal the details of the formation, development and transformation of the complex martensite/twin structure.

Based on the results, the deformation mechanisms involved in the macroscopic mechanical behavior are discussed by considering the heterogeneous nature of the microstructural evolution.

Experimental

100g ingots of Ti-12Mo-xZr were fabricated by vacuum arc-melting using high purity Ti, Mo and Zr elements. These alloys were melted at least five times under a high purity argon atmosphere.

The ingots were annealed in vacuum tubular furnace under 2.0 × 10 -4 Pa at 1173K for 1.8 ks, followed by water-quenching. The quenched ingots were cold-rolled to 0.5mm sheets, corresponding to reduction rate of 95%. The cold-rolled sheets were solution-treated at 1173K under high vacuum (2.0 × 10 -4 Pa) for 1.8 ks and then water-quenched to retain the fully recrystallized β microstructure. Tensile tests were conducted by using an INSTRON 5982 machine at room temperature with an initial strain rate of 10 -3 s -1 . Phase characterization was carried out by X-ray diffraction (XRD) using Bruker D8 ADVANCE with CuKα radiation operating at 40 kV and 40 mA. Deformation microstructures were analyzed by using in-situ tensile tests coupled with electron backscatter diffraction (EBSD) scans, performed on a Proxima 100-Micromecha machine in a field emission gun scanning electron microscope (FEG-SEM ZEISS LEO1530) operating at 20 kV and a NORDIF EBSD camera. A JEOL 2100plus transmission electron microscope (TEM) operating at 200 kV was used to carry out fine-scale microstructure analysis. Samples for TEM observations were prepared by a twin-jet polishing technique, using a solution of 4% perchloric acid in methanol, held at about 250 K.

Results

Characterization of solution treated (ST) states

Fig. 1 shows the XRD profiles, optical micrographs and selected-area electron diffraction (SAED) patterns of the ST samples of Ti-12Mo-3Zr, Ti-12Mo-6Zr, and Ti-12Mo-10Zr alloys, respectively. As shown in Fig. 1(a), the presence of the β phase is confirmed in the three alloys without other additional diffraction peaks. It can be noticed that the peak position of the β phase (2θ: degree) slightly shifts towards lower angles with the increasing of Zr content, suggesting a small expansion of the β phase lattice parameter consistent with the bigger size of Zr atoms. Fig. 1(b) displays that the microstructure of ST samples consists in typical recrystallized β grains (with a grain size in the range 50-80 μm). Moreover, besides the β phase, the nanometer athermal ω (ωath) phase is actually detected by TEM technique, as displayed in Fig. 1(c). It can be observed that the reflections of ωath phase become weaker and more diffuse with the increased Zr content, when normalized to the β matrix reflection intensities, suggesting that the ωath precipitation could probably be lower in volume fraction. Pang et al. [24] has observed that Zr addition in Ti-24Nb-[0-8]Zr (at.%) alloys could suppress the formation of ω phase. In contrast to the results of Pang et al [24] , the reflections of ωath phase is still visible in Ti-12Mo-10Zr alloy.

Tensile properties

Fig. 2 shows the true strain -true stress curves of the present alloys at ST state. The stress level of the plastic flow curves clearly increases with Zr addition, highlighting the notable solid solution strengthening contribution of Zr element in the Ti-12Mo-xZr system. In particular, the yield stress of Ti-12Mo-10Zr alloy approaches about 740 MPa, which is notably high for a strain-transformable Ti-alloy. Likewise, the uniform elongation of all three alloys is higher than ε = 0.24, displaying large ductility. It is worth noting that the plateau-like platform can be observed at the very beginning of plastic deformation in Ti-12Mo-3Zr alloy (marked in Fig. 2). Sun et al. [6] suggested that the appearance of the plateau is related to the strain-induced martensitic transformation. Aiming at clarifying the evolution of strain-hardening behavior, the corresponding strain-hardening rate (dσ/dε) curves of the three alloys are shown in inset of Fig. 2. The three alloys exhibit different types of evolution of the strain-hardening rate. The strain-hardening rate curves are non-monotonic, similar to those observed in TRIP/TWIP β Ti-alloys [6][7][8][9]14] . It can be observed that Ti-12Mo-6Zr displays a noticeable difference with respect to the other two alloys, regarding the presence of a stable strain-hardening rate at about 1750 MPa in the range of ε=0.05-0.13. The underlying mechanisms are investigated in detail by in-situ EBSD and TEM techniques.

Deformed microstructure characterization

XRD profiles

Fig. 3 shows the XRD profile of Ti-12Mo-xZr (x=3, 6, 10) fractured samples. The residual SIM α″ phase can be detected in the Ti-12Mo-3Zr and Ti-12Mo-6Zr alloys. The appearance of residual α″ is related to the occurrence of SIM α″ in the plastic deformation stage in Ti-12Mo-3Zr and Ti-12Mo-6Zr alloys, partially preserved in the alloys after unloading. Such similar phenomenon is also observed in Ti-12Mo alloy [6,25] and Ti-9Mo-6W alloy [7] . In contrast, only β phase is detected in Ti-12Mo-10Zr alloy without SIM α″, suggesting either SIM α″ did not occur or has completely reversed to β phase during unloading [11] .

In-situ EBSD analysis

In order to clarify the activation and the evolution of the TRIP/TWIP microstructure of the present alloys, in-situ EBSD mapping under tension is performed in the early deformation stage. Fig. 4 shows the EBSD analysis of Ti-12Mo-3Zr sample taken from the same region at different strain. The inverse pole figure (IPF) map at ε=0.02 in Fig. 4(a) shows that plate-like deformation bands with different orientations (coded by color contrast) appear in the β grain, upon loading. The green plate-like deformation bands are identified to be SIM α″, as indicated in the corresponding α″ phase maps (Fig. 4(b)). For further clarifications, as shown in Figure 4c, the orientation relationship between green SIM α″ deformation bands and β matrix has been analyzed and found to be in agreement with the classic crystallographic orientation relationship {211}β //{110}α″ , suggesting that SIM α″ could likely originate directly from the β matrix. It can be seen from Fig. 4(b) that the quantity and thickness of SIM α″ bands (inside Grain A) increase upon loading from ε=0.02 to 0.04.

After unloading from ε=0.04, most of the SIM α″ bands are preserved in β matrix. Similar to the observations reported in superelastic Ti-27at.%Nb [26] and TRIP/TWIP Ti-V-Cr-Al [14] , internal twinning of SIM α″ as it is observed that the color contrast of some α″ bands become blue upon loading (ε=0.04 loading) and then transformed back to green during unload (ε=0.04 unloading). The α″ twinning mode is identified as {130} <310>α″ martensite twinning by crystallographic stereo projection analysis in Fig. 4(c). Fig. 5 shows in-situ EBSD mapping of Ti-12Mo-6Zr sample. Both {332} <113> twinning (hereafter called as 332T) and SIM α″ can be observed in the deformed sample. The IPF maps in Fig. 5(a) shows that plate-like deformation bands with two different color contrasts (purple and blue) appear in the β matrix. It can be seen from Fig. 5(c) that the first deformation band appears at loading to ε=0.02. A quick multiplication of the bands happens when increasing the tensile strain to ε=0.04. Most of the bands can be indexed as β phase respecting strictly the 332T relationship with respect to the β matrix, confirmed by the stereo projection analysis (Fig. 5(c)). Two variants of 332T can be observed in the grain traced by in-situ EBSD, belonging to two different 332 planes of the matrix (the two common poles circled in Fig. 5(c)). As shown in Fig. 5(a), the two variants of 332T keep unchanged in their crystal orientations and thicknesses after unloading (ε=0.04 unloading), which suggests that the mechanical twinning between the β matrix and 332T is irreversible. In addition, SIM α″ is also formed in the untwinned area of the β matrix, shown in green in the Image Quality (IQ) map in Fig. 5(b). The orientation relationship of SIM α″ with β matrix is {211}β //{110}α″ according to the analysis of stereographic projection pole in Fig. 5(c), indicating that SIM α″ should originate directly from the β matrix. Similar to the observation in unloaded Ti-12Mo-3Zr alloy, the length and thickness of the SIM α″ decrease after unloading from ε = 0.04, meaning that the martensitic transformation β → α″ should be partially reversible.

In a similar way, in-situ EBSD observations are also performed in the Ti-12Mo-10Zr sample, as shown in Fig. 6. The multiplication of deformation bands, identified as 332T by stereographic projection (Fig. 6(c)), can be observed in the β matrix in the IPF map (Fig. 6(a)) during tensile loading. The SIM α″ is not detected during the in-situ loading-unloading (Fig. 6(b)), being in coincidence with the absence of SIM α″ in the post-mortem XRD profile (Fig. 3). It suggests that the SIM α″ is fully suppressed in Ti-12Mo-10Zr sample.

TEM observations

Fig. 7(a) and (b) show the morphology of SIM α" in the deformed Ti-12Mo-3Zr sample by TEM bright-field (BF) and dark-field (DF) images, respectively. It can be seen that a set of SIM α" plates with different orientations appears in the β matrix. The SIM α" plates are, as well, observed in the Ti-12Mo-6Zr sample. Fig. 7(c) and (d) shows a set of parallel SIM α" plates (width < 100 nm) occur in the β matrix between two 332T bands, which is consistent with the in-situ EBSD observations (Fig. 5). In contrast, no SIM α" is observed at TEM length scale in Ti-12Mo-10Zr sample, complementary to the EBSD mapping resolution limited at micrometric scale. The dislocation slip bands are also observed in the BF images of Ti-12Mo-6Zr and Ti-12Mo-10Zr, which can be seen as band-like contrasts in the untwinned β matrix (marked in Fig. 7(c) and (e), respectively). The contrast is due to the inclination of the thin slip bands with respect to the imaging zone axis, resulting in parallel bands cut by the TEM foil surfaces. It is noted in Fig. 7(e) that secondary deformation bands take place inside the primary 332T bands, which is identified by SAED (inset of Fig. 7(f)) as secondary 332T (2 nd 332T). The trace of the 2 nd 332T plane is almost parallel to the slip bands. It is thought that such coincidence could be related to the accommodation of the stress-strain concentration at primary 332T boundary, when assuming that this twin boundary acts as strong barrier to the intersecting dislocation slip. Aiming at clarify the situation, the Ti-12Mo-6Zr and Ti-12Mo-10Zr foils are re-tilted to <111> zone axis to study the dislocation pile-ups in the β matrix near 332T boundary.

Evidence is presented in Fig. 8, displaying the situation of the dislocation pile-up at a 332T twin boundary in an intersecting slip band. An array of different sets of dislocation lines is observed aligned towards the 332T interface. A gradient of the inter-distance between every two lines can be seen decreasing when approaching the twin boundary. The burgers vectors of the dislocation lines is analyzed by using the three g.b=0 extinction conditions around [111]β zone axis (Fig. 8(b)-(d)). It can be seen a set of dislocation lines is composed by a group of screw dislocations attributed to three Burgers vectors b=a/2<111> (labeled by b1, b2 and b3 in the Fig. 8(a)). The b1 dislocations are aligned and piled up well against the 332T interface (Fig. 8(c)), tangled by the dislocations of the other two burgers vectors, showing complex interactions between piled-up dislocations lines and the gliding dislocations on other slip planes. The additional TEM observation is also performed at SIM α" boundary in the Ti-12Mo-6Zr foil. It seems that a dislocation transmission from incoming slip band to the interior of SIM α" occurs at the intersection of the two bands. (Supplementary Fig. S1).

Discussion

Based on the in-situ EBSD analysis and TEM observations, Fig. 9 shows a global schematic illustration about the evolutions of the deformation microstructure as a function of Zr concentrations. The initial structures of the three alloys at ST state are the equiaxed β grains. And the SIM α″and/or twins could be activated after yielding. As described in Fig. 9(a), the SIM α″ appears just after yielding in Ti-12Mo-3Zr sample, and then {130}<310>α″ twinning forms inside the SIM α″. However, the 332T is not activated all along the loading and unloading process. The mechanical stability in the metastable Ti-alloys is close related to the deformation products, since the required fault energy in activating the deformation modes follows the order: SIM α″ < 332T < 112T < dislocation slip [27] . Therefore, the absence of 332T suggests that the critical resolved shear stress (CRSS) of SIM α″ is lower than that of 332T, due to the "low" beta stability of Ti-12Mo-3Zr.

By increasing Zr content to 6%, the 332T is activated prior to SIM α″, presenting a combined TWIP+TRIP deformation mechanism. Further increasing of Zr content to 10%, the SIM α″ is fully suppressed during the loading-unloading process. It can be seen that the fine-tuning of deformation mechanism, from single TRIP via combined TRIP+TWIP to single TWIP, is achieved by progressively increasing the Zr content in Ti-12Mo-xZr alloy. Nevertheless, it has to be noted that the Moeq estimation of Zr fails to predict the mechanical stability of Ti-12Mo-xZr. According to the activation of deformation mechanisms at room temperature, the mechanical stability order of the three alloys at ST state should be proposed as: Ti-12Mo-3Zr < Ti-12Mo ≈ Ti-12Mo-6Zr < Ti-12Mo-10Zr. The result suggests that the Zr effect on the chemical stability in metastable Ti-alloys is not linear as proposed in the Moeq estimation, especially at the low concentration range Zr % < 6%. It has been reported that the Zr addition can achieve dual effects in Ti-Nb based alloys.

One is the compositional effect, namely, Zr addition could decrease Ms temperature, so that suppress the martensitic transformation from the β phase to SIM α″ [28] . The other one is the suppression effect of ω formation for SIM α″ [24] . Al-Zain et al. [29] has observed that the ωath phase can suppress the SIM α″ in Ti-Nb-Mo alloys, implying that the weaker ω precipitation could increase Ms temperature. Therefore, the mechanical stability of Ti-12Mo-xZr alloys could depend on the dual effects of Zr addition. The Zr content up to 3% in Ti-12Mo-xZr suggests that the effect of Zr on the suppression of ωath phase is stronger than the compositional effect. The addition of Zr leads to a lower mechanical stability in Ti-12Mo-3Zr than that in Ti-12Mo. With the increase of Zr content (from 3% to 10%), the effect of suppressing ω phase become weaker and the compositional effect becomes dominant. Therefore, the mechanical stability of Ti-12Mo-xZr (x = 3, 6, 10) increases with the Zr addition. The phenomena are similar with the Sn content on Ti-Nb-Mo-based alloys [30,31] and the Ta content on Ti-Nb based alloys [32] .

Regarding the mechanical behavior, the three alloys display remarkable differences in their yield stress, uniform elongation and strain-hardenability. The yield stress could be described to the synergetic co-operations among the grain-size strengthening, the Zr solution hardening and the CRSS of SIM α″ and/or 332T. The average grain-sizes of the three alloys are measured to be about 100 µm, implying that the grain-size provides an almost equal contribution to the yield stress. The effect of the Zr solution hardening is proportional to its content in the multi-component Ti-alloys [33] .

In addition, regarding the fault energy order mentioned in the above discussion section, a reasonable presumption is the CRSS of SIM α″ should be lower than that of 332T. Thus, plausible explanation can be suggested about the evolution of the yield stress in both Ti-12Mo-3Zr and Ti-12Mo-10Zr alloys from 520MPa to 740MPa. However, it seems that the evolution of the yield stress is less than expected in Ti-12Mo-3Zr and Ti-12Mo-6Zr alloys (from 520MPa to only 590MPa). In other words, the evolution of yield stress seems to be non-linear. The phenomenon probably originates from two-fold factors: 1) the weakening of ωath precipitation hardening effect, due to the ωath phase decrease with the increasing Zr content in Ti-12Mo-xZr alloys; 2) the formation of SIM α″ in Ti-12Mo-6Zr also could be responsible for the limited yield stress. The SIM α″, although observed with low area fraction in some grains (Fig. 5), could still be active enough in other grains orientated favorably to SIM α″, leading to a heterogeneous yielding process. The grains favored SIM α″ could initiate martensitic transformation before the twin nucleation and growth in the grains favored 332T, resulting in an earlier yielding than expected. The similar mechanism has been also reported in the TRIP/TWIP Ti-12Mo [6] and Ti-V-Cr-Al [14] alloys.

Additionally, the three alloys with large uniform elongation show an individual evolution of strain-hardening rate. The large uniform elongation corresponds to enough strain-hardenability, suggesting delay the occurrence of necking [34] . The Ti-12Mo-6Zr alloy shows the highest strain-hardening rate when compared with Ti-12Mo-3Zr and Ti-12Mo-10Zr. The difference in strain-hardening rates could result from the various evolution of deformation microstructure. The deformation microstructure evolution in Ti-12Mo-xZr alloys indicates deformation pathways between the two deformation modes (TRIP and TWIP) as a function of the progressive increase of Zr addition. In other words, the tuning of the deformation mechanisms could be realized by fine-tuning Zr content. Regarding the comparison of the strain-hardening behavior among the TWIP alloys, the back-stress due to the dislocation pile-ups at boundary/slip band intersections can be considered as the key factor responsible for the kinematic hardening [10] . By TEM investigations, the occurrence of dislocation pile-ups is confirmed near the matrix/332T and matrix/SIM α″ boundaries (shown in Fig. 8 andS1). Usually the dislocation pile-ups in these interfaces could cause the back-stress fields. And the stress fields provide a long-range strengthening effect throughout the plastic deformation. Qualitatively, the back-stress is usually evaluated to be linearly proportional to the inverse of the mean inter-boundary spacing, i.e. the distance between two adjacent boundaries of the deformation bands (SIM α″ and/or 332T) [35] . Since the thickness of the bands observed in the three alloys are quite similar, the inverse of the mean inter-boundary spacing is therefore governed by the transformed or twinned fraction of the alloy at a given strain. The kinetic of the nucleation and growth of the deformation bands has been studied by Olson and Cohen [36] for SIM α″ nucleation, and the same model has been adopted by Bouaziz [37] for twinning. The kinetic model has been applied successfully to describe the stress-strain behavior in single TWIP Ti-alloys [10] . Indeed, The strain-hardening behavior of Ti-12Mo-10Zr is consistent with the literature modelling results [10] since the deformation mechanism operated is single TWIP, compatible with the simple shear band assumption used in Olson and Cohen kinetic model. However, in the case of Ti-12Mo-6Zr with mixed TRIP/TWIP mechanisms, the same kinetic TWIP model is not well supported by the experimental results. The Ti-12Mo-6Zr presents a much higher strain-hardening rate than the other two alloys after the yield point (Fig. 2). It is thought that the difference of strain-hardening rate could be related to the localized SIM α″formation in the untwinned β matrix, which is not considered in the current kinetic model. It can be seen from the EBSD and TEM images that short SIM α″ bands are likely to form in the untwinned β matrix between two adjacent twin boundaries. Their formation could be due to the accommodation of the highly localized stress or strain misfits near the twin boundaries and grain boundaries. Similar situations have been reported in Ti-V-Cr-Al alloy [14] and Ti-Mo-W alloy [7] . Then it is quite obvious evident that these local SIM α″ bands further increase the total volume fraction of deformation products than predicted by Olson and Cohen model, and therefore further reduce the inter-spacing of the deformation-induced boundaries in the β matrix (illustrated schematically in Fig. 9). This special fact in the microstructure evolution combing TRIP and TWIP may be one of the reasons to explain the high strain-hardening rate of Ti-12Mo-6Zr when comparing to Ti-12Mo-10Zr. As for Ti-12Mo-3Zr, the strain-hardening behavior at early deformation stage (ε<0.05) is very similar with that of Ti-12Mo-10Zr. The finding is corroborated by the similar microstructure evolutions observed in the two alloys. It might suggest that the Olson and Cohen's kinetic could be valid to describe the nucleation and growth of SIM α″ at the early stage of the plastic deformation in TRIP Ti-alloys. However, in the late stage, the strain-hardening rate increases rapidly to the level close to the strain-hardening rate of Ti-12Mo-6Zr, which means the increasing strain-hardening rate can be due to other reasons. The exact mechanism requires further investigations.

Conclusions

In this work, the influence of Zr element on microstructure evolution, mechanical properties and deformation mechanism of the Ti-12Mo-xZr (x=3, 6, 10) alloys are investigated. The following main results are obtained.

(1) The transition from TRIP effect to TWIP effect in the Ti-12Mo-xZr alloys can be realized by fine-tuning the addition of Zr element.

(2) All the three alloys display an excellent strength-ductility. But the yield stress of Ti-12Mo-10Zr alloy which is much higher than that of Ti-12Mo-6Zr and Ti-12Mo-3Zr, coming from the absence of SIM during deformation, as well as solid solution strengthen of Zr element.

(3) Ti-12Mo-6Zr alloys display a higher strain-hardening rate when compared to Ti-12Mo-3Zr and Ti-12Mo-10Zr, which resulted from a mixture of both TRIP and TWIP mechanisms.

Figure list : Fig. 1 .

 :1 Figure list:

Fig. 2 .

 2 Fig. 2. True strain-stress curves of the ST Ti-12Mo-xZr (x= 3, 6, 10) samples, and the corresponding

Fig. 3 .

 3 Fig. 3. XRD profiles of the ST sample of Ti-12Mo-xZr (x= 3, 6, 10) alloys after deformation until

Fig. 4 .

 4 Fig. 4. EBSD analysis of the Ti-12Mo-3Zr deformed sample taken from the same region at strain

Fig. 5 .

 5 Fig. 5. EBSD analysis of the Ti-12Mo-6Zr (Type A) deformed sample taken from the same region

Fig. 6 .

 6 Fig. 6. EBSD analysis of the Ti-12Mo-10Zr deformed sample taken from the same region at strain

Fig. 8 .

 8 Fig. 8. TEM micrographs of Ti-12Mo-6Zr at ε=0.05: (a) BF image of an array of dislocations piled

Fig. 9 .

 9 Fig. 9. Schematic formation sequence of 332T and SIM α″ during loading and unloading process in

Acknowledgments

The work was supported by National Natural Science foundation of China (Grant No. 51601216 and 51901193); China Postdoctoral Science Foundation (Grant No. 2018M632414); Fund of State Key Lab of Advanced Metals and Materials, University of Science and Technology Beijing (Grant No.2019-ZD03), Fundamental Research Funds for the Central Universities (Grant No. 2017XKQY009), Funds of Industry-University-Research Cooperation in Jiangsu Province (Grand No.BY2018075) and Key Research and Development Program of Shaanxi (Grant No. 2019GY-151). Mr. Bingnan QIAN is sponsored by China Scholarship Council.


Jiawei Cao

Yiming Chen
email: chenym@ysu.edu.cn

Yuanhui Wang

Gang Cheng

Thierry Barrière

Lei Wang

Numerical analysis of fractional viscoelastic column based on shifted Chebyshev wavelet function

Keywords: Fractional constitutive model, Viscoelastic column, Shifted Chebyshev wavelet function, Convergence analysis, Numerical solution

An innovative numerical procedure for solving the viscoelastic column problem based on fractional rheological models, directly in the time domain, is investigated. Firstly, the governing equation is established according to the fractional constitutive relation. Secondly, the resulting equation is transformed into algebraic equation and solved by using the shifted Chebyshev wavelet function. Furthermore, the convergence analysis and the retained numerical benchmarks are carried out to validate the performance of the proposed method. A small value of the absolute error between numerical and accurate solution is obtained. Finally, the dynamic analysis of viscoelastic beam-column problems is investigated with different crosssection shape (circular and square) under various loading conditions (axial compressive force and harmonic load). The displacement, strain and stress of the viscoelastic column at different time and position are determined. The deformation and stress of the viscoelastic column of different materials under the same loading condition are compared. The results in the paper show the highly accuracy and efficiency of the proposed numerical algorithm in the dynamical stability analysis of the viscoelastic column.

Introduction

Non-integer or fractional viscoelastic constitutive models attract more and more attention [START_REF] Machado | Fractional order generalized information[END_REF]. These models are widely used to describe the viscoelastic behaviour of materials with the development of fractional order differential operators. The memory properties of viscoelastic materials are described better with these operators. Xu et al. [START_REF] Xu | Vibration analysis of complex fractional viscoelastic beam structures by the wave method[END_REF] used the wave method to analyze the structure of fractional-order viscoelastic beams, and discussed the effect of fractional order on the structure during the analysis. Mendiguren et al. [START_REF] Mendiguren | A generalised fractional derivative model to represent elastoplastic behaviour of metals[END_REF] proposed a one-dimensional fractional-order elastoplastic model. They analyzed the special cases of first, second, and third terms, which fitted with experimental data, and compared with the classic Hollomon and Ramberg-Osgood models. Sherief et al. [START_REF] Sherief | Effect of variable thermal conductivity on a half-space under the fractional order theory of thermoelasticity[END_REF] studied the variable thermal conductivity of material with half-space method in the context of fractionalorder thermoelasticity theory. The effect of fractional derivative parameters on the behaviour of the solution was also investigated. Therefore, fractional order is more suitable than integer order for solving problems in real life. Liaskos et al. [START_REF] Konstantinos | Implicit analytic solutions for a nonlinear fractional partial differential beam equation[END_REF] derived the implicit analytical solutions for a nonlinear fractional governing equation of viscoelastic Euler-Bernoulli beam. Lin et al. [START_REF] Liu | Fractional boundary layer flow and heat transfer over a stretching sheet with variable thickness[END_REF] established the governing equations related to the left and right Riemann-Liouville fractional derivatives. A comprehensive study was effectuated on the flow and heat transfer of the space fractional boundary layer with a variable thickness stretching sheet. Their results showed considerable importance of fractional calculus in the numerical modelling of viscoelastic behaviour of materials.

The column structure may exhibit complex dynamic behaviour due to different dynamic loads and working environments in building structures and aircraft interior products. The study focus on the dynamic stability of viscoelastic column is paid much attention [START_REF] Ling | Response and stability of SDOF viscoelastic system under wideband noise excitations[END_REF]. In recent years, the analysis of the stability of this structure using Lyapunov functions is well developed [START_REF] Aman | Stochastic meshfree method for elastic buckling analysis of columns[END_REF][START_REF] Jin | An analysis of dynamic stability of an elastic column[END_REF].

Viscoelastic materials are of great significance in vibration system design [START_REF] Jones | Handbook of viscoelastic vibration damping[END_REF]. Viscoelastic damping exhibited in glassy materials. Yang et al. [START_REF] Yang | Numerical analysis of static and dynamic stabilities of viscoelastic columns[END_REF] proposed an adaptive segmentation algorithm for time to describe the static and dynamic stability of the viscoelastic column. Deng et al. [START_REF] Deng | Stochastic stability of a fractional viscoelastic column under bounded noise excitation[END_REF] used the stochastic averaging method to investigate the stability of viscoelastic beam-column under random axial loading. A new set of mean variables was proposed to approximate the response of the dynamical system. The Lyapunov exponent was retained to resolve the eigenvalue problem. A highorder stochastic averaging method was applied to investigate the stochastic stability of a fractional order viscoelastic column with axially random force loading [START_REF] Deng | Higher-order stochastic averaging to study stability of a fractional viscoelastic column[END_REF]. Amit et al. [START_REF] Yadav | Non-linear dynamic instability analysis of mono-symmetric thin walled columns with various boundary conditions[END_REF] studied the non-linear flexural-torsional vibration and instability analysis of the thin-walled column with different cross-sections under cyclic axial loading. The effect of eccentricity, crosssectional shape and boundary condition on the response of the frequency amplitude and instability region was analysed. Floris et al. [START_REF] Floris | Stochastic stability of a viscoelastic column axially loaded by a white noise force[END_REF] examined the stability of a hinged-hinged viscoelastic column under random axial forces of Gaussian white noise. The response of the differential equation was obtained by using Itô's differential law. The necessary and sufficient numerical conditions for stability of the viscoelastic column were that the coefficient matrix were negative real eigenvalues and complex eigenvalues with negative real part. Leung et al. [START_REF] Leung | Parametric bifurcation of a viscoelastic column subject to axial harmonic force and time-delayed control[END_REF] used the Runge-Kutta method to research the bifurcation behaviour of viscoelastic column at different fractional orders and material properties. The stability of viscoelastic column was illustrated by the curves of amplitude and frequency response. Zhang et al. [START_REF] Zhang | An efficient and accurate method for modeling nonlinear fractional viscoelastic biomaterials[END_REF] developed a fractional nonlinear viscoelastic model to describe the large deformation of biomaterials. The fractional viscoelastic constitutive model was integrated in a finite element framework. The proposed method showed high computational efficiency and better convergence rate.

The dynamic analysis of viscoelastic column is not limited to the establishment of the material behaviour equations, but also to solving the fractional differential equations with an efficient numerical algorithm. The conventional approaches are based on the multi-scale [START_REF] Pu | A fractional partial differential equation based multiscale denoising model for texture image[END_REF], galerkin method [START_REF] Wang | Split-step spectral Galerkin method for the two-dimensional nonlinear space-fractional Schrödinger equation[END_REF], finite element method [START_REF] Yang | Finite element method for nonlinear Riesz space fractional diffusion equations on irregular domains[END_REF]. Betancur-Herrera and Munoz-Galeano [START_REF] David | A numerical method for solving Caputo's and Riemann-Liouville's fractional differential equations which includes multi-order fractional derivatives and variable coefficients[END_REF] proposed an innovative numerical method to solve the fractional derivative equations with Riemann-Liouville's or Caputo's form. The main contribution of the method was transforming the fractional derivatives to recurrence equations, which allowed to solve the fractional equations as a system of algebraic equations. Chen et al. [START_REF] Chen | A tunable finite difference method for fractional differential equations with nonsmooth solutions[END_REF] proposed a correction method to obtain the numerical solutions of the fractional differential equations. The improved Lubich's method with few correction terms demonstrated high accuracy compared with other presented methods. Many time domain problems are converted into frequency-domain by using Laplace transform. The solution in the time domain is determined by using the inverse of Laplace transform in the frequency-domain. The inconvenient of this method is the complexity of Laplace transform and its inverse transform. Qi and Guo [START_REF] Qi | Transient fractional heat conduction with generalized Cattaneo model[END_REF] applied the fractional calculus approach to find the numerical solutions of an initial-boundary value problem. The Laplace transform was used to present the solution under series forms of H-function to facilitate the numerical computation. The wavelet basis function is proposed to solve fractional differential equations directly in time domain. Its advantages are high-accuracy, fewer iterations and less computation time [START_REF] Chen | Numerical study of a class of variable order nonlinear fractional differential equation in terms of Bernstein polynomials[END_REF]. Legendre [START_REF] Yuttanan | Legendre wavelets approach for numerical solutions of distributed order fractional differential equations[END_REF], Bernstein [START_REF] Wang | Numerical algorithm for time-fractional Sawada-Kotera equation and Ito equation with Bernstein polynomials[END_REF][START_REF] Bourne | The computation of the degree of an approximate greatest common divisor of two Bernstein polynomials[END_REF] and Chebyshev polynomials [START_REF] Sakran | Numerical solutions of integral and integro-differential equations using Chebyshev polynomials of the third kind[END_REF][START_REF] Foucart | Computation of Chebyshev polynomials for union of intervals[END_REF][START_REF] Sweilam | Second kind shifted Chebyshev polynomials for solving space fractional order diffusion equation[END_REF] are often employed to construct wavelet algorithms. Chen et al. [START_REF] Chen | Numerical algorithm to solve system of nonlinear fractional differential equations based on wavelets method and the error analysis[END_REF] employed wavelet method for solving the non-linear fractional differential equations. The absolute error between the numerical and the algebraic solution was estimated. The shifted Chebyshev polynomials with variable coefficients was proposed to solve the generalized fractional pantograph equations by Wang et al. [START_REF] Wang | Numerical algorithm to solve generalized fractional pantograph equations with variable coefficients based on shifted Chebyshev polynomials[END_REF]. Yu et al. [START_REF] Yu | A numerical method for solving fractional-order viscoelastic Euler-Bernoulli beams[END_REF] applied Quasi-Legendre polynomials to solve the governing equation of the viscoelastic Euler-Bernoulli beam. Hassani et al. [START_REF] Hassani | Generalized shifted Chebyshev polynomials: Solving a general class of nonlinear variable order fractional PDE[END_REF] used a novel class of basis functions based on the shifted Chebyshev polynomials to solve the nonlinear variable order fractional derivative equations. The accuracy and efficiency of the method was confirmed by the convergence analysis and several numerical examples.

Although many researches on the stability of viscoelastic column in various aspects were effectuated, less method has been proposed to solve directly the displacement, stress and strain of viscoelastic column with different materials in the time domain. In this paper, the shifted Chebyshev wavelet function is proposed to perform the dynamic analysis of viscoelastic column directly in the time domain. The material properties are described with a fractional viscoelastic model. The governing equation of column with this fractional viscoelastic model is established. An efficient numerical algorithm based on the shifted Chebyshev wavelet functions is proposed to solve the governing equations. The performance and efficiency of the developed algorithm are investigated. The numerical solutions concerning the displacement, stress and strain of the column have been obtained under axial compressive force and harmonic load conditions. The displacement solutions of viscoelastic columns with different cross-sections and materials under the same load are compared. It provides theoretical basis and numerical approach for the application of viscoelastic materials in engineering.

This paper is structured as follows, section 2 is the introduction of fractional order and viscoelastic material parameters. Section 3 shows the governing equation of the viscoelastic column. Section 4 presents the proposed algorithm. The convergence of the algorithm is analyzed in section 5. The displacement, strain and stress of the viscoelastic column with different materials under different loading conditions are obtained and compared in section 6. The research work is concluded in section 7.

Preliminary

Definition of fractional derivative

In fractional calculus, the different operators called Riemann-Liouville, Caputo and Grünwald-Letnikov operators are respectively defined in [START_REF] Yi | Legendre wavelets method for the numerical solution of fractional integro-differential equations with weakly singular kernel[END_REF], [START_REF] Taghvafard | Phase and anti-phase synchronization of fractional order chaotic systems via active control[END_REF] and [START_REF] Diethelm | The Analysis of Fractional Differential Equations[END_REF]. The definition and properties of the fractional derivatives of Caputo are introduced as follows.

Definition 1 Fractional Caputo derivative c D α t f (t) = d m f (t) dt m , α = m ∈ N + 1 Γ(m-α) t a f (m) (τ ) (t-τ ) α-m+1 dτ, 0 ≤ m -1 < α < m (1)
where 0 < α < 1 is fractional derivative order, c D α t is Caputo fractional derivative operator, f (t) is a continuous function, m is positive integer, Γ (•) is Gamma function, and Γ (z) = ∞ 0 e -t t z-1 dt. Based on the above definition, the following equations are obtained c D α t t m = 0, m ∈ N 0 and m < α

Γ(m+1) Γ(m+1-α) t m-α , m ∈ N 0 and m ≥ α or m / ∈ N 0 and m > α (2)
The Caputo fractional derivative of constant is zero.

c D α t C = 0 ( 3
)
where C is constant. The Caputo fractional derivative is a linear operator.

c D α t (C(f (x))) = C c D α t f (x) (4)

Constitutive equation of linear viscoelasticity

Many materials used in mechanical engineering exhibit elastic and viscous behaviours simultaneously. Lots of physical models have been proposed to describe their viscoelastic characteristic. This viscoelastic relationship between stress and strain could be described more completely and accurately by the fractional order constitutive equation compared to integer order constitutive equation.

The linear viscoelastic properties of the materials can be modelled by different combinations of elastic and viscous elements. The fractional constitutive viscoelastic models, such as Four-parameter Zener model (FZM) [START_REF] Pritz | Analysis of four-parameter fractional derivative model of real solid materials[END_REF], Poynting-Thomson model (PTM) [START_REF] Xu | Creep constitutive models for viscoelastic materials based on fractional derivatives[END_REF], Fractional Maxwell model (FMM) [START_REF] Shen | Unsteady MHD flow and heat transfer of fractional Maxwell viscoelastic nanofluid with Cattaneo heat flux and different particle shapes[END_REF] and Fractional Kelvin-Voigt model (FKVM) [START_REF] Ren | The fractional Kelvin-Voigt model for Rayleigh surface waves in viscoelastic FGM infinite half space[END_REF] have been used to describe the viscoelasticity of the materials. Liu et al. [START_REF] Liu | Boundary layer flow of fractional Maxwell fluid over a stretching sheet with variable thickness[END_REF] employed the time fractional Maxwell constitutive equation to investigate the boundary layer flow. The fractional order operator exhibits better memory characteristic in describing the flow in a highly non-homogeneous transport process. The classical Kelvin-Voigt model is represented by a purely viscous dashpot and a purely elastic string connected in parallel. In fractional Kelvin-Voigt model, the purely viscous dashpot is replaced by a fractional element, as shown in Eq. ( 5). When the fractional order α is equal to 1, the model reduces to the classical Kelvin-Voigt model. The FKVM requires fewer parameters to be identified by the experimental tests. Farno et al. [START_REF] Farno | Comparison between classical Kelvin-Voigt and fractional derivative Kelvin-Voigt models in prediction of linear viscoelastic behaviour of waste activated sludge[END_REF] showed that FKVM can better explain the performance of viscoelastic materials than classic Kelvin-Voigt model. In this paper, FKVM is retained to describe the constitutive relationship between stress and strain of viscoelastic material. FKVM is represented in Figure . 1 and its constitutive equation is as follows:

σ (x, t) = [E 0 + η • c D α t ] ε (x, t) (5)
where

E 0 is elastic modulus, η is viscosity, σ (x, t) is stress, ε (x, t) is strain.

FKVM parameters

The material physical properties can be characterized by the dynamic mechanical analysis. The viscoelastic parameters in the fractional behaviour law, such as fractional order, elastic modulus, can be identified according to the regression analysis. Leung et al. [START_REF] Leung | Parametric bifurcation of a viscoelastic column subject to axial harmonic force and time-delayed control[END_REF] found that the stability domain of the viscoelastic column increases with time-delays when α = 0.8, η = 0.05 by comparing the results with different values of α and η. The numerical results obtained by the direct integration based on the Runge-Kutta 4 th order method are approximately equivalent to the analytical solution. The description of geometrical and material properties [START_REF] Bellum | Investigation on modulus of elasticity of fly ash-ground granulated blast furnace slag blended geopolymer concrete[END_REF][START_REF] Cheng | Physical modelling, numerical simulation and experimental investigation of microfluidic devices with amorphous thermoplastic polymers using a hot embossing process[END_REF] are shown in Table .  1.

Development of fractional governing equation of viscoelastic column

When the axial compressive load F (t) is applied on the viscoelastic column, the transverse displacement appears, as shown in Figure . 2, where the transverse displacement is noted as v (x, t), x and z are axial and transverse coordinate respectively.

The fractional dynamic equation of the column [START_REF] Deng | Higher-order stochastic averaging to study stability of a fractional viscoelastic column[END_REF] is given as

∂ 2 M (x, t) ∂x 2 = ρA ∂ 2 v (x, t) ∂t 2 + β 0 ∂v (x, t) ∂t + F (t) ∂ 2 v (x, t) ∂x 2 , ( 6
)
where density is ρ, damping coefficient is β 0 , cross-sectional area is A, and the bending moment is M (x, t). The bending moment is

M (x, t) = A σ (x, t)zdA. ( 7
)
When the column exhibits the small deformation, the strain of the column is expressed as

ε (x, t) = - ∂ 2 v (x, t) ∂x 2 z (8)
Substituting Eq. ( 5) into Eq. ( 7), it can get

M (x, t) = A [E 0 + η• c D α t ]ε (x, t) zdA (9)
Based on Eq. ( 8) and ( 9), the bending moment of viscoelastic column is

M (x, t) = -E 0 I ∂ 2 v (x, t) ∂x 2 -ηI c D α t ∂ 2 v (x, t) ∂x 2 , (10)
where

I = A z 2 dA is the moment of inertia. Then ∂ 2 M (x, t) ∂x 2 = -E 0 I ∂ 4 v (x, t) ∂x 4 -ηI c D α t ∂ 4 v (x, t) ∂x 4 (11)
Therefore, the fractional governing equation of viscoelastic column is obtained as

ρA ∂ 2 v (x, t) ∂t 2 +β 0 ∂v (x, t) ∂t +E 0 I ∂ 4 v (x, t) ∂x 4 +ηI c D α t ∂ 4 v (x, t) ∂x 4 +F (t) ∂ 2 v (x, t) ∂x 2 = 0. ( 12
)

Shifted Chebyshev wavelet function algorithm

A numerical algorithm is proposed to obtain the solutions of the fractional governing equations. The operator matrices of integer and fractional order are deduced according to the shifted Chebyshev wavelet function.

Shifted Chebyshev wavelet function

The Chebyshev polynomials defined above [-1, 1] are derived. The expression of Chebyshev polynomials is obtained by the following recursive formula

T i+1 (z) = 2zT i (z) -T i-1 (z) , i = 1, 2 • • • ( 13
)
where T is a function expression and

T 0 (z) = 1, T 1 (z) = z. The conversion variable z = 2x L -1 is introduced to obtain the shifted Chebyshev wavelet function in the interval of [0, L]. T i+1 2x L -1 = 2 2x L -1 T i 2x L -1 -T i-1 2x L -1 , i = 1, 2 • • • (14)
Let

T i 2x L -1 = G i (x)
, the shifted Chebyshev wavelet function becomes:

G i+1 (x) = 2 2x L -1 G i (x) -G i-1 (x) , i = 1, 2 • • • (15)
where

G 0 (x) = 1, G 1 (x) = 2x L -1.
The specific expression of the shifted Chebyshev wavelet function of order i is

G i (x) = i i k=0 (-1) i-k (i + k -1)!2 2k (i -k)! (2k)!L k x k , i = 1, 2 • • • (16)
where

G i (0) = (-1) i , G i (L) = 1.
The shifted Chebyshev wavelet function satisfies the following orthogonality relation

L 0 G j (x) G k (x) v L (x)dx = h k ( 17
)
where

v L (x) = 1 √ Lx-x 2 , h k = b k 2 π, k = j 0, k = j , b 0 = 2, b k = 1, k ≥ 1.
Then a series of shifted Chebyshev wavelet function matrix Φ (x) can be written as

Φ (x) = [G 0 (x) , G 1 (x) , • • • , G n (x)] T = A n Z (x) (18)
where

Z (x) = [1, x, • • • , x n ] T , A n = [p ij ] n i,j=0 , p ij =    1, i = 0, j = 0 0, i < j 2 2 L p i-1,j-1 -p i-1,j -p i-2
,j , i ≥ j, i, j = 0 where A n is full rank and invertible.

Approximation of unknown function

For any continuous function of one variable v (x) ∈ L 2 ([0, L]), it can be approximated by the shifted Chebyshev wavelet function.

v (x) = ∞ i=0 c i G i (x) (19)
where c i is the coefficient of shifted Chebyshev wavelet function. The orthogonality of the shifted Chebyshev wavelet function is used to get

c i = 1 h i L 0 v (x) G i (x) v L (x)dx, i = 0, 1, 2, • • •
A finite number of terms are truncated, the function could be written as

v (x) ≈ n i=0 c i G i (x) = C T Φ (x) , ( 20
)
where n is the number of terms of the Chebyshev wavelet function,

C T = [c i ] n i=0 , Φ (x) = [G 0 (x) , G 1 (x) , • • • , G n (x)] T . Any continuous function of two variables v (x, t) ∈ L 2 ([0, L] × [0, T ]) based on the shifted Chebyshev wavelet function can be approximated v (x, t) = ∝ i=0 ∝ j=0 v ij G i (x) G j (t) , (21)
where

v ij = 1 h i h j L 0 T 0 v (x, t) G i (x) G j (t)v L (x) v T (t) dtdx, i, j = 0, 1, 2 • • • .
A finite number of terms are truncated, the function could be obtained

v (x, t) ≈ n i=0 n j=0 v ij G i (x) G j (t) = Φ T (x) V Φ (t) (22)
where

Φ (x) = [G 0 (x) , G 1 (x) , • • • , G n (x)] T ; Φ (t) = [G 0 (t) , G 1 (t) , • • • , G n (t)] T , V = [v ij ] n i,j=0 is the coefficient matrix. If v (x, t) / ∈ L 2 ([0, L] × [0, T ]
) but belongs to the continuous function space, according to the function approximation theory, the effect of the function approximation will be better. Therefore, the shifted Chebyshev wavelet function may be applied in multivariate fractional calculus.

Determination of operator matrix 4.3.1. Integer order operator matrix

Definition 2 If there is a matrix P x such that Φ (x) = P x Φ (x), then P x is called the first order operator matrix of shifted Chebyshev wavelet function.

Φ (x) = (A X Z (x)) = A X Z (x) = A X      1 x . . . (x n )      = A X      0 1 . . . nx n-1      = A X M Z (x) (23)
where A X is obtained by replacing n by X in Eq. [START_REF] Pu | A fractional partial differential equation based multiscale denoising model for texture image[END_REF], and

M = [m ij ] n i,j=0 , m ij = 0, i = j + 1 i, i = j + 1 Because of Z (x) = A -1 X Φ (x)
, then the Eq. ( 23) can be written as

Φ (x) = A X M A -1 X Φ (x) = P x Φ (x) (24)
where

P x = A X M A -1
X is a matrix of first order differential operator for shifted Chebyshev wavelet function.

Definition 3 If there is a matrix G x such that Φ (x) = G x Φ (x), then G x is called the second order operator matrix of shifted Chebyshev wavelet function.

Φ (x) = (Φ (x)) = A X M A -1 X Φ (x) = A X M A -1 X Φ (x) = A X M A -1 X A X M A -1 X Φ (x) = (P x ) 2 Φ (x) = G x Φ (x) ( 25
)
where

G x = (P x ) 2 = A X M A -1 X
2 is a matrix of second order differential operator for shifted Chebyshev wavelet function.

Based on the derivative of first and second order differential operator matrices, the integer order differential operator matrix of the shifted Chebyshev wavelet function can be expressed as

P m x = A X M A -1 X m (26)
where m ∈ Z.

The integer order differential of Φ (x) can be expressed as

d m Φ (x) dx m = P m x Φ (x) = A X M A -1 X m Φ (x) (27)
The partial differential equation terms in Eq. ( 12) could be obtained:

∂ 2 v (x, t) ∂x 2 ≈ ∂ 2 Φ T (x) V Φ (t) ∂x 2 = ∂ 2 Φ T (x) ∂x 2 V Φ (t) = Φ T (x) P x T 2 V Φ (t) = Φ T (x) A X M A -1 X 2 V Φ (t) (28) ∂ 4 v (x, t) ∂x 4 ≈ ∂ 4 Φ T (x) V Φ (t) ∂x 4 = ∂ 4 Φ T (x) ∂x 4 V Φ (t) = Φ T (x) P x T 4 V Φ (t) = Φ T (x) A X M A -1 X 4 V Φ (t) (29) ∂v (x, t) ∂t ≈ ∂ Φ T (x) V Φ (t) ∂t = Φ (x) V ∂Φ (t) ∂t = Φ T (x) V P t Φ (t) = Φ T (x) V A T M A -1 T Φ (t) (30)
∂ 2 v (x, t) ∂t 2 ≈ ∂ 2 Φ T (x) V Φ (t) ∂t 2 = Φ (x) V ∂ 2 Φ (t) ∂t 2 = Φ T (x) V (P t ) 2 Φ (t) = Φ T (x) V A T M A -1 T 2 Φ (t) (31)
4.3.2. Fractional order operator matrix Definition 4 If there is a matrix Q t , such that c D α t Φ (t) = Q t Φ (t), then Q t is called the fractional order operator matrix of shifted Chebyshev wavelet function.

Based on the definition of the Caputo fractional derivative

c D α t t n = Γ (n + 1) Γ (n + 1 -α) t n-α
so fractional derivative of Φ (t) can be written as

c D α t Φ (t) = c D α t (A T Z (t)) = A T c D α t Z (t) = A T c D α t      1 t . . . t n      = A T      0 Γ(2) Γ(2-α) t 1-α . . . Γ(n+1) Γ(n+1-α) t n-α      = A T EZ (t) (32)
where A T is obtained by replacing n by T in Eq. [START_REF] Pu | A fractional partial differential equation based multiscale denoising model for texture image[END_REF], and

E = [e ij ] n i,j=0 , e ij = Γ(i) Γ(i+1) t -α , i = j, i = 1 0, otherwise Because of Z (t) = A -1
T Φ (t), then the Eq. ( 32) can be written as

c D α t Φ (t) = A T EA -1 T Φ (t) = Q t Φ (t) (33)
where

Q t = A T EA -1
T is the fractional differential operator matrix of shifted Chebyshev wavelet function.

The fractional partial differential equation term in Eq. ( 12) could be obtained:

c D α t ∂ 4 v (x, t) ∂x 4 ≈ c D α t Φ T (x) A X M A -1 X 4 V Φ (t) = Φ T (x) A X M A -1 X 4 V c D α t Φ (t) = Φ T (x) A X M A -1 X 4 V A T EA -1 T Φ (t) (34)
4.4. Discretisation of the governing equation of viscoelastic column Eq. ( 12) could be converted into the following form:

ρAΦ T (x) V A T M A -1 T 2 Φ (t) + β 0 Φ T (x) V A T M A -1 T Φ (t) +E 0 IΦ T (x) A X M A -1 X 4 V Φ (t) + IηΦ T (x) A X M A -1 X 4 V A T EA -1 T Φ (t) +F (t) Φ T (x) A X M A -1 X 2 V Φ (t) = 0 (35)
Based on the collection method, the reasonable match points

x i = 2i-1 2(n+1) L, i = 1, 2, • • • , n, t j = 2j-1 2(n+1) T, j = 1, 2, • • •
, n are used to discretize the variable (x, t) to (x i , t j ). Eq. ( 35) is transformed into a set of algebraic equations.

The coefficient

v ij (i = 0, 1, 2, • • • , n; j = 0, 1, 2, • • • , n
) is determined by using Matlab platform and least square method. The numerical solution of the fractional derivative equation can be obtained.

The proposed algorithm can be summarized as shown in Table . 2.

Table 2: Numerical algorithm to resolve the fractional governing equation.

Input: n, α, E 0 , β 0 , η, l, ρ, A, I, F (t) Output: 1. Function approximation v (x, t) ≈ Φ T (x) V Φ (t) 2.
Derive the operator matrix of Φ(x) and Φ(t) in Eq. 24 and 33.

3.

Substitute Eq. 28 -31 and Eq. 34 into Eq. 12.

4.

Eq. 12 could be converted into Eq. 35 5.

Let

x i = 2i-1 2(n+1) L, i = 1, 2, • • • , n, t j = 2j-1 2(n+1) T, j = 1, 2, • • • , n 6.
Solve algebraic equations with MATLAB mathematical software 7.

Obtain the solution of the initial equation v (x, t)

Convergence analysis and validation

In this section, the convergence analysis and numerical example are given to prove the validity and accuracy of the proposed algorithm in solving fractional order equations. The fractional equation solved in section 5.2 is a mathematical example with known algebraic solution. The obtained numerical solutions are compared with the algebraic solution to confirm the efficiency of the proposed algorithm.

Convergence analysis

In the domain Λ = [0, L] × [0, T], the norm of any continuous function

v (x, t) is v (x, t) = sup (x,t)∈Λ |v (x, t)| (36)
The absolute error is defined as

|e n (x, t) | = |v (x, t) -v n (x, t) | (37)
where v (x, t) and v n (x, t) are the algebraic and numerical solutions of fractional governing equation of viscoelastic column.

Theorem 1 If v (x, t) ∈ C 3 (Λ), the absolute error bound is e n (x, t) = v (x, t) -v n (x, t) ≤ N h 3 = O h 3 (38) Proof 1 Let e n,ij (x, t) = v (x, t) -v n (x, t) , (x, t) ∈ Λ n 0, (x, t) ∈ Λ -Λ n (39)
where

Λ n = {(x, t) |ih ≤ x < (i + 2) h, jh ≤ t < (j + 2) h, i, j = 0, 2, • • • , n -2}, v n (x, t) is the quadratic polynomials interpolation function of v (x, t) on Λ n , then e n,ij (x, t) = 1 6
∂ 3 v (ξ 1,i , t) ∂x 3 i+2 i =i (x -x i ) + 1 6 ∂ 3 v (x, ζ 1,j ) ∂t 3 j+2 j =j (t -t j ) - 1 36
∂ 6 v (ξ 2,i , ζ 2,j ) ∂x 3 t 3 i+2 i =i (x -x i ) j+2 j =j (t -t j ) ( 40
)
where x, ξ 1,i , ξ 2,i ∈ [x i , x i+2 ) , t, ζ 1,j , ζ 2,j ∈ [t j , t j+2 ) , i, j = 0, 2, • • • , n -2.

Thus the expression can get

e n,ij (x, t) ≤ 1 6

∂ 3 v (ξ 1,i , t) ∂x 3 i+2 i =i (x -x i ) + 1 6
∂ 3 v (x, ζ 1,j ) ∂t 3 j+2 j =j (t -t j ) + 1 36
∂ 6 v (ξ 2,i , ζ 2,j ) ∂x 3 t 3 i+2 i =i (x -x i ) j+2 j =j (t -t j ) ( 41
)
where i+2 i =i (x -x i ) = sup x∈[x i ,x i+2 ) i+2 i =i (x -x i ) .
Based on the maximum value of

i+2 i =i (x -x i ) in x = i + 1 - √ 3 3
h, the following expression could be obtained

i+2 i =i (x -x i ) ≤ 2 √ 3h 3 9 , ∀x ∈ [x i , x i+2 ) ( 42
) j+2 j =j (t -t j ) ≤ 2 √ 3h 3 9 , ∀t ∈ [t j , t j+2 ) ( 43
)
Substituting Eq. ( 42) and ( 43) into the Eq. ( 41) to obtain:

e n (x, t) ≤ √ 3h 3 27
∂ 3 v (x, t) ∂x 3 + ∂ 3 v (x, t) ∂t 3 + √ 3h 3 27
∂ 6 v (x, t) ∂x 3 ∂t 3 = N h 3 = O(h 3 ) (44)
As conclusion, the theorem is proved.

Based on the proof of the theorem, the accuracy of the proposed algorithm is confirmed and it can used to solve the fractional order equations effectively.

Application of developed algorithm for solving mathematical problem

The objective of this section is to investigate the accuracy of the proposed numerical algorithm.

The order of convergence is calculated to describe the effectiveness of the proposed numerical algorithm. The order of convergence q is expressed as follows:

q = log e n max e n-1 max log hn h n-1 (45)
where e n max and e n-1 max are the maximum absolute error,

h n = π 2n , h n-1 = π 2(n-1) .
The following equation is considered as a numerical example:

50 ∂ 2 v (x, t) ∂t 2 + 1.6 ∂v (x, t) ∂t + 0.15 ∂ 4 v (x, t) ∂x 4 + 0.15D α t ∂ 4 v (x, t) ∂x 4 +0.1F (x, t) ∂ 2 v (x, t) ∂x 2 = 0 ( 46
)
The boundary and initial conditions are:

   v (0, t) = v (5, t) = 0 v (0, t) = v (5, t) = 0 v (x, 0) = v (x, 0) = 0 ( 47
)
where α = 0.89, F (x, t) =

x 2 (5-x) 2 (100+3.2t)+3.6(t 2 + Γ(3) Γ(2.11) t 1.11 ) (6x-1.2x 2 -5)t 2 .
The algebraic solution of the Eq. ( 46) is v (x, t) = x 2 (5-x) 2 t 2

145

. The absolute error, defined in Eq. [START_REF] Diethelm | The Analysis of Fractional Differential Equations[END_REF], between the numerical and algebraic solutions decreases gradually with the number of terms of shifted Chebyshev wavelet functions. The order of convergence is improved with the increase of the number of terms. The efficiency and accuracy of the proposed algorithm based on the shifted Chebyshev wavelet function is demonstrated. Therefore, n = 6 is retained in the numerical algorithm for the following solutions. It is proved as a highly useful method to solve the fractional governing equation of the viscoelastic column under various load conditions.

Results and discussions

The actual fractional governing equation of the viscoelastic column is solved by using the proposed algorithm.

The displacement, stress and strain of the viscoelastic columns with different cross-section shape and materials under various load conditions are investigated in this section. The square and circular cross-section parameters are shown in Table . 5.

The boundary and initial conditions are as following: The column with the square cross-section is selected in this study. The metal material parameters of the viscoelastic properties are defined in Table .  1. When the axial compression force 1000 N is applied on the column and the numerical solutions of the displacement of the column at various position and time are shown in Figure . 4. The displacement of the two ends of a viscoelastic column is always zero and is not affected by time, which is coherent with the boundary conditions. In other locations, the displacement of the viscoelastic column gradually increases with time and the maximum value is obtained at t = 1 s. The maximum transverse displacement of the viscoelastic column is obtained at The displacement of the fractional viscoelastic column under the periodic load was analyzed by Li et al. [START_REF] Li | Dynamical stability of viscoelastic column with fractional derivative constitutive relation[END_REF]. Therefore, the proposed algorithm was limited to a specific displacement solution. The algorithm used in this paper exhibits wider application and higher efficiency.

v (0, t) = v (l, t) = 0, ∂v (0, t) ∂x = ∂v (l, t) ∂x = 0 (48)
0 1 1 5 v(x,t) [m] ×10 -12 2 4 t [s] 0.5 3 x [m]

Effect of the column cross-section on the numerical solution of governing equation

The effect of cross-section parameters on the transverse displacement of the viscoelastic column is investigated in this section. Two viscoelastic columns with different cross-section shape (square and circular) are retained in this study. Their cross-sectional areas are the same and the cross-section parameters are summarized in Table . 5. The transverse displacement of the viscoelastic column with circular cross-section shape under the axial and8, the evolution of the transverse displacement of the column with square and circular crosssection is similar. However, the transverse displacement of the column with circular cross-section is smaller than that with square cross-section, when the columns exhibit the same cross-sectional area and are under the same axial compressive load condition. square cross-section and subjected the same applied compressive force 1000 N. The proposed algorithm based on the shifted Chebyshev wavelet function is used to calculate the transverse displacement of the viscoelastic column. The maximum value of the transverse displacement of concrete column is approximately 1.913×10 -10 m, however, the maximum value of the transverse displacement of polymer column is approximately 8.352 × 10 -4 m. The maximum value of transverse displacement of polymer column is larger than that of the concrete column. The stress and strain curves of the viscoelastic columns with polymer and concrete material in function of position and time are shown in Figures. 9 and 10. The stress and strain are calculated according to Eqs. ( 5) and [START_REF] Aman | Stochastic meshfree method for elastic buckling analysis of columns[END_REF]. The trend of the stress and strain curves is similar with different materials. Therefore, the values of stress and strain of polymer column are extremely larger than those of the concrete column. The maximum value of the strain in polymer column is approximately 5 × 10 -4 , when the maximum value of the strain in concrete column is approximately 1 × 10 -10 . The maximum value of the stress in polymer column is approximately 1900 Pa, when the maximum value of the stress in concrete column is approximately 2.5 × 10 -3 Pa. This is mainly due to the smaller value of elastic modulus of the polymeric material.

0 1 5 5 ×10 -8 v(x,t) [m] 4 t [s] 0.5 3 x [m]

Conclusions

In this paper, a numerical algorithm based on the shifted Chebyshev wavelet function is developed to solve the fractional governing equation of the viscoelastic column directly in the time domain. The fractional derivative model is used to describe the viscoelastic behaviour of the column with different materials. The transverse displacement, stress and strain of the viscoelastic column are obtained under various load conditions: axial compressive force and harmonic load. 1. Base on the analysis of the mathematical example, it is proved that the proposed algorithm exhibits high accuracy in solving the fractional order equations. The accuracy of the algorithm improves with the number of terms retained in the shifted Chebyshev wavelet function.

2. The transverse displacement of viscoelastic column under axial compressive force and harmonic load are obtained. The displacement increases with the applied axial compression force. The magnitude and period of the harmonic load affects the deformation of the viscoelastic column.

3. The transverse displacement of the viscoelastic column with circular cross-section shape is smaller than that with square cross-section, when the columns exhibit the same cross-sectional area and are under the same axial compressive force. 4. The polymer column exhibits the larger displacement than the concrete column, due to its lower elastic modulus. The proposed algorithm demonstrates good sensitivity of material in the numerical simulation of the deformation of the viscoelastic column.

Figure 1 :

 1 Figure 1: Schematic representation of FKVM.

Figure 2 :

 2 Figure 2: The displacement of viscoelastic column.

Figure 3 :

 3 Figure 3: The comparison of the numerical and algebraic solutions at different matching point when n = 5. The numerical solution of the equation is obtained with different number of terms of the shifted Chebyshev wavelet function (n = 4, n = 5 and n = 6) based on the proposed algorithm at (x i , t j ) ∈ [0, 5] × [0, 1]. The numerical solution is highly consistent with the algebraic solution, as shown in Figure. 3 (n = 5). The absolute error values between the numerical and algebraic solutions at different matching points are listed in Table.3. The maximum absolute error value and order of convergence of different n are shown in Table.4.The absolute error, defined in Eq.[START_REF] Diethelm | The Analysis of Fractional Differential Equations[END_REF], between the numerical and algebraic solutions decreases gradually with the number of terms of shifted

  3. The maximum absolute error value and order of convergence of different n are shown in Table.

  4.

  168 × 10 -5 71.9530 5 6.557 × 10 -7 136.6654 6 1.162 × 10 -8 199.0116

Figure 4 :

 4 Figure 4: The displacement solution of viscoelastic column under axial compressive force: (a) F = 1000 N, (b) at different times.

Figure 5 :

 5 Figure 5: Evolution of (a) stress (Pa) and (b) strain of viscoelastic column under the axial compressive force.

Figure 6 :

 6 Figure 6: The numerical solution of displacement under different load conditions: (a) F = 1500 N, (b) F = 2000 N.

Figure 7 :

 7 Figure 7: Transverse displacement of viscoelastic column under different harmonic loads: (a) F = 1000 cos (0.05πt), (b) F = 1000 cos (0.1πt), (c) F = 1000 cos (0.15πt), (d) F = 1500 cos (0.05πt)

6. 3 .

 3 Effect of material property on the numerical solution of governing equationThe deformation of the viscoelastic columns with two different materials (concrete and polymer) are compared in this part. The material parameters are summarized in Table.

1 .Figure 8 :

 18 Figure 8: Displacement of a viscoelastic column with circular cross-section at different times.

Figure 9 :

 9 Figure 9: Stress and strain of concrete viscoelastic column.

Figure 10 :

 10 Figure 10: Stress and strain of polymer viscoelastic column.

Table 1 :

 1 Description of geometrical and material properties of the viscoelastic column.

	Metal	5	7750	600	21
	Concrete	5	4000	600	22680
	Polymer	5	1300	600	3162.4

l(m) ρ(kg/m 3 ) β 0 (1) E 0 (MPa)

Table 3 :

 3 Evolution of absolute error using the developed method in function of several values of n.

	(x, t)	Algebraic solution	n = 4	Absolute error value n = 5	n = 6
	(0, 0) (0, 0.5) (0, 1) (2.5, 0) (2.5, 0.5) (2.5, 1) (4, 0) (4, 0.5) (4, 1)	0 0 0 0 0.06735 0.2694 0 0.02759 0.1103	1.182 × 10 -8 7.567 × 10 -10 1.284 × 10 -11 3.966 × 10 -9 2.259 × 10 -10 1.19 × 10 -12 1.02 × 10 -8 6.515 × 10 -10 1.988 × 10 -11 3.766 × 10 -9 5.163 × 10 -10 1.258 × 10 -11 5.856 × 10 -6 3.29 × 10 -7 5.841 × 10 -9 1.168 × 10 -5 6.557 × 10 -7 1.162 × 10 -8 4.726 × 10 -9 6.988 × 10 -10 1.872 × 10 -11 2.236 × 10 -6 1.251 × 10 -7 2.225 × 10 -9 4.462 × 10 -6 2.498 × 10 -7 4.436 × 10 -9

Table 4 :

 4 The maximum absolute error value and order of convergence for different n.

Table 5 :

 5 Square and circular cross-section parameters of viscoelastic column.

		d(m)	A(m 2 )	I(m 4 )
	square cross-section circular cross-section 0.1128 0.1	0.01 0.01	8.3 × 10 -6 8.03 × 10 -6
	v (x, 0) =	∂v (x, 0) ∂t	= 0.	(49)
	6.1. Effect of the load condition on the numerical solution of governing
	equation
	6.1.1. Axial compressive force

Acknowledgements

This work is supported by the Natural Science Foundation of Hebei Province (A2017203100) in China and the EIPHI Graduate School (contract ANR-17-EURE-0002).