Abstract
This paper is devoted to a general methodological study on sensitivity analysis and robust optimization for a planar crank-slider mechanism in presence of joint clearances and random parameters and investigate the effects of parameter uncertainty on optimization results when joint clearance sizes are constantly changing due to wear. The first-order sensitivity analysis based on the response surface proxy model is performed. Then, a multiobjective robust optimization algorithm based on sensitivity analysis is carried out to reduce the undesirable effects of joint clearances and random parameters. In the algorithm, a multiobjective robust optimization model derived from the mean and variance of the objective function is constructed. Here, the objective function is defined based on the consideration of reducing the contact force generated at all clearance joints. Additionally, in order to balance computational accuracy and efficiency in the multiobjective robust optimization process, high-precision Kriging agent models are established. The optimum values of design variables are determined by combining Monte Carlo sampling and multiobjective particle swarm optimization method. By combining the Baumgarte approach with Lankarani–Nikravesh contact force model and Coulomb friction model, the dynamic equations of the planar multibody system with clearance joints are established. The uniform probability distribution is applied for characterizing random parameters. Simulation results show that the influence of design variable variations on the objective function changes in relation to the joint clearance size, but their relative influence degree on the objective function will not vary with the size of joint clearances. Moreover, the optimal solution selected on the Pareto front will affect the average levels and peak fluctuations of the dynamic responses in multibody systems.
1. Introduction
Due to manufacturing tolerance, material irregularity, and joint clearances, uncertainties always exist in any engineering system. Obviously, this will lead to robustness of optimization design in engineering design problems [1]. Moreover, the joint clearances continue to increase as wear progresses, which will seriously reduce the performance of the mechanical system and is also an important factor affecting the optimization results in multibody systems. Therefore, it is absolutely necessary to conduct robust optimization of multibody system considering the unceasing change of joint clearances and the randomness of design parameters.
Nowadays, more and more attention has been paid to the topics of modeling joints with clearance in multibody system dynamics. Based on the permanent contact condition, Earles and Wu [2] firstly introduced a model in which the clearance simulated a virtual link without mass connecting the bearing center to the journal center. Lankarani and Nikravesh [3] later established a contact force model based on the Hertz contact law. In this model, a hysteretic damping term was added to represent the energy consumption during impact. This model has been widely used by many researchers to this day. Flores and his research groups [4–6] used the contact force approach and a modified coulomb’s friction law to systematically study the dynamic behaviors of a rigid crank-slider mechanism with one clearance joint in plane and space. Based on the general trend of the Hertz contact model, Machado et al. [7] made a comprehensive and systematic comparison between several compliant contact force models, such as Lankarani–Nikravesh contact force model, Hunt–Crossley contact force model, Herbert–McWhannell contact force model, and provided information on the application scope and accuracy of these models in different contact scenarios. Along with the deepening of research, more and more scholars began to pay attention to the dynamic model with multiple clearance joints, to study the effect of the number of joint clearances on dynamic behavior of mechanical system. Flores and Lankarani [8, 9] investigated the dynamic behavior of a planar crank-slider mechanism with multiple joint clearances. Later, Bai and Sun [10] numerically investigated a planar four-bar mechanical system including multiple revolute clearance joints. The main results of their work were drawn that the interaction between multiple clearance joints is very strong, and the position and size of the joint clearances have a great influence on the dynamic responses of multibody mechanical systems. Not only that, Tan et al. [11] theoretically and experimentally investigated the dynamic behavior of a crank-slider mechanism with two revolute clearance joints. They found that the presence of the two clearance joints in the mechanism resulted in a mutual coupling region in the drawn contact figures. More recently, Wang et al. [12] synthesized to study the dynamic responses of a planer slider-crank mechanism with revolute clearance joint. They designed an experimental test rig to monitor the dynamic responses of the mechanism system by using a linear voltage differential transducer (LVDT) and an accelerometer. Li et al. [13] proposed an approach for modeling a clearance revolute joint with a constantly updating wear profile in a multibody system and studied the wear characteristics and dynamic response of a clearance revolute joint by numerical simulation and experimental tests.
Considering the inevitability of joint clearances, a growing number of researchers began to study the component flexibility, joint lubrication, and design optimization in multibody system from the perspective of reducing the influences of joint clearances on dynamic behavior. Taking the crank-slider mechanism as demonstrative application example, Schwab et al. [14] modeled the mechanical system in the case of the link flexibility and the lubrication in the clearance joint and analyzed the effects of the link flexibility and the joint lubrication on contact force peak values. Khemili and Romdhane [15] carried out an investigation on the dynamic responses of a slider-crank mechanism in the presence of flexible link. The response of the contact force generated at the clearance joint indicates that, due to the cushioning action of connecting rod flexibility, the response peak values were weakened. Erkaya and Uzmay [16] implemented a simulation of kinematics and dynamics for a planar four-bar mechanism considering link flexibility and joint clearance. They found that the response peaks in flexible mechanisms were smaller than those in rigid mechanisms. Furthermore, Sun et al. [17] introduced harmonic drive into the crank-slider mechanism with one clearance joint and compared the suspension effect of harmonic drive and connecting rod flexibility by the peak value of driving moment and contact force. In their study, the modeling of flexible components was based on the absolute nodal coordinate formulation (ANCF) that had also been extensively and thoroughly documented by Tian et al. [18–20]. In addition, Tian et al. [19] also developed an approach to establish the dynamic model of spatial multibody system with flexible components. In order to reduce friction, increase damping, and improve the stability of relative motion between journal and bearing, joint clearances are usually filled with some lubricant fluid. Tian et al. [20] established the dynamic model of the lubricated cylindrical joint of spatial multibody systems using hydrodynamic theory. The estimation of the maximum contact force indicates that, under the condition of lubrication, the peak value of contact force was greatly weakened. In recent years, there seems to be a growing interest in the problem of optimizing the multibody systems with joint clearances to reduce the effect of joint clearances on dynamic behavior. Erkaya and Uzmay [21–23] utilized a genetic algorithm to adjust the optimum values of optimized link in four-bar mechanisms for alleviating the unwanted vibration effect caused primarily by joint clearance. They also conducted kinematic optimization for the four-bar and slider-crank mechanisms with joint clearances, to decrease the effect of clearance on path error. Then, Erkaya [24] presented a novel modeling and optimization approach to decrease the undesired effects of joint clearances on a walking mechanism. Both Mukras [25] and Li et al. [26] optimized the wear performance of a planar crank-slider mechanism with two clearance joints. The difference between their studies rested with the calculation of the collision force at clearance joints, which is based on the contact force model and the elastic foundation model, respectively. Varedi et al. [27] carried out a dynamic optimization for a planar crank-slider mechanism with one clearance joint. Following, Varedi-Koulaei et al. [28] implemented simultaneously dynamic and kinematic optimization for a planar 3-RRR parallel manipulator with clearance joint. In their papers, an optimization algorithm based on particle swarm optimization (PSO) was proposed to track the optimum parameters of components in mechanisms.
It is important to note that the abovementioned researches on optimal design of multibody system is based on the assumption that design variations, such as structural and material parameters, physical dimensions of parts, and joint clearance sizes, are perfect; that is, deterministic optimization design does not consider the effect of these variations, and as a result, the design solution may be very sensitive to these variations. Based on reliability sensitivity analysis, Gao et al. [29] presented a general method for reliability optimization of a planar slider-crank mechanism with two clearance joints by considering the mass distribution of each component in the mechanism as random variables so as to reduce the effects of clearance variances on dynamic response. Huang and Zhang [30] applied an improved Taguchi method to investigate the robust tolerance design of a four-bar function generation mechanism with joint clearances. Through plotting sensitivity diagrams, they analyzed the influence of parameter changes on performance function of this mechanism. In order to improve analysis accuracy, Luo et al. [31] modeled a new robust design model for the four-bar function generating mechanism and proposed a new mechanism synthesis methodology in which clearance and dimension are handled as interval variables and truncated random variables, respectively. Batou et al. [32] investigated the robust dynamical design of a constrained multibody system, in which the single-objective robust optimization model was derived from the mean and variance of the random performance vector, and the Monte Carlo (MC) method was used as stochastic solver. Considering the sensitivity of design parameters’ uncertainties to deterministic optimization results, Bilel et al. [33] performed a multiobjective robust design optimization for a sewing machine system. In order to reduce the effects of joint clearance and parameters’ uncertainties on kinematic error and improve kinematic accuracy of mechanism system, Sun et al. [34] presented a novel multiobjective robust optimization methodology for a slider-crank mechanism with joint clearance. In this approach, the confidence region method (CRM) was intended for deriving the robust optimization model, and the robust design of the multibody system was performed using the Kriging surrogate model, MC simulation method, and multiobjective particle swarm optimization (MOPSO) method.
In the past studies of optimal design of multibody systems with clearance joints and random parameters, the sensitivity analysis, helping to evaluate trade-offs between available design variables, had often been ignored. However, it is essential for design optimization and optimal control. Currently, the direct differentiation and adjoint variable methods are two main techniques to perform sensitivity analysis of the multibody system [35]. Carrying out sensitivity analysis for such mechanical systems need to calculate derivative of its dynamic equation which is usually expressed as a highly nonlinear ordinary differential equations or differential-algebraic equations. Obviously, this requires the dynamic equation must be differentiable. In addition, in practical engineering, there may be situations, such as discontinuous load, complex derivative calculation, or lack of derivative conditions for the objective function. At this time, if the sensitivity calculation is carried out by the direct differential method or adjoint variable method, there will be great difficulties, even are not feasible.
The primary purpose of the present work is to present a general sensitivity analysis method and a multiobjective robust optimization methodology for sensitivity analysis and robust optimization of a constrained mechanical system with joint clearances and random parameters. The first-order sensitivity analysis based on the response surface proxy model is performed to find design parameters that have a significant impact on the objective function. Here, the mass, the mass center, and the moment of inertia of the moving components are taken as the design variables, while the journal center locus of all clearance joints is defined as the objective function. In order to verify the accuracy of sensitivity calculation, a method to determine the convergence interval of design variables is proposed, in which the sensitivity does not change even if the interval length of design variables is still shortened while keeping the midpoint of the interval constant. Based on sensitivity analysis, a double-loop multiobjective robust optimization model derived from the mean and variance of the objective function is constructed. In order to solve the multiobjective robust optimization problem, the MOPSO placed in the outer cycle is applied to minimize the optimization objectives. In the inner loop, the Kriging agent model and MC method are used to determine the optimization objectives for robust optimization design. Assuming existence of the ever-changing clearance size due to wear in the revolute joints, a planar crank-slider mechanism with two clearance joints is used as a typical example to verify the efficiency and optimality of the proposed approaches.
2. Dynamic Modeling of Multibody System with Revolute Joints in presence of Clearance
There are three modes of relative motion between the journal and the bearing, namely, contact mode, free-flight mode, and impact mode. Figure 1(a) shows the plane contact mode, in which a normal contact force and a tangential contact force between journal and bearing will be generated. Obviously, when no contact exists between them, the forces will vanish [28], i.e.,where is the relative penetration depth, as illustrated in Figure 1(b), and is given asin which is the clearance size; and are radii of the bearing and the journal, respectively, as depicted in Figure 1(a); and is the absolute value of eccentricity vector that connects and , as shown in Figure 1(b). It should be highlighted that is the unit vector that defines the normal direction of the plane of collision, and vector is obtained by rotating the vector in the counterclockwise direction by 90 degrees, as shown in Figures 1(a) and 1(b). Besides, the value of can be expressed as , where and represent the horizontal and vertical displacement of the journal center relative to the bearing center, respectively, which introduce two additional DOFs to the system.
(a)
(b)
Considering both the elastic and damping effects during the impact between journal and bearing, the contact force model proposed by Lankarani and Nikravesh is utilized in this work. The normal reaction force at the revolute joint in the presence of clearance is as follows [3, 4]:where denotes the elastic effect and denotes the damping effect during collision between clearance joints. is a constant whose value depends on material properties of the contact surfaces. In this article, the value of is 1.5 [29]. is the stiffness coefficient depending on the physical properties of the contact surfaces and their geometry and represents the damping coefficient. Their expressions are given as follows [3]:where and are the initial impact velocity and the restitution coefficient. Here, and are related to the elastic modulus and Poisson’s ratio of the bearing and journal, respectively, and are defined as follows [3]:
Substituting the values of and into equation (3) yields [3]and thus, the normal contact force vector is expressed as follows [4]:
It is noteworthy that the damping effect is related not only to the energy dissipated in the impact process but also to the friction dissipation effect on the contact surface. In this paper, the basic Coulomb friction model is adopted, which is expressed as follows [4]:where is the friction coefficient obtained from experiments. This paper takes .
The friction force vector is expressed as follows [4]:
The motion equation for a constrained mechanical system with joint clearance based on the Baumgarte approach is expressed as [36]where is the generalized mass matrix, is the generalized acceleration vector, is the Lagrange multipliers vector, and and are the constraint equation and constraint Jacobi matrix. . and are stabilization parameters. This paper takes [8]. represents the generalized force vector composed with elastic force and generalized external force. In this paper, the motion equation solved by using the fourth-order adaptive Runge–Kutta method.
3. Optimization Process
3.1. Optimization Problem Definition
Dynamic analysis indicates that the dynamic performance of mechanism system, in the presence of joint clearance, will change significantly if the mass of the links is redistributed [27]. Therefore, the mass, the mass center, and the moment of inertia of the crank and the connecting rod as well as the mass of the slider in the planar crank-slider mechanism displayed in Figure 2 are taken as all the design variables here, i.e.,where the subscripts 1, 2, and 3 of these design variables correspond to the component number, as shown in Figure 2.
Studies in literatures [8, 9, 25, 26] have shown that the collision strength and frequency between revolute joints of mechanisms with two clearance joints are much higher than those with only one clearance joint, which means that the dynamic performance of mechanisms is greatly affected by the clearance joints. In order to improve the dynamic performance of the mechanism, the direct strategy is to reduce or eliminate the impact forces generated at clearance joints. This can be guaranteed by keeping the journal and the bearing always in contact; i.e., the center locus of the journal relative to the bearing should be as close as to a circle. In the algorithm presented here, integral steps in a rotating period after stable operation are extracted. Therefore, the objective function in this paper is defined aswhere denotes the time interval of interest. and are the coordinates of the journal center with respect to the bearing center, if which and refer to the numerical integration step and the clearance joint, respectively. is the amplitude of the joint clearance, and is the number of clearance joints. is a vector containing design variables that have great influence on the objective function and .
In order to obtain the design variable vector and minimize the number of design variables as much as possible, sensitivity analysis is needed. Sensitivity, i.e., derivative information, reflects the influence of design variables or parameters’ changes on objective function or constraint function. Mathematically, it can be understood that if a function is differentiable, its first-order sensitivity can be expressed aswhere represents the design variable in . represents the sensitivity value of design variable at point in the design interval. The larger the absolute value of is, the more significant the influence of design variable on the objective function is. Moreover, the negative indicates the objective function decreases with the increase in design variable , while the positive indicates the objective function increases with increasing design variable .
In the actual mechanical system, design tolerance, fabrication, and assembly errors are unavoidable. Moreover, joint clearance and structural and material parameters are often imperfect, which seriously affect the dynamic behaviors of the mechanical system and always trigger uncertainties in the systems. In this paper, we consider that the information of uncertain variables is complete, and we use the probability method to describe uncertain variables. Under this circumstance, the objective function considering the randomness of parameters can be written aswhere is a vector containing some design variables. is a vector containing random parameters. .
If random parameters belong to probability space , respectively, then vector belongs to probability space , in which is a set of possible values of , is used to define the probability of elements in set , , and is a set of constraints contained in set . is defined as [34]
In particular, when variable is independent of each other,
More generally, the probability space expressed by probability density function iswhere .
3.2. Multiobjective Robust Optimization Definition
The concept of robust design was put forward by Dr. Taguchi in 1970s [37] and was introduced into product design. Later, Gabrel et al. [38] and Yao et al. [39] gave a comprehensive review of robust optimization theory. Robustness refers to the insensitivity of system performance to uncertainties, while robust optimization refers to seeking the optimal feasible solution of design variables in the case of unavoidable fluctuations and making the design objective function insensitive to the fluctuations of design variables, so as to minimize the adverse impact of uncertainty on system performance. The schematic diagram of objective robustness is shown in Figure 3.
As shown in Figure 3, assuming that point A is the optimal solution obtained by deterministic optimization and point B is the optimal solution obtained by robust optimization, point A obviously outperforms point B in terms of absolute value. However, when design variables change in the range of centering on point A and point B, respectively, the corresponding variations of function values are marked as and . From the graph, we can see , which indicates that, in the same range of design variables, point A is more sensitive to the change in design variables than point B. In this way, even if the objective function at point A is smaller than that at point B, it is not robust enough. When the variables change, the degree of deviation of the objective function from the optimization results may be very large, which may seriously cause the system to fail to meet working requirements.
The conventional optimization model considering an engineering design problem is described aswhere and are vectors of design variables and random parameters, respectively, as indicated in equation (15). is the constraint, and and are vectors of the upper and lower bounds of design variables, respectively.
In robust design, design variables and design parameters could all be the contributing sources of variations. Consequently, the expected system performance is a random function. At this time, both its mean value and variance are expected to be minimized at the same time. The general form of the objective function in equation (19) becomes
In this case, the multiobjective robust optimization model can be expressed as follows:
The main drawback of the optimization results obtained by equation (21) is that no compromise is recommended. In fact, the minimum of the first objective function is not necessary for the minimum of the second objective function . However, a Pareto front could draw all possible compromises between the two objectives. These compromises are some points where one of the two objectives is improved and the other will be deteriorated. Therefore, additional criteria need to be introduced to select the optimum on the Pareto front. In this paper, the solution with minimum optimization objective is defined as solution, while the solution with minimum optimization objective is defined as solution. The optimization results are expressed by the average values of and solution after four times of multiobjective robust optimization.
4. Multiobjective Robust Optimization Strategy Integrating Sensitivity Analysis, Kriging Surrogate Model, Response Surface Model (RSM), MC Method, and MOPSO
In order to solve the multiobjective robust optimization model in equation (21), a double-loop multiobjective optimization strategy is presented, which combines sensitivity analysis, Kriging surrogate model, RSM, MC method, and MOPSO. The outer loop is to search design variables to minimize the objective function, while the inner loop is to determine the objective function. The procedures of the multiobjective robust optimization strategy are shown in Figure 4.
4.1. Outer Loop Multiobjective Particle Swarm Optimization
In the multiobjective robust optimization strategy, MOPSO is applied in the outer loop to search design variables to minimize the optimization objectives. PSO is an evolutionary computation technology, invented by Kennedy and Eberhart [40]. It is similar to the genetic algorithm and is an iterative optimization tool. However, without crossover and mutation used by the genetic algorithm (GA), the PSO searches the optimal solution in the solution space according to the optimal particle. This algorithm is very suitable for solving nonconvex nonlinear optimization problems and is fast convergence and uses less computational time.
In PSO, each particle is a potential solution. The position and velocity of each particle in -dimensional space are, respectively, expressed as a vector. Suppose that the position and velocity of the particle in the -dimensional space are, respectively, and . In each iteration, the particle updates itself by tracking two optimal solutions. The first is the optimal solution found by the particle itself; i.e., the best local position , and the second is the optimal solution currently found by the whole population, i.e., the global best position . At this time, the particle updates its speed and new position through the following formula:where is the number of particles; is the dimension of each particle; represents the iteration; and and are learning factors that represent the ability to self-summarize and learn from excellent individuals in group, which are usually between 0 and 4. In this paper, they are taken as 1.8. and are independent random numbers between 0 and 1. denotes the inertia weight factor. In view of the premature performance and the oscillation near the global optimal solution in the later stage of the PSO algorithm, this paper adopts the linear decreasing method to determine . Let the inertia weight decrease linearly from its maximum to minimum . The formula of change with iteration of the algorithm is as follows:where is the maximum iteration step. This paper takes and as 1.2 and 0.1, respectively.
The steps of the outer loop multiobjective particle swarm optimization are as follows [34]:(1)Choose the number of particles, , and the maximum iteration step .(2)Initialize particle position and velocity.(3)Through inner loop, evaluate the mean and variance of objective function and then get the fitness of particles.(4)Estimate whether particles have been recycled. If no, go back to the previous step. If yes, do the next calculation.(5)Compare and substitute with the best individual pareto optimal solutions.(6)Compare and substitute with the global pareto optimal solutions.(7)Judge if the maximum iteration has been reached. If yes, stop and output global pareto optimal solutions. If no, do the next calculation.(8)Go to the next iteration, .(9)Determine individual and global best position.(10)Update particle position and velocity by equation (22), then return to Step (3).
4.2. Inner Loop Sensitivity Calculation and Random Analysis
In the inner loop, in order to balance computational accuracy and efficiency of random analysis, Kriging-based Monte Carlo method is applied for system performance random analysis to obtain the optimization objectives. In addition, in order to improve the efficiency of optimization, RSM-based sensitivity analysis of the multibody system is carried out at first. Equation (14) has clearly shown that the premise of sensitivity calculation is to get the explicit expression of objective function, then get derivative with respect to each design variable, and finally obtain the sensitivity value of each design variable at a certain point.
The steps of sensitivity calculation and random analysis are as follows:(1)Determine the convergence interval of each design variable for sensitivity calculation by continuously shortening the interval length of design variables while keeping the midpoint of the interval constant. The first step is to establish a high-precision polynomial function based on the response surface model between system performance and design variables within the global design interval of design variables, then to take each design variable derivatives of the polynomial function, and finally to calculate the sensitivity value at midpoint of the interval for every design variable. In the second step, the length of the design interval is reduced by increasing the lower limit of the design variable interval obtained in the first step by , decreasing the upper limit by , and keeping the midpoint of the interval unchanged, to obtain a new design interval. Here, is the variation of the upper and lower limits of the interval. Then, according to the sensitivity calculation method in the first step, the sensitivity is calculated again based on the changed design variable interval. In the third step, the sensitivity values calculated in the first two steps are compared. If the two sensitivity values are the same, the design variable interval obtained in the second step is the final determined convergence interval of design variables. Otherwise, according to the method of obtaining the design variable interval in the second step, a new design variable interval is obtained to calculate the sensitivity again and then the sensitivity calculated in the two times is compared. Repeat the above steps until the sensitivity values calculated in the previous two steps are the same, so as to determine the convergence interval of design variables. The overview of the response surface methodology for sensitivity analysis is as follows. RSM originated from fitting the input-output relationship of experimental data in physical experiments [41]. Because of random experimental error , the response value obtained by experiment is different from the expected value . Assume that is a vector consisting of design variables. The relationship between them is as follows: The approximating function is typically a low-order polynomial and can be decomposed as follows: where is the coefficient vector solved by the least square method and . If is a -order polynomial, the dimension of coefficient vector is . In this paper, in order to build the high-precision RSM, a group of initial sample points for design variable vectors was extracted by the Latin hypercube sampling (LHS) method [42], the number of which was 45. Then, they were brought into the dynamic equation of multibody system equation (11) to evaluate the performance values of impact force generated at all clearance joints, as expressed by equation (13), and 45 performance values were eventually obtained. With the initial set of sample points, the RSM which is initially assumed to be 3-order polynomial is constructed. At this time, the relative maximum absolute error and the complex correlativity coefficient , as follows, representing the local and the global fitting accuracy, respectively, are used to test the fitting accuracy of the response surface agent model: where is the real response value estimated by the dynamics equation, is the approximate response value evaluated by RSM, and are the mean and standard deviation of the test samples, is the number of the test samples, and this paper takes . If is close to 1 and is close to 0, it indicates that the RSM built has a high fitting accuracy. Otherwise, it represents the overall fitting accuracy of the surrogate model is poor, and it is necessary to resample or change the order and item of polynomial until meeting requirements.(2)Construct the response surface model between system performance and design variables in the convergence interval.(3)Obtain the derivative of every design variable at the midpoint of the convergence interval.(4)Compare the magnitude of sensitivity value and identify design variables that have greater impact on system performance.(5)Construct high-precision Kriging agent model between system performance and design variables in global interval. Here, the Kriging agent model is based on the DACE toolbox in MATLAB established by Lophaven [43]. In order to construct a high-precision Kriging agent model in global interval of design variables for optimum design of multibody system in this paper, first and foremost, it is necessary to determine the subintervals with high fitting accuracy of each design variable. Then, the Kriging agent models in the subinterval of each design variable are established. Assuming that the global interval lengths of each design variable are and the lengths of the determined subintervals are , then the Kriging agent model in the global interval of design variables will be composed of Kriging agent model built in the subintervals. In the process of particle swarm optimization, when a particle in the global design interval is given, the corresponding elements in the particle representing each design variable judge their own subinterval belonging at this time, and then the response value based on the Kriging proxy model in this subinterval is calculated. Here, it should be noted that equations (26) and (27) are also used to test the accuracy of the Kriging proxy model to determine whether a subinterval has a high fitting accuracy. If the fitting accuracy in the subinterval is not high, the original interval will continue to be reduced or resampled until a convergent subinterval with high fitting accuracy is obtained.(6)Choose the number of sample points of random parameters, .(7)For random parameters in each particle from the outer loop, perform MC sampling according to its distribution. For the nonrandom parameters in each particle, spread its sample points into identical values.(8)Introduce the samples of design parameters into the constructed high-precision Kriging model in turn to compute system performance .(9)Estimate whether sample points have been recycled. If no, go back to step (8). If yes, do the next calculation.(10)Calculate the fitness function equation (20).
5. Case Study
In this section, the crank-slider mechanism is used to demonstrate the efficiency of the proposed method in this paper [26]. The kinematic model of the mechanism with two clearance joints marked as O12 and O23 is depicted in Figure 2. Moreover, the component 1 (crank) and the component 2 (connecting rod) are assumed rigid, and the joint between component 1 and component 4 (frame) is assumed rigid and without any clearance. Initially, the crank and the connecting rod are collinear, and thus, the journal and bearing centers are coincident. Table 1 [26] gives the geometric and inertia properties of the original mechanism in which the mass center of the crank and the connecting rod are in the middle of the links. The material properties associated with journals and bearings and the radii of the journals and the bearings are listed in Table 2 [25]. From Table 2, it can be inferred that, initially, the size of joint clearances is 0.05 mm. Table 3 [27] gives the constraints of design variables and the changing range of joint clearances. Moreover, the problem is solved using the MOPSO algorithm, which leads to the optimal solution in Table 4. The crank is driven with a constant velocity equal to 5000 rpm.
6. Results and Discussion
In order to clarify the sensitivity of design variables shown in Table 3 with the objective function , it is necessary to establish the functional relationship between the objective function and design variables. In the process of multiobjective robust optimization, determining the minimum value of and for any set of sample point (, , , , , , and ) requires a complete optimizing process, which could involve hundreds or even thousands of cycles, and each cycle also contains several thousand numerical integration steps for solving the dynamic equation of multibody system. But beyond that, if uncertainties of design variables and the changing joint clearance size due to some errors which may appear during fabrication and assembly of the mechanical system components are taken into account, the number of sample points will increase sharply; then, the optimization process is quite computationally intensive and time-consuming during constantly evaluating the objective function. Therefore, the surrogate-based optimization techniques, Kriging agent model and response surface model with high fitting accuracy, are developed and utilized in this research to perform the optimization processing.
In this illustrative example, in order to describe the influence of design variables on the objective function with the change in joint clearances more intuitively, four high-accuracy Kriging proxy models are constructed. The surface plots between objective function and design variables as well as clearance magnitudes are shown in Figures 5(a)–5(d). The fitting accuracy of these four Kriging proxy models is measured by comparing with the dynamic equation and the surrogate model to solve the objective function value , as described in Figures 6(a)–6(d), respectively. Furthermore, it is shown that, in Figure 6, the global fitting accuracy tends to be 1 and the local fitting accuracy tends to be 0, which indicates that the four Kriging agent models have high fitting accuracy and can accurately characterize the impact of design variables on the objective function.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
From these surface plots in Figures 5(a)–5(d), it can be seen clearly that the larger the clearance size , the more significant the objective function value representing the collision strength between the journal and the bearing at the clearance joints. The reason for this phenomenon is that increasing clearance size will increase the nonperiodicity of the mechanism dynamic behavior, resulting in more severe collision between the journal and the bearing. From Figures 5(a) and 5(b), it is clear that as the mass of the connecting rod and the slider increases, the objective function value grows, which can also be concluded from the sensitivity curve shown in Figures 5(a) and 5(b). However, it is also clear that the steady changing trend of the objective function value is obtained with an increase in the mass of the crank , even if the size of joint clearance is increasing. From Figures 5(c) and 5(d), it can be concluded that the effects of the mass center and the moment of inertia of the connecting rod and on the objective function are more obvious than that of the crank and . Moreover, the mass center of the connecting rod is negatively correlated with the objective function , and the moment of inertia of the connecting rod is positively correlated with the objective function .
Figures 5(a)–5(d) reveal the sensitivity changing at midpoint of global interval of the design variables with the joint clearance size. They are investigated via four cubic polynomial response surface model built. Figures 7(a)–7(d), corresponding to Figures 5(a)–5(d), respectively, show the variation of fitting accuracy of the four response surface models with joint clearance magnitude within different variation intervals of design variables. From Figures 7(a) and 7(b), it is clear that high fitting accuracy can be obtained in the global design interval of the design variables , , and , while Figures 7(c) and 7(d) indicate that only when the convergence intervals of design intervals , , , and are determined by reducing their global interval length, the higher fitting accuracy can be obtained, and the local fitting accuracy is still not very high. The main reason can be explained as the influence of these design variables , , , and on the objective function varying considerably with the clearance joint size. Moreover, the lower local fitting accuracy is also reflected in Figure 5(d) whose surfaces are not smooth enough.
(a)
(b)
(c)
(d)
Based on the above sensitivity analysis, one can infer that the design variables having great influence on the objective function are the mass of the connecting rod , the mass of the slider , the mass center, and the moment of inertia of the connecting rod and . Therefore, this paper ultimately determines to use these four parameters as design variables to carry out the following optimization design, and the change in joint clearance is taken into account in the whole optimization process.
With the Kriging surrogate model, the slider-crank mechanism with two clearance joints is optimally designed under two cases: (1) uncertainty of the four design variables , , , and are ignored, while clearance size follows uniform distribution U(0.3, 1) mm. (2) The clearance size and the four design variables , , , and are considered as uniform random variables, which are given by U(0.3, 1) mm, U(, ) kg, U(, ) kg, U(, ) m, and U(, ) . , , , and are design values of the connecting rod mass, the slider mass, the mass center, and the moment of inertia of the connecting rod, respectively. The fitting accuracy of the Kriging proxy model which describes the relationship between objective function and variables , , , , and is shown in Figure 8. The Pareto fronts from a certain optimization of the multiobjective function for the two cases are shown in Figure 9. The final optimization results are given in Table 4. In addition, it must be pointed out that the multiobjective robust optimization process of the multibody system in presence of joint clearances with changing size only takes about two seconds.
The optimized parameters in Table 4 are introduced into the dynamic equation of the observed mechanism system with two clearance joints. It is assumed that the size of the two clearance joints has changed to 0.3 mm as wear progresses, and thereafter, the aperiodicity of dynamic responses becomes more and more complex. Figure 10 illustrates the comparison of the dynamic response results for and solutions, which are quantified by plotting the contact force 1 produced on O12, contact force 2 existing in O23, crank reaction moment acting on the crank, and slider acceleration. These dynamic response curves are plotted against those from the initial design, i.e., before optimization, and reported for two full rotation cycle of the crank after steady state.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Figures 10(a)–10(d) show the comparison of dynamic behavior curves from solution of deterministic optimization and uncertain optimization, which clearly reveals that, compared with before optimization, the fluctuation peaks of contact force 1, contact force 2, and the crank reaction moment are reduced sharply for these two optimization results. Moreover, except for the acceleration response of the slider, the obtained uncertain optimization results can reduce the peak fluctuation of the dynamic response more effectively.
Figures 10(e)–10(h) illustrate the comparison of dynamic behavior curves from deterministic optimization and uncertain optimization for solution. It is easy to see that the peak fluctuation of dynamic responses, including the acceleration response of the slider, is significantly reduced under the condition of solution obtained by uncertain optimization.
Comparing Figures 10(a)–10(d) with Figures 10(e)–10(h), it can be found that the dynamic behavior curves for solution obtained by uncertain optimization have particularly much more peak values than that for the solution obtained by uncertain optimization. Compared with the peak value of contact force 2 curve shown in Figure 10(f), the peak value of the contact force 2 curve shown in Figure 10(b) is especially prominent. Moreover, the same phenomenon can be found by comparing Figure 10(a) with Figure 10(e); that is, the curve of contact force 1 obtained under uncertainty optimization solution has more pronounced peaks than that obtained under uncertainty optimization solution . In addition, these two impulse contact forces would propagate through the rigid component 2 and component 1 of the mechanism and, ultimately, have a direct impact on the acceleration response of the slider and the crank reaction moment, respectively. Therefore, the curves of the slider acceleration and the crank reaction moment change more smoothly under uncertain optimization solution , as shown in Figures 10(g) and 10(h), compared with that under uncertain optimization solution , as shown in Figures 10(d) and 10(c), respectively.
In addition, the mean and variance of the objective function and are also reflected in the dynamic response curve of the system when uncertain optimization is carried out. From Figures 10(a) and 10(e), it can be observed that the mean value of contact force 1 obtained by uncertainty optimization solution is larger than that obtained by solution of uncertainty optimization. However, the peak fluctuation of contact force 1 obtained by uncertainty optimization solution is smaller than that obtained by uncertainty optimization solution. Therefore, the uncertainties of design variables can affect optimal results of the mechanism, and the optimization results representing the mean and variance of the optimization objective will have an impact on the average level and peak fluctuation of dynamic response, respectively. The larger the variance of the optimization objective, the more obvious the peak value of the response curve produced by the corresponding optimization result. However, the larger the average value of the optimization objective, the higher the average level of the response curve generated by the corresponding optimization result.
7. Conclusion
In this paper, a general approach of sensitivity analysis and multiobjective robust optimization design for the multibody system in the presence of varying clearance size of revolute joints is presented to investigate the influence of parameter uncertainties on optimization results when clearance size of revolute joints is constantly changing. Here, the change in joint clearance size considered in the robust optimization design is to take into account the uncertainty of wear between journal and bearing. The first-order sensitivity based on the response surface proxy model is used for sensitivity analysis to minimize the number of design variables as much as possible and thus to reduce uncertainties. In order to verify the accuracy of the calculated sensitivity, it relies on the fact that when the interval length of the design variable is shortened and the midpoint of the interval is constant, the sensitivity value in the convergence interval is unchanged. The robust optimization design is considered as the double-loop process. In the inner loop, high precision Kriging surrogate models are constructed, and the Monte Carlo method is applied for random analysis of system performance to obtain the objective function. In the outer loop, the MOPSO is proposed to search design variables to minimize the objective function.
The numerical simulation results of sensitivity analysis show that when considering constant change of the joint clearance sizes, design variables, including the mass of the connecting rod , the mass of the slider , the mass center, and the moment of inertia of the connecting rod and , have great influence on the objective function . Moreover, the sensitivity of each design variable varies evidently in relation to the joint clearance size, but their relative influence on the objective function will not change with the size of joint clearance. The results of multiobjective robust optimization show that joint clearances and parameter uncertainty obviously affect optimal results of the mechanism. The robust optimal design presented here offers better dynamical performances over the deterministic optimal design and original design. In addition, the optimal solutions and , which represent the minimum optimization objective and , respectively, selected at the Pareto front also have a great influence on the dynamic behavior of multibody system. These two solutions and affect the average levels and peak fluctuation degree of the dynamic response, respectively. The average level of dynamic response from solution is lower than that from solution . However, the fluctuation degree of dynamic response from solution is higher than that from solution .
Nomenclature
| : | Vector containing all design variables |
| : | Vector containing some design variables |
| : | Vector containing random parameters |
| B: | Possible values set |
| : | Constraint set |
| : | Probability density function |
| : | Objective function |
| : | First-order sensitivity |
| : | Mean and variance of objective function |
| : | Constraint function |
| : | Solution with minimum |
| : | Solution with minimum |
| : | Position of the particle |
| : | Velocity of the particle |
| : | Best local position |
| : | Global best position |
| : | The number of particles |
| : | Dimension of each particle |
| : | The iteration |
| : | Learning factors |
| : | Independent random numbers |
| : | Inertia weight factor |
| : | Response value obtained by experiment |
| : | Random experimental error |
| : | Coefficient vector of approximating function |
| : | Unit vectors defining normal collision direction |
| : | Unit vectors defining tangential collision direction |
| : | Radius |
| : | Penetration depth |
| : | Restitution coefficient |
| : | Stabilization parameters |
| : | Friction coefficient |
| : | Young’s modulus |
| : | Poisson’s ratio |
| : | Normal contact force and its vector |
| : | Friction force and its vector |
| : | Generalized contact force vector |
| : | Restitution coefficient |
| : | Penetration velocity |
| : | Initial penetration velocity |
| : | Determination coefficient |
| : | Root mean square error |
| : | Global mass matrix |
| : | Constraint vector |
| : | Constraint Jacobian matrix |
| : | Vector of generalized coordinates |
| : | Vector of applied loads |
| : | Vector of Lagrange multipliers. |
Data Availability
The Matlab data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by a project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions and the Initial Scientific Research Fund by University of Jinan (Grant no. 160100323).