Metal -ceramic functionally graded materials (FGMs) have been extensively used in aerospace engineering where high strength and excellent heat insulation materials are desired . In this paper , load distribution in threads of the Thermal Protection System used bolted joint made up of porous ZrO 2/(ZrO 2+Ni) FGMs is investigated by ABAQUS codes . The bolted joint is subjected to reentry heating corresponding to th e Access to Space Vehicle. Effects of bolt -nut parameters including thread tooth profile, thread pitch, and modulus ratio of bolt to nut on load distribution in threads are analyzed in detail. It is found that uneven load distribution in threads occurs at elevated temperature, which mainly focus es on the first two threads closest to the nut bearing surface , with the first thread carr ying 74% of the total load . Bolt-nut parameters have great effects on load distribution in threads, with trapezoidal thread, extra fine thread and greater modulus ratio of bolt to nut leading to more evenly distributed load . Further studies show that nut shape has significant effect s on load distribution in thread s, the optim ized nut is designed to make the maximum load bearing r atio of the thread decrease to 30 .21%, and thus the service reliability of the bolted joint is greatly improved. Keywords: Functionally graded composite joints ; Thermomechanical; Finite element modeling; Load distribution * Corresponding author. Tel.: +86 010 627 57417 Email address: wbzhou@pku.edu.cn 2 1. Introduction Thermal Protect ion System (TPS) has been widely used as the thermal protecting layer which is usually covered on the surface of the aircraft . Because of the complex structure and shape of the aircraft , a unibody design and manufacture of TPS is infeasible , therefore it is inevitable to connect lots of thermal protection components . Consequently , heat short s resulting from the gap and joints between thermal protection components cannot be ignored [1]. Heat short s may lead the interior part of the fuselage structure to an extreme high temperature , severely damagi ng the heat insulation property of the TPS. During the process of service, the TPS has to withstand not only thermal loads, but also mechanical loads as well as harsh chemical environments, repeatedly without failure [2]. Therefore, in the research of the TPS, the analysis of heat short s effect and the use of high temperature resistant joints with high strength are indispensable. Many researche s have been conducted on metallic thermal protection system, including the design and characterization of thermal insulation material/ structure [3 -8], however, rare studi es have been done on heat shorts blocking joints . In addition , most researches focused on metal joints [ 9-15], which had quite poor heat insulation properties. Some studies were carried out on C/SiC joints [1 6-18], yet the com plicated fabrication approaches for geometrically complex composites such as bolted joints have certain limitations in terms of size and shape, which have greatly l imited their wider applica tions. The studies for the thermodynamic properties of these joints which are quite necessary for joints used on TPS are less . Functionally graded materials (FGMs) are ideal candidates for the applications requiring multifunctional performance. For example , the metal -ceramic FGMs can be designed to take advantages of the mechanical strength of metals and the heat and corrosion resistance of ceramics [ 19-20]. In our previous work, we fabricated and experimentally investigated the microstructure and therm omechanical properties of the porous ZrO 2/(ZrO 2+Ni) ceramics and found that they had a low thermal conductivity and high strength [21-22]. Therefore, the ceramic -metal FGMs can be a potential candidate for heat short s blocking joints. Most of the bolted joint s have uneven load distribution in threads , which are destined to affect the mechanical property and service reliability of the joints. Studies have show n that failure in threads are the most common failure mode of the bolted joint when subjected to axial tensile load [16, 23]. This failure mode makes the strength of the bolted joint become less fully utilized. So the load distribution in threads is an important f actor that cannot be neglected . However, few studies have been reported to date on the load dist ribution in threads of the bolted joint under thermomechanical circumstances. Therefore, i t is of great importance to study the load distribution in threads of porous ZrO 2/(ZrO 2+Ni) functionally graded bolted joint (FGBJ) subjected to thermomechanical load ing, as well as the factors influencing the load distribut ion. 3 Finite element analysis (FEA) has been widely used to analyze the thermal and mechanical properties of the joints [24-26], since the method through experimental study to obtain the thermomechan ical properties of the metal -ceramic functionally graded bolted joint (FGBJ) is accurate but time -consuming, so the theoretical simulation and optimization would be rather important for the successful implementation of metal -ceramic functionally graded joi nting technology in a wide variety of high temperature application. In this paper, firstly a thermomechanical coupled finite -element (FE) model is constructed to study the load distribution in threads of the metal -ceramic functionally graded bolted joint (FGBJ) , whose material systems (porous ZrO 2/(ZrO 2+Ni) FGMs) are fabricated by cold isostatic pressing and pressureless sintering (CIP -PLS) . Microstructures and mechanical properties of the FGMs are experimental ly studied in our previous work [ 21- 22]. Then load distribution in threads of the FGBJ with different bolt -nut parameters including thread tooth profile, thread pitch, and modulus ratio of bolt to nut are invest igated, thus the optimized bolt -nut parameters are proposed for the FGBJ. Finally the optim al structural design of the bolted joint under thermomechanical circumstances is done to reduce the uneven ly distribut ed load , and thus the service reliability of the bolted joint is also improved. The findings will highlight the optimization of the bolt-nut parameters and the structure of the bolted joint in a wide variety of high temperature application . 2. Preparation of the bolted joint In our previous work [ 21], we have fabricated porous ZrO 2/(ZrO 2+Ni) functionally graded samples, as shown in Fig. 1a. The height of the bearing load layer (layer a , 12% porosity ), transition layer (layer b , 15% porosity ) and middle layer (layer c , 30% porosity ) are la=10mm, lb=2mm, and lc=10mm in the axial direction , respectively . These samples are to be further manufactu red into the functionally graded bolted joints . Fig. 1b shows t he schematic of the bolted joint. The bolt head , bolt thread and nut all co nsist of the bearing load layer (ZrO 2+30 vol. %Ni ) to ensure the mechanical strength and service reliability . The bolt shank without threads is mainly comp rised of the middle layer (porous ZrO 2), which fulfills the function of heat resistance. The transition layer ( ZrO 2+15 vol. %Ni ) on the bolt shank is designed to reduce the residual stress es which are generated inevitab ly due to the mismatch in coefficients of thermal expansion between different layers. Also, t he transition layer can serve the purpose of cutting the fillet. According to the standard handbook of fastening and joining [27], the height and diameter of the b olt head are k=5.5mm and s=12.8mm , respectively. The diameter of the bolt shank is d=7.8mm, which is also the major diameter of the external thread. The height and diameter of the nut are m=6.8mm and s=13mm, respectively. The thread pitch is P=1.25mm. The number of engaged threads on the nut is thus calculated to be five. The root radius of the thread is taken to be L/6, where L is the thread height. 4 The bolted joint with the preload is called high -strength bolt, which usually need s only 2~3 threads out of the nut when the bolt and nut are fully engaged. Therefore, we suppose the number of threads on bolt shank is eight, and the height of the bolt shank with threads is accordingly to be b=10mm. 3. Theoretical simulation method 3.1. Finite element model ing Based on ABAQUS code s, a two-dimensional ( 2D) axisymmetric finite element model is establish ed using the advantage of symmetry, as shown in Fig. 2, to study the load distribution in threads of the bolted joint subjected to thermomechanical loading . The error caused by neglecting the helix angle and modeling the bolt and nut axisymmetrically has been ignored in this study , since t he helix angle is quite small (less than 2 deg) on most bolts and thus the three -dimensional effect can be neglected . Actually , using 2D approximation for such a problem was also reported in the earlier experimental photoelastic studies [28]. This model comprises 48147 elements and 48976 nodes. The finite element mesh at t he crest and root of the thread is refined using a very small el ement length . To improve the convergence , nodes between two mutually contacted threads are coincident. The thermomechanical response of the bolted joint is acquired with two steps. First ly, heat transfer is simulated to obtain the time -dependent temperatur e distribution of the bolted joint subjected to the applied thermal loads and boundary conditions. Then , the load distribution in threads caused by the transient temperature distributions coupled with the external mechanical load are determined. The soluti ons are facilitated by using the same mesh for both the heat transfer and thermomechanical problems. The DCAX 4 elements are adopte d for the heat transfer analysis and the CAX 4 elements for the therm omechanical analysis. Table 1 -2 show the material propert ies used in the finite element modeling . It should be pointed out that, the material properties of ZrO 2 and Ni are only experimental ly tested at some discrete temperature points, other parameters that have not been measured will be given by interpolation i n the simulations. The uncertainty in the material properties inevitably occurs, the scatter in the data here is below 5% . Each composition ’s material constants , such as Young’s modulus, thermal conductivity, coeffic ient of thermal expansion (CTE), specifi c heat and Poisson’s ratio can be calculated as follows [ 29-30]: 0 0 0 0 0 0 0 0 0 0 2/3 2/3 00 1/3 1/3 2/3 2/3 008 (1 )(1 )(23 8 ) 8(1 )(23 8 ) (5 )(37 8 ) ,, [(1 ) ] (1 ) [(1 ) ]g gEPEP CC PP P P P P                              (1) 5 where P is the porosity, g is the thermal conductivity of air , 0E, 0C, 0, 0and 0 are all dependent on the ceramic volume fraction and can be expressed from the Vegard's rule below [ 31] 0 03( )[1 ]3 ( )(1 )ce ce m m m ce ce ce m ce mM M V M V V         (2) where M represents the material property parameters ( E , , C,  and  ), ceV donates the ceramic volume fraction, and subscript s ‘ce’ and ‘m’ represent ceramic and metal, respectively. The material model adopte d is based on the tension experiment o f the porous ZrO 2/(ZrO 2+Ni) FGMs. Fig. 3 shows t he tensile load-displacement relationship. It should be pointed out that, the stress and strain acquired by experimental test are only the nominal stress and nominal strain, which should be transformed into t he true stress and true stain used in ABAQUS. The relationship between the nominal stress (strain) and true stress (strain) can be expressed as (1 ) ln(1 )ture nom nom ture nom       (3) where nom , nom are the nominal stres s and nominal stain, ture , ture are the true stress and true stain, respectively. By Eq. (3), the relationship between plastic strain and true stress which is input into ABAQUS is given by (1 )ln(1 )nom nom pl true el nomE        (4) where Eis the Young’s modulus and pl , el are the plastic strain and elastic stain, respectively . 3.2. Simulation pro cedure Load distribution in threads of the functionally gra ded composite joint subjected to thermomechanical loading is acquired with two steps. Firstly, the heat transfer simulation is conducted to acquire the time - dependent temperature field of the bolted joint . The initial temperature of the bolted joint is 25oC (i.e., room temperature), the ambient air temperature of the top surface and bottom surface of the bolted joint are set as 1200oC and 25oC, respectively. Table 2 shows t he temperature dependent convection and radiation coefficient s. The total loading tim e is 2000s. Then a coupled thermal -mechanical approach based on ABAQUS nonlinear finite element codes is adopted after the heat transfer process . Contact is defined between all the contacted surfaces . The friction coefficient of bolted joints usually range s from 0.1 ~0.2 [32], and thus the same friction coefficient 0.15 is used here between 6 all the contacted surfaces . The axis of symmetry of the bolt is restrained in the radial direction. The upper surface of the nut contacted with the plate is restrained in the horizontal y direction, while the lower surface of the nut is free. An initial preload of F0=1600N (0-2003s) is applied to the bolt through the bolt load method in ABAQUS. It should be noted that this preload set by this method can automatically vary with the length of the bolt when the thermomechanical process starts. A constant surface load of F=1000N (3-2003s) is applied to the top surface of the bolt head . It should be pointed out that, the load s are applied gradually (little load s are applied initially) and the constraints are release d stepwise ( extra constraints are need ed initially) to realize the convergence . The thermal load (3-2003s) is applied by importing the transient temperature distribution calculated f rom the heat transfer process . Therm omechanical coupling behavior of the bolted joint in 3 -2003s is numerically simulated by applying the mechanical load and thermal load simultaneously, which is called transient thermodynamic behavior. 4. Results and Discussion 4.1. Load distribution in thr eads Fig. 4a shows the axial load distribution in engaged threads of the bolted joint in 0 -2003s . As shown in. Fig. 4a, the load distribution at the end of the thermomechanical process is 743.4N , 257.0N for thread 1 and thread 2, respectively. The other t hreads almost carry no load s. Here, thread 1 is the engaged thread which is closest to the nut bearing surface, thus thread 5 is the engaged thread which is closest to the free surface of the nut. All the loads in threads add up to 1000.4N, which is very c lose to the applied surface load ( F=1000N). This also means that the preload relaxation of the bolted joint happens at elevated temperature, the preload gradually decreases to a near zero value at the end of the thermomechanical process. Dividing the load distribution i n threads by the corresponding external load ( F=1000N), the load bearing ratio of thread 1 -5 is calculated as 0.74, 0.26, 0, 0, and 0, respectively. Therefore, we conclude that uneven load distribution in threads of the bolted joint occurs at elevated temperature, and the load distribution in threads under thermomechanical circumstances mainly focuses on the first two threads which are closest to the nut bearing surface, with the first thread carrying almost 74% of the total load. To better il lustrate s this phenomenon, we also plot the load distribution in threads of the bolted joint only subjected to mechanical loading. Fig. 4b show s the partial enlarged drawing of Fig. 4a in 0-100s. It should be noted that the bolted joint is only subjected t o mechanical loading ( preload of F0=1600N) in 0-3s. The results show that the load distribution at 3s is 571.3N , 305.8N , 222.8N , 206.8N and 293.1N for thread 1 to thread 5, respectively. All the loads in threads add up to 1599.8N, which is very close to th e applied surface load 7 (F0=1600N). Dividing the load distribution i n threads by the corresponding load ( F0=1600N), the load bearing ratio of thread 1 -5 only subjected to mechanical loading is calculated as 0.36, 0.19, 0.14, 0.13, and 0.18 respectively. From Fig. 4b, we can see that the load bearing ratio of threads under thermomechanical circumstantces is more uneven than that of the mechanical loading alone. Also, we may figure out how the load distribution behaves at the initial stage of the thermomechani cal process. When the thermomechanical process starts, the loads in threads all increase firstly (3 -10s) although the preload decreases, then all gradually decrease (3-36s) before they are reallocated and reach steady at the coupled loading of the preload, surface load (F0=1600N) and thermal load. 4.2. Effects of bolt -nut parameters on load distribution in threads In this section, t he effects of bolt -nut parameters including thread tooth profile, thread pitch, and modulus ratio of bolt to nut on the load di stribution in threads and service reliability of the bolted joint subjected to thermomechanical loading are analyzed in detail. Load distribution in threads of the bolted joint with d ifferent thread tooth profiles such as square thread , trapezoidal thread and triangular thread are investigated respectively and are shown in Fig. 5 . For each thread tooth profile, the bolted joint is subjected to the same thermomechanical loading as mentioned before. The thread pitch is fixed at P=1.25mm, and the modulus ratio of bolt to nut is kept to be Eb/En=1. Here, Eb and En are the modulus of the bolt and nut, respectively. As shown in Fig. 5, load bearing ratio of thread 1 is the biggest for triangular thread, followed by square thread and trapezoidal thread. Therefore, thread tooth profile has certain effects on the load distribution i n threads of the bolted joint, trapezoidal thread leads to the most uniform load distribution in threads, with square thread and triangular thread coming next, respectively. But the effects are relatively limited, more effective parameters still need to be investigated to better reduce the uneven load distribution in threads. Fig.6 shows l oad distribution in threads of the bolted joint with different thread pitch P, including coarse thread (P=1.25mm), fine thread ( P=1mm) and extra fine thread ( P=0.75mm). The diameter of the bolt shank is fixed at d=7.8mm, then the corresponding values of P/d are 0.16 , 0.128 and 0.096 , respectively. Trapezoidal threads are used here, and the modulus ratio of b olt to nut remains to be Eb/En=1. For each thread pitch P, the bolted joint is subjected to the same thermomechanical loading as stated before. As shown in Fig. 6, with the decrease of the thread pitch, load bearing ratio of thread 1 gradually decreases. However, t he second biggest load bearing ratio (thread 2) remains almost unchange d, which means the decrease of load bearing ratio of thread 1 has mainly been transferred to the last several threads (thread 3 -5) which carry almost no load before, 8 thus makin g the load bearing ratio in threads more uniform. Therefore, thread pitch has great effects on the load distribution is threads of the bolted joint, extra fine thread leads to the most uniform load distribution in threads. This is because with the decrease of the thread pitch, thread root becomes thinner and its stiffness will decrease correspondingly, which makes the thread easier to suffer from the bending defo rmation . Thus the deformation compatibility between the threads and bolt -nut body will also incr ease, making the load distribution in threads of the bolted joint become more even . Load distribution in threads of the bolted joint with modulus ratio of bolt to nut varying from Eb/En=1:5 to Eb/En=5:1 are investigated respectively and are illustrated in Fig. 7 a-c. For each modulus ratio of bolt to nut, the bolted joint is subjected to the same thermomechanical loading as previously mentioned. Trapezoidal threads are used here . In Fig. 7a, the thread pitch is fixed at P=1.25mm. As shown in Fig. 7 a, with th e increase of modulus ratio of bolt to nut, load bearing ratio of thread 1 gradually decreases, while other threads except thread 2 all gradually increase correspondingly. Load bearing ratio of thread 2 firstly increases when modulus ratio of bolt to nut i ncreases from Eb/En=1:5 to Eb/En=1:3, then remains almost unchanged when modulus ratio of bolt to nut keeps increasing. Also, we can see from Fig. 7 a that with the increase of modulus ratio of bolt to nut, especially when modulus of the bolt is greater tha n that of the nut, unevenly distributed load in threads of the bolted joint can be improved significantly. This is because when the nut is softer than the bolt, it is possible for nut threads to distort a small amount to conform to the expanding contours o f the bolt threads when tensioned. It has been reported before that softer nut leads to more even load distribut ion in threads of the bolted joint only subjected to mechanical loading [33], here our results confirm that it also applies to the bolted joint under thermomechanical circumstances. Therefore, modulus ratio of bolt to nut has remarkable effects on the load distribution in threads of the bolted joint, greater modulus ratio of bolt to nut leads to more uniform load distribution in threads, thus the service reliability of the bolted joint is also improved. We also study the cases where the thread pitch is slightly perturbed, especially for the first two thread pitches. Fig. 7b shows the case where thread pitch of thread 2 is chang ed to be extra fine t hread ( P=0.75mm), while the others remain as P=1.25mm. Fig. 7c shows the case where thread pitch of thread 1 is changed to be P=0.75mm, while the others keep to be P=1.25mm. As shown in Fig. 7b -c, with the increase of modulus ratio of bolt to nut, unevenly distributed load in threads can be improved greatly. This conclusion is in acco rdance with the case where the thread pitch is fixed at P=1.25mm in Fig. 7a . Moreover , by c ompar ing Fig. 7b -c with Fig.7a, we can see that when thread pitch es of thread 1 and t hread 2 are slightly perturbed (decreased), the unevenly distributed load in threads are both improved for the two cases , and the improvement for the case 9 where thread pitch of thread 1 is perturbed is more pronounced than the case where thread pitch of th read 2 is perturbed. This is also consistent with our result in Fig. 6 that with the decrease of the thread pitch, load distribution in threads become s more even . Here from Fig. 7b -c we further conclude that the thread pitch of thread 1 is the domi nant fac tor for improving the unevenly distributed load through changing the thread pitch. 4.3. Optimization of the structure of the bolted joint Unevenly distributed load in threads of the functionally graded bolted joint subjected to thermomechanical loading can greatly degrade its high-temperature mechanical performance. Research has show n that by designing the nut, uneven load distribution in threads can be significantly improved [34], but whether it is effective at elevated temperature still needs to be studie d. Inspired by previou s research, an internal groove and an external groove are constructed on the nut and studied under transient thermomechanical circumstance , as shown in Fig. 8-11. The one without a groove is also plotted as comparisons. Fig. 8 shows t he structure of the bolted joint with an internal groove , where b and l are the width of the nut after and before grooving . Firstly we study the effect of b/l on load bearing ratios of threads, and t hree different cases are simulated for the functionally g raded bolted joint with an internal groove , (i) b/l=1/4, (ii) b/l=2/4, (iii) b/l=3/4. The depth of the internal groove is fixed at a=P/4, where P=1.25mm is the thread pitch . The load distribution in threads of the bolted joint is illustrated in Fig. 9 a. As shown in Fig. 9a, the maximum load bearing ratio of threads still occurs at the first thread, and gradually decreases as the decrease of b/l, while the load bearing ratio of the last three threads all increase with the decrease of b/l, thus making the loa d bearing ratio in threads more uniform. Then, we study the effect of a on load bearing ratios of threads . Here a is restricted to be not greater than 0.5 P to ensure the fully engagement of the first thread. Three different cases are simulated, (i) a/P=1/4, (ii) a/P=1/3, (iii) a/P=1/2. The width of the nut after grooving is fixed at b/l=1/4. Fig. 9b shows that the maximum load bearing ratio of threads gradually increases with the increase of a/P, while the load bearing ratio of the last three threads all de crease as the increase of a/P, thus making the load bearing ratio in threads more uneven. When b/l=1/4, a/P=1/4, the maximum load bearing ratio of threads decreases to 41 .70% for the first thread . Compared with the maximum load bearing ratio of the bolted joint without the groove, namely 74% as calculated before, the maximum load bearing ratio of threads of the bolted joint decreases by 4 3.65%. Therefore, the internal groove has great effects on the load distribution is threads of the bolted joint, by decre asing b/l and a/P, uneven load distribution in threads can be greatly improved. The structure of the bolted joint with an external is shown in Fig. 10. Four different cases are simulated for the functionally graded bolted joint with an external groove , (i) a/h=1/4, b/l=2/3, (ii) a/h=1/4, b/l=5/6, (iii) 10 a/h=1/2, b/l=2/3, (iv) a/h=1/2, b/l=5/6. Fig. 11 illustrates the load distribution in threads of the bolted joint with an external groove. It is clear from Fig. 1 1 that the maximum load bearing ratio of threa ds of the bolted joint with the external groove all have a decrease for the four cases , leading to more evenly distributed load in threads. When a/h=1/4, b/l=5/6, the maximum load bearing ratio of threads decreases to 30.21% for the first thread . Compared with the bolted joint without the groove, the maximum load bearing ratio of threads of the bolted joint decreases by 59.18%. Comparing the internal groove with the external groove, the nut with an external groove ( a/h=1/4, b/l=5/6) acts best performance to improve the unevenly distributed load in thread s and achieves the best service reliability of the bolted joints. Therefore, the shape of the nut has a great effect on the stress concentration of the bolted joint. B oth the internal and external groove can make the loading position of the nut stretches outward along the radial direction, which makes the nut bear the bending moment in an anticlockwise direction ( zM , based on the coordinate system shown in Fig. 2). Hence, the contact tightness between the first thread of the bolt and the nut decreases, while the contact tightness between the last thread of the bolt and the nut increases, resulting in the decrease of the maximum load bearing ratio of threads of the bolted joint which occurs at the first thread. It should be noted that the uncertainty in the material properties will have an impact on the optimization of the structure of the bolted joint, since the optim ization is quantitative analysis and sensitive. Further investigatement is needed to quantify the effects of these uncertainties on the optimization of the structure of the bolted joint. 5. Conclusions Thermodynamic behavior and load distribution in thread s of the Thermal Protection System (TPS) used bolted joints made up of porous ZrO 2/(ZrO 2+Ni) FGMs are numerically analyzed in detail by ABAQUS codes. From the numerical simulations some conclusion s can be drawn: (1) Unevenly distributed load in threads of th e bolted joint occurs at elevated temperature, load distribution in threads mainly focuses on the first two threads which are closest to the nut bearing surface, with the first thread carrying almost 7 4% of the total load . (2) Bolt-nut parameters have great effects on load distribution in threads, with trapezoidal thread, extra fine thread and greater modulus ratio of bolt to nut leading to more evenly distributed load in threads and better service reliability of the bolted joint. 11 (3) The shape of the nut has a significant effect on the load distribution in threads, the nut with an external groove ( a/h=1/4, b/l=5/6) can make the load bearing ratio of the first thread decrease to 30.21%, and thus the service reliability of the bolted joint is greatly improved. Acknowledgments The authors are grateful for the support by the National Natural Science Foundation of China (Grant nos. 11090330, 11090331, 11102003 and 11472038 ), the National Basic Research Program of China (973 Program, 2010CB83270 1, 2011CB610303), and t he Foundation of the Author of National Excellent Doctoral Dissertation of China (No. 201029) .