Analytical Modeling of Riveting Squeezing Force Considering Non-Uniform Deformation of Rivets in Aeronautical Structures

Analytical modeling of the squeezing force for aircraft wings and fuselage panels in the existing literature usually assumes uniform deformation of the rivets, while in reality, the deformation of the rivets is non-uniform. To achieve high-quality squeezing force modeling, this paper introduces Coulomb’s friction and four critical adjustments to the original equation: the non-uniform rivet/sheet interference along the sheet’s hole axial ordinate; the barreling effect when calculating the driven head’s volume; the spring-back of the driven head’s dimensions; the modified height of the driven head; and the modified sheet-hole expanded diameter considering the convex structure of the driven head. The calculated values of the proposed new model demonstrate an improved level of accuracy, particularly under squeeze ratios commonly encountered in the aerospace industry.


Introduction
Aircrafts are complex systems composed of various structures, differing in materials, shapes, and dimensions [1][2][3].Panels, created by riveting skin and stringers together, are widely used in aircraft [4,5], as illustrated in Figure 1a,b.In the riveting process, the squeezing force is one of the main setup parameters used to control the joining quality, mainly in automated riveting machines [6,7].During riveting, a rivet undergoes nonuniform deformations, such as the driven head becoming barrel shaped driven head or the shank becoming bell shaped, as shown in Figure 1c.The squeezing force (F sq ) is directly related to the driven head's dimensions (D i and H i ) and the non-uniform expanded rivet shank dimension (d i (z), including the interference between the rivet and the sheet hole) are depicted in Figure 1c.
Riveting involves five stages, detailed in Figure 1d, that directly influence the squeezing force, and are essential for aircrafts' structural strength and safety [2,6].Initially, the rivet is placed centrally in the hole between the bucking bar and punch.In stage 2, under the squeezing force, the rivet shank experiences elastic deformation, retaining a gap between the rivet shank and the sheet hole.By stage 3, plastic deformation is generated until the rivet comes into contact with the inner wall of the hole.This process involves axial compression and radial expansion of the rivet shank due to punch pressure.During stage 4, after the shank comes into contact with the inner wall of the hole, the material flow of the rivet shank inside the hole is radially constrained by the inner wall of the hole.As the squeezing force increases, the rivet shank outside the hole gradually forms a barrel-shaped driven head through friction.Simultaneously, the inside shank deforms from compression, with material from the driven head flowing into the hole.This results in non-uniform expansion of the inner wall of the hole, forming a bell-shaped rivet shank and a convex structure near the driven head.In stage 5, after maintaining the squeezing force for a certain period of time, the bucking bar and punch separate from the rivet, and the squeezing force is removed accordingly.Spring-back occurs in both the driven head and shank.
compression, with material from the driven head flowing into the hole.This results in non-uniform expansion of the inner wall of the hole, forming a bell-shaped rivet shank and a convex structure near the driven head.In stage 5, after maintaining the squeezing force for a certain period of time, the bucking bar and punch separate from the rivet, and the squeezing force is removed accordingly.Spring-back occurs in both the driven head and shank.In riveting, the dimensions of the driven head determine the connection strength [3][4][5], while the interference between the rivet and the sheet hole affects the fatigue performance.Furthermore, the squeezing force primarily causes diametral expansion of the hole, leading to riveting-induced deformation or material growth [6][7][8][9][10].Concurrently, with the aerospace industry's widespread use of automated and robot riveting, squeezing force has emerged as a crucial parameter for riveting quality.Müller [10] was the first scholar to systematically study the riveting process, discovering the relationship between squeezing force and the dimensions of the driven head.Schijve [11], through an analytical review of Müller's experimental data, discerned a latent correlation between the ratio of the driven head's diameter (D) to its height (H) and the squeezing force applied.De Rijck [12,13], based on Müller and Schijve's research, proposed a widely used model for calculating riveting squeezing force (as shown in Equations ( 1) and ( 2)) based on the assumption of the constant volume of the rivet shank.
( ) [ ] Figueira et al. [14,15] combined the previous studies by Müller [10], Schijve [11] and De Rijck [12,13] and introduced two factors (riveted grip thickness and friction among rivet, sheet, and punch tool) into the revised algebraic model shown in Equation (3), which is often more accurate under the medium squeeze ratio.In riveting, the dimensions of the driven head determine the connection strength [3][4][5], while the interference between the rivet and the sheet hole affects the fatigue performance.Furthermore, the squeezing force primarily causes diametral expansion of the hole, leading to riveting-induced deformation or material growth [6][7][8][9][10].Concurrently, with the aerospace industry's widespread use of automated and robot riveting, squeezing force has emerged as a crucial parameter for riveting quality.Müller [10] was the first scholar to systematically study the riveting process, discovering the relationship between squeezing force and the dimensions of the driven head.Schijve [11], through an analytical review of Müller's experimental data, discerned a latent correlation between the ratio of the driven head's diameter (D) to its height (H) and the squeezing force applied.De Rijck [12,13], based on Müller and Schijve's research, proposed a widely used model for calculating riveting squeezing force (as shown in Equations ( 1) and ( 2)) based on the assumption of the constant volume of the rivet shank.
Figueira et al. [14,15] combined the previous studies by Müller [10], Schijve [11] and De Rijck [12,13] and introduced two factors (riveted grip thickness and friction among rivet, sheet, and punch tool) into the revised algebraic model shown in Equation (3), which is often more accurate under the medium squeeze ratio.
Materials 2024, 17, 2756 3 of 27 Cheraghi [16] employed an adaptive mesh scheme, designed to counteract the nonlinear characteristics of rivet materials, to simulate the riveting process.This approach significantly reduces the amount of element distortion that occurs during the driven head's formation, thereby enhancing the accuracy of deformation estimations in the critical interface region between the sheets and the driven head.Chang et al. [17] accomplished a detailed study of riveting-induced deformation on wing panels generated by the rivet squeezing process and proposed a numerical simulation methodology based on FEM.The results obtained from the tests on a single curvature wing panel showed that there was good agreement between the experimental and simulated results.Chang et al. [18] developed a model for calculating the riveting squeezing force, which was constructed based on an analysis of the limit stress created when a thick-walled cylinder enters the plastic state.They also provided a method to solve model parameters based on the principle of a constant volume of the rivet shank.Furthermore, through numerical simulation, they validated the rationality of the division of stages in the riveting process.Zhang et al. [19] focuses on mathematically modeling and simulating the riveting process to enhance the precision of deformation analysis in thin-walled sheet metal parts, employing an ABAQUSbased calculation to analyze the stress and deformation of the model.Zeng et al. [20] derived a model of the riveting process to calculate squeezing force, and improved the model's accuracy by taking into account the flow of rivet material into plate holes.Their model incorporates a reduction in the volume fraction of the driven rivet head and is achieved through a reduction factor determined by numerical simulation.The models of the squeezing force in the existing literature usually assume uniform deformation of the rivet, while in reality, the deformation of the rivet is non-uniform, including the deformation of the driven head and the rivet shank within the hole.This leads to an overestimation or underestimation of the squeezing force because the calculated volume and the diameter of driven head are larger, and the diameter of the rivet shank is smaller.
This study aims to propose a revised model for calculating squeezing force using the principal stress method, addressing the rivet's non-uniform deformation.It redefines four parameter sets: the maximum expanded diameter of the sheet hole (d 2a ), the driven head's equivalent diameter accounting for barrel-shape (D eq ), the dimensions of the driven head after the spring-back (D 3 and H 3 ), and the modified dimensions of the expanded diameters of the driven head and sheet hole considering the convex structure of the driven head (H 4 and d 2b ).And the purpose of the research is to identify the most critical parameter and modify the calculation equation to minimize the deviation from the experimental results.

Squeeze Force Modeling Considering the Non-Uniform Deformation of the Rivet
Figueira et al. [15] demonstrated that Equations (1) and ( 2) have an overall good agreement with experimental or FEM results.There is only a slight tendency to obtain underestimated values for high squeeze ratios (D/D 0 ), which are generally above a ratio of 1.6; and in median ratios, from 1.3 up to 1.6, there is a slight tendency to obtain overestimated values, although that may not be exactly true in all cases (see Figure 2).To address this issue, Figueira et al. introduced Equation (3), which applies a strain correction for the rivet material at the hole boundaries and takes two factors (riveted grip thickness and friction among rivet, sheet, and punch tools) into consideration.However, experimental findings suggest that this correction, while helpful, is often insufficient, as the overestimation of squeezing force still occurs in many instances.Specimens A2 and A4 have the same thickness (t 1 = t 2 ≈ 2 mm) and use the same rivets (NAS1097D5-6), making them comparable.thickness (t1 = t2 ≈ 2 mm) and use the same rivets (NAS1097D5-6), making them comparable.

Figure 2.
The comparison of the results obtained by Figueira et al. [15] and the results obtained by De Rijck [12] (specimens A2 and A4).
The observed deviations can primarily be attributed to the fact that Equations ( 1)-( 3) are based on the assumption of uniform deformation in rivets, as shown in Figure 3a.Contrary to this assumption, the deformation of the driven head and the shank within the rivet hole is actually non-uniform due to the effects of friction and the constraints imposed by the rivet hole, as shown in Figure 3b.Consequently, the maximum diameter of the driven head (D) used in Equations ( 1)-( 3) inevitably leads to an overestimation of the squeezing force because the calculated volume of the driven head is larger than the actual barrel-shaped volume.Additionally, the integration area of stress is also larger than the actual one.Moreover, the uniform diameter (d2) of the rivet shank used in the Equation ( 2) is smaller than d2a (as shown in Figure 3a).This results in an overestimation of the size of the barrel-shaped portion of the driven head (zone Ⅱ), given by (D − d2) > (D − d2a), thereby leading to an overestimation of the calculated squeezing force.This is because the stress value ( iii σ ) in zone Ⅱ is higher than that ( ii σ ) in zone Ⅰ due to a greater effect of strain hardening outside the shank/hole interface area [15].
Furthermore, the spring-back of the driven head, which results in an increase in diameter and a decreased height post-deformation, contributes to a reduction in the squeezing force, a factor not accounted for in Equations ( 1)-( 3).The non-uniformity of  [15] and the results obtained by De Rijck [12] (specimens A2 and A4).
The observed deviations can primarily be attributed to the fact that Equations ( 1)-( 3) are based on the assumption of uniform deformation in rivets, as shown in Figure 3a.Contrary to this assumption, the deformation of the driven head and the shank within the rivet hole is actually non-uniform due to the effects of friction and the constraints imposed by the rivet hole, as shown in Figure 3b.Consequently, the maximum diameter of the driven head (D) used in Equations ( 1)-( 3) inevitably leads to an overestimation of the squeezing force because the calculated volume of the driven head is larger than the actual barrel-shaped volume.Additionally, the integration area of stress is also larger than the actual one.Moreover, the uniform diameter (d 2 ) of the rivet shank used in the Equation ( 2) is smaller than d 2a (as shown in Figure 3a).This results in an overestimation of the size of the barrel-shaped portion of the driven head (zone II), given by (D − d 2 ) > (D − d 2a ), thereby leading to an overestimation of the calculated squeezing force.This is because the stress value (σ iii ) in zone II is higher than that (σ ii ) in zone I due to a greater effect of strain hardening outside the shank/hole interface area [15].
which can make the calculation of squeezing force more precise.So, our study integrates the aforementioned four factors into a revised analytical model aimed at mitigating computational errors, irrespective of high, medium, or low squeeze ratios.In this section, four research tasks are delineated for execution: (1) the determination of rivet volume at various stages; (2) the squeezing force derivation using the principal stress method; (3) the calculation of the four new parameters; (4) the establishment of a new calculation model.

Determination of Rivet Volume at Various Stages
We assume that the volume of the rivet remains constant during stages 1 to 5 of the riveting process, as depicted in Figures 3-6.The key dimensions of the rivet, such as the height and diameter of the driven head, directly impact the calculation of squeezing force, which can be obtained by calculating the volume of the rivet shank at each stage (throughout the entire calculation process, the riveting structure is consistently assumed to be symmetrical).
The modeling process requires the adoption of a polar coordinate system (z, r, θ), positioned on the bottom of manufactured head, as shown in Figure 4. Within this system, the z-axis corresponds to the rivet shank, and its origin is located on the upper surface of the sheet 1.The model initially assumes an axisymmetric geometry for the whole model Furthermore, the spring-back of the driven head, which results in an increase in diameter and a decreased height post-deformation, contributes to a reduction in the squeezing force, a factor not accounted for in Equations ( 1)-(3).The non-uniformity of interference leads to maximum interference on one side of the contact between the rivet hole and the driven head.In this region, the actual strain is also smaller; hence, the predicted squeezing force should be correspondingly smaller.And the spring-back of the driven head results in an increase in height and a decrease in diameter, while the material of the driven head inserted into the rivet hole leads to the formation of a convex structure near the driven head.These are factors that are not accounted for in Equations ( 1)-(3) which can make the calculation of squeezing force more precise.So, our study integrates the aforementioned four factors into a revised analytical model aimed at mitigating computational errors, irrespective of high, medium, or low squeeze ratios.In this section, four research tasks are delineated for execution: (1) the determination of rivet volume at various stages; (2) the squeezing force derivation using the principal stress method; (3) the calculation of the four new parameters; (4) the establishment of a new calculation model.

Determination of Rivet Volume at Various Stages
We assume that the volume of the rivet remains constant during stages 1 to 5 of the riveting process, as depicted in Figures 3-6.The key dimensions of the rivet, such as the height and diameter of the driven head, directly impact the calculation of squeezing force, which can be obtained by calculating the volume of the rivet shank at each stage (throughout the entire calculation process, the riveting structure is consistently assumed to be symmetrical).
(rivet and the sheets), the volume of rivet is assumed to be constant and the thicknesses of the sheets (t1 and t2) are assumed to be kept constant.The subscript of the rivet volume V is standardized according to the following criteria: Arabic numerals are employed to denote the different stages of the riveting process, whereas Roman numerals are used to indicate the various components of the rivet.In riveting stage 1, the volume of rivet shank V0 is given by Between riveting stage 2 and riveting stage 3, under the influence of the squeezing force, the diameter of the rivet shank expands until it matches the diameter of the hole [10,11], denoted as D1.Given that the volumes in both the second and third stages are equivalent, denoted as V1, a solution is derived by concurrently solving Equations ( 4) and ( 5), leading to Equation (6).Equation ( 6) enables the calculation of the length of the rivet shank in stage 3, designated as L1, Ⅰ represents the volume within the hole and Ⅱ represents the volume of the driven head under punching pressure, as shown in Figure 5.
In stage 4, as illustrated in Figure 6, upon contact between the rivet and the hole wall, frictional resistance occurs at the part where the bottom of the rivet shank comes into contact with the punch, the portion of the rivet shank outside the hole gradually assumes a "barrel-like shape", while the rivet shank inside the hole develops a non-uniform deformation, and assumes a "bell-like shape".The volume of the rivet is denoted as V2 in stage 4, and considering that di(z) is the non-uniform diameter of the rivet shank in the hole along the z axis, the volume of the rivet shank is VⅠ, and the volume of driven head is VⅡ (which will be calculated in Section 2.3.1),then it can be found out that  And the diameter of the rivet shank di(z) is no longer constant; it is slightly larger on the driven head side at any coordinate z.In Stage 5, when the final dimensions of the driven head (D and H) are reached, the squeezing force is withdrawn.Following a brief period of spring-back in both the rivet shank and the sheet hole diameters, alterations in the dimensions of the rivet ensue, Ⅲ represents the volume within the hole and Ⅳ represents the volume of the driven head undergoes spring-back, as shown in Figure 7.The volume of the rivet is denoted as V3 in stage 5, the non-uniform dimensions of rivet shank diameter is still denoted as di(z), and the non-uniform driven head diameter considering the barrel-shape, Di(z); V3 can be expressed using Equation ( 9).And the diameter of the rivet shank di(z) is no longer constant; it is slightly larger on the driven head side at any coordinate z.In Stage 5, when the final dimensions of the driven head (D and H) are reached, the squeezing force is withdrawn.Following a brief period of spring-back in both the rivet shank and the sheet hole diameters, alterations in the dimensions of the rivet ensue, Ⅲ represents the volume within the hole and Ⅳ represents the volume of the driven head undergoes spring-back, as shown in Figure 7.The volume of the rivet is denoted as V3 in stage 5, the non-uniform dimensions of rivet shank diameter is still denoted as di(z), and the non-uniform driven head diameter considering the barrel-shape, Di(z); V3 can be expressed using Equation ( 9).The modeling process requires the adoption of a polar coordinate system (z, r, θ), positioned on the bottom of manufactured head, as shown in Figure 4. Within this system, the z-axis corresponds to the rivet shank, and its origin is located on the upper surface of the sheet 1.The model initially assumes an axisymmetric geometry for the whole model (rivet and the sheets), the volume of rivet is assumed to be constant and the thicknesses of the sheets (t 1 and t 2 ) are assumed to be kept constant.The subscript of the rivet volume V is standardized according to the following criteria: Arabic numerals are employed to denote the different stages of the riveting process, whereas Roman numerals are used to indicate the various components of the rivet.
In riveting stage 1, the volume of rivet shank V 0 is given by Between riveting stage 2 and riveting stage 3, under the influence of the squeezing force, the diameter of the rivet shank expands until it matches the diameter of the hole [10,11], denoted as D 1 .Given that the volumes in both the second and third stages are equivalent, denoted as V 1 , a solution is derived by concurrently solving Equations ( 4) and ( 5), leading to Equation (6).Equation ( 6) enables the calculation of the length of the rivet shank in stage 3, designated as L 1 , I represents the volume within the hole and II represents the volume of the driven head under punching pressure, as shown in Figure 5.
In stage 4, as illustrated in Figure 6, upon contact between the rivet and the hole wall, frictional resistance occurs at the part where the bottom of the rivet shank comes into contact with the punch, the portion of the rivet shank outside the hole gradually assumes a "barrel-like shape", while the rivet shank inside the hole develops a non-uniform deformation, and assumes a "bell-like shape".The volume of the rivet is denoted as V 2 in stage 4, and considering that d i (z) is the non-uniform diameter of the rivet shank in the hole along the z axis, the volume of the rivet shank is V I , and the volume of driven head is V II (which will be calculated in Section 2.3.1),then it can be found out that And the diameter of the rivet shank d i (z) is no longer constant; it is slightly larger on the driven head side at any coordinate z.
In Stage 5, when the final dimensions of the driven head (D and H) are reached, the squeezing force is withdrawn.Following a brief period of spring-back in both the rivet shank and the sheet hole diameters, alterations in the dimensions of the rivet ensue, III represents the volume within the hole and IV represents the volume of the driven head undergoes spring-back, as shown in Figure 7.The volume of the rivet is denoted as V 3 in stage 5, the non-uniform dimensions of rivet shank diameter is still denoted as d i (z), and the non-uniform driven head diameter considering the barrel-shape, D i (z); V 3 can be expressed using Equation (9).
Meanwhile, it is clear that:

Riveting Squeezing Force Modeling
In stage 4, just before the squeezing force removal, the squeezing force Fsq can be calculated by: As shown in Figure 8, it is clear that the deformation and hardening levels are different in Zone Ⅰ and Zone Ⅱ [21].Based on the assumption of the uniform deformation of the rivet as shown Figure 8a, and introducing Coulomb's friction condition [22,23] and the Hollomon relationship [24], the power law expression for the flow curve (e.g., Equation ( 17)) can be used as an approximation to predict the plastic stress-strain behavior under riveting conditions.Figueira [15] provided an adapted equation for each respective Zone's strain condition (σi and σiii in Zones Ib and II, respectively, as shown in Figure 8a): In Figueira's study [15], it is assumed that the interference fit during the riveting process is uniform.According to this assumption, the diameter of the rivet shank within the hole remains constant, and therefore, as shown in Figure 3a, d 2 is commonly considered to be equivalent to the diameter of the rivet shank when it is pressed into the hole, based on the principle of volume constancy.And, as shown in Figure 3, we have Therefore, a further derivation can be obtained: Meanwhile, it is clear that:

Riveting Squeezing Force Modeling
In stage 4, just before the squeezing force removal, the squeezing force F sq can be calculated by: As shown in Figure 8, it is clear that the deformation and hardening levels are different in Zone I and Zone II [21].Based on the assumption of the uniform deformation of the rivet as shown Figure 8a, and introducing Coulomb's friction condition [22,23] and the Hollomon relationship [24], the power law expression for the flow curve (e.g., Equation ( 17)) can be used as an approximation to predict the plastic stress-strain behavior under riveting conditions.Figueira [15] provided an adapted equation for each respective Zone's strain condition (σ i and σ iii in Zones I b and II, respectively, as shown in Figure 8a): Based on the assumption of the non-uniform deformation of the rivet, as shown Figure 8b, Equation (17) needs to be rewritten as: In which Based on the assumption of the non-uniform deformation of the rivet, as shown Figure 8b, Equation (17) needs to be rewritten as: In which ε sq lb represent the effective plastic strain of the portion of the driven head with a diameter of d2b, d2a, d2 under different interference distributions; ( ) ε Ⅱ sq represents the effective plastic strain of the remaining portion of the driven head.
Refering to Figure 8b, the squeezing force can be calculated by By simultaneously combining Equations ( 16)-( 20), the overall expression Equation ( 21) can be derived in which  Refering to Figure 8b, the squeezing force can be calculated by By simultaneously combining Equations ( 16)-( 20), the overall expression Equation ( 21) can be derived in which Considering the convex structure of the driven head shown in Figure 8c, the derivation process is the same as Equation (19), leading us to another model: in which There are some parameters in Equations ( 21) and ( 22) that need to be determined, including: d 2a and D eq , d 2b and ∆H, the amount of the spring-back ∆ sbr and ∆ sbz .
To determine the abovementioned unknown parameters, there are four aspects of the research that need to be addressed in Section 2.3: (1) consider the influence of nonuniform interference fit; (2) obtain the equivalent diameter of driven head; (3) calculate the dimensions of driven head after the spring-back of the rivet; (4) modify the dimensions of driven head, and the expanded sheet hole diameter considering the convex structure of the driven head.This allows us to revise the calculation model for squeezing force, thereby enhancing the overall accuracy of the squeezing force calculation.During the formation of the rivet head, as the riveting punch presses down, the rivet material outside the hole undergoes axial compression and radial expansion.This process is accompanied by friction (f ) between the contact surfaces, resulting in a barrel-shaped contour on the edges of the compressed rivet head [24][25][26][27][28].
As shown in Figure 9, in the barrel-shaped driven head, the key parameters include the maximum diameter D 2 (which can be measured), the height of the driven head H 2 (typically measurable), and the minimum diameter d 3 (which needs to be calculated).increment for the ith pressing; and Δzi is the axial displacement increment for the ith pressing.The specific calculation and schematic process are illustrated in Figure 11.

Final state
The iterations for i 1、i and i+1 The process of the formation of a rivet driven head is similar to the upsetting process of a cylindrical workpiece, as shown in Figure 9. Therefore, under the same dimensional conditions, we equate the solution for d 3 in the riveting process to the solution for d 3 in the forging process: To obtain d ′ 3 , Sun [29] used the energy method to theoretically model the upsetting process of a cylindrical workpiece.They quantitatively derived the barrel-shaped contour of the deformed cylindrical metal material after upsetting through displacement increment expressions, as indicated in Equations ( 23)- (28).The initial calculation conditions are the coordinates (r 0 , z 0 ) of any arbitrary point A on the cylindrical metal material and the final pressing amount ∆H.The entire pressing process is divided into n computational steps, as shown in Figure 10.

Initial state
The iterations for i−1、i and i+1 Using A (r i−1 , z i−1 ) as the input, the coordinates A (r i , z i ) of any point on the shape's contour after the ith pressing are determined.By iterating this process, the radial displacement ∆r and axial displacement ∆z at any point on the barrel-shaped contour can be obtained, ultimately yielding the barrel-shaped contour after the upsetting deformation.In the formula, B is an undetermined coefficient related to friction; H 0 is the initial height of the cylindrical metal material; n is the number of iterations; H i−1 is the calculated height of the cylindrical metal material after the (i − 1)th pressing; ∆r i is the radial displacement increment for the ith pressing; and ∆z i is the axial displacement increment for the ith pressing.The specific calculation and schematic process are illustrated in Figure 11.

Final state
The iterations for i−1、i and i+1 In these equations, B represents an undetermined coefficient that is associated with friction.
where r 0 defines the location on the upper surface where the maximum frictional shear stress (τ max ) occurs [30,31], which refers to the absolute distance between this point and the center of the upper surface.
By solving Equations ( 25)-( 32) simultaneously, we can determine the values of r n and z n .Then, we can use Equations ( 24) and ( 26) obtain the diameter d 3 of the driven head.
As shown in Figure 12, to obtain an accurate solution to V II , we assume that d 3 = d 5 .Liang et al. [32] simplified the driven head riveting process into a free upsetting procedure.Based on this simplification, a mathematical model describing the stress field of the driven head during the deformation process was proposed, followed by derivation of the driven head contour equation.Through expansion of this equation using the Taylor series, it was concluded that the external contour of the driven head model could be effectively fitted with a univariate quadratic equation. .
By performing an integral calculation on the outer contour of the driven head, we derived the volume of driven head V Ⅳ : considering Equations ( 31)-( 33), V Ⅳ becomes: 3 ) 60 Then, according to the principle of constant volume, the equivalent diameter of the driven head can be calculated at the same height:  In the ZOr coordinate system, by defining three key points, A, B, and C, on the outer contour of the driven head as shown in Figure 12, and by substituting the coordinate values of these points into Equation (33), the following results were obtained: By performing an integral calculation on the outer contour of the driven head, we derived the volume of driven head V IV : considering Equations ( 31)-( 33), V IV becomes: Materials 2024, 17, 2756 13 of 27 Then, according to the principle of constant volume, the equivalent diameter of the driven head can be calculated at the same height: Rans and Li [33,34] used a finite element simulation to derive the typical distribution trends of rivet interference fit.They found that the radial expansion of the hole edge on sheet 1 remains relatively constant and easy to obtain: (39) So, the total volume of the rivet shank (V I I I ) is assumed to consist of a cylinder (V I I I a ) and a frustum of a cone (V I I I b ), as shown in Figure 13.

The Expanded Hole Diameter d2a Considering the Non-Uniform Interference Fit
Rans and Li [33,34] used a finite element simulation to derive the typical distribution trends of rivet interference fit.They found that the radial expansion of the hole edge on sheet 1 remains relatively constant and easy to obtain: So, the total volume of the rivet shank ( V Ⅲ ) is assumed to consist of a cylinder ( a V Ⅲ ) and a frustum of a cone ( b V Ⅲ ), as shown in Figure 13.
By combining Equations ( 37) and (40) through Equation ( 43), we can obtain d2a: (44) . d2a at the interface between the driven head and the lower surface of sheet 2.

The Influence of the Material Inserted into the Hole on the Dimensions of the Driven Head
In Equation (45), it is assumed that the volume of the rivet shank inside the hole remains unchanged [20]; however, in contrast, Equation (46) acknowledges that a portion of the driven head material is still squeezed into the hole, forming a convex structure with a height of ΔH, as depicted in Figure 14.
When the rivet is in contact with the hole wall, the volume of the material outside the hole is equal to the volume after the formation of the upsetting head; as such:

The Influence of the Material Inserted into the Hole on the Dimensions of the Driven Head
In Equation (45), it is assumed that the volume of the rivet shank inside the hole remains unchanged [20]; however, in contrast, Equation (46) acknowledges that a portion of the driven head material is still squeezed into the hole, forming a convex structure with a height of ∆H, as depicted in Figure 14. 42 The solutions for ΔH and d2b can be obtained through iterative calculations.After the release of the squeezing force, the interference fit within the plate and the driven head will release their stored elastic potential energy, thereby undergoing springback, leading to a certain degree of volume deformation within the material, and resulting in a change in the dimensions of the driven head (D and H) (as shown in Figure 15).In stage 5, the portion of the material flowing into the rivet hole will undergo spring-back and return to the outside of the rivet hole, resulting in an increase in the volume of the driven head during this stage.In the illustration, the red area represents the diameter and height of the driven head in riveting stage 4, while the blue area represents the diameter and height of the rivet after the spring-back.The calculation of the average stress can be achieved by solving Equation (48).When the rivet is in contact with the hole wall, the volume of the material outside the hole is equal to the volume after the formation of the upsetting head; as such: The solutions for ∆H and d 2b can be obtained through iterative calculations.

The Spring-Back and Influence of Rivet after the Pressure Has Been Removed
During riveting stage 5, after the squeezing force has been removed, spring-back, ∆ sbr and ∆ sbz , occurs in the driven head and shank.
After the release of the squeezing force, the interference fit within the plate and the driven head will release their stored elastic potential energy, thereby undergoing springback, leading to a certain degree of volume deformation within the material, and resulting in a change in the dimensions of the driven head (D and H) (as shown in Figure 15).In stage 5, the portion of the material flowing into the rivet hole will undergo spring-back and return to the outside of the rivet hole, resulting in an increase in the volume of the driven head during this stage.In the illustration, the red area represents the diameter and height of the driven head in riveting stage 4, while the blue area represents the diameter and height of the rivet after the spring-back.The calculation of the average stress can be achieved by solving Equation (48).
Therefore, the driven head size considering the spring-back is Therefore, the driven head size considering the spring-back is

Calculation Results and Discussion
The riveting dimensions and material parameters used for analytical calculations in this paper are obtained from the following sources: the rivet parameters are derived from the relevant parameters of specimen a2 in Müller's experimental sample [10], while the sheet material parameters are taken from MMPDS 07 [35] (as shown in Table 1).Additionally, the riveting dimensions and material parameters used for the calculations also come from De Rijck's experimental samples A1, A2, A3, A4, and A13 [12] (as shown in Table 2).
Through experimentation, the following investigations will be carried out: (1) Importance analysis of the proposed parameters including d3, Deq, d2a, and SP in squeezing force calculation, aiming to identify the factor exerting the most significant influence on squeezing force calculation.(2) Considering the non-uniform deformation of the rivet, comparisons need to be made among the calculation results obtained from Equation ( 21), the experimental results obtained by Müller [10] and De Rijck [12], and the calculation results obtained from Equation (3) (as shown in Table 3).(3) Considering the convex structure of the driven head, comparisons need to be made among the calculation results obtained from Equation ( 22), the experimental results obtained by Müller [10] and De Rijck [12], and the calculation results obtained from Equation (3) (as shown in Table 3).

Calculation Results and Discussion
The riveting dimensions and material parameters used for analytical calculations in this paper are obtained from the following sources: the rivet parameters are derived from the relevant parameters of specimen a2 in Müller's experimental sample [10], while the sheet material parameters are taken from MMPDS 07 [35] (as shown in Table 1).Additionally, the riveting dimensions and material parameters used for the calculations also come from De Rijck's experimental samples A1, A2, A3, A4, and A13 [12] (as shown in Table 2).
Through experimentation, the following investigations will be carried out: (1) Importance analysis of the proposed parameters including d 3 , D eq , d 2a , and S P in squeezing force calculation, aiming to identify the factor exerting the most significant influence on squeezing force calculation.(2) Considering the non-uniform deformation of the rivet, comparisons need to be made among the calculation results obtained from Equation ( 21), the experimental results obtained by Müller [10] and De Rijck [12], and the calculation results obtained from Equation (3) (as shown in Table 3).(3) Considering the convex structure of the driven head, comparisons need to be made among the calculation results obtained from Equation ( 22), the experimental results obtained by Müller [10] and De Rijck [12], and the calculation results obtained from Equation (3) (as shown in Table 3).[12] and De Rijck et al. [13] experimentally measured data or assumed by default.

Data Value
Coulomb  3. Comparative of Müller's [10] and De Rijck's [13] experimental data and the analytical calculations based on the parameters listed in Tables 1 and 2.

Analysis of Parameters' Importance
To analyze the importance of d 3 , D eq , d 2a , and S P in the squeezing force calculation equation, we selected specimens A2 and A4 from Table 2. Firstly, we extracted the riveting squeezing force measurement data for specimens A2 and A4 from Table 3.Then, we extracted the calculated values obtained through Equation (3), depicted, respectively, by red dashed lines and blue dashed lines in Figures 16 and 17.
In Equation ( 3), by replacing D with d 3 (obtained from Equation ( 24)), as shown in Equation (55), the calculated riveting squeezing force could be plotted as red dots; using D eq (obtained from Equation (38)) instead of D, as shown in Equation (56), the resulting riveting squeezing force was plotted as blue triangles; substituting d 2a (obtained from Equation (44)) for d 2 , as shown in Equation (57), the resulting riveting squeezing force could be plotted as black circles; replacing H with H 3 (obtained from Equation ( 53)) and D with D 3 (obtained from Equation (54)), as shown in Equation (58), the resulting riveting squeezing force was plotted as mauve squares.

Considering Non-Uniform Deformation of the Rivet
Taking into account the non-uniform deformation of the rivet, the calculated squeezing force values obtained from Equation (21) are shown in Table 4. Figure 18 is plotted to compare Equations ( 1), ( 3) and (21), and the experimental data, aiming to analyze the accuracy of the predicted squeezing force calculated from Equation (21), as shown in Table 4.
The overall trend in the squeezing force obtained from Equation ( 21) is consistent with the experimental values, but in most cases, it slightly underestimates the experimental results.The standard deviations for specimens a2, A1, A2, A3, A4, and A13 are 0.53%, 2.24%, 3.73%, 17.66%, 7.6%, and 11.16%, respectively.In Figure 18a, for specimen A2, within the medium squeeze ratio range (1.4 ≤ D/D0 ≤ 1.7), the predicted values from Equation ( 21) are more accurate compared to Equations ( 1) and ( 3), with a deviation of less than 5.28%.For specimen A4 (see Figure 18e), within the medium and high squeeze ratio range (1.39 ≤ D/D0), the accuracy of predictions from Equation ( 21) is higher compared to Equations ( 1) and (3), with a deviation from the experimental values of 7.03%.For specimens A1 (see Figure 18b) and A13 (see Figure 18f), overall, the predictions from Equation (21) are lower than the experimental values, but the deviation is still within 13.02%.Particularly for specimen A2 (see Figure 18c) at low to medium squeeze ratios (D/D0 < 1.7), and specimen A3 (see Figure 19d) at medium and high squeeze ratios (1.4 ≤ D/D0), the accuracy of predictions from Equation ( 21) is higher compared to Equations ( 1) and (3).The maximum deviation between predicted and experimental values is less than 4%, with deviations of 3.59% and 2.97%, respectively.In Figures 16 17, we must find out the major factors influencing the accuracy of the calculation (MFAC), and the major factors influencing the calculation results (MIFC).
Comparing with the experimental data and the calculated results, that the following can be seen: (1) For A2 and A4, in the majority of cases, substituting d 3 or D eq with D, or d 2a with d 2 , or applying S P to adjust the dimensions D and H of the driven head, resulting in a squeezing force calculation that is smaller than the computed value.(2) In Figure 16, d 3 emerges as the major influencing factor for the calculation results (MIFC), while D eq is the major factor influencing the accuracy of the calculation (MFAC).(3) In Figure 17, d 3 remains the MIFC, whereas d 2 appears to be the MFAC.(4) Additionally, it is easy to observe from Figures 16 and 17 that, relative to the experimental values, the riveting force calculated using D eq is more accurate than that obtained using d 3 .Therefore, in the following calculation process, D eq is chosen as the calculation region for σ iii .

Considering Non-Uniform Deformation of the Rivet
Taking into account the non-uniform deformation of the rivet, the calculated squeezing force values obtained from Equation ( 21) are shown in Table 4. Figure 18 is plotted to compare Equations ( 1), ( 3) and ( 21), and the experimental data, aiming to analyze the accuracy of the predicted squeezing force calculated from Equation ( 21), as shown in Table 4.This is the case because Equations ( 1) and (3) use the maximum rivet head diameter D, which results in an overestimation due to the larger volume of the rivet head.In contrast, the predictions from Equation (21) adequately consider the influence of the barrelshaped of the rivet head on the calculation results, resulting in an improved accuracy.The overall trend in the squeezing force obtained from Equation ( 21) is consistent with the experimental values, but in most cases, it slightly underestimates the experimental results.The standard deviations for specimens a2, A1, A2, A3, A4, and A13 are 0.53%, 2.24%, 3.73%, 17.66%, 7.6%, and 11.16%, respectively.In Figure 18a, for specimen A2, within the medium squeeze ratio range (1.4 ≤ D/D 0 ≤ 1.7), the predicted values from Equation ( 21) are more accurate compared to Equations ( 1) and ( 3), with a deviation of less than 5.28%.For specimen A4 (see Figure 18e), within the medium and high squeeze ratio range (1.39 ≤ D/D 0 ), the accuracy of predictions from Equation ( 21) is higher compared to Equations ( 1) and ( 3), with a deviation from the experimental values of 7.03%.For specimens A1 (see Figure 18b) and A13 (see Figure 18f), overall, the predictions from Equation (21) are lower than the experimental values, but the deviation is still within 13.02%.Particularly for specimen A2 (see Figure 18c) at low to medium squeeze ratios (D/D 0 < 1.7), and specimen A3 (see Figure 19d) at medium and high squeeze ratios (1.4 ≤ D/D 0 ), the accuracy of predictions from Equation ( 21) is higher compared to Equations ( 1) and (3).The maximum deviation between predicted and experimental values is less than 4%, with deviations of 3.59% and 2.97%, respectively.

Considering the effect of the Material Inserted into the Hole on the Height of the Driven Head
Considering the convex structure of the driven head, the squeezing force values calculated using Equation ( 22) are shown in Table 4. Subsequently, Figure 19 was be plotted  This is the case because Equations ( 1) and (3) use the maximum rivet head diameter D, which results in an overestimation due to the larger volume of the rivet head.In contrast, the predictions from Equation ( 21) adequately consider the influence of the barrel-shaped of the rivet head on the calculation results, resulting in an improved accuracy.

Considering the Effect of the Material Inserted into the Hole on the Height of the Driven Head
Considering the convex structure of the driven head, the squeezing force values calculated using Equation ( 22) are shown in Table 4. Subsequently, Figure 19 was be plotted to compare Equations ( 1), ( 3) and ( 22) with the experimental data, aiming to analyze the accuracy of the squeezing force predicted by Equation (22) as shown in Table 4.The overall trend in the squeeze force obtained from Equation ( 22) is consistent with the experimental values.The standard deviations for specimens a2, A1, A2, A3, A4, and A13 are 8.46%, 3.32%, 1.32%, 6.87%, 4.27%, and 1.34%, respectively.For specimen A1 (see Figure 19b), under medium squeeze ratios (1.3 ≤ D/D 0 ≤ 1.62), Equation (22) shows a significant improvement in the overall accuracy, with deviations less than 3.77%.For specimen A13 (see Figure 19f), the prediction accuracy of Equation ( 22) has also improved, with deviations below 5.45% under medium-to-high squeeze ratios (1.4 ≤ D/D 0 ≤ 1.74).For specimen a2 (see Figure 19a), under medium-to-high squeeze ratios (1.27 ≤ D/D 0 ), the deviation between the calculated values from Equation (22) and the experimental values is within 5.98%.For specimen A2 (see Figure 19e), at low squeeze ratios (D/D 0 ≤ 1.3), the accuracy of Equation ( 22) is higher compared to that of Equations ( 1) and ( 3), with deviations from the experimental values of 13.49%.
This improvement is attributed to considering the material inserted into the rivet hole when correcting the driven head size, resulting in a closer match to experimental results.

Conclusions
A new analytical model for calculating the riveting squeezing force has been proposed, taking into account the non-uniform deformation of the rivet, rather than the uniform deformation that was previously used in the literature, in which four sets of parameters are redefined, including maximum expanded sheet hole diameter (d 2a ); the equivalent diameter of the driven head considering the barrel shape (D eq ); the dimensions of driven head after spring-back (D 3 and H 3 ); the modified height of driven head (H 4 ); and the modified expanded sheet hole diameter (d 2b ) after considering the convex structure of the driven head.
The analysis leads to the following points: (1) The introduction of both d 2a and D eq can reduce the calculated values of the riveting force, while the introduction of spring-back can slightly increase the calculated values of the riveting force compared to the previous model, as shown in Equation ( 3).
(2) The new analytical models for calculating riveting squeezing force, represented by Equations ( 21) and (22), which take into account the non-uniform deformation of the rivet, significantly reduced the deviation between the calculated squeezing force and experimental data, especially in the medium squeeze ratio range (1.3 ≤ D/D 0 ≤ 1.6).The previous formulation (Equations ( 1) and ( 3)) generally overestimates the required squeezing force due to an overestimation of the true strain of the rivet material within the hole vicinity.Equations ( 21) and ( 22) significantly improved the prediction accuracy, and the deviation from the experimental values is reduced to within 10% for some specimens (while the deviation for Equations ( 1) and (3) sometimes exceed 20%), this is because the non-uniform deformation of the rivet is considered fully.(3) After considering the convex structure of the driven head, the calculation accuracy is significantly improved, and most importantly, the consistency of the predicted squeezing force is improved, with a standard deviation of 8.46%.However, Equation ( 22) had a standard deviation of 17.66%.(4) For other types of sheet materials, such as CFRP, and different fasteners like headless rivets, the proposed model requires modifications to maintain accuracy.And, the current study did not account for the friction between the rivet shank and the hole.
To further the research in this area, the proposed model must take into account the friction that occurs between the rivet shank and the holes.This improvement would help reduce the tendency to underestimate the riveting forces in the low squeeze ratio range.

Figure 1 .
Figure 1.Typical riveted aircraft structure and the typical riveting process.

Figure 1 .
Figure 1.Typical riveted aircraft structure and the typical riveting process.

Figure 2 .
Figure 2. The comparison of the results obtained by Figueira et al.[15] and the results obtained by De Rijck[12] (specimens A2 and A4).

Figure 3 .
Figure 3. Four factors affecting the non-uniform deformation of the rivet.

Figure 3 .
Figure 3. Four factors affecting the non-uniform deformation of the rivet.

Figure 4 .
Figure 4. Stage 1: initial dimensions of the rivet and the polar reference coordinate system.

Figure 4 .
Figure 4. Stage 1: initial dimensions of the rivet and the polar reference coordinate system.

Figure 5 .
Figure 5. Stage 2-stage 3: the rivet shank touches the sheet hole wall under the squeezing force.

Figure 6 .
Figure 6.Stage 4: formation of the driven head under punching pressure. sq

FFigure 5 .
Figure 5. Stage 2-stage 3: the rivet shank touches the sheet hole wall under the squeezing force.

Figure 5 .
Figure 5. Stage 2-stage 3: the rivet shank touches the sheet hole wall under the squeezing force.

Figure 6 .
Figure 6.Stage 4: formation of the driven head under punching pressure. sq

Figure 7 .
Figure 7. Stage 5: the squeezing force is removed; the driven head undergoes spring-back.

Figure 7 .
Figure 7. Stage 5: the squeezing force is removed; the driven head undergoes spring-back.

Figure 8 .
Figure 8. Axial stresses at rivet heads and the interface stress σ i between Zones (I a ) and (I b ) and σ iii (II) (only the axial stresses).

2. 3 .
The Determination of the Unknown Parameters 2.3.1.Estimation of D eq Considering the Barrel-Shaped Effect

Figure 10 .
Figure 10.The contour change diagram of the cylindrical workpiece during upsetting process.

Figure 11 .Figure 9 .
Figure 11.Schematic diagram of iterative solution method for contours of the cylindrical workpiece.

Materials 2024 ,
17, x FOR PEER REVIEW 11 of 29increment for the ith pressing; and Δzi is the axial displacement increment for the ith pressing.The specific calculation and schematic process are illustrated in Figure11.

Figure 10 .
Figure 10.The contour change diagram of the cylindrical workpiece during upsetting process.

Figure 10 .
Figure 10.The contour change diagram of the cylindrical workpiece during upsetting process.

Figure 10 .
Figure 10.The contour change diagram of the cylindrical workpiece during upsetting process.

Figure 11 .Figure 11 .
Figure 11.Schematic diagram of iterative solution method for contours of the cylindrical workpiece.

Figure 12 .
Figure 12.Equating the barrel-shaped volume of the driven head to the volume of a cylinder (VIV = VV).

Figure 12 .
Figure 12.Equating the barrel-shaped volume of the driven head to the volume of a cylinder (V IV = V V ).

Figure 13 .
Figure 13.d 2a at the interface between the driven head and the lower surface of sheet 2.

Figure 14 .
Figure 14.The diameter and height of the convex structure of the driven head material entering the hole.2.3.4.The Spring-Back and Influence of Rivet after the Pressure has been Removed During riveting stage 5, after the squeezing force has been removed, spring-back, sbr

Figure 14 .
Figure 14.The diameter and height of the convex structure of the driven head material entering the hole.

Figure 15 .
Figure 15.Fifth stage of riveting: the driven head final dimensions (D3 and H3) after the driven head undergoes the spring-back.

Figure 15 .
Figure 15.Fifth stage of riveting: the driven head final dimensions (D 3 and H 3 ) after the driven head undergoes the spring-back.

Figure 16 .
Figure 16.The individual effects of each factor on the squeezing force of specimen A2.

Figure 16 . 29 Figure 16 .
Figure 16.The individual effects of each factor on the squeezing force of specimen A2.

Figure 17 .
Figure 17.The individual effects of each factor on the squeezing force of specimen A4.

Figure 17 .
Figure 17.The individual effects of each factor on the squeezing force of specimen A4.

Figure 18 .
Figure 18.The calculated results based on the combined factors of non-uniform interference d2a and equivalent driven head diameter Deq, and the spring-back.

Figure 18 .
Figure 18.The calculated results based on the combined factors of non-uniform interference d 2a and equivalent driven head diameter D eq , and the spring-back.

Materials 2024 , 29 Figure 19 .
Figure 19.The calculated results based on the combined factors of non-uniform interference d2b, the height of convex structure ΔH, and the spring-back.

Figure 19 .
Figure 19.The calculated results based on the combined factors of non-uniform interference d 2b , the height of convex structure ∆H, and the spring-back.

Table 2 .
Rivet material parameters used in analytical calculations based on De Rijck's

Table 4 .
Results calculated by four improving factors.
L 3 the formed length of rivet shank after elastic recovery d 1 the diameter of original sheet-hole (before rivet shank expansion) d 2 the diameter of expanded sheet-hole considering the uniform deformation d 2a the diameter of the bottom of expanded sheet-hole considering the non-uniform deformation d 2b the diameter of the bottom of expanded sheet-hole considering convex structure of the