MINIMIZING CONTACT STRESSES IN AN ELASTIC RING BY RESPONSE SURFACE OPTIMIZATION

(1)

MINIMIZING CONTACT STRESSES IN AN

ELASTIC RING BY RESPONSE SURFACE

OPTIMIZATION

Asim Rashid

THESIS WORK 2010

(2)

MINIMIZING CONTACT STRESSES IN AN

ELASTIC RING BY RESPONSE SURFACE

OPTIMIZATION

Asim Rashid

This thesis work is performed at Jönköping Institute of Technology within the subject area Product Development and Materials Engineering. The work is a part of the master’s degree.

The author is responsible for the given opinions, conclusions and results. Supervisor: Niclas Strömberg

Credit points: 30 ECTS credits Date:

(3)

Abstract

In this thesis an optimization routine is implemented to minimize the contact pressure. The routine is based on the successive response surface methodology and is applied to the problem of minimizing contact pressure between a snap ring and its housing. Design of Experiment techniques are used to decide the combination of design variables. A Python script is developed that performs the computer simulations in Abaqus. Maximum contact pressure information generated by Abaqus is imported in Matlab. Combination of design variables and results from Abaqus are used to build second order response surfaces by using least squares method. Then quadratic programming is used to determine the optimum of the response surface. This process is repeated according to the rules of successive response surface methodology until convergence is achieved. Two different geometries, polynomial and super elliptical curves are considered as the prospective candidates for improved design of snap ring. The results show that a polynomial curve is more suitable than a super elliptical curve. A polynomial curve gives a lower value of maximum contact pressure and better distribution over the whole length as compared to the super elliptical design. The optimization has been compared with a commercial software and shows promising results.

Key Words

Successive Response Surface methodology, SRSM, Contact Pressure, Contact Stress, Snap Ring, Optimization

(4)

1 Introduction

Snap rings are widely used in mechanical assemblies, both internal and external. When an internal snap ring is placed in position in the housing it exerts the pressure on the walls of housing. During service they suffer rapid wear at the open ends, due to higher contact pressure at these locations. Conventionally a circular arc with constant radius is used as base curve for snap ring geometry. In this thesis, super elliptical and polynomial curves are investigated as alternative to circular arc. Contact pressure calculations are combined with successive response surface methodology (SRSM) to determine the parameters of base curves for minimizing the contact pressure and hence wear.

1.1 Background

Response surface methodology (RSM) was first introduced in Box and Wilson (1951). Since then it has been widely used for optimization problems. Successive Response surface methodology (SRSM) is the extension of RSM, Where the design space is gradually reduced. This gradual reduction in design space is presented in Stander and Craig (2002). This technique has been implemented and improved by Gustafsson and Strömberg (2008).

The SRSM has then been used, for example, for shape optimization of castings (Gustafsson & Strömberg, 2008) and optimization of sheet metal forming processes (Jansson, 2005)

Previous work to determine the optimal shape of a snap ring was done by Strömberg (2007). In that work a polynomial curve was considered and an in-house contact software was used. Neural networks with SRSM were used to find the optimal geometry.

1.2 Purpose and aims

The aim is to determine the equations for super elliptical and polynomial curves that minimize the maximum contact pressure and compare the results for both geometries. The total length of all new proposed designs of snap ring must be equal to the length of original design with circular arc. Mathematically it can be expressed as

minmax

s.t. (1)

Here represents the collection of all design variables and represents the length of the original snap ring. Figure 1 shows geometry of the original snap ring in free state.

(6)

Figure 1: Original Geometry of Snap Ring

Least squares method will be used to fit the second order response surface and optimization problem is solved by using Quadratic Programming. The solution obtained by Quadratic Programming is also compared with optimization performed in commercial software, HyperStudy, which is a part of the software package HyperWorks (HyperWorks, 2009). In HyperStudy, method of moving least squares will be chosen to fit the response surface.

1.3 Delimits

Only second order models will be used to build a response surface. Higher order models or the techniques like neural networks or Kriging model will not be considered for this work.

1.4 Outline

The thesis is organized in several parts. In the current introductory section, the project background, aims and scope have been described.

In section 2, a more thorough theoretical discussion on the main ideas and concepts behind the Response Surface Methodology, Successive Response Surface Methodology, Contact Mechanics and different types of base curves considered here, are discussed.

Section 3 gives the details how the project was implemented. Which software were used and how the optimization routine works.

In Section 4, the results achieved in this work are presented. In section 5, the conclusion of the performed work is summarized. In appendix, dimensions of the original snap ring and housing are given.

(7)

2 Theoretical background

This section will provide information about the theoretical knowledge on which the implementation of the project will be based.

2.1 Response Surface Methodology

“Response Surface Methodology (RSM) is a collection of statistical and mathematical techniques useful for developing, improving and optimizing processes.” (Myers & Montgomery, 2002, p. 1)

Different techniques are used for optimization of a function. But in most practical problems the relation between response and design variables is not known. In such situations, several experiments or computer simulations are performed and relation (function) between independent variables and response is determined. Once the relation has been determined then the optimum of the original problem can be estimated by determining the optimum of this function.

Basic Concepts

Before going into the details of how response surface methodology is implemented, first a few basic concepts will be explained.

Response

There is always some criterion or measure to judge the performance of a product or system. This criterion is called response.

Design Variable

The variables that affect the response are called design variables, independent variables or factors.

Response Surface

The algebraic expression that describes the response as a function of the independent variables is known as response surface. In other words, response surface describes the relationship between the response and the independent variables. Figure 2, shows graphically the relationship between response and independent variables x ( ₁,x₂) for the following function.

(8)

Figure 2: Response Surface -6 -4 -2 0 2 4 6 -5 0 5 -8000 -6000 -4000 -2000 0 2000 x₁ x₂

The graphical display is a very useful tool that elaborates the relationship between response and variables in a very convenient way. This graphical display of the response as a surface has led to the terminology of response surface.

If a dependent variable y is a function of k independent variables then mathematically it is expressed as

y f x , x , x … … … . . x (2) Such a relation is called the functional relation. For given values of independent variables, the function f indicates the corresponding value of y. Figure 3 shows functional relation graphically.

Independent Variable Dependent Variable

Figure 3: Example of Functional Relation

But the function is not always known for most practical problems. In such situations, several experiments or computer simulations are performed and relation between independent variables and response is determined. In this situation it is called statistical relation and expressed as.

(9)

y f x , x , x … … … . . x ε (3) Here ε is the random error term. Figure 4 shows statistical relation along with responses graphically.

Independent Variable Response

Figure 4: Example of Statistical Relation

A statistical relation, unlike a functional relation, is not perfect. This means that the curve (in case of single independent variable) does not, in general, passes through the data points. Instead it is a close fit to the data points.

Implementation of RSM

Implementation of RSM is a three step process. These steps are following 1) Performing Experiments

2) Building Response Surface

3) Optimization of Response Surface

2.1.1 Performing Experiments

To determine the statistical relation between independent variables and response, it is always required to perform sufficient number of experiments or computer simulations. Computer simulations have high computational cost for complex problems. Therefore it is desired to know about a problem much, with only a few simulations. DoE (Design of Experiments) is the method to choose the combinations of design variables to perform experiments or simulations. This is also called designed experiments or experimental design. DoE give maximum information using minimum runs and hence minimum cost.

(10)

Basic Terminology

Before going into the details, first basic terminology will be explained.

Variable Domain

The set containing all the values between lower and upper limit of a variable is called variable domain.

Level

Level is defined as the value of a design variable considered for an experimental run. For two level designs, the levels are taken as the lower and upper limit values of a design variable. For three level designs, middle value of the variable domain, in addition to the lower and upper limit values are considered as the levels of variable. This is elaborated in Figure 5.

Two Levels

Three Levels

Figure 5: Levels of a Variable

Design Space

Design space is composed of all possible combinations of values from the domains of variables involved in an experiment. Figure 6 elaborates the design space graphically.

Max

Coded Variable

It is customary that the design variables are scaled such that the lower and upper limit of each variable is given by -1 and 1respectively. Let ξ represent a variable in natural units and x represent the variable in coded units, then a variable can be transformed by the following formula

Variable 1

Min Max

Min

Design Space Variable 2

(11)

/2

/2 (4)

Design Point

One run of an experiment or computer simulation can be represented by a single point in Design space. It is required to specify the levels of all variables to define a design point. Figure 7 elaborates the design point graphically

Design Point Min Max Min Max Variable 1 Variable 2 Design Space

Figure 7: Design Point for two variables Different Experimental Designs

Different experimental designs are used depending upon the area of application. Only few will be presented here.

2.1.1.1 Factorial Experiments

An experimental design is called factorial design if all the design variables involved are considered only for a specified number of levels. It contains all possible level combinations of design variables. Such a design consists of Lk design points, here L denotes number of levels and k denotes the number of variables. The disadvantage of factorial experiments is that the number of experiments increases rapidly with the increase in number of design variables. Figure 8 shows a 2k

and 3k

designs.

2k Factorial Design

A design with lower and upper limit combinations of all design variables is called a full factorial design in two levels.

3k Factorial Design

A design with middle, lower and upper limit combinations of all design variables is called a full factorial design in three levels.

(12)

X2

X1 X2

X1

Figure 8: 2k_{and 3}K_{Factorial Designs}

2.1.1.2 Central Composite Design

This is obtained by adding one center point and a set of star points along the axis to 2k

design. The star points are located at a distance from the origin. If | | 1, then the design is called circumscribed and if | | 1 , then design is called inscribed as elaborated in Figure 9.

-1.5 -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 1.5 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Circumscribed Inscribed Figure 9 : Central Composite Designs

2.1.1.3 Box-Behnken design

First all possible combination of variables pairs are identified. Then variables in these pairs are combined according to 22 factorial design and remaining variables are set to zero in each design point. Then a center point is added. More than one center points can be added if required. Figure 10 shows a Box-Behnken design for three variables.

(13)

1 2 3 1 1 0 1 1 0 1 1 0 1 1 0 1 0 1 1 0 1 1 0 1 1 0 1 0 1 1 0 1 1 0 1 1 0 1 1 0 0 0

Figure 10: Box-Behnken Design for three variables

Choosing an Experimental Design

Experimental design chosen for a given problem should minimize the number of runs without losing quality. Following are the general guidelines for choosing experimental design.

1) 2k

factorial design is appropriate for a first order response surface.

2) Central composite and Box-Behnken designs are appropriate for a second order response surface.

2.1.2 Building Response Surface

After performing the desired number of experiments or computer simulations, next step is to determine the relation that describes the response as a function of the independent variables. The relation is an approximate model and statistical techniques are used to build it. This involves choosing a polynomial of suitable order. A first order model, for two variabl is es

y β β x β x ε (5)

In this equation β’s are unknown constants and need to be determined. β’s are called regression coefficient . A se o d or er mo el, for s c n d d two a v riables is

y β β x β x β x β x β x x ε (6)

In generalized form first and second order polynomials can be written as

y β β x K ε y β β x K β x K β x x K ε (7)

Here K is the number of variables.

-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1

(14)

Selection of a suitable order of polynomial depends upon the complexity of the problem. When there is significant curvature in the responses, then first order model is inadequate. In such case second order model is reasonable choice. Figure 11 elaborates when higher order polynomials should be chosen.

Third Order Second Order

First Order

Figure 11: Curvature and Order Of Polynomial

The second order model is usually preferred as it is very flexible and can model quite complex processes. It is generally a good estimation of the true response surface.

Estimation of the Regression Coefficients

To estimate the values of regression coefficients in the polynomials, it is required to perform a minimum number of computer simulations. This number of computer simulations depends upon the number of variables and order of polynomial. Table 1 lists the minimum number of simulations required (Strömberg, 2009).

Order of model Minimum Number of Simu lations (Nmin)

First K 1

Second K 1 K 2

2

Table 1: Minimum Number of DoE

Following three methods are used to estimate the regression coefficients. 1. Least Squares

2. Weighted Least Squares 3. Moving Least Squares

(15)

2.1.2.1 Le tas Squares

Suppose that computer simulations are performed to compute the responses for at least n N experimental points. Let x denote the ith

level of variable x and the responses be designated as y , y , y … … . y . Then the first order model in terms of the indivi ud al re ponses can be described as s

y β β x β x β x ε , i=1,2,3…….n (8) This equation can be expressed in terms of all responses in matrix form as

(9) Where , 1 ₁₁ ₁₂ ₁ 1 21 22 2 1 ₁ ₂ , β β1 β , ε ε₂ ε

The method of least squares estimates the β’s so that the sum of the squares of the errors, ε, are minimized. Graphically it can be explained as minimizing the sum of the squares of the vertical deviations of responses from estimated model as shown in Figure 12.

Estimated Model Response Value

Figure 12: Deviation of Data Error in the matrix form can be expressed as

(10)

Then sum of the squares of the errors can be expressed as

L ε (11)

(16)

According to method of least squares, the estimator of β is the b that minimizes the sum of squares of errors. his im lies that T p

∂L ∂ (12) b b₁ b Where

The estimator b can be determined from the following equation

(13)

So the fitted model is

Ŷ Xb (14)

2.1.2.2 Weighted Least Squares

Residual can be defined as the difference between the Response value yi (obtained by computer simulation) and the corresponding fitted value ŷi. Mathematically it

can be expressed as

ŷ (15)

Magnitude of residual is represented by the vertical distance of the response value y_i from the corresponding point on the fitted model as shown in Figure 13.

ŷi yi Ŷ Xb Response Independent Variable Figure 13: Illustration of Residual

Residuals are used to evaluate whether a fitted model is appropriate for the response data at hand. Figure 14 shows scatter plot of responses and fitted model.

(17)

It can be seen that the residuals increase with the increase in independent variable. Higher the value of the independent variable, the more spread out the residuals are.

Ŷ X

Figure 14: Scatter Plot and Fitted Model

Responses with small spread are more reliable than those with large spreads. Weighted Least Squares (WLS) technique is used to take this relative reliability of responses into consideration. Responses with small spread are given more importance or weightage to determine the estimators of regression coefficients. In Weighted Least Squares technique, each response is given a weight and this weight is incorporated into the criterion for determining the estimators of regression coefficients.

L w ε (16)

If error variance (σ ) is known, then weight can be expressed as (Kutner, Nachtsheim & Neter, 2004)

w 1 ₍₁₇₎

The normal equation can be exp ssre ed as follows.

(18)

And the estimator b can be determined from fo low l ing equation

(19)

W is the matrix containing the weights for all responses.

0 0

(18)

2.1.2.3 Moving Least Squares

Moving least squares (MLS) method is a generalization of weighted least squares technique. Once a response surface has been built, response values can be computed for a new set of variables without running a computer simulation of original model. Let be a point where fitted model is evaluated and a DoE sampling point. Then the distance eb tween these two points can be written as

(20)

In MLS the weight is the function o tf his δ

w (21)

The main difference between WLS and MLS is that the weights associated with individual experimental points remain constant in WLS but they are a function of δ in MLS. The weight associated with a particular experimental point , decays as the point moves away

Examples of some commonly used f un ict ons are (Strömberg, 2009)

ω ,

ω δ exp δ

2α 0

(22)

Because the weights are a function of , the estimators of regression coefficients are also dependent on . So

(23)

2.1.3 Optimization of Response Surface

Optimization can be defined as the process of achieving optimal level or state of something. Optimal state can be either the maximum or minimum of a function. Once a response surface has been built, response values can be computed for a new set of variables without running a computer simulation of original model. Optimum value for the process can also be determined by determining the optimum of response surface.

If a second order model is used for a response surface, then we can easily compute the Stationary Point of a Response surface by differentiating the fitted model and equating it to zero. But additional information is needed to decide the nature of the stationary point i.e. if it is a point of maximum response, minimum response or saddle point. In case of saddle point, response is neither maximum nor minimum. Figure 15, Figure 16, and Figure 17 show these cases for two variables.

(19)

Figure 15: Stationary Point is Maximum Response -5 0 5 -5 0 5 -800 -600 -400 -200 0 200 x₁ x₂ Figure 16: Stationary Point is Minimum Response -5 0 5 -5 0 5 -200 0 200 400 600 800 x 1 x₂ Figure 17: Stationary Point is Saddle Point -6 -4 -2 0 2 4 6 -5 0 5 -400 -200 0 200 400 600 x₁ x₂

Second order polynomial fo k variables can be expressed as r

y β β x β x β x β x β x

(20)

Then the fitted model in matrix notation can be written as (Myers & Montgomery, 2002). ŷ (25) Where x x2 x , b b₂ b , /2 /2 /2 .

Stationary Point can be determined by differentiating the fitted model and equating it to zero.

ŷ 0 So

By denoting the stationary point as s

2 0 (26)

x

(27)

The predicted response at the s att ionary point is

ŷ

(28)

Let λ1, λ1,…….. λk be the eigenvalues of B. The signs of λ’s determine the nature

f stationary point. (Myers & Montgomery, 2002) o

1. If λ1, λ1,…….. λk are all negative, the stationary point is a pint of maximum response.

2. If λ1, λ1,…….. λk are all positive, the stationary point is a pint of minimum response.

(21)

2.2 Successive Response Surface Methodology

As the response surface is an approximation of the true response of a product or process, the optimum point determined from a response surface is not necessarily the optimum of the actual process. So to determine the optimum with sufficient accuracy, additional DoE is performed and a new optimum is determined. The design space for this new DoE is centered around the previous optimum and smaller than initial design space. In this way successive DoE are performed around the optimums within a design space which is reduced in each successive DoE. This procedure is repeated several times until the following condition is satisfied (Strömberg, 2009).

(29) Here x denotes the optimum design point determined at kth

iteration and y denotes the response calculated at this point. is a small number to check the convergence.

This procedure of performing successive DoE and determining new improved optimums is referred as Successive Response Surface Methodology (SRSM). The varying design space is usually called Region of Interest (RoI). The principle of SRSM is presented in Figure 18.

SRSM is an iterative method. Initially a DoE for a specified design space is performed, response surface is built and optimum point is determined from this surface. Then a new DoE is performed for a Design space that is smaller than initial design space and centered on the optimum point. This procedure narrows down the search for optimum point. The centering of new RoI around the optimum and its gradual reduction in size is usually called Panning and Zooming technique. This technique is very well presented by Gustafsson and Strömberg (2008).

(22)

RoI 1 RoI 2 RoI 4 Opt 1 Opt 2 Opt 3 Opt 4 Variable 2 RoI 3 Variable 1

(23)

2.3 Contact Mechanics

When two solid bodies touch each other at one or more points, localized stresses are developed in the body and deformation takes place.

Signorini’s Contact Conditions

Signorini’s contact conditions are very important while studying the contact phenomenon. These conditions are as follows

1. A node on elastic body cannot penetrate the obstacle. 2. Normal contact force cannot be negative.

3. Normal contact force should be zero if the node is not in contact and greater than zero if contact establishes.

To express these conditions mathematically, consider that an elastic body is deformed by an external force F and the deformation is restricted by a rigid body, which is fixed, as shown in Figure 19

Rigid Body

Figure 19: An elastic body obstructed by a rigid body Elastic Body

A set of nodes are defined on the boundary of elastic body where the two bodies are expected to contact each other. For each node a normal direction is determined. Now let denote the displacement of node and denote the initial normal distance between the node and rigid body. Then the first Signorini’s contact condition can be expressed a s.

. 0 (30)

Let denote the normal contact force for a node, and then second condition can be expressed as

(24)

0 (31) And the third condition can be x e pre sed as s

. 0 (32)

For all contact nodes, these equations can be expressed as o (33) In these equations • N de •

notes the matrix containing for all contact nodes fo

denotes the matrix containing r all nodes • denotes the matrix containing for all nodes • denotes the matrix containing for all nodes • o denotes the Hadamard product.

The last two conditions can also be derived by minimizing the potential energy in the elastic body while considering the no-penetration of rigid body as constraint. Mathematically it can be expressed as

min Π ₍₃₄₎

Here the potential energy can be describe as d

Π (35)

The solution to this problem is given by the following KKT conditions (Strömberg, 2009).

, , o

(36) Here is a column vector containing Lagrange multiplier for each contact constraint. These are interpreted as the contact forces in normal direction. Normal distance can be represented as . . Signorini’s contact conditions are illustrated in Figure 20

(25)

Contact Clearance

Contact Pressure

0 0

Figure 20: Signorini’s Contact Conditions

2.3.1 Contact Algorithms

Following algorithms are usually used to solve the contact problems. 1. Penalty Method

2. Uzawa’s method 3. Newton’s method 4. Smooth approximation

Penalty Method

In penalty formulation, Signorini’s contact conditions are replaced by following (Strömberg, 2009).

(37) Here the notation used can be explain by ed

| | 2

r is a penalty coefficient. By letting ∞ penalty model resembles Signorini’s contact conditions. But this does not work in practice. With this method the contact force is proportional to the penetration distance, so it always implies a certain penetration. Penalty method is illustrated in Figure 21

(26)

Uzawa’s Method

This method is described as follows (Strömberg, 2009). 1. Let r > 0 and a small number. Set 0 and 2. Solve

3. Update the contact forces according to

4. If then stop otherwise let 1 and go to step 2

Newton’s Method

This method is described as follows (Strömberg, 2009). ,

First introduce

(38) Following search criterion can be derived by aylo T r expansion.

(39) Then Newton’s approach can be described as

1. Set 0 and pick a small . Guess an initial state 2. Determine

3. Then find

4. If then stop otherwise let 1 and go to step 2

Smooth Approximation

There is a drawback in the Newton method that Projection is not smooth. This is overcome in Smooth Approximation by replacing Signorini’s contact conditions by a differentiable relation.

(27)

2.4 Super Elliptical Curve

Super elliptical curve can be described by the following equation. The symbols used are illustrated in Figure 22

1 (40)

This can also be expressed as

r

Rx Y

Ry

X

Figure 22: Illustration of Symbols used in Super Elliptical curve equation The super elliptical curve equ tioa n can a so be ex resl p sed as

(41)

Length of Curve

The length of curve between X = A and X= B can be found from

1 (42)

Y

X can be found be taking derivative of equatio 1) n

dX (4 dX dX (43) By chain rule dX dY . dY dX n dY dX (44) So dY dX / ₍₄₅₎

(28)

Instantaneous Radius

Instantaneous radius, r can be dervied from super elliptical curve equation by expressing X and Y in terms of r. So super elliptic l cura ve can be written as

1 (46)

And it can be written as

(47)

2.5 Polynomial Curve

A curve expresed by a polynomial is called polynomial curve. A Polynomial curve where radius is the function of angle and order is four can be described by the following equation. The the symbols used are illustr t d in Fa e igure 23

r α C C α C α C α C α

r

Here C , C , C , C and C are coefficients.

(48)

α

γ

r(α)

Figure 23: Illustration of Symbols used in Polynomial curve equation

Length of Curve

The length of can be found from

(29)

3 Implementation

First of all existing geometry of the part is analyzed to see the contact pressure distribution along the length of outer edge of the ring. For this purpose commercial FE software Abaqus (Abaqus, 2007) is used. Because the geometry of the ring is symmetric about horizontal axis, as shown in Figure 24, only the half of ring will be analyzed in Abaqus. It will reduce the computation time significantly.

Figure 24: Symmetric Geometry of Snap Ring

The cross-section of the ring is a full circle. It can be represented as a semi circle by using symmetry conditions. The cad model of snap ring modeled in Abaqus is shown in Figure 25.

Figure 25: Snap Ring representation in Abaqus

(30)

The housing will be modeled as a rigid body, because it is much stiffer than the snap ring. It is represented by a single surface as shown in Figure 26

Figure 26: Representation of Housing as a Surface

Augmented Lagrange method is selected to solve the contact problem between housing and ring. First order elements are used for meshing the ring as contact algorithm may not work properly with higher order elements. Figure 27 shows the contact pressure graph obtained by simulating in the Abaqus

(31)

For this simulation 17164 elements were generated for the snap ring and maximum contact pressure of 107 MPa was observed as shown in the graph. True representation of the parts in 3D is computationally expensive. So it was thought it could be time saving if 2D representation can be sufficient to represent the problem. Figure 28 shows the contact pressure graph obtained by simulating in the Abaqus for 2D representation of the problem.

Figure 28: Contact pressure graph for 2D representation

For this simulation 1864 elements were generated for the snap ring and maximum contact pressure of 24.6 MPa was observed as shown in the graph.

By comparing the contact pressure graphs for 3D and 2D representations, it is obvious that although scale of magnitudes are different but the overall profile of the contact pressure is almost same of both represenatations.

APython script is developed for this problem by recording macros in the Abaqus. A multi-point spline is used to build the base curve of snap ring, while recording macros. The spline curve is flexible enough to represent any geometry required for this project. The python script is modified so that it can read a text file containing coordinates of points of the base curve and build the spline according to this coordinate data. This script also performs the analysis in Abaqus and writes maximum contact pressure information for each simulation in a text file.

A circular arc is used as base curve in the original design of the snap ring. By looking at the contact pressure graph, it is apparent that a curve with constant radius is not sufficient for obtaining the uniform contact pressure distribuition. To obtain the uniform contact pressure distribuition, two type of curves were considered i.e. these curves were used as base curve for creating the geometry of snap ring. These curves are as follows

1. Super Elliptical 2. Polynomial

(32)

3.1 Super Elliptical Curve

Figure 29 shows the snap ring in compressed state (deformed shape) after simulating in Abaqus. A Gap can be seen between the snap ring and housing in the right quarter. But, on the left quarter there is no gap between the ring and housing.

Figure 29: Compressed state showing the gap

The gap shown in Figure 29 implies that for super elliptical base curve, different parameters should be used for the right and left quarters of base curve. This can be expressed graphically as shown in Figure 30.

Rx2 Ry Rx1 XL X Y

(33)

Here is the x-coordinate of end point of the curve. So the equation for left and right quarters can be written as

1

(50)

So optimization problem can be expressed as follows.

min a s.t. 1 m x 1 (51)

Here represents the length of original snap ring. For given values of , , , and , is computed from the constraint in (51). Then the coordinates for finite number of points at equal angular distances are computed by using equations in (50).

3.2 Polynomial Curve

In the polynomial considered here, radius is the function of angle and order is three. The equation can be written as follows. The symbols used are illustrated in Figure 23.

r α C C α C α C α (52)

So optimization problem can be expressed as follows. min max

s. t. C C α C α C α

For given values of C , C , C and C , is computed from the constraint in (53)

(53). The radius at finite number of points at equal angular distances is computed by using equation in (52). Then x and y coordinates are calculated by the radius and angle information at these points.

Matlab coding was developed to generate the coordinates of base curves and export them as text file. To evaluate the integrals in (51) and (53), built-in command “quad” was used.

(34)

3.3 Optimization Routine

Optimization routine is written in Matlab. It is used to determine the parameters in the curve equations for optimum geometry and can be described as follows. Figure 31 shows the routine as a flowchart.

1. A suitable start point and initial range are specified.

2. Central Composite Design is used to decide the combination of design variables. For these set of design variables, coordinates of curves are generated in Matlab, while length constraint is implemented. These coordinates are then exported in a text file.

3. A python script reads the text file containing the coordinates and performs the analysis in Abaqus. Then the script writes maximum contact pressure information for each simulation in a text file.

4. The text file containing maximum pressure information is imported in Matlab. This information combined with the DoE variable combinations are used to build the quadratic response surface. Least squares method is used to fit the response surface. The optimum of this response surface is determined by using the quadratic programming.

5. If 1, where k is the iteration number, check the convergence according to the following criteria (Gustafsson & Strömberg, 2008).

Here denotes the optimum design point determined at kth iteration and RSM denotes the contact pressure calculated at this point by the response surface. FEA denotes contact pressure calculated by FEA simulation, at kth iteration, in Abaqus. If and then the optimum has converged and is the optimum point. If or then solution has not converged and further iteration is needed. Here is a small number to check the convergence.

6. Determine new RoI by Pan and Zoom technique as presented by Gustafsson and Strömberg (2008) and go to step 2.

(35)

In this routine, tasks of building of response surface and computation of optimum are performed by developing the code in Matlab. This optimization is also compared with commercial software, HyperStudy. In that case, step 4 is modified such as DoE variable combinations and maximum contact pressure information is exported as a text file in a format that is readable to HyperStudy. HyperStudy is launched in batch mode from Matlab. Moving Least Squares method is chosen with second order approximation to perform optimization. After optimization, the optimum design point and optimum response are imported into Matlab. Other steps are same in both cases.

Set Start Point and initial Range

Setup a Design Of Experiments

Solve with FEA

Build Response Surface and Optimize

Yes Check Convergence

Pan and Zoom No

Stop

Figure 31: Flowchart for Optimization Routine

Matlab provides the built-in functions to generate the DoE and solve quadratic programming problems. “ccdesign” function was used to generate central composite design and “quadprog” function was used to solve the quadratic programming problem.

(36)

4 Results

In this section, results obtained from the successive response surface optimization are described. It is important to mention that sometimes the results of iterations, even after convergence, have been presented for comparison purposes.

4.1 Super Elliptical Base Curve

It can be concluded from the Figure 29 that for uniform contact pressure distribution.

2

2 (54)

Therefore it is supposed that just by finding suitable values for these two variables, the target of uniform contact pressure with lower maximum magnitude can be reached. So an optimization study is performed for these two design variables and all others are kept constant.

4.1.1 2 Design Variables

Following are chosen as the starting point for this successive response surface optimization.

2.0

2.2 (55)

Figure 32 shows the graph for the maximum contact pressure at the location of optimum design point ( ) which is determined from the response surface at every iteration by quadratic programming. Solid line shows the contact pressure as computed from the response surface and dotted line shows the contact pressure as computed by the actual simulation in Abaqus at . During initial iterations, the contact pressure computed from response surface differs much from the actual simulation results in Abaqus. It can be expected due to a large design space in the beginning. But as the design space (or RoI) continues to reduce, the results predicted by the response surface come closer to the actual simulation results.

(37)

0 2 4 6 8 10 12 14 16 18 1 2 3 4 5 6 7 8 9 10 11 12 13 Maximum Contact Pressure[MPa ] No. of Iterations Response Surface Simulation Result

Figure 32: Maximum Contact Pressure at the optimum design point for 2 variables. Solid line shows results as computed by Quadratic Programming from response surface. Dotted line shows the results obtained by simulation in Abaqus at optimum design point.

At iteration 13, convergence was achieved, at minimum contact pressure of 8.5 MPa with the following magnitu s de of design variables.

1.994609

2.108754 (56)

Figure 33 shows the graph of contact pressure for this optimum design point.

(38)

Figure 34 shows the snap ring in compressed state (deformed shape) for this optimum design point. A gap can be seen between the snap ring and housing.

Figure 34: Compressed state showing the gap

Although the maximum contact pressure has decreased to 8.5 MPa but still it is not distributed over the whole ring. Therefore it was decided to run the optimization study for three design variables, where the third variable will be the Ry.

It was supposed that by increasing the value of Ry, gap can be eliminated or reduced. It was also supposed that new optimum for and will be near the optimums obtained while optimizing for 2 variables. So following are chosen as the starting point for this successive response surface optimization.

2.0 2.1 Ry 0.1

(57)

Figure 35 shows the graph for the maximum contact pressure at the location of optimum design point ( ) at every iteration. During initial iterations, contact pressure computed from the response surface differs much from the actual simulation results in Abaqus. But as the design space reduces, the results predicted by the response surface come closer to the actual simulation results, but this behavior does not persist for further iterations. The difference keeps on changing. It can be concluded that the relation between the response and design variables is too nonlinear, both at global and local levels, that second order response surface cannot approximate it with sufficient accuracy.

(39)

0 2 4 6 8 10 12 14 16 18 1 2 3 4 5 6 7 8 9 10 11 12 Maximum Contact Pressure[MPa] No. of Iterations Response Surface Simulation Result

As the response surface is unable to represent the true response, so it is necessary to look at the results of DoE at individual iterations. By looking at the results of DoE at iteration 8, lowest response of 6.57 MPa was observed for the following information of design variables.

1.991253 2.097422 0.099572

(58)

Figure 36 shows the graph for this design point. This design point is different from the optimum point calculated by the response surface.

(40)

Although maximum contact pressure has decreased to 6.57 MPa but it is not distributed over the whole ring and gap, between ring and housing, is still there. Therefore it was decided to run the optimization study for five design variables. The knowledge gained from the optimization studies for 2 and 3 variables will be utilized to decide the start point.

It was supposed that by decreasing Rx₁ and Rx₂, gap can be eliminated or reduced. It was also supposed that new optimum for , and Ry will be near the optimums obtained while optimizing for 3 variables. So following are chosen as starting point for this successive response surface optimization.

2.0 2.1 0.0995 0.0985 0.0985 (59)

Figure 37 shows the graph for the maximum contact pressure at the location of optimum design point ( ) at every iteration. The contact pressure computed from response surface differs from the actual simulation results in Abaqus, throughout all the iterations. The difference between both keeps on changing, furthermore contact pressure computed from response surface shows negative values that is unrealistic. It can be concluded that the relation between the response and design variables becomes more and more nonlinear with the increase in design variables. ‐6 ‐4 ‐2 0 2 4 6 8 10 12 1 2 3 4 5 6 7 8 Maximum Contact Pressure[MPa ] No. of Iterations Response Surface Simulation Result

(41)

As the response surface is unable to represent the true response, so it is necessary to look at the results of DoE at individual iterations. By looking at the results of DoE at iteration 5, minimum contact pressure of 4.51 MPa was found for following information of design variables.

1.986683 2.079620 0.096846 0.097973 0.098760 (60)

Figure 38 shows the graph for this design point. The optimum design point which is computed from response surface has the response of 5.5 MPa.

Figure 39 shows the graph for the design point shown in (60) when simulated for 3D representation, here the maximum contact pressure is 17.66 MPa.

Figure 38: Contact Pressure graph for 2D representation when 5 variables are considered for optimization

(42)

Figure 39: Contact Pressure graph for 3D representation, when 5 variables are considered for optimization

Although maximum contact pressure has decreased significantly as compared to original design but still it is not distributed over the whole ring. Figure 40 shows the snap ring in compressed state (deformed shape) for this optimum design point. A gap can be seen between the snap ring and housing. Although the gap is considerably smaller than previous designs but it is still there.

(43)

In general, FEA softwares are considered more reliable to compute the von Mises stresses as compared to the contact stresses. So another optimization study was performed with the objective of minimizing the von Mises stress instead of contact stress to see if it can propose a better design of snap ring. Figure 41 shows the result of this optimization study for the start point shown in (59). It was observed that the optimum designs proposed by minimizing of maximum von Mises stress problem are different from those of minimizing of maximum contact pressure problem. After five iterations, significant gaps between the ring and housing begin to appear. ‐800 ‐600 ‐400 ‐200 0 200 400 1 2 3 4 5 6 7 8 Maximum Von Mises Stress[MPa] No. of Iterations Response Surface Simulation Result

Figure 41: Maximum von Mises Stress at the optimum design point for 5 variables. Solid line shows results as computed by Quadratic Programming from response surface. Dotted line shows the results obtained by simulation in Abaqus at optimum design point.

4.1.4 Comparison with HyperStudy

Figure 42 shows the result for the start point shown in (59), while HyperStudy is used for optimization. Convergence is achieved at 5th

iteration and optimums computed during iterations 3, 4 and 5 from response surface confirm with the actual simulation results in Abaqus. Minimum pressure of 6.7 MPa is obtained at convergence.

(44)

‐6 ‐4 ‐2 0 2 4 6 8 10 12 14 1 2 3 4 5 6 7 8 Maximum Contact Pressure[MPa ] No. of Iterations Response Surface Simulation Result

Figure 42: Maximum Contact Pressure at the optimum design point for 5 variables. Solid line shows results as computed by HyperStudy. Dotted line shows the results obtained by simulation in Abaqus at optimum design point.

By looking at the results of DoE at iteration 8, minimum contact pressure of 4.93 MPa was found for the following information of design variables.

1.998974 2.098207 0.098602 0.098604 0.098761 (61)

By comparing the graphs in Figure 42 and Figure 37, it can be concluded that HyperStudy performed better to attain convergence but the optimum found is not the real minimum.

Figure 43 shows the result of another optimization study for the start point shown in (59), when HyperStudy is used for optimization but instead of maximum contact pressure, maximum von Mises stress is minimized. It was observed that the optimum designs proposed by minimizing of maximum von Mises stress and maximum contact pressure problem are different from the first iteration.

(45)

175 180 185 190 195 200 205 210 215 220 1 2 3 4 5 6 7 8 Maximum Von Mises Stress[MPa] No. of Iterations Response Surface Simulation Result

Figure 43: Maximum von Mises stress at the optimum design point. Solid line shows results as computed by HyperStudy. Dotted line shows the results obtained by simulation in Abaqus at optimum design point.

(46)

4.2 Polynomial Base Curve

To determine the important variables in polynomial curve equation for response surface study, few simulations were performed where all coefficients were fixed and only was changed in the polynomial equation. It was found that does not have much influence on the contact pressure distribution over the snap ring although it changes the maximum contact pressure. So it was decided to keep its magnitude fixed at the original circular radius of the snap ring and consider the other three coefficients i.e. , and as design variables during the successive response surface optimization.

As compared to the super elliptical curve, it was difficult to guess the start point for design variables. So it was decided to run the optimization study over a large design space to explore it. During the investigation of the results of this study, it was found that there are some designs which give lower contact pressure but are not feasible i.e. they had large gaps between snap ring and housing in the compressed state. So another optimization study with the start point which would generate designs of snap ring that are most probably feasible is performed. The following was chosen as starting point for this successive response surface optimization.

100/10 250/10 200/10

(62)

Figure 44 shows graph for the maximum contact pressure at the location of optimum design point ( ) at every iteration. It is evident that maximum contact pressure is decreasing rapidly in every iteration. For the initial five iterations, maximum contact pressure computed from the response surface and actual simulation results are nearly same but later they differ. It is interesting to note that when design space is large, quadratic response surface can calculate the response value with sufficient accuracy but as the response value reaches optimum, quadratic response surface cannot predict the response accurately. It can be concluded that when the response is near the optimum, the shape of true response surface becomes highly non-linear that a second order response surface approximation is not sufficient to represent the true response surface.

(47)

‐5 0 5 10 15 20 25 30 35 1 2 3 4 5 6 7 8 9 10 11 Maximum Contact Pressure [MPa] No. of Iterations Response Surface Simulation Result

Figure 44: Maximum Contact Pressure at the optimum design point. Solid line shows results as computed by Quadratic Programming from response surface. Dotted line shows the results obtained by simulation in Abaqus at optimum design point.

During the study of DoE results at iteration 6 and 7 it was found that design points with the lower responses are scattered arbitrarily in the design space. It might be due to the presence of multiple stationary points in the design space. While looking at the DoE results for iteration 7, the following design point was found with lowest response value of 1.74 MPa.

0/10 750/10

55/10

(63)

Figure 45 shows the graph of contact pressure for this design point.

(48)

Although maximum contact pressure has reduced significantly for this design point as compared to the original design with circular arc as base curve, but still the pressure is not very well distributed over the whole ring, especially at ends. This requirement is difficult to be expressed in the form of a constraint and hence visual inspection of contact pressure graph is necessary. Therefore contact pressure graphs for some other design points with lower contact pressure during iterations 6, 7 and 8 were drawn. But none has pressure distributed over the whole length of the snap ring. To fulfill this requirement, some new design points were guessed by looking at the results for the design points where contact pressure was low. By little experimentation, the following design point was found with maximum contact pressure of 2.06 MPa and better distribution over the whole ring.

1 0/10 605/10

30/10

(64)

Figure 46 shows the graph for this new design point.

Figure 46: Contact Pressure graph

It was thought that a better design point with the lower maximum contact pressure and better distribution can be found in the vicinity of this design point. So to explore the design space around this design point, another successive response surface optimization was performed with this design point as starting point. The range about this start point was set to be very small. After a few iterations, following point was found with maximum contact pressure of 1.46 MPa and better distribution o rve th whole ring. e

8.053846/10 605.939638/10 30.641689/10

(49)

Figure 47 shows the graph for this new improved design point. Figure 48 shows the graph for same point when simulated for 3D representation; here the maximum contact pressure is 5.673 MPa.

Figure 47: Contact Pressure graph for 2D representation

Figure 48: Contact Pressure graph for 3D representation

Figure 49 shows the result of another optimization study for the start point shown in (62) but instead of maximum contact pressure, maximum von Mises stress is minimized. It was observed that the optimum designs proposed by minimizing of maximum von Mises stress and maximum contact pressure problem are same until sixth iteration, and after that they begin to differ. After 10th iteration, unfeasible designs with large gap between the ring and housing appear.

(50)

0 50 100 150 200 250 300 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Maximum Von Mises Stress[MPa] No. of Iterations Response Surface Simulation Result

Figure 49: Maximum von Mises stress at the optimum design point. Solid line shows results as computed by Quadratic Programming from response surface. Dotted line shows the results obtained by simulation in Abaqus at optimum design point.

4.2.1 Comparison with HyperStudy

Figure 50 shows the result for the start point shown in (62) when HyperStudy is used for optimization. It is evident that maximum contact pressure is decreasing rapidly in each iteration. By comparing the results with those shown in Figure 44, it is apparent that results obtained by using HyperStudy are almost same as for quadratic programming. By looking at the results of DoE at iteration 7, minimum contact pressure of 2.3 MPa was found for following values of design variables.

44.4114/10 830/10

75/10

(51)

‐5 0 5 10 15 20 25 30 35 1 2 3 4 5 6 7 8 9 10 11 Maximum Contact Pressure[MPa ] No. of Iterations Response Surface Simulation Result

Figure 50: Maximum Contact Pressure at the optimum design point. Solid line shows results as computed by HyperStudy. Dotted line shows the results obtained by simulation in Abaqus at optimum design point.

Figure 51 shows the result of another optimization study for the start point shown in (62) when HyperStudy is used for optimization but instead of maximum contact pressure, maximum von Mises stress is minimized. It was observed that the optimum designs proposed by minimizing of maximum von Mises stress and maximum contact pressure problem are different from the first iteration.

0 50 100 150 200 250 300 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Maximum Von Mises Stress[MPa] No. of Iterations Response Surface Simulation Result

Figure 51: Maximum von Mises stress at the optimum design point. Solid line shows results as computed by HyperStudy. Dotted line shows the results obtained by simulation in Abaqus at optimum design point.

(52)

5 Conclusion

After the optimization studies it was concluded that although the super elliptical base curve has reduced the peak contact pressure significantly, almost by 6 times as compared to the original design with circular arc, but it cannot give the geometry which has contact pressure distributed all over the snap ring.

The polynomial curve showed better results as compared to the super elliptical curve. It has reduced the maximum contact pressure almost by 19 times and contact pressure is distributed all over the snap ring. But the solution is not as good as the solution presented in Strömberg (2007). This might be due to the finite element treatment of the contact problem. In Strömberg (2007), an in-house software was utilized but in this work Abaqus has been used.

For the super elliptical curve, the optimization routine did not work well both with quadratic programming and HyperStudy. It can be concluded that the relation between the response and design variables is too nonlinear, both at global and local levels, that second order response surface cannot approximate it with sufficient accuracy.

For the polynomial curve, the optimization routine worked well both with quadratic programming and HyperStudy globally. But locally it did not give accurate computations. It might be because the shape of true response surface locally becomes highly non-linear that a second order response surface approximation is not sufficient to represent it.

A third order response surface or more advanced like Kriging or Neural Networks might be needed to represent the true response accurately.

(53)

6 References

Abaqus (2007). Version 6.7-4 [Computer software]. ABAQUS Inc., Providence, RI, USA.

Box, G.E.P. & Wilson K.B. (1951). On the experimental attainment of optimum conditions. Journal of the Royal Statistical Society, 13(1), 1–45.

Gustafsson, E. & Strömberg, N. (2008). Shape optimization of castings by using successive response surface methodology. Structural and Multidisciplinary

Optimization, 35(1), 11-28.

Hyperworks (2009). Version 10.0 [Computer software]. Altair Engineering Inc., Troy, MI, USA.

Jansson, T. (2005) Optimization of sheet metal forming processes (Ph.D. Thesis). Department of Mechanical Engineering, Linköping University, Sweden. Kutner, M.H., Nachtsheim, C.J. & Neter, J. (2004). Applied Linear Regression

Models. London: McGraw-Hill Education – Europe.

Myers, R.H. & Montgomery, D.C. (2002). Response Surface Methodology: Process

and Product Optimization Using Designed Experiments. New York: John

Wiley & Sons Inc.

Stander, N. & Craig, K.J. (2002) On the robustness of simple domain reduction scheme for simulation-based optimization. Engineering Computations, 19(4), 431–450.

Strömberg, N. (2007). What is the Optimal Shape of a Snap Ring? Proceedings of

the 6th International Conference on Multibody Systems, Nonlinear Dynamics, and Control. ASME, Las Vegas.

Strömberg, N. (2009, fall semester). Optimization Driven Design, Class Lecture. Jönköping University, Sweden.

(54)

7 Appendix

Dimensions of Original Snap Ring and Housing

Figure 52: Free State of Snap Ring

MINIMIZING CONTACT STRESSES IN AN ELASTIC RING BY RESPONSE SURFACE OPTIMIZATION

MINIMIZING CONTACT STRESSES IN AN

ELASTIC RING BY RESPONSE SURFACE

OPTIMIZATION

Asim Rashid

THESIS WORK 2010

MINIMIZING CONTACT STRESSES IN AN

ELASTIC RING BY RESPONSE SURFACE

OPTIMIZATION

Abstract

Table of Contents

1 Introduction

1.1 Background

1.2 Purpose and aims

1.3 Delimits

1.4 Outline

2 Theoretical background

2.1 Response Surface Methodology

2.2 Successive Response Surface Methodology

2.3 Contact Mechanics

2.4 Super Elliptical Curve

2.5 Polynomial Curve

3 Implementation

3.1 Super Elliptical Curve

3.2 Polynomial Curve

3.3 Optimization Routine

4 Results

4.1 Super Elliptical Base Curve

4.2 Polynomial Base Curve

5 Conclusion

6 References

7 Appendix