Abstract

It is essential to consider microscopic contact deformation on the joint surface when the dynamic properties of flange connection structure were characterized. In the study, an equivalent modeling method of virtual material based on microscopic contact experiment was proposed. First, the macroscopic model of flange connection structure was established and assembly process was simulated. Several subregions were divided according to the pressure distribution of joint surface. In each subregion, the pressure distribution was approximately uniform. Second, the microtopography of joint surface was measured experimentally, and the fractal parameters were calculated by power spectral density function method. At the same time, the contact mechanics experiment by industrial CT equipment was conducted to accurately measure the actual contact state of joint surface, and the mathematical model of contact pressures and ratios of actual contact area to nominal contact area was established. Next, the contact stiffness of joint surface was calculated based on fractal contact theory and the actual contact area measured by industrial CT equipment. Finally, the parameters of equivalent virtual material were calculated by weighted average of each subregion. In addition, a modal experiment of flange connection structure was conducted to compare with the simulation result. The result shows that the average relative error of natural frequency was 4.69% between experiment and simulation. It verified the accuracy and reliability of the equivalent modeling method of virtual material based on contact mechanics experiment.

Nomenclature

     
  • |$\lambda $|

    The ratio of actual contact to nominal contact area

  •  
  • |$\gamma $|

    Shear strain

  •  
  • v

    Poisson’s ratio

  •  
  • μ

    Coefficient of friction (N/m2)

  •  
  • |$\rho $|

    Density (kg/m3)

  •  
  • f

    Contact pressure (MPa)

  •  
  • h

    Tickness (mm)

  •  
  • s

    Sampling frequency (Hz)

  •  
  • D

    Fractal dimension

  •  
  • G

    Fractal roughness (m)

  •  
  • N

    Spectral length (1/m)

  •  
  • P

    Power spectral density (μm3)

  •  
  • |$\omega $|

    Spatial angular frequency (1/μm)

  •  
  • |$\bar F$|

    Average contact pressure (N)

  •  
  • |${A_n}$|

    Nominal contact area (mm2)

  •  
  • |$E,{\rm{\,\,}}{E_n}$|

    Elastic modulus (Pa)

  •  
  • |${E_x}$|

    Elastic modulus in the x-axis direction (Pa)

  •  
  • |${E_y}$|

    Elastic modulus in the y-axis direction (Pa)

  •  
  • |${E_z}$|

    Elastic modulus in the z-axis direction (Pa)

  •  
  • |${G_t}$|

    Shear modulus (Pa)

  •  
  • |${K_n}$|

    Normal contact stiffness (N/mm)

  •  
  • |$K_t^{\rm{*}}$|

    Dimensionless parameter of normal contact stiffness

  •  
  • |${K_t}$|

    Tangential contact stiffness (N/mm)

  •  
  • |$K_t^{\rm{*}}$|

    Dimensionless parameter of tangential contact stiffness

Highlights
  • The contact mechanics experiment at microscopic scale was first introduced using industrial CT equipment.

  • Developed subregion method to handle nonuniform pressure distribution on joint surfaces.

  • Improved equivalent modeling method of virtual material by integrating microscopic contact mechanics and subregion method.

  • Achieved 4.69% average relative error in natural frequency compared with modal experiment.

1. Introduction

The flange joint structure is widely used in pipeline systems, vehicle transmissions and other fields, due to its advantages of simple structure, low cost, repeatable disassembly, and high joint strength (Lv et al., 2024; Lee et al., 2021). As the requirement of product performance continues to be increased and service environment becomes increasingly complex and harsh, the quality of flange connection structure becomes more and more important (Worswick & Finn, 2000; Luan et al., 2012). The quality depends not only on the design and material selection, but also significantly on the contact performance of flange joint surface. The joint surface is defined to be the contact area between flanges. According to statistics, about 60–80% of the total stiffness and more than 90% of the total damping in mechanical structures come from joint surface (Bruno et al., 2005). Therefore, it is essential to study the contact performance of joint surface in the flange connection structure.

The machined surfaces are commonly rough, consisting of many micro-asperities (Blanco-Lorenzo et al., 2023; Jang, 2023; Wang et al., 2021). The actual contact areas are the wave peak parts of these micro-asperities (Xiao et al., 2015; Yu et al., 2022), resulting in the actual contact area of joint surface being far smaller than the nominal contact area (Miao & Huang, 2014; Yastrebov et al., 2017). If the influence of microscopic contact is ignored, there will be a large deviation when analyzing the contact performance of flange connection structure. However, the contact behavior of rough surface mainly contributes on the microscopic scale, and it is difficult to accurately construct a macro-micro coupling model. Even if the coupling modeling is built successfully, a huge cost will be paid for the process of solution and calculation. Therefore, it is crucial to equivalently model the contact behavior of joint surface in flange connection structure considering microscopic contact. Common equivalent modeling methods include spring-damping element method (Li & Hao, 2016; Wei et al., 2020), thin layer element method (Ahmadian et al., 2006; Beaudoin & Behdinan, 2019; Zhao et al., 2018), the virtual material method (Liu et al., 2022; Iranzad & Ahmadian, 2012; Ye et al., 2016), and among others. The methods of spring-damping element and thin layer element are simple in modeling. It is difficult to characterize the mechanical properties of joint surface considering the contact behavior of rough surfaces.

In contrast, virtual material method applies an imaginary material to be directly coupled and fixedly connected with the materials on both sides of joint surface. Then, the corresponding material parameters are assigned to the virtual material layer, according to the pressure and physical characteristics of contact surface. In this way, the rough surface in the microscopic scale is transformed into smooth surface in the macroscopic scale, thus the equivalent simulation of dynamic characteristics is achieved. The virtual material method has good computational speed and accuracy, fully considering the contact characteristics of microscopic topography. It has been applied in the prediction of dynamic characteristics of various types of equipment. The costs of modeling and calculating large equipment can be greatly reduced by simplifying joints such as bolts, bearings and ball screws. However, the virtual material layer is used to directly bond the connected components to form a whole without the presence of fasteners. The method is only used for the prediction of dynamic properties and it is difficult to analyze the mechanical properties.

For the equivalent modeling method of virtual material, it is significant to accurately identify the virtual material parameters of joint surface. Most scholars combined the Hertz contact theory with fractal theory for the identification of virtual material parameters (Lei et al., 2022; Li et al., 2022; Shi & Zhang, 2022; Yu et al., 2021). Tian calculated the elastoplastic contact deformation of asperities considering the fractal characteristics of rough surfaces. A method for identifying the virtual material parameters of joint surface was proposed considering the interaction of tangential and normal contacts (Tian et al., 2011). Based on the model proposed by Tian, Zhao further investigated the influence of uneven contact stress distribution (Zhao et al., 2016). Liao derived the relationship between contact stiffness and normal contact load based on the Hertz contact theory and fractal geometry theory, and proposed a more effective method to identify virtual material parameters (Liao et al., 2016).

Some scholars also adopted other theoretical methods to identify the parameters of virtual material model. For example, Shi constructed the dynamic stiffness relationship of joint surface by GW contact theory, and deduced the material properties according to the stiffness relationship (Shi & Zhang, 2019). However, this method only calculated the normal stiffness and did not consider the tangential stiffness of joint surface. Han established a contact mechanics model by Persson contact theory, and calculated the virtual material parameters of the plate structure (Han et al., 2019). This method directly calculated the elastic parameters of virtual material layer based on the parameters of micro-asperities, without calculating the contact stiffness of joint surface. It causes that the accuracy is not as good as that of fractal contact theory.

In the process of equivalent modeling for virtual material, it is crucial to accurately calculate the actual contact area of joint surface. However, most scholars adopted theoretical analytical methods to calculate the actual contact area, which might deviate significantly from actual engineering. Recently, Shen covered pigments on the joint surface and applied pressure to mark the contact area. Subsequently, a three-dimensional microscopic system was used to measure the size of marked area, which was regarded as the actual contact area (Shen et al., 2023). However, this measurement method has relatively poor accuracy. To the best of the authors’ knowledge, no research has been conducted to accurately measure the actual contact area of non-transparent joint surface by experiment, and the approximate value may cause unreliable and inaccurate equivalent modeling for virtual material. To solve the above issue, this study introduced contact mechanics experiment at microscopic scale to accurately measure the actual contact area, and proposed an equivalent modeling method of virtual material combining contact mechanics experiments and fractal contact theory.

The contributions of the study are as follows:

  1. The contact mechanics experiment at microscopic scale was first introduced using industrial CT equipment. Specialized fixtures and test specimen were designed and machined. An actual contact area of flange connection structure was measured accurately, which would provide the foundation for calculating contact stiffness and the parameters of virtual material.

  2. Considering the nonuniform pressure distribution on the joint surface, the idea of subregion was proposed. At each subregion, the pressure distribution was almost uniform. The equivalent parameters of virtual material were calculated accurately by weighted average of each subregion.

  3. The traditional virtual material method was improved by bringing in contact mechanics experiment at microscopic scale and the idea of subregion. The results of modal experiment and simulation indicated that the average relative error of natural frequency was only 4.69%. It fully validated the accuracy and reliability of the equivalent modeling of virtual material method.

The organization of the paper is as follows. In the first section, the actual assembly process of flange connection structure was simulated by finite-element method. The pressure distribution on the joint surface was analyzed, based on which the subregion division was achieved. Second, contact mechanics experiment based on industrial CT equipment was conducted. The ratios of the actual contact to nominal contact area under different pressures were accurately measured, and the nonlinear relationship between them was obtained through function fitting. In the third section, the fractal parameters were accurately obtained by power spectral density function method after measuring the actual microtopography of joint surface. The actual contact stiffness of joint surface was calculated by fractal contact theory. Next, all the parameters of virtual material were calculated. Finally, a modal testing experiment was designed, and the experimental result was compared with the modal simulation result to verify the accuracy of the equivalent modeling method of virtual material proposed in this study.

2. Subregion Division of Flange Connection Structure

The flange connection structure analyzed in this paper is shown in Figure 1, which consists of two flanges, six bolts and nuts. The upper and lower flanges come into contact under the preload generated by tightening the bolts and nuts (the contact area is called joint surface). The detailed structural dimensions of flange are shown in Figure 2, and the material is aluminum alloy. The specification of bolt and nut is M10 × 1.5, and the material is structural steel. In the assembly process, there is obvious phenomenon of uneven pressure distribution on the joint surface after tightening bolts and nuts. If the average pressure distribution on the joint surface is assumed instead of the uneven pressure distribution, the equivalent modeling accuracy of virtual material will be decreased. Therefore, the finite-element analysis (FEA) method is applied to simulate the assembly process of flange connection structure. The actual pressure distribution on the joint surface is extracted, and based on it, the subregions with the same pressure are divided.

Geometrical model of flange connection structure.
Figure 1:

Geometrical model of flange connection structure.

Structural dimensions of flange.
Figure 2:

Structural dimensions of flange.

2.1. Finite-element modeling

The finite-element model of flange connection structure was established using software Hypermesh 12.0, as shown in Figure 3. The finite-element model contained a total of 9,36,132 hexahedral meshes. For two flanges, an automatic meshing method was used to obtain uniform and dense hexahedral meshes. For bolts and nuts, the accurate modeling method of thread mesh proposed by Fukuoka and Nomura was employed (Fukuoka & Nomura, 2008). The divided hexahedral meshes were well distributed, with high-mesh quality and computational accuracy.

Finite-element model of flange connection structure.
Figure 3:

Finite-element model of flange connection structure.

2.2. Subregion division

First, the parameters of finite-element model were set. The type of all the meshes was Solid185, which are mainly applied to three-dimensional solid structures and sensitive to deformation and stress. In terms of material properties, bolts and nuts were structural steel material with a density of 7850 kg/m3, Young’s modulus of 206 GPa, and Poisson’s ratio of 0.3, respectively. The flanges were aluminum alloy material with a density of 2700 kg/m3, Young’s modulus of 69 GPa and Poisson’s ratio of 0.33, respectively. In terms of contact interfaces, the model contained thread surfaces, joint surface, and among others. The friction coefficients of all contact interfaces were set to 0.15, which was a typical value adopted by other researchers (Izumi et al., 2005). In terms of boundary conditions, fixed constraints were applied to the bottom of flange and bolt. Meanwhile, the flange sides were also constrained in the x- and y-directions to prevent sliding during assembly (the coordinate system was shown in Figure 2). The same tightening angle was applied sequentially to the six nuts to simulate the assembly process. The bolt size was M10 and the actual preload generated in the engineering application was approximately 20 kN, resulting in a rotation angle of approximately 28°.

Subsequently, the established finite-element model of flange connection structure was imported into the commercial software Ansys Workbench 19.2 for assembly simulation solution. During the solving process, the global stiffness matrix of all the nodes used to represent the relationship between the displacements and forces was derived. The virtual work method or variation method was be applied to solve the global stiffness matrix. After the simulation was completed, the mechanical behavior of material was presented, and the pressure distribution of joint surface was extracted, as shown in Figure 4. Subregions on the joint surface were divided according to the pressure distribution pattern in the figure. The joint surface was divided into seven subregions, and the pressure distribution at each subregion was almost the same. Subregions 1–7 were numbered according to the distance from bolt axis. The subregion close to the bolt axis was subregion 1, and the subregion away from the bolt axis was subregion 7. Table 1 shows the average contact pressure |$\bar F$| and nominal contact area An for each subregion, which will be used for the calculation of virtual material parameters.

Pressure distribution and subregion division of joint surface.
Figure 4:

Pressure distribution and subregion division of joint surface.

Table 1:

Average contact pressures and nominal contact areas of subregions.

Subregion number1234567
|$\bar F$|/MPa81.29171.93662.58253.22843.87428.84312.114
An /mm2591.122429.212239.970262.687283.6041196.2868308.250
Subregion number1234567
|$\bar F$|/MPa81.29171.93662.58253.22843.87428.84312.114
An /mm2591.122429.212239.970262.687283.6041196.2868308.250
Table 1:

Average contact pressures and nominal contact areas of subregions.

Subregion number1234567
|$\bar F$|/MPa81.29171.93662.58253.22843.87428.84312.114
An /mm2591.122429.212239.970262.687283.6041196.2868308.250
Subregion number1234567
|$\bar F$|/MPa81.29171.93662.58253.22843.87428.84312.114
An /mm2591.122429.212239.970262.687283.6041196.2868308.250

3. Microtopography Measurement and Fractal Characterization

In order to accurately calculate the actual contact stiffness of joint surface, the microtopography of joint surface was measured and fractally characterized. Super View W1 series optical 3D surface profiler was applied for measuring the surface topography. Figure 5 shows the equipment, which consists of an interferometry head, support column, lens group, removable platform, shock-absorbing base, joystick, and computer. The experimental specimen is shown in Figure 6A. It was placed on the platform under the lens during measurement. The relative position was adjusted by joystick until the lens and testing surface achieved a focusing distance and the interference fringes appeared on the video window of computer. The microtopography of experimental specimen could be obtained by automatic measurement method. Based on the above operations, the three-dimensional (3D) microtopography of joint surface with 5 × 5 mm area was measured, as shown in Figure 6B.

Optical 3D surface profiler.
Figure 5:

Optical 3D surface profiler.

Experimental specimen and measured 3D microtopography. (A) Experimental specimen. (B) 3D microtopography.
Figure 6:

Experimental specimen and measured 3D microtopography. (A) Experimental specimen. (B) 3D microtopography.

After the microtopography on the joint surface was obtained, the method of power spectral density function was used to calculate fractional parameters. This method described the actual spatial distribution of surface profile by calculating the spectral composition of surface profile waveform. The detailed procedure was as follows. First, the cross-sectional point cloud coordinates were intercepted as the input signal, and the input signal was divided into multiple overlapping segments. Next, a window function was applied to each segment of signal to reduce spectral leakage. The discrete Fourier transform was employed to each segment of signal to obtain spectrum,

(1)

where xn is the nth coordinate sample of input signal, N is the length of segment, and k is the frequency index. The power spectrum of each segment was calculated by squaring the spectrum of each segment. For real signal, the power spectral density was given by

(2)

The spatial angular frequency ω was calculated based on the sampling frequency s and the spectral length N,

(3)

The measured point cloud of the actual surface profile could be converted into power spectral density P and spatial angular frequency ω through Equations 2 and 3. The power spectral density function was plotted by taking the logarithm of both ω and P as the horizontal and vertical coordinates, respectively, as shown in Figure 7. The data points in the figure were fitted linearly to obtain the slope of the fitted line kPSD = –2.41263 and the intercept bPSD = –0.74056. The logarithmic expression of the power spectral density function was known from the reference (Pan et al., 2017),

(4)

where D is the fractal dimension of rough surface, and the physical significance of D is the extent of space filling by rough surface, with larger values corresponding to denser fillings; G is the fractal roughness of rough surface, and a rougher surface is characterized by larger G values; P is power spectral density; ω is the spatial angular frequency; and γ represents the spatial frequency of random profile (usually taking the value of 1.5). Subsequently, Equation 4 was regarded as a first-order linear function, with the lgP(ω) as the ordinate and the lgω as the abscissa. The expressions for the corresponding slope and intercept were obtained:

(5)

where γ, kPSD and bPSD were known. Therefore, D = 2.2937 and G = 4.6196 × 10−10 m could be calculated.

Power spectral density and function fitting.
Figure 7:

Power spectral density and function fitting.

4. Measurement of Contact Area Based on Industrial CT Equipment

For the flange joint surface, the actual contact area is much smaller than the nominal contact area. According to the fractal contact theory, the actual contact area of joint surface is directly related to the contact stiffness. This will significantly affect the equivalent modeling accuracy of virtual material. In this section, industrial CT equipment (industrial computerized tomography equipment) was employed to measure the actual contact states of joint surfaces in the microscopic scale. First, the experimental specimens and fixtures were designed and machined according to the requirements of equipment. Next, the ratio of actual contact areas (ratio of actual contact area to nominal contact area) on the joint surface at different pressures were obtained from the measurements by industrial CT equipment. Finally, a fitting function was selected to describe the relationship between the ratio of actual contact to nominal contact area and contact pressure.

4.1. Aluminum alloy specimen and fixture

Industrial CT equipment can only detect the specimen with small size, thus the joint surfaces of flange connection structure cannot be measured directly (Zhang et al., 2020). Herein, two small-sized aluminum alloy specimens were designed for contact mechanics experiment in the microscopic scale. In order to ensure the consistency of microscopic contact characteristics with the joint surface of flange connection structure, the machining method of both specimens was the same with that of the flanges. Two cylindrical aluminum alloy specimens were machined as shown in Figure 8. The upper circular surface was the contact surface with a diameter of 3 mm. The height of the cylinder was 10 mm. During the experiment, a fixture required to be designed to make the specimens fixed and a certain pressure applied to the specimens. The design scheme of fixture was shown in Figure 8. The experimental fixture consisted of a fastening bolt, end caps, loading bars, tube, blocks, pressure sensor, and base. The specimens were fixed in the center of the tube. During the experiment, the tested samples was scanned by X-rays after penetrating the tube. In order to obtain the better experimental results, the tube was made of plastic material. This material had low density and little influence on the penetration ability of X-rays. A certain pressure could be applied to the specimens through the fastening bolt. The pressure value could be read through sensor on the lower part. The physical pictures of experimental fixture and pressure sensor were shown in Figure 9.

Details of aluminum alloy specimens and experimental fixture.
Figure 8:

Details of aluminum alloy specimens and experimental fixture.

Experimental fixture and pressure sensor.
Figure 9:

Experimental fixture and pressure sensor.

4.2. Contact mechanics experiment by industrial CT equipment

Contact mechanics experiment based on industrial CT equipment was carried out. The nanoVoxel-5000 series dual-ray source CT system was used for the experiment. The spatial resolution of CT system is 500 nm, and the measurement accuracy is relatively high. The detailed procedures of contact mechanics experiment were as follows. First, the fastening bolt was disassembled, and the two specimens were installed into the tube of experimental fixture. Second, the fastening bolt was installed and the contact pair of specimens was adjusted for the convenience of testing. Finally, the experimental fixture was fixed on the stage, and a pressure was applied to the contact interface by tightening the fastening bolt. The assembly of the experimental fixture and specimen was shown in Figure 10.

Experimental equipment and fixture installation.
Figure 10:

Experimental equipment and fixture installation.

During the measurement, X-rays were emitted from the ray source under the control of computer. After passing through the measured area and experiencing partial energy loss, X-rays were recorded on the detector. The stage was repeatedly rotated at a small angle, moved up and down, in order to achieve a comprehensive measurement of the contact area. Finally, the computer processed the signal recorded by all detectors to reconstruct a complete three-dimensional model of contact area. The overall measurement process was shown in Figure 11.

Schematic of X-ray computed tomography.
Figure 11:

Schematic of X-ray computed tomography.

The schematic of contact state distribution under 20 MPa pressure was obtained. The 3D contact region could be projected along the axis direction to obtain a 2D diagram, as shown in Figure 12. Figures 12A and B shows the 3D and 2D views, respectively. Generally, the regions with voids along the axis direction (blue regions) were considered as non-contact regions and other regions (black regions) were considered as contact regions, as shown in Figure 13. In the figure, the blue regions were the non-contact area, and the black regions were the actual contact area. The ratio of the black area to the total circular contact area was 0.1915, thereby the ratio of actual contact to nominal contact area under 20 MPa pressure was 0.1915.

Measurement result of contact area. (A) Three-dimensional view. (B) Two-dimensional view.
Figure 12:

Measurement result of contact area. (A) Three-dimensional view. (B) Two-dimensional view.

Schematic of actual contact state.
Figure 13:

Schematic of actual contact state.

4.3. Relationship between the ratio of actual contact to nominal contact area and contact pressure

In order to explore the actual contact states of joint surface under different contact pressures, the fastening bolt was further tightened to increase the contact pressure. The contact pressure was changed to 20, 40, 60, and 80 MPa. The contact states and the ratios of actual contact to nominal contact area under different pressures were shown in Table 2.

Table 2:

Contact states and the ratios of actual contact to nominal contact area under different pressures.

graphic
graphic
Table 2:

Contact states and the ratios of actual contact to nominal contact area under different pressures.

graphic
graphic

The data in Table 2 was the ratios of actual contact to nominal area under four specific pressures. A third-order polynomial function λ(f) could be used to equivalently represent the relationship between the contact pressure and the ratio of actual contact to nominal contact area,

(6)

where λ is the ratio of actual contact to nominal contact area; f is the contact pressure; a0, a1, a2, and a3 are polynomial coefficients, a= 0.11243, a= 0.00468, a2 = −4.0286 × 10−5, and a= 2.10493 × 10−7. The distribution between contact pressure and the ratio of actual contact to nominal contact area is shown in Figure 14. In the figure, the scattered points were the experimental results, and the curve was the fitting function. This function would be used to calculate the actual contact area in the virtual material method.

Distribution between contact pressure and the ratio of actual contact to nominal contact area.
Figure 14:

Distribution between contact pressure and the ratio of actual contact to nominal contact area.

5. Calculation of Virtual Material Parameters

The equivalent modeling approach of virtual material provides the better description of the microtopography, contact stiffness, and mechanical behavior of contact area. The equivalent virtual material layer is inserted between two contacted components, replacing the joint surface, as shown in Figure 15. The parameters of equivalent virtual material are calculated through theoretical modeling, including equivalent thickness, elastic modulus, shear modulus, Poisson's ratio, and density.

Equivalent modeling approach of virtual material.
Figure 15:

Equivalent modeling approach of virtual material.

5.1. Equivalent thickness

According to the reference (Zhao et al., 2016), the thickness of virtual material layer could be regarded as the sum of the peak height of micro-asperities and the thickness of deformation area during the machining process. In this study, the sum of the above machining parameters of the aluminum alloy component is approximately 0.5 mm, thus the thickness of the virtual material layer is h = 1 mm.

5.2. Equivalent elastic modulus

The elastic modulus is defined as the strain under unit stress. When the equivalent thickness of the virtual material layer h is relatively small, it is assumed that the average stress of virtual material layer is the ratio of the external load to the contact area. Based on it, the strain of virtual material layer is the ratio of the change in thickness to itself. Therefore, the equivalent elastic modulus En of virtual material can be obtained as follows:

(7)

where Kn is the normal contact stiffness of joint surface. The normal contact stiffness of joint surface is provided by the normal contact stiffness of each micro-asperity. According to the fractal contact theory, the dimensionless parameter of normal contact stiffness can be obtained as follows (Yang et al., 2020),

(8)

where λ is the ratio of actual contact to nominal contact area, which is calculated by Equation 6; D is the fractal dimension, which is calculated by Equation 5; |$g( D ) = {(2 - D)^{D/2}} \cdot {D^{( {2 - D} )/2}} \cdot {(D - 1)^{ - 1}}$|⁠; |$a_c^* = {G^{*2}} \cdot {[H/(2E)]^{ - 2/(D - 1)}}$|⁠; H is the hardness of soft material; G* is the dimensionless parameter of characteristic length scale, |${G^*} = G/\sqrt {{A_n}} $|⁠; An is the nominal contact area; G is the fractal roughness. The relationship between the dimensionless normal contact stiffness and the actual normal contact stiffness of joint surface is as follows:

(9)

where E is the equivalent elastic modulus of two flanges, |${E^{ - 1}} = ( {1 - \nu _1^2} )E_1^{ - 1} + ( {1 - \nu _2^2} )E_2^{ - 1}$|⁠; ν1, ν2 and E1, E2 are the Poissons ratios and elastic moduli of the upper and lower flange materials, respectively. In the study, they are the material parameters of aluminum alloy. The normal contact stiffness of joint surface can be calculated by substituting Equations 6 and 8 into Equation 9.

The elastic modulus of joint surface in the x-axis direction is related to the tightness of contact surface and the elastic modulus of each flange. The dimensionless contact area λ reflects the tightness of contact surface. The elastic modulus Ex of virtual material layer in the x-axis direction is determined by the following formula:

(10)

In the study, the virtual material layer is regarded to be isotropic. Therefore, the elastic modulus of the virtual material layer in the y-axis direction Ey is equal to Ex.

Finally, the weight of subregions is considered (Wang et al., 2021). It is assumed that the joint surface is divided into n subregions in total. According to the calculation method for elastic modulus of uniform composite materials, the average elastic modulus is calculated based on each subregion. It is the weighted average of the equivalent elastic modulus of each subregion En and its corresponding nominal contact area An,

(11)

According to the above method, the equivalent elastic modulus |${\bar E_n}$| = 6.0066 × 1013 Pa; |${\bar E_x}$| = |${\bar E_y}$| = 1.4562 × 1010 Pa had been calculated.

5.3. Equivalent shear modulus

The equivalent shear modulus of virtual material is the ratio of the shear stress τ to the shear strain γ. When the shear strain γ is relatively small, the shear strain γ and the tangential deformation δt satisfy γ ≈ tanγ = δt/h. The average shear stress of virtual material is the ratio of the tangential load to the nominal contact area. Therefore, the equivalent shear modulus Gt of virtual material can be obtained as follows:

(12)

where Kt is the tangential contact stiffness of joint surface. The tangential contact stiffness of joint surface is contributed by the tangential contact stiffness of micro-asperities on the joint surface. The dimensionless parameter of tangential contact stiffness is given

(13)

where T and P are the tangential and normal loads, respectively; µ is the coefficient of friction; v is the comprehensive Poisson's ratio of the contact pair, which is the Poisson’s ratio of aluminum alloy in the study. The relationship between the dimensionless tangential contact stiffness and the actual tangential contact stiffness is as follows:

(14)

where |$\bar G$| is the comprehensive shear modulus of the contact pair, |$\bar{G} = E/( {2 + 2\nu } )$|⁠. The tangential contact stiffness of joint surface can be calculated by substituting Equations 6 and 13 into Equation 14.

Similarly, the weight of subregions is considered, it is assumed that the joint surface is divided into n subregions in total. The average shear modulus is the weighted average of the equivalent shear modulus of each subregion Gt and its corresponding nominal contact area An,

(15)

According to the above method, the equivalent shear modulus |${\bar G_t}$| = 2.4075 × 109 Pa had been calculated.

5.4. Equivalent density

The equivalent density of virtual material layer is the ratio of the equivalent mass to the equivalent volume,

(16)

where V is the volume of contact part; V1 and V2 are the volumes of the micro-asperities in the contact area, respectively; m1 and m2 are the masses of micro-asperities in the contact area, respectively; h1 and h2 are the thicknesses of the micro-asperities in the contact area, respectively, and h1 = h= 0.5 mm; ρ1 and ρ2 are the densities of the upper and lower flanges, respectively, which is the density of aluminum alloy material in the study.

Similarly, the weights of subregions are considered. It is assumed that the joint surface is divided into n subregions in total. The average density is the weighted average of the equivalent density of each subregion ρ and its corresponding nominal contact area An,

(17)

According to the above method, the equivalent density |$\bar \rho {\rm{\,\,}}$|= 561.6801 kg/m3 had been calculated.

5.5. Equivalent Poisson’s ratio

There are a large number of gaps between the contacted micro-asperities. According to the fractal contact theory, when the joint surface is subjected to normal loads, normal and tangential deformations will be generated on the contacted micro-asperities. Subsequently, the tangential deformation under the normal loads will fill the gaps between each micro-asperity. Therefore, the overall tangential deformation of joint surface is zero. It is considered that the equivalent Poisson’s ratio of virtual material is zero.

6. Experimental Validation

In order to verify the accuracy and reliability of equivalent modeling method, a finite-element model was established according to the virtual material parameters calculated in Section 5. The modal and natural frequency were obtained through finite-element simulation method. At the same time, an experimental apparatus was built to measure the mode and natural frequency of the testing specimen. By comparing the results of simulation and experiment, the reliability and accuracy of the proposed equivalent modeling method were verified.

6.1. Modal simulation analysis

The finite-element model including the virtual material layer was established by using software Hypermesh 12.0, as shown in Figure 16. Compared with the model in Section 2.1, a thin virtual material layer was added between the two flanges. The material parameters of bolts, nuts and flanges have been defined in Section 2. And the material parameters of the virtual material layer have been calculated in Section 5. As a comparison, the microscopic contact of the joint surface is not considered. The virtual material parameters of the joint surface were calculated by the contact stiffness calculation method proposed in reference (Li et al., 2020). The virtual material parameters calculated by the above two methods were shown in Table 3. The model was imported into the software Ansys Workbench 19.2 for modal simulation. All contact types were set as bonding contact. Since the free mode was required to be analyzed, there was no need to apply constraints. After the above parameters were set, the modal simulation was done in the software. The first four order vibration modes and natural frequencies obtained from the simulation were shown in Figure 17.

Finite-element model of flange connection structure with virtual material layer.
Figure 16:

Finite-element model of flange connection structure with virtual material layer.

The first four order vibration modes and natural frequencies.
Figure 17:

The first four order vibration modes and natural frequencies.

Table 3:

The parameters of virtual material.

Parameters of virtual materialMethod proposed in the paperNo consideration of micro contact (control group)
|${\bar E_z}$|6.0066 × 1013 Pa9.0832 × 108 Pa
|${\bar E_x}$|1.4562 × 1010 Pa9.0832 × 108 Pa
|${\bar G_t}$|2.4075 × 109 Pa2.4075 × 109 Pa
|$\upsilon $|00
|$\bar \rho $|561.6801 kg/m32750 kg/m3
Parameters of virtual materialMethod proposed in the paperNo consideration of micro contact (control group)
|${\bar E_z}$|6.0066 × 1013 Pa9.0832 × 108 Pa
|${\bar E_x}$|1.4562 × 1010 Pa9.0832 × 108 Pa
|${\bar G_t}$|2.4075 × 109 Pa2.4075 × 109 Pa
|$\upsilon $|00
|$\bar \rho $|561.6801 kg/m32750 kg/m3
Table 3:

The parameters of virtual material.

Parameters of virtual materialMethod proposed in the paperNo consideration of micro contact (control group)
|${\bar E_z}$|6.0066 × 1013 Pa9.0832 × 108 Pa
|${\bar E_x}$|1.4562 × 1010 Pa9.0832 × 108 Pa
|${\bar G_t}$|2.4075 × 109 Pa2.4075 × 109 Pa
|$\upsilon $|00
|$\bar \rho $|561.6801 kg/m32750 kg/m3
Parameters of virtual materialMethod proposed in the paperNo consideration of micro contact (control group)
|${\bar E_z}$|6.0066 × 1013 Pa9.0832 × 108 Pa
|${\bar E_x}$|1.4562 × 1010 Pa9.0832 × 108 Pa
|${\bar G_t}$|2.4075 × 109 Pa2.4075 × 109 Pa
|$\upsilon $|00
|$\bar \rho $|561.6801 kg/m32750 kg/m3

6.2. Experimental validation of modal testing

A modal testing apparatus was built to measure the mode and natural frequency of the real flange connection structure, as shown in Figure 18. The testing apparatus consists of a force hammer, a torque wrench, acceleration sensors, a dynamic signal testing and analysis system and a computer. The sensitivity of the force hammer was 1.059 mV/N, and the sensitivity of the acceleration sensors was 10.29 mV/g.

Schematic of modal testing apparatus.
Figure 18:

Schematic of modal testing apparatus.

The experimental process of modal testing was as follows. First, a tightening torque of 37.4 Nm was applied to the nuts by torque wrench. A preload of 20 kN was generated in the bolts to achieve the assembly of flange connection structure. An elastic rope was used to suspend the experimental component to simulate an approximate free state, as shown in Figure 19. Subsequently, a force hammer was used to strike the specific positions of the experimental component. The acceleration sensor was employed to obtain the vibration responses at different excitation points. Based on this, the dynamic signal testing and analysis system calculated the frequency–response function relationships between the corresponding excitation points and vibration pickup points. Finally, the amplitude–frequency curve and phase–frequency curve of the frequency–response function were depicted. The low order natural frequencies of the experimental component were obtained. During the modal testing, there was a certain mass in the sensor itself. In order to reduce the influence of the additional mass on the testing result, the method of vibration pickup by single point and excitation by multi points was adopted for the modal testing. That is to say, the position of the sensor remained unchanged during the testing process. The low order natural frequencies obtained from the experiment were compared with the simulation results, as shown in Table 4. It can be seen that the relative errors of each order natural frequency obtained from the simulation and the experiment are all less than 5% (except that the error of the first order natural frequency is relatively high). The average relative error of the first four order natural frequencies is 4.69% while that of the control group is 14.19%. It indicates that the prediction accuracy of our method is much higher than that of the control group. It validates the accuracy and reliability of the equivalent modeling method of virtual material proposed in the study.

Testing site of elastic rope suspension.
Figure 19:

Testing site of elastic rope suspension.

Table 4:

Comparison of experiment and simulation results.

Natural frequencies order1234
Experimental natural frequency (Hz)57.19884.262103.724155.973
Simulated natural frequency (Hz)50.26781.394105.637158.165
Relative error–12.11%–3.40%1.84%1.41%
Simulated natural frequency of control group (Hz)42.38872.18695.362144.726
Relative error–25.89%–14.33%–8.06%–8.49%
Natural frequencies order1234
Experimental natural frequency (Hz)57.19884.262103.724155.973
Simulated natural frequency (Hz)50.26781.394105.637158.165
Relative error–12.11%–3.40%1.84%1.41%
Simulated natural frequency of control group (Hz)42.38872.18695.362144.726
Relative error–25.89%–14.33%–8.06%–8.49%
Table 4:

Comparison of experiment and simulation results.

Natural frequencies order1234
Experimental natural frequency (Hz)57.19884.262103.724155.973
Simulated natural frequency (Hz)50.26781.394105.637158.165
Relative error–12.11%–3.40%1.84%1.41%
Simulated natural frequency of control group (Hz)42.38872.18695.362144.726
Relative error–25.89%–14.33%–8.06%–8.49%
Natural frequencies order1234
Experimental natural frequency (Hz)57.19884.262103.724155.973
Simulated natural frequency (Hz)50.26781.394105.637158.165
Relative error–12.11%–3.40%1.84%1.41%
Simulated natural frequency of control group (Hz)42.38872.18695.362144.726
Relative error–25.89%–14.33%–8.06%–8.49%

7. Conclusions

This paper proposed an equivalent modeling method of virtual material for flange connection structure based on contact mechanics experiment and fractal contact theory. The conclusions were as follows:

  1. Industrial CT experiment was conducted to measure the actual contact area under different pressures. A third-order polynomial function was applied to nonlinear fitting on the contact pressure and the ratio of actual contact area to nominal contact area. The fitting function could accurately calculate the ratio of actual contact area to nominal contact area under different pressures.

  2. The idea of subregion was proposed. Seven subregions were divided and pressure distribution at each subregion was approximately the same based on the result of FEA. The actual microtopography at each subregion was measured experimentally and the fractal parameters were obtained. Fractal contact theory was applied to calculate actual contact stiffness of each subregion. The parameters of equivalent virtual material were calculated by weighted average of each subregion, which had higher accuracy.

  3. A modal testing experiment for the flange connection structure was conducted. A finite-element model containing virtual material was established. The virtual material parameters were obtained according to the method proposed in this study, and a modal simulation was also performed. The result indicated that the average relative error of the low-order natural frequencies obtained from the experiment and the simulation was only 4.69%, which validated the accuracy and reliability of the equivalent modeling of virtual material method proposed in the study.

In summary, our method can achieve the accurate equivalent modeling for virtual material by fully considering the actual contact area and pressure distribution of joint surface, and microscopic contact deformation. The experimental results of modal testing validated the high prediction accuracy. However, the equipment of industrial CT is expensive, and the measurement procedures are complicated. Specialized fixtures and test specimen require to be designed and machined. The experimental costs are high. In order to apply our method in engineering, in the further work, neural network can be employed to predict the actual contact area in the conditions of different input parameters, including surface morphology, pressure distribution, material properties, and among others.

Conflicts of Interest

The authors declare no conflict of interest.

Author Contributions

Xingjie Wang: Conceptualization, Formal analysis, Writing—original draft. Hao Gong: Funding acquisition, Writing—review & editing, Methodology. Jianhua Liu: Investigation, Supervision. Penghao Zhao: Software, Visualization. Rongquan Zhu: Data curation, Resources.

Funding

This work was support by National Key R&D Program of China (No. 2022YFB3304200, and No. 2022YFB3403800), National Natural Science Foundation of China (grant numbers U22A20203, 52105503, and 51935003), and Hebei Natural Science Foundation (E2023105059).

References

Ahmadian
 
H.
,
Mottershead
 
J. E.
,
James
 
S.
,
Friswell
 
M. I.
,
Reece
 
C. A.
(
2006
).
Modelling and updating of large surface-to-surface joints in the AWE-MACE structure
.
Mechanical Systems and Signal Processing
,
20
,
868
880
.. .

Beaudoin
 
M. A.
,
Behdinan
 
K.
(
2019
).
Analytical lump model for the nonlinear dynamic response of bolted flanges in aero-engine casings
.
Mechanical Systems and Signal Processing
,
115
,
14
28
.. .

Blanco-Lorenzo
 
J.
,
Liu
 
S.
,
Santamaria
 
J.
,
Meehan
 
P. A.
,
Vadillo
 
E. G.
(
2023
).
Frictional contact analysis in a spherical roller bearing
.
Journal of Computational Design and Engineering
,
10
,
139
159
.. .

Bruno
 
D.
,
Greco
 
F.
,
Lonetti
 
P.
(
2005
).
A 3D delamination modelling technique based on plate and interface theories for laminated structures
.
European Journal of Mechanics-A/Solids
,
24
,
127
149
.. .

Fukuoka
 
T.
,
Nomura
 
M.
(
2008
).
Proposition of helical thread modeling with accurate geometry and finite element analysis
.
Journal of Pressure Vessel Technology
,
130
,
011204
. .

Han
 
R.
,
Li
 
G.
,
Gong
 
J.
,
Zhang
 
M.
,
Zhang
 
K.
(
2019
).
Equivalent method of joint interface based on persson contact theory: Virtual material method
.
Materials
,
12
,
3150
. .

Iranzad
 
M.
,
Ahmadian
 
H.
(
2012
).
Identification of nonlinear bolted lap joint models
.
Computers & Structures
,
96
,
1
8
.. .

Izumi
 
S.
,
Yokoyama
 
T.
,
Iwasaki
 
A.
,
Sakai
 
S.
(
2005
).
Three-dimensional finite element analysis of tightening and loosening mechanism of threaded fastener
.
Engineering Failure Analysis
,
12
,
604
615
.. .

Jang
 
S.
(
2023
).
Deterministic surface roughness effects on elastic material contact with shear thinning fluid media
.
Journal of Computational Design and Engineering
,
10
,
2312
2331
.. .

Lee
 
I. D.
,
Lee
 
I.
,
Han
 
S.
(
2021
).
3D reconstruction of as-built model of plant piping system from point clouds and port information
.
Journal of Computational Design and Engineering
,
8
,
195
209
.. .

Lei
 
S.
,
Mao
 
K.
,
Tian
 
W.
,
Li
 
L.
(
2022
).
A three-dimensional transition interface model for bolt joint
.
Machines
,
10
,
511
. .

Li
 
L.
,
Wang
 
J.
,
Shi
 
X.
,
Chen
 
Z.
,
Cai
 
A.
(
2020
).
Dynamic properties of bolted joints characterized by functionally gradient virtual materials method
.
Journal of Vibration Engineering
,
33
,
525
532
.. (in Chinese)​ ​​ ​​​​.

Li
 
Y.
,
Hao
 
Z.
(
2016
).
A six-parameter Iwan model and its application
.
Mechanical Systems and Signal Processing
,
68
,
354
365
.. .

Li
 
Z.
,
Zhang
 
Y.
,
Li
 
C.
,
Xu
 
M.
,
Dai
 
W.
,
Liu
 
Z.
(
2022
).
Modeling and nonlinear dynamic analysis of bolt joints considering fractal surfaces
.
Nonlinear Dynamics
,
108
,
1071
1099
.. .

Liao
 
J.
,
Zhang
 
J.
,
Feng
 
P.
,
Yu
 
D.
,
Wu
 
Z.
(
2016
).
Interface contact pressure-based virtual gradient material model for the dynamic analysis of the bolted joint in machine tools
.
Journal of Mechanical Science and Technology
,
30
,
4511
4521
.. .

Liu
 
X.
,
Sun
 
W.
,
Liu
 
H.
,
Du
 
D.
,
Ma
 
H.
(
2022
).
Semi-analytical modeling and analysis of nonlinear vibration of bolted thin plate based on virtual material method
.
Nonlinear Dynamics
,
108
,
1247
1268
.. .

Luan
 
Y.
,
Guan
 
Z. Q.
,
Cheng
 
G. D.
,
Liu
 
S.
(
2012
).
A simplified nonlinear dynamic model for the analysis of pipe structures with bolted flange joints
.
Journal of Sound and Vibration
,
331
,
325
344
.. .

Lv
 
T.
,
Chen
 
Z.
,
Wang
 
D.
,
Du
 
X.
(
2024
).
Parametric design and modeling method of carbon fiber reinforcement plastic-laminated components applicable for multi-material vehicle body development
.
Journal of Computational Design and Engineering
,
11
,
261
287
.. .

Miao
 
X.
,
Huang
 
X.
(
2014
).
A complete contact model of a fractal rough surface
.
Wear
,
309
,
146
151
.. .

Pan
 
W.
,
Li
 
X.
,
Wang
 
L.
,
Guo
 
N.
,
Mu
 
J.
(
2017
).
A normal contact stiffness fractal prediction model of dry-friction rough surface and experimental verification
.
European Journal of Mechanics-A/Solids
,
66
,
94
102
.. .

Shen
 
J. X.
,
Chen
 
Y.
,
Xu
 
P.
,
Yu
 
Y. H.
(
2023
).
Virtual material simulation of contact properties of steel-BFPC joint surface based on discrete principle
.
Engineering Mechanics
,
40
,
205
214
.
(in Chinese)
.

Shi
 
K.
,
Zhang
 
G.
(
2019
).
A parameterized model of fixed joint interface based on virtual material
.
Journal of Mechanical Science and Technology
,
33
,
5209
5217
.. .

Shi
 
W.
,
Zhang
 
Z.
(
2022
).
An improved contact parameter model with elastoplastic behavior for bolted joint interfaces
.
Composite Structures
,
300
,
116178
. .

Tian
 
H.
,
Li
 
B.
,
Liu
 
H.
,
Mao
 
K.
,
Peng
 
F.
,
Huang
 
X.
(
2011
).
A new method of virtual material hypothesis-based dynamic modeling on fixed joint interface in machine tools
.
International Journal of Machine Tools and Manufacture
,
51
,
239
249
.. .

Wang
 
R.
,
Liu
 
J.
,
Zhang
 
F.
,
Ding
 
X.
(
2021
).
An approach to evaluate the sealing performance of sealing structures based on multiscale contact analyses
.
Journal of Computational Design and Engineering
,
8
,
1433
1445
.. .

Wei
 
C.
,
Zhu
 
H.
,
Lang
 
S.
(
2020
).
A modified complete normal contact stiffness model of a fractal surface considering contact friction
.
Fractals
,
28
,
2050081
. .

Worswick
 
M. J.
,
Finn
 
M. J.
(
2000
).
The numerical simulation of stretch flange forming
.
International Journal of Plasticity
,
16
,
701
720
.. .

Xiao
 
H.
,
Shao
 
Y.
,
Brennan
 
M. J.
(
2015
).
On the contact stiffness and nonlinear vibration of an elastic body with a rough surface in contact with a rigid flat surface
.
European Journal of Mechanics-A/Solids
,
49
,
321
328
.. .

Yang
 
H.
,
Che
 
X.
,
Yang
 
C.
(
2020
).
Investigation of normal and tangential contact stiffness considering surface asperity interaction
.
Industrial Lubrication and Tribology
,
72
,
379
388
.. .

Yastrebov
 
V. A.
,
Anciaux
 
G.
,
Molinari
 
J. F.
(
2017
).
On the accurate computation of the true contact-area in mechanical contact of random rough surfaces
.
Tribology International
,
114
,
161
171
.. .

Ye
 
H.
,
Huang
 
Y.
,
Li
 
P.
,
Li
 
Y.
,
Bai
 
L.
(
2016
).
Virtual material parameter acquisition based on the basic characteristics of the bolt joint interfaces
.
Tribology International
,
95
,
109
117
.. .

Yu
 
X.
,
Sun
 
Y.
,
Li
 
H.
,
Wu
 
S.
(
2022
).
An improved meshing stiffness calculation algorithm for gear pair involving fractal contact stiffness based on dynamic contact force
.
European Journal of Mechanics-A/Solids
,
94
,
104595
. .

Yu
 
Y.
,
Cheng
 
H.
,
Liang
 
B.
,
Di
 
Z. H. A. O.
,
Junshan
 
H. U.
,
Zhang
 
K.
(
2021
).
A novel virtual material layer model for predicting natural frequencies of composite bolted joints
.
Chinese Journal of Aeronautics
,
34
,
101
111
.. .

Zhang
 
F.
,
Liu
 
J.
,
Ding
 
X.
,
Yang
 
Z.
(
2020
).
A discussion on the capability of X-ray computed tomography for contact mechanics investigations
.
Tribology International
,
145
,
106167
. .

Zhao
 
G.
,
Xiong
 
Z.
,
Jin
 
X.
,
Hou
 
L.
,
Gao
 
W.
(
2018
).
Prediction of contact stiffness in bolted interface with natural frequency experiment and FE analysis
.
Tribology International
,
127
,
157
164
.. .

Zhao
 
Y.
,
Yang
 
C.
,
Cai
 
L.
,
Shi
 
W.
,
Liu
 
Z.
(
2016
).
Surface contact stress-based nonlinear virtual material method for dynamic analysis of bolted joint of machine tool
.
Precision Engineering
,
43
,
230
240
.. .

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact [email protected]