# Inverse Identification of Virtual Material Parameters Using Surface Response Methodology

Authors: Fangyuan Cui; Dengxin Hua; Pengyang Li; Shikun Lu; Yan Li; Lingfei Kong
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The virtual material model is now widely applied for modeling the dynamical performance of assembled structures since it can effectively represent the complicated contact behavior of joint interfaces despite being relatively simple to create. In this study, a virtual material model is adopted for modeling the dominant physics of a bolted joint subject to a set of pretightening conditions. The unknown virtual material parameters are acquired by an inverse identification procedure that uses the surface response methodology. The greatest advantage of this approach is the ease with which it acquires the joint parameters without taking apart a built-up structure to do special measurements on each separated component. Intricate theoretical calculations can also be avoided when this method is used. This study addresses the responses of virtual material parameters under different pretightening considerations. Predictions based on the identified virtual material parameters are compared with the corresponding results obtained using the analytical method. The correlation between the two sets of results at all preload levels is promising, which indicates the successful identification of the virtual material parameters.


---

## Body

## 1. Introduction

Most mechanical structures have to be connected together using various types of joints. It is well known that the structural discontinuities introduced by joints usually induce a considerable decrease in stiffness and increase in damping locally and thus cause changes in dynamic characteristics like the natural frequencies of the full-body model of a structure [1]. Therefore, the proper modeling and accurate determination of the relevant joint parameters are critical to the vibration control and performance prediction in the successful design of a mechanical structure.Research in these areas can be primarily divided into two streams: fundamental investigations at the microscopic level and engineering applications at the macroscopic level. Studies of the former type have resulted in the development of several contact models for microsurface asperities. One of the earliest models, based on the classical Hertz contact theory, is the Greenwood and Williamson (GW) model [2]. The GW model is actually a pure elastic contact model based on one asperity without considering the elastic-plastic deformation. Several studies were undertaken to overcome this drawback and a series of extensive asperity contact models describing elastic-plastic and complete plastic deformation were established [3–5]. However, few researchers have concerned themselves with applying these micromechanical models to the analysis and prediction of the behavior of actual mechanical structures with joint interfaces.On the other hand, studies belonging to the latter group often rely on simulations for the modeling of joint interfaces, especially those used for engineering applications. At this time, finite element analysis (FEA) technique including modeling and simulation is one of the most commonly used methods for studying the dynamic characteristics of mechanical structures [6]. Several common methods in FEA involve the use of beam element [7], Iwan’s model [8, 9], and spring-damper system [10, 11]. Among these methods, the use of spring-damper system for describing the characteristics of joint interface is the preferred approach so far. The spring-damper-based model is usually constructed by a quasi-static method and can visually represent the static characteristics of the joints. However, the stiffness coefficient used in such models is the static contact stiffness which is independent of the dynamic properties of the structure. Moreover, the application of this method is limited by the difficulties in determining accurately the number of connection nodes and the distribution of their positions, which heavily influence the results of the FE simulation.For several decades, researchers have been attempting to create predictive FE models for jointed interfaces. The emergence of the thin-layer interface concept [13] provided an alternative way for modeling and capturing the dominant physics of the joint interface regardless of the above shortcomings. In this method, a fixed joint interface between substructures can be modeled as a thin layer of elastic-plastic material, which replicates the performance of the contact interface. By tuning the effective material properties in this narrow band, the final updated FE model will represent its physical counterpart more closely [14]. Such thin-layer elements can be applied in two ways, i.e., for linear performance analyses [15, 16] and for nonlinear characteristic analyses [17]. Satisfactory results were reported when this approach was applied to the study of different types of mechanical joints such as spot welds joints [18], bolted lap joints [19, 20], and flanged joints [16]. Besides, this approach is capable of modeling the nonlinear behavior of joint interface parameters, for example, the damping ratios [21].The greatest challenge in modeling the thin layer is the determination of the unknown parameters for the added material layer. Generally speaking, direct and indirect approaches may be applied to extract these unknown parameters [22, 23]. The first one determines the parameters directly in an analytic manner. The difficulty of this task lies in identifying what exactly to model, since current knowledge about the intrinsic complexities of the behavior of the joint is severely limited [24, 25]. On the other hand, the indirect identification approach based on modal data is an inverse method to compensate for the incompleteness of the assumptions in the modeling process. Thus, the second method has been receiving increasing attraction and a significant amount of literature is available. These studies are mainly concerned with parameters such as Young’s and shear moduli [16, 18]. However, the contact behavior of joints were also found to be governed by the geometrical properties of the joint such as the thickness of the thin layer [26]. Iranzad and Ahmadian [27] have assumed the interface layer to have a thickness equal to the distance between the neutral axes of the two mating 2D lap beams. Other authors [14] have assumed the thickness of the layer to be zero to avoid the introduction of extra mass to the full-body model of the assembly.Based on the literature reviewed above, it can be pointed out that there is still no clear strategy for determining the optimum thickness of the thin-layer elements. The thickness of the thin layer is, in fact, governed by many other factors such as the surface roughness, surface medium film, geometrical dimensions of the joints, and preloading conditions. Besides, not much research effort has been focused on validating the effectiveness of the results that have been obtained.Adopting a similar strategy for the thin-layer theory, Tian [28] proposed a virtual material hypothesis-based model for the first time and gave a set of analytical expressions for the virtual material parameters by applying the classical Hertz contact theory together with the fractal theory. This parametric model is applicable to any fixed mechanical joint interface and can provide good correlation with test results [12]. Nevertheless, using analytical methods for acquiring virtual material parameters is usually a difficult process, because, apart from the intrinsic complexity of the analytical method, incorrectly chosen parameter values may produce results that are totally wrong. Accordingly, one should seek a more reliable procedure to identify these joint parameters.Motivated by [29], a virtual material model was developed for modeling the contact behavior of bolted joints. The joint parameters were identified using an inverse approach based on the surface response method. The aim of this study is to provide an accurate and simple procedure for the identification of virtual material parameters in the predictive model of bolted structures under various pretightening conditions. In addition, detailed analyses of the relationships between the objective function used in the optimization procedure and the updating parameters were conducted. The predictions using the analytical method were used as reference models to verify the success of the identified values.The contents of this paper are organized as follows. A brief introduction to the virtual material theory along with its limitations is provided in Section2. In Section 3, the methodology used in this study is presented in detail. This is followed by details of the FE model and the parameters of the test specimen. The results of all the case studies are presented in Section 5. Sections 6 and 7 show all the summaries and conclusions.


## 2. Virtual Material Theory

In Figure1, the picture on the left presents the schematic of an assembly consisting of a flexible joint interface. It is observed that, at the microscopic level, the joint interface is composed of a large number of asperities formed by peaks and troughs whose radii of curvature vary widely. Due to the intermittent contact of the interfacial asperities, the flexible joint interface will exhibit a stiffness softening phenomenon. According to the virtual material theory, adding a layer of virtual material to the assembled structures can generate a relatively simple model without any complicated joint interface. The picture on the right shows the simplified model based on the virtual material theory. The so-called virtual material is modeled as a new component rigidly bonded to the two mating substructures. Two sets of material properties are assigned for the FE model of the assembly: one pertaining to the new component to resemble the contact stiffness of the joint interface and the other for the rest of the components to simulate the stiffness of the substructures. However, the parameters for this new virtual material component are unknown and need to be identified. The key concern in this theory is the determination of the four virtual material parameters, namely, Young’s modulus, Poisson’s ratio, thickness, and density.Figure 1
Schematic diagram of an assembly consisting of a flexible joint interface (a) and its simplified model with a virtual material layer (b).
(a) (b)


### 2.1. Theoretical Formulations of Virtual Material Model

The unknown physical and geometrical parameters that must be identified are presented in(1)E=EPn,E1,E2,ν1,ν2,Ra1,Ra2(2)ν=νPn,E1,E2,ν1,ν2,Ra1,Ra2(3)ρ=ρρ1,ρ2(4)h=hh1,h2in whichE, ν,ρ, andhare Young’s modulus, Poisson’s ratio, density, and thickness of the virtual material, respectively; the first term Pn in (1)-(2) refers to the normal load applied on the joint interface; E1 and E2 indicate Young’s modulus of substructures 1 and 2, respectively; similarly, ν1, ν2, Ra1, Ra2,ρ1, andρ2 are the respective Poisson ratios, surface roughness values, and densities. The remaining parameters, i.e., h1 and h2, denote the thickness of asperity layers in parts 1 and 2.By adopting the Hertz contact theory and the fractal theory, Tian [28] stated that Young’s modulusE for the virtual material is as follows:(5)E=E∗E′=2Dψ1-0.5D3π2G∗1-D2-DDψ0.5D-1Ar∗0.5Dac∗-0.5-2-DDψ0.5D-1Ar∗-0.5E′where E∗ is the dimensionless Young modulus of the equivalent virtual material, E′ is the equivalent Young modulus of the two mating substructures,D, ψ, andG are the fractal parameters determined by surface roughness values Ra1 and Ra2, ac∗ is the dimensionless critical area of transformation from elastic to plastic deformation, and Ar∗ is the dimensionless real contact area, which is closely related to the normal load Pn.The final form of Poisson’s ratio of the virtual material is given by(6)ν=1+ν′E∗Gx∗-1where ν′ is the equivalent Poisson ratio of the two mating substructures and Gx∗ is the dimensionless shear modulus of the equivalent virtual material.Assume the sum of the thickness of the asperity layers of the two mating substructures areh1 and h2, masses m1 and m2, and volumes V1 and V2, respectively. If Ac denotes the contact area of the joint interface, then the density of the virtual material layer can be given by(7)ρ=m1+m2V1+V2=ρ1Ach1+ρ2Ach2Ach1+h2=ρ1h1+ρ2h2h1+h2It is worth noting thath1 and h2 (in (7)) are microscopic parameters which cannot be measured accurately.The thickness of the virtual material layer is usually set at a constant value of 1 mm in this model.


### 2.2. Problem Statement

Due to the effect of machining methods used, the surface microtopography of mechanical joints is uneven to varying extents. Figure2 presents a representative visualization of the topography of a machined surface. As can be seen, the structure of the microsurface is a spatial region composed of gas, liquid, and solid particles. The entire microstructure, from the outermost layer to the innermost, can be divided into the following: a film of adsorbed gas-liquid molecules, an oxide film, a Beilby layer, a layer of considerable deformation, and a layer of slight deformation. The virtual material method assumes that the dominant performance parameters of the joint interface, like the contact stiffness, are related only to these asperity layers. During the modeling of joints, therefore, the thickness of virtual materialh is defined as the sum of the thicknesses of the asperity layers in the two mating substructures. According to literature [28], when the surface roughness value of the contact surface is 0.8 um, the approximate value of the thicknessh is 1 mm. So, the virtual material is assigned as a thickness of 1 mm during modeling.Figure 2
Microstructure of asperities on a typical machined surface.However, different machining methods and lubrication conditions will give rise to diverse surface roughness which further influence the thickness of the asperity layers. Small loads usually cause infinitesimal deformations in the asperities while severe bulk deformation will occur under large normal loads. Chang [30] proposed that the bulk substrate contributes significantly to the overall deformation in addition to the effect of asperity deformation. Hence, metal surfaces that have been subjected to different processing conditions usually have asperity layers of different thicknesses. For example, in the case of a grinding surface with a surface roughness Ra = 0.8um, which is the same as that used in literature [28], the total measured thickness of the asperity layers is approximately 0.26 mm [31].While deriving the theoretical formulas, Tian [28] did not take the bulk deformation into consideration. Moreover, the asperity layers mentioned above differ in their physical, chemical, and mechanical properties, which means that the densities of the virtual material may also differ to a certain extent. Thus, the available theoretical expressions that do not consider the bulk deformation or density variation may not fully describe the contact behavior of the joint interface in a satisfactory manner. In other words, if the two aspects are also considered, the calculation results will be more representative of the effects of surface microphenomenon on the dynamic characteristics of the structure.


## 2.1. Theoretical Formulations of Virtual Material Model

The unknown physical and geometrical parameters that must be identified are presented in(1)E=EPn,E1,E2,ν1,ν2,Ra1,Ra2(2)ν=νPn,E1,E2,ν1,ν2,Ra1,Ra2(3)ρ=ρρ1,ρ2(4)h=hh1,h2in whichE, ν,ρ, andhare Young’s modulus, Poisson’s ratio, density, and thickness of the virtual material, respectively; the first term Pn in (1)-(2) refers to the normal load applied on the joint interface; E1 and E2 indicate Young’s modulus of substructures 1 and 2, respectively; similarly, ν1, ν2, Ra1, Ra2,ρ1, andρ2 are the respective Poisson ratios, surface roughness values, and densities. The remaining parameters, i.e., h1 and h2, denote the thickness of asperity layers in parts 1 and 2.By adopting the Hertz contact theory and the fractal theory, Tian [28] stated that Young’s modulusE for the virtual material is as follows:(5)E=E∗E′=2Dψ1-0.5D3π2G∗1-D2-DDψ0.5D-1Ar∗0.5Dac∗-0.5-2-DDψ0.5D-1Ar∗-0.5E′where E∗ is the dimensionless Young modulus of the equivalent virtual material, E′ is the equivalent Young modulus of the two mating substructures,D, ψ, andG are the fractal parameters determined by surface roughness values Ra1 and Ra2, ac∗ is the dimensionless critical area of transformation from elastic to plastic deformation, and Ar∗ is the dimensionless real contact area, which is closely related to the normal load Pn.The final form of Poisson’s ratio of the virtual material is given by(6)ν=1+ν′E∗Gx∗-1where ν′ is the equivalent Poisson ratio of the two mating substructures and Gx∗ is the dimensionless shear modulus of the equivalent virtual material.Assume the sum of the thickness of the asperity layers of the two mating substructures areh1 and h2, masses m1 and m2, and volumes V1 and V2, respectively. If Ac denotes the contact area of the joint interface, then the density of the virtual material layer can be given by(7)ρ=m1+m2V1+V2=ρ1Ach1+ρ2Ach2Ach1+h2=ρ1h1+ρ2h2h1+h2It is worth noting thath1 and h2 (in (7)) are microscopic parameters which cannot be measured accurately.The thickness of the virtual material layer is usually set at a constant value of 1 mm in this model.


## 2.2. Problem Statement

Due to the effect of machining methods used, the surface microtopography of mechanical joints is uneven to varying extents. Figure2 presents a representative visualization of the topography of a machined surface. As can be seen, the structure of the microsurface is a spatial region composed of gas, liquid, and solid particles. The entire microstructure, from the outermost layer to the innermost, can be divided into the following: a film of adsorbed gas-liquid molecules, an oxide film, a Beilby layer, a layer of considerable deformation, and a layer of slight deformation. The virtual material method assumes that the dominant performance parameters of the joint interface, like the contact stiffness, are related only to these asperity layers. During the modeling of joints, therefore, the thickness of virtual materialh is defined as the sum of the thicknesses of the asperity layers in the two mating substructures. According to literature [28], when the surface roughness value of the contact surface is 0.8 um, the approximate value of the thicknessh is 1 mm. So, the virtual material is assigned as a thickness of 1 mm during modeling.Figure 2
Microstructure of asperities on a typical machined surface.However, different machining methods and lubrication conditions will give rise to diverse surface roughness which further influence the thickness of the asperity layers. Small loads usually cause infinitesimal deformations in the asperities while severe bulk deformation will occur under large normal loads. Chang [30] proposed that the bulk substrate contributes significantly to the overall deformation in addition to the effect of asperity deformation. Hence, metal surfaces that have been subjected to different processing conditions usually have asperity layers of different thicknesses. For example, in the case of a grinding surface with a surface roughness Ra = 0.8um, which is the same as that used in literature [28], the total measured thickness of the asperity layers is approximately 0.26 mm [31].While deriving the theoretical formulas, Tian [28] did not take the bulk deformation into consideration. Moreover, the asperity layers mentioned above differ in their physical, chemical, and mechanical properties, which means that the densities of the virtual material may also differ to a certain extent. Thus, the available theoretical expressions that do not consider the bulk deformation or density variation may not fully describe the contact behavior of the joint interface in a satisfactory manner. In other words, if the two aspects are also considered, the calculation results will be more representative of the effects of surface microphenomenon on the dynamic characteristics of the structure.


## 3. Methodology

Figure3 presents the forward and inverse approaches to access the virtual material parameters. In the forward parameter recognition method, it is desirable to first prepare specimens for each of the separated components and measure their surface roughness. Then, based on the applied preload as well as the material properties of the two mating substructures, the parameters of the unknown virtual material layer can be calculated based on (5) to (7). The final step is the incorporation of these theoretically calculated virtual material parameters into the FEA software to obtain the theoretical predictions. In order to verify its effectiveness, a modal test is performed on the whole assembly to compare and validate these theoretical predictions.Figure 3
The forward and inverse methods to access the virtual material parameters.The key issues in the forward parameter acquisition method is the determination of fractal parameters such as the fractal dimensionDand the fractal roughness parameterG which are essential for the theoretical calculation of the virtual material parameters. The determination of these two basic parameters requires special measurements on each separated component. In addition, it was found that these two parameters are also dependent on the means of detection and the resolution of the measuring instrument. This means that once there is any change in the joint interface, all these parameters need to be redefined. Thus, the correct determination of these values is of great importance, since they significantly influence the results of the final calculation. These inconveniences restrict the usage of the direct forward determination method for joint parameters in engineering application since the large number of parameters involved necessitates extensive computation. Hence, an alternative way for the determination of joint parameters would be of immense use.At the present time, the inverse identification approach based on modal data is the preferred one and much literature has been published in this connection [32, 33]. Therefore, this method of identification is introduced as an alternative method for the determination of the virtual material parameters. In comparison with the direct calculation method, the inverse process is much easier and more universal where the assessment of joint parameters is concerned. The main advantage of this procedure is the ease of acquisition of virtual material parameters without taking apart a built-up structure to perform special measurements on each separated part, and a large amount of computational effort can be saved.In this method, the only requirement for the identification of joint parameters is the modal data of the whole assembly, which is easy to access. The experimental modal data are used in the optimization procedure with an aim to develop a more reliable FE model. The virtual material is assumed to be linear isotropic and, thus, can be characterized using the four joint parameters namely, Young’s modulus, Poisson’s ratio, thickness, and density. These joint parameters for the virtual material model can be modified so that the influence of the deviation between the FE predictions and experimental observations on the objective function is minimized. The optimization of the objective function is carried out by the use of the surface response optimization module available in ANSYS 15.0.In this study, the difference between the eigenvalues obtained experimentally and numerically is taken as the objective function for the optimization:(8)objfun=min∑i6ωiFE-ωiEXPωiEXP2where ωiFE and ωiEXP are, respectively, the analytically determined and experimentally measured natural frequencies. The first six measured natural frequencies will be used to identify the unknown parameters of the virtual material model by tuning the parameters present in the FE model so that a good agreement between the predictions of the updated FE model and the measured data is achieved.Minimization of the objective function (8) is carried out using an iterative linear eigensensitivity method with the following expression [16, 18]:(9)ωiFE2-ωiEXP2=ϕiFETΔK-ωiFE2ΔMϕiFEϕiFETMϕiFEwhere ϕiFE is the ith analytical mode shape;M is the mass matrix; and ΔK and ΔM are the variations in the stiffness and mass matrices.


## 4. FE Model

The application of the proposed parameter acquisition method for the joint is demonstrated by a test specimen with bolted joints. The geometric characteristics of the experimental specimen [12] are shown in Figure 4. The two mating substructures are coupled using two M12 inner-hexagonal bolts. The contact between the two substructures is modeled as a dry contact. The joint interface has a rectangular area of 80 × 30 mm2. For reasons of simplicity, the bolts and bolt holes were ignored in the FE model used in this study.Figure 4
Geometry of the specimen (dimensions in mm).The FE model of the specimen was developed using ANSYS 15.0 software to predict its dynamical behavior. Figure5 shows the configuration of the elements in each subzone for the meshed FE model. A thin layer of virtual material elements with an initial thickness of 1 mm was introduced for modeling the interface of the bolted joint, shown in subzone 3. The dynamical behavior of the FE model is assumed to be controlled by this narrow region whose material and geometrical property parameters can be modified.Figure 5
Arrangement of the elements in each subzone.Since the number and position of the nodes and elements of the meshed model are involved in the FE solution, the element quality has a direct effect on the accuracy of the dynamical response predictions of the model. The geometry of this test structure presents difficulties in using 3D hexahedral elements, which are generally more accurate than tetrahedral finite elements, for meshing. For this reason, the 3D mesh is constructed by dividing the whole structure into five segments as shown in Figure5, with bonded contacts to make a continuous mesh. The sweep method is adopted for subzone 2, subzone 3, and subzone 4 while the multizone method is employed for subzone 1 and subzone 5. The mesh metric provided by ANSYS 15.0 helps in quantifying the mesh quality and evaluating the mesh model. The average value of the mesh metric ranges from 0 to 1. The higher the average value of the mesh metric, the better the mesh quality. As can be seen from Table 1, the average value of mesh metric (0.99813) for the model used is very close to 1, which indicates the mesh refinement is sufficient to ignore the influence of mesh quality on the simulation results.Table 1
Statistics for mesh results.
Nodes Elements Mesh metric Min Max Average Standard deviation 145961 33104 0.35667 0.99999 0.99813 0.03459The main emphasis in the present study is placed on the validation of the effectiveness of the proposed virtual material parameter acquisition method. So, the theoretical virtual material parameters and the experimental work reported in literature [12] were considered for the confirmation of the validity and accuracy of the method proposed in this paper.The material properties for the two mating substructures are presented in Table2. Based on the parameters given in this table and the corresponding theoretical formulations in Section 2.1, the parameters of the virtual material under different preload conditions can be theoretically calculated, as shown in Table 3.Table 2
Material properties for the two mating substructures [12].
Parameter Material E (GPa) ν ρ (kg/m3) R a (μm) σ y (MPa) H (MPa) Substructure 1 HT250 116 0.27 7340 3.2 240 700 Substructure 2 HT250 116 0.27 7340 3.2 240 700Table 3
Theoretical parameters of virtual material at three different preload conditions [12].
TorqueT (N•m) E (GPa) ν ρ (kg/m3) 30 0.631 0.22 7340 60 0.736 0.24 7340 90 1.38 0.27 7340In order to investigate the effects of the bolt pretightening torque on the joint contact parameters, the assembled structure was tested at three different levels of bolt pretightening torque, i.e., 30 N·m, 60 N·m, and 90 N·m. The tightening torque was imposed by an electronic torque wrench (SATA 96525) with a precision of ±2%. Figure 6 shows the experimental setup used. It can be seen that the test setup includes the test specimen, an impulse hammer, seven accelerometers, the LMS Test. Lab software, and a computer. The excitation was done by the impulse hammer, which also measured the exerted force through the force sensor mounted on the hammer tip. Five impacts were applied at each excitation point to ensure reliability. The response signals of the structure were acquired using the accelerometers. Data processing and analysis were handled by the LMS Test. Lab software. The measured natural frequencies of the test specimen under three conditions were tabulated in Table 4.Table 4
Measured natural frequencies at three different preload conditions [12].
TorqueT (N•m) 1st 2nd 3rd 4th 5th 6th 30 354.54 377.89 620.31 1600.53 1610.66 1760.12 60 363.04 400.75 649.05 1611.49 1626.97 1801.2 90 390.44 489.49 733.72 1662.74 1705.49 2026.8Figure 6
Experimental setup used [12].In the next section, the parameters of the virtual material under three different preload cases will be identified based on the experimental results presented in Table4.


## 5. Case Study

### 5.1. Identification of Joint Parameters at Pretightening Torque of 30 N·m

Initially, the virtual material parameters at a relatively small preload will be identified using the modal parameters pertaining to the 30 N·m pretightening torque in Table4.In order to include the effects of thickness and density on the predictions from the FE model, four design variables were selected in this paper, namely, Young’s modulus, Poisson’s ratio, density, and thickness of the virtual material. As explained in Section3, the deviations between the measured and predicted natural frequencies were assumed to be related to those four parameters.Table5 shows the four design variables and their permissible ranges of variation used in the optimization procedure. At the beginning of the procedure, the material properties shown in Table 2 were assigned to those design variables as their initial values. Note that all the case studies have identical initial values. The choice of permissible range of variation for Young’s modulus was made based on observations [4, 29] indicating that it decreased sharply by about 3-4 orders of magnitude compared with its initial values, after the optimization procedure was conducted. Poisson’s ratio changed within the range of 0.2 to 0.3. The lower and upper bounds of the rest of the parameters were set to be 0.9 and 1.01 times the initial values, respectively, according to the default setting in the surface response optimization procedure in ANSYS 15.0.Table 5
Design variables and their ranges (at 30 N·m pretightening torque).
Design variables E (GPa) ν ρ (kg/m3) h(mm) Initial value 116 0.27 7340 1 Lower bound 0.6 0.2 6606 0.9 Upper bound 0.7 0.3 8074 1.1Figure7 presents the variations of the objective function with any pair of the four design variables. It can be seen that the objective function experiences changes between 0.017 and 0.045. Also, it can be observed that there are optimal values for Young’s modulus and thickness (see Figures 7(a) and 7(e)). With the increase in Poisson’s ratio, the value of the objective function will also increase. A higher value of Poisson’s ratio introduces larger errors between the predictions made by the FE model and the experiment (see Figures 7(b) and 7(d)). In contrast to the other three updated parameters, the density shows minor differences as the objective function changes (see Figures 7(a), 7(d), and 7(e)), which seems reasonable since the contact type assigned in this case study is dry friction.Figure 7
Changes in the objective function with multiple design variables at the pretightening torque of 30 N·m.
(a) (b) (c) (d) (e) (f)The convergence history of the objective function is shown in Figure8. There are 100 design points to be evaluated in each iteration. This figure indicates that, after seven iterations, the objective function converges to the minimum target values. So, the total number of design points considered during the minimization process is 700.Figure 8
Convergence history of objective function along the optimization process.The virtual material parameter values identified in the final optimization stage are listed in Table6. It is observed that, after the optimization, the reduction of Young’s modulus in the joint region is about 185 times lower than its initial value. Similarly, there is a 26% decrease in Poisson’s ratio. In contrast, the density and thickness experience a very small diminution (4% and 7%, respectively) compared with their initial values. These results indicate that Young’s modulus is more significant for the stiffness softening behavior of joints than the other design variables.Table 6
Identified parameters of the virtual material (at 30 N·m pretightening torque).
Parameter E (GPa) ν ρ (kg/m3) h (mm) Initial value 116 0.27 7340 1 Final value 0.626 0.2 7042 0.935Comparisons of natural frequencies from the test (measured) and the analytical method (reference) together with the updated FE model (updated) at a pretightening torque of 30 N•m are shown in Figure9. Comparing the three sets of results, one notices that the predictions from the updated model using the proposed method agree quite well with theoretically obtained results. In addition, the sum of the percentage errors of all the modes obtained from the updated model (Error 2) is a little smaller than the sum in the reference model (Error 1). This indicates that the updating process can facilitate good agreement between the updated FE model and the analytical predictions.Figure 9
Comparisons of natural frequencies from the test (measured), the analytical model (reference), and the updated FE model (updated) at a pretightening torque of 30 N·m.


### 5.2. Identification of Joint Parameters at Pretightening Torque of 60 N·m

In this case, the experimentally measured natural frequencies for a medium pretightening torque of 60 N·m, shown in Table4, were taken as the target values for the parameter identification of the virtual material.Table7 presents the initial values and permissible ranges of variation for the relevant design variables, whose selection criteria are the same as explained in Section 5.1.Table 7
Design variables and their ranges (at 60 N·m pretightening torque).
Design variables E (GPa) ν ρ (kg/m3) h (mm) Initial value 116 0.27 7340 1 Lower bound 0.7 0.2 6606 0.9 Upper bound 0.8 0.3 8074 1.1Variations of the objective function with any pair of parameters are shown in Figure10. It can be observed that the objective function experiences a small variation from 0.013 to 0.035. Also, there are optimal values for Young’s modulus and thickness (see Figures 10(a) and 10(e)). But the optimal Young modulus in this case is a little higher than that in the previous case. This is because an increase in the pretightening torque will lead to an increase in the contact stiffness, which results in a higher Young’s modulus for the virtual material. Poisson’s ratio shows a similar trend to that in the previous case; i.e., one will obtain the minimum value for the objective function at a smaller value of Poisson’s ratio (see Figures 10(b) and 10(d)). The effect of density on the objective function during the optimization process is not obvious (see Figures 10(a), 10(d), and 10(e)).Figure 10
Changes in objective function with multiple design variables at the pretightening torque of 60 N·m.
(a) (b) (c) (d) (e) (f)Figure11 shows the variation in the objective function when the four selected design variables vary simultaneously. Obviously, convergence of the objective function occurs after seven iterations with 100 design points being assessed in every iteration.Figure 11
Convergence history of the objective function during optimization.The final virtual material parameters obtained are listed in Table8. It is interesting to note that Young’s modulus confirms the fact that a higher preload gives rise to a higher contact stiffness. Young’s modulus in this case has decreased from the initial value of 116 GPa by three orders of magnitude to 0.798 GPa. The value of Poisson’s ratio decreased from 0.27 to 0.2. However, the changes in density (7340 kg/m3 to 7091 kg/m3) and thickness (1 mm to 1.026 mm) are only marginal.Table 8
Identified parameters of the virtual material (at 60 N·m pretightening torque).
Parameter E (GPa) ν ρ (kg/m3) h (mm) Initial value 116 0.27 7340 1 Final value 0.798 0.2 7091 1.026Figure12 shows a bar chart of the experimental natural frequencies (measured), analytical predictions from the reference model (reference), and updated FE model (updated) at a pretightening torque of 60 N·m. In this figure, the comparison results between the test natural frequencies and the analytical model are denoted by Error 1. Similarly, Error 2 refers to the difference between the natural frequencies of the updated FE model and the measured values. As can be seen, Error 2 values for all the models lie within ±10%. Moreover, the summation of the percentage errors in the updated model is a little smaller than that in the reference model, which indicates the successful identification of the virtual material parameters.Figure 12
Comparisons of natural frequencies from the test (measured), the analytical model (reference) and the updated FE model (updated) at a pretightening torque of 60 N·m.


### 5.3. Identification of Joint Parameters at Pretightening Torque of 90 N·m

In order to explore the relationships between the high preload and virtual material parameters, the case of a pretightening torque of 90 N·m was considered. Again, the test results corresponding to 90 N·m, presented in Table4, were chosen as the target model for the simulation.Table9 shows the design variables and their ranges in this case. The upper and lower limits for Poisson’s ratio, density, and thickness are the same as those in the previous cases while Young’s modulus ranges from a lower value of 1.0 GPa to an upper value of 2.0 GPa.Table 9
Design variables and their ranges (at 90 N·m pretightening torque).
Design variables E (GPa) ν ρ (kg/m3) h (mm) Initial value 116 0.27 7340 1 Lower bound 1.0 0.2 6606 0.9 Upper bound 2.0 0.3 8074 1.1The relationships between the objective function and any given pair of the updating parameters are shown in Figure13. The value of the objective function ranges from 0.006 to 0.036. Figures 13(a)–13(c) indicate that the optimal Young modulus is achieved at the vicinity of 1.25 GPa. By tuning Poisson’s ratio between 0.2 and 0.3, the objective function first decreased and then increased. Thus, the optimal Poisson ratio close to the value of 0.025 (see Figures 13(d) and 13(f)) corresponds to the minimum value of the objective function. It can be observed from Figures 13(a) and 13(d) that the change in density induces only a weak response in the objective function. As for Figures 13(c), 13(e), and 13(f), one can notice that the greater the thickness, the smaller the objective function.Figure 13
Changes in objective function with multiple design variables at the pretightening torque of 90 N·m.
(a) (b) (c) (d) (e) (f)The optimal parameters for the virtual material are identified by the minimization of the objective function. Figure14 presents the minimization process of the objective function. The figure shows that the result converges after eight iterations where 100 design points are evaluated per iteration.Figure 14
Convergence history of the objective function along the optimization process.The values of the virtual material parameters, obtained by minimization of the objective function, are given in Table10. It can be clearly observed that Young’s modulus shows a noticeable change from the initial value of 116 GPa to the final value of 1.246 GPa. The initial value of Young’s modulus reduces by almost two orders of magnitude. In this case, we obtained an optimal value (0.226) for Poisson’s ratio. The identified value of density (7193.1) is also smaller than its initial value (7340) as found in previous cases while the identified value of thickness (1.077) is a little larger than the initial value of 1 mm.Table 10
Identified parameters of the virtual material (at 90 N·m pretightening torque).
Parameter E (GPa) ν ρ (kg/m3) h (mm) Initial value 116 0.27 7340 1 Final value 1.246 0.226 7193.1 1.077The percentage errors in the natural frequencies after updating the FE model with the identified virtual material parameters (Error 2) are presented in Figure15. Contrary to the Error 1 values in the reference model, the Error 2 values obtained using the identified virtual material parameters are smaller, which means the new method agrees better with the experimental results.Figure 15
Comparisons of natural frequencies from the test (measured), the analytical model (reference) and the updated FE model (updated) at 90 N·m pretightening torque.


## 5.1. Identification of Joint Parameters at Pretightening Torque of 30 N·m

Initially, the virtual material parameters at a relatively small preload will be identified using the modal parameters pertaining to the 30 N·m pretightening torque in Table4.In order to include the effects of thickness and density on the predictions from the FE model, four design variables were selected in this paper, namely, Young’s modulus, Poisson’s ratio, density, and thickness of the virtual material. As explained in Section3, the deviations between the measured and predicted natural frequencies were assumed to be related to those four parameters.Table5 shows the four design variables and their permissible ranges of variation used in the optimization procedure. At the beginning of the procedure, the material properties shown in Table 2 were assigned to those design variables as their initial values. Note that all the case studies have identical initial values. The choice of permissible range of variation for Young’s modulus was made based on observations [4, 29] indicating that it decreased sharply by about 3-4 orders of magnitude compared with its initial values, after the optimization procedure was conducted. Poisson’s ratio changed within the range of 0.2 to 0.3. The lower and upper bounds of the rest of the parameters were set to be 0.9 and 1.01 times the initial values, respectively, according to the default setting in the surface response optimization procedure in ANSYS 15.0.Table 5
Design variables and their ranges (at 30 N·m pretightening torque).
Design variables E (GPa) ν ρ (kg/m3) h(mm) Initial value 116 0.27 7340 1 Lower bound 0.6 0.2 6606 0.9 Upper bound 0.7 0.3 8074 1.1Figure7 presents the variations of the objective function with any pair of the four design variables. It can be seen that the objective function experiences changes between 0.017 and 0.045. Also, it can be observed that there are optimal values for Young’s modulus and thickness (see Figures 7(a) and 7(e)). With the increase in Poisson’s ratio, the value of the objective function will also increase. A higher value of Poisson’s ratio introduces larger errors between the predictions made by the FE model and the experiment (see Figures 7(b) and 7(d)). In contrast to the other three updated parameters, the density shows minor differences as the objective function changes (see Figures 7(a), 7(d), and 7(e)), which seems reasonable since the contact type assigned in this case study is dry friction.Figure 7
Changes in the objective function with multiple design variables at the pretightening torque of 30 N·m.
(a) (b) (c) (d) (e) (f)The convergence history of the objective function is shown in Figure8. There are 100 design points to be evaluated in each iteration. This figure indicates that, after seven iterations, the objective function converges to the minimum target values. So, the total number of design points considered during the minimization process is 700.Figure 8
Convergence history of objective function along the optimization process.The virtual material parameter values identified in the final optimization stage are listed in Table6. It is observed that, after the optimization, the reduction of Young’s modulus in the joint region is about 185 times lower than its initial value. Similarly, there is a 26% decrease in Poisson’s ratio. In contrast, the density and thickness experience a very small diminution (4% and 7%, respectively) compared with their initial values. These results indicate that Young’s modulus is more significant for the stiffness softening behavior of joints than the other design variables.Table 6
Identified parameters of the virtual material (at 30 N·m pretightening torque).
Parameter E (GPa) ν ρ (kg/m3) h (mm) Initial value 116 0.27 7340 1 Final value 0.626 0.2 7042 0.935Comparisons of natural frequencies from the test (measured) and the analytical method (reference) together with the updated FE model (updated) at a pretightening torque of 30 N•m are shown in Figure9. Comparing the three sets of results, one notices that the predictions from the updated model using the proposed method agree quite well with theoretically obtained results. In addition, the sum of the percentage errors of all the modes obtained from the updated model (Error 2) is a little smaller than the sum in the reference model (Error 1). This indicates that the updating process can facilitate good agreement between the updated FE model and the analytical predictions.Figure 9
Comparisons of natural frequencies from the test (measured), the analytical model (reference), and the updated FE model (updated) at a pretightening torque of 30 N·m.


## 5.2. Identification of Joint Parameters at Pretightening Torque of 60 N·m

In this case, the experimentally measured natural frequencies for a medium pretightening torque of 60 N·m, shown in Table4, were taken as the target values for the parameter identification of the virtual material.Table7 presents the initial values and permissible ranges of variation for the relevant design variables, whose selection criteria are the same as explained in Section 5.1.Table 7
Design variables and their ranges (at 60 N·m pretightening torque).
Design variables E (GPa) ν ρ (kg/m3) h (mm) Initial value 116 0.27 7340 1 Lower bound 0.7 0.2 6606 0.9 Upper bound 0.8 0.3 8074 1.1Variations of the objective function with any pair of parameters are shown in Figure10. It can be observed that the objective function experiences a small variation from 0.013 to 0.035. Also, there are optimal values for Young’s modulus and thickness (see Figures 10(a) and 10(e)). But the optimal Young modulus in this case is a little higher than that in the previous case. This is because an increase in the pretightening torque will lead to an increase in the contact stiffness, which results in a higher Young’s modulus for the virtual material. Poisson’s ratio shows a similar trend to that in the previous case; i.e., one will obtain the minimum value for the objective function at a smaller value of Poisson’s ratio (see Figures 10(b) and 10(d)). The effect of density on the objective function during the optimization process is not obvious (see Figures 10(a), 10(d), and 10(e)).Figure 10
Changes in objective function with multiple design variables at the pretightening torque of 60 N·m.
(a) (b) (c) (d) (e) (f)Figure11 shows the variation in the objective function when the four selected design variables vary simultaneously. Obviously, convergence of the objective function occurs after seven iterations with 100 design points being assessed in every iteration.Figure 11
Convergence history of the objective function during optimization.The final virtual material parameters obtained are listed in Table8. It is interesting to note that Young’s modulus confirms the fact that a higher preload gives rise to a higher contact stiffness. Young’s modulus in this case has decreased from the initial value of 116 GPa by three orders of magnitude to 0.798 GPa. The value of Poisson’s ratio decreased from 0.27 to 0.2. However, the changes in density (7340 kg/m3 to 7091 kg/m3) and thickness (1 mm to 1.026 mm) are only marginal.Table 8
Identified parameters of the virtual material (at 60 N·m pretightening torque).
Parameter E (GPa) ν ρ (kg/m3) h (mm) Initial value 116 0.27 7340 1 Final value 0.798 0.2 7091 1.026Figure12 shows a bar chart of the experimental natural frequencies (measured), analytical predictions from the reference model (reference), and updated FE model (updated) at a pretightening torque of 60 N·m. In this figure, the comparison results between the test natural frequencies and the analytical model are denoted by Error 1. Similarly, Error 2 refers to the difference between the natural frequencies of the updated FE model and the measured values. As can be seen, Error 2 values for all the models lie within ±10%. Moreover, the summation of the percentage errors in the updated model is a little smaller than that in the reference model, which indicates the successful identification of the virtual material parameters.Figure 12
Comparisons of natural frequencies from the test (measured), the analytical model (reference) and the updated FE model (updated) at a pretightening torque of 60 N·m.


## 5.3. Identification of Joint Parameters at Pretightening Torque of 90 N·m

In order to explore the relationships between the high preload and virtual material parameters, the case of a pretightening torque of 90 N·m was considered. Again, the test results corresponding to 90 N·m, presented in Table4, were chosen as the target model for the simulation.Table9 shows the design variables and their ranges in this case. The upper and lower limits for Poisson’s ratio, density, and thickness are the same as those in the previous cases while Young’s modulus ranges from a lower value of 1.0 GPa to an upper value of 2.0 GPa.Table 9
Design variables and their ranges (at 90 N·m pretightening torque).
Design variables E (GPa) ν ρ (kg/m3) h (mm) Initial value 116 0.27 7340 1 Lower bound 1.0 0.2 6606 0.9 Upper bound 2.0 0.3 8074 1.1The relationships between the objective function and any given pair of the updating parameters are shown in Figure13. The value of the objective function ranges from 0.006 to 0.036. Figures 13(a)–13(c) indicate that the optimal Young modulus is achieved at the vicinity of 1.25 GPa. By tuning Poisson’s ratio between 0.2 and 0.3, the objective function first decreased and then increased. Thus, the optimal Poisson ratio close to the value of 0.025 (see Figures 13(d) and 13(f)) corresponds to the minimum value of the objective function. It can be observed from Figures 13(a) and 13(d) that the change in density induces only a weak response in the objective function. As for Figures 13(c), 13(e), and 13(f), one can notice that the greater the thickness, the smaller the objective function.Figure 13
Changes in objective function with multiple design variables at the pretightening torque of 90 N·m.
(a) (b) (c) (d) (e) (f)The optimal parameters for the virtual material are identified by the minimization of the objective function. Figure14 presents the minimization process of the objective function. The figure shows that the result converges after eight iterations where 100 design points are evaluated per iteration.Figure 14
Convergence history of the objective function along the optimization process.The values of the virtual material parameters, obtained by minimization of the objective function, are given in Table10. It can be clearly observed that Young’s modulus shows a noticeable change from the initial value of 116 GPa to the final value of 1.246 GPa. The initial value of Young’s modulus reduces by almost two orders of magnitude. In this case, we obtained an optimal value (0.226) for Poisson’s ratio. The identified value of density (7193.1) is also smaller than its initial value (7340) as found in previous cases while the identified value of thickness (1.077) is a little larger than the initial value of 1 mm.Table 10
Identified parameters of the virtual material (at 90 N·m pretightening torque).
Parameter E (GPa) ν ρ (kg/m3) h (mm) Initial value 116 0.27 7340 1 Final value 1.246 0.226 7193.1 1.077The percentage errors in the natural frequencies after updating the FE model with the identified virtual material parameters (Error 2) are presented in Figure15. Contrary to the Error 1 values in the reference model, the Error 2 values obtained using the identified virtual material parameters are smaller, which means the new method agrees better with the experimental results.Figure 15
Comparisons of natural frequencies from the test (measured), the analytical model (reference) and the updated FE model (updated) at 90 N·m pretightening torque.


## 6. Summary of Results

In this section the results obtained from all case studies in Section5 are summarized. The identified parameters of the virtual material, listed in Tables 6, 8, and 10, are presented in Figure 16.Figure 16
Variation of the identified virtual material parameters with pretightening torques.As is evident from this line chart, the general trends of the identified parameters indicate the following.(1) The identified Young modulus (E) increases with the increase in bolt pretightening torque, which seems reasonable since, from a microscopic perspective, a higher normal load on the joint interface imposes a larger real contact area among the asperities. Thus, a greater contact stiffness will be obtained on the joint interface.(2) The identified Poisson ratio (ν) of the virtual material shows a gentle upward trend when the bolt tension is varied, which is similar as those theoretical ones in Professor Tian’s papers [12, 28]. A similar interpretation can be made that an increase in the normal load will lead to an increase in the real contact area on the joint interface. Accordingly, the dimensionless Young modulus of the equivalent virtual material E∗ becomes bigger. Equation (6) shows that a bigger E∗ gives rise to a largerν.(3) A further observation reveals that at different bolt pretightening torques, the density of virtual material element varies with a moderate slope. This is indicative of the relatively weaker influence of this design variable on the modal parameters compared with the two material properties mentioned previously.(4) It should be noted that the quality of simulation of the bolted joints depends on geometrical attributes like the thickness of virtual material (h) aside from material properties. The thickness increases slightly with the increase of pretightening torque. In fact, there is an optimum thickness pertaining to the virtual material which is determined by the normal preload.


## 7. Conclusion

Application of virtual material models to represent the face-to-face contact in joints is gaining much popularity nowadays. The methods for the forward calculation of virtual material parameters usually require a large amount of computation, thus limiting their application in engineering. In this work, the virtual material parameters of an FE model were determined by an inverse identification approach. This new approach was demonstrated on an assembly structure having bolted joints, under different preload conditions.The identified results suggest that considering the influence of the normal load on the microstructure of the joint interface may provide more reliable predictions. Also, the relationships between the pretightening torque and the dominant physics of the joints such as stiffness attenuation and damping augmentation can be depicted by using these effective material properties and geometrical attributes parameters.Comparisons of the predictions from the reference models (Forward calculation) and the updated models (Inverse identification) demonstrate that the proposed parameter acquisition method is capable of representing the dominant physics of joints with an acceptable amount of accuracy. Moreover, the comparison results show that when the identified parameters are used, the predictions of the updated model are in better agreement with experiments. Furthermore, parameter identification approach is simple enough to be easily incorporated into existing commercial FE software, which supports the proposed methodology. This study provides a methodology that is useful for macroscale modeling and dynamic characteristic analyses such as the eigenvalue problem of a joint interface.



---

*Source: 7678219-2018-07-18.xml*