Maximum in-plane stress of a plane with a through-hole is given by
$$ \sigma = \frac{Kp}{(D-d)t},$$where $t$ is the thickness of the plate, $p$ is the pressure applied, and $K$ is the stress concentration factor. $K$ is given by
$$K = 1.11 +1.11(\frac{d}{D})^{-0.18}.$$We must find the hole size that minimizes the $\sigma$. We shall make the notational simplification $\frac{d}{D} = x.$ We can note from the problem formulation that we cannot write $\sigma$ strictly in terms of $x$; we must assume a value for $D$. We shall thence assume $D$, as well as $p$, and $t$ are all equal to 1. Given our assumptions we can we the stress equation as
$$\sigma (1 - x) = K.$$This can be further simplified in terms of an explicit objective funtion as
$$\sigma = \frac{1.11 + 1.11 x^{-0.18}}{1 - x}.$$The objective function is plotted below.
In [1]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from __future__ import division
%config InlineBackend.figure_formats=['svg']
%matplotlib inline
plt.rc('pdf',fonttype=3) # for proper subsetting of fonts
plt.rc('axes',linewidth=0.5) # thin axes; the default for lines is 1pt
al = np.linspace( 0.05, 0.15, 500)
plt.plot(al, (1.11 + 1.11*al**(-0.18))/(1 - al), 'k')
plt.axis([0.05, 0.15, 3.09,3.17])
plt.title("Objective Function")
plt.ylabel("f")
plt.xlabel("x/D")
plt.show()
We can see (at least qualitatively), from the plot of the objective funtion that on the interval $0.05 < \alpha < 0.15$ the optimum value lies somewhere between 0.09 and 0.10, and the function evaluated in that range has an average value of about 3.10.
We shall proceed to find the minimum of the objective function by using both an equal interval search algorithem and a polynomial approximation.
In [18]:
#Equal Interval Search
#Erin Schmidt
#Adapted, with significant modification, from Arora et al.'s APOLLO
#implementation found in "Introduction to Optimum Design" 1st Ed. (1989).
import numpy as np
def func(al, count): #the objective function
count = count + 1
f = (1.11 + 1.11*al**(-0.18))/(1 - al)
return f, count
def mini(au, al, count): #evaluates f at the minimum (or optimum) stationary point
alpha = (au + al)*0.5
(f, count) = func(alpha, count)
return f, alpha, count
def equal(delta, epsilon, count, al):
(f, count) = func(al, count)
fl = f #function value at lower bound
#delta = 0.01 #step-size
#au = 0.15 #alpha upper bound
while True:
aa = delta
(f, count) = func(aa, count)
fa = f
if fa > fl:
delta = delta * 0.1
else:
break
while True:
au = aa + delta
(f, count) = func(au, count)
fu = f
if fa > fu:
al = aa
aa = au
fl = fa
fa = fu
else:
break
while True:
if (au - al) > epsilon: #compares interval size to convergence criteria
delta = delta * 0.1
aa = al #intermediate alpha
fa = fl #intermediate alpha function value
while True:
au = aa + delta
(f, count) = func(au, count)
fu = f
if fa > fu:
al = aa
aa = au
fl = fa
fa = fu
continue
else:
break
continue
else:
(f, alpha, count) = mini(au, al, count)
return f, alpha, count
#run the program
delta = 0.01
epsilon = 1E-3
count = 0
al = 0.01 # alpha lower bound
(f, alpha, count) = equal(delta, epsilon, count, al)
print('The minimum is at {:.4f}'.format(alpha))
print('The function value at the minimum = {:.4f}'.format(f))
print('Total number of function calls = {}'.format(count))
In [22]:
# Polynomial approximation (4-point cubic)
# -Erin Schmidt
import numpy as np
from math import sqrt
# make an array with random values between 0.05 and 0.15 with 4 entries
x = (0.05 + np.random.sample(4)*0.15)
# make an array of function values at the 4 points of x
def f(x): # the objective function
return (1.11 + 1.11*x**(-0.18))/(1 - x)
f_array = []
i = 0
while i <= len(x) - 1:
f_array.append(f(x[i]))
i += 1
# use the equations from Vanderplaats 1984 to solve coefficients
q1 = x[2]**3 * (x[1] - x[0]) - x[1]**3 * (x[2] - x[0]) + x[0]**3 * (x[2] - x[1])
q2 = x[3]**3 * (x[1] - x[0]) - x[1]**3 * (x[3] - x[0]) + x[0]**3 * (x[3] - x[1])
q3 = (x[2] - x[1]) * (x[1] - x[0]) * (x[2] - x[0])
q4 = (x[3] - x[1]) * (x[1] - x[0]) * (x[3] - x[0])
q5 = f_array[2] * (x[1] - x[0]) - f_array[1] * (x[2] - x[0]) + f_array[0] * (x[2] - x[1])
q6 = f_array[3] * (x[1] - x[0]) - f_array[1] * (x[3] - x[0]) + f_array[0] * (x[3] - x[1])
a3 = (q3*q6 - q4*q5)/(q2*q3 - q1*q4)
a2 = (q5 - a3*q1)/q3
a1 = (f_array[1] - f_array[0])/(x[1] - x[0]) - \
a3*(x[1]**3 - x[0]**3)/(x[1] - x[0]) - a2*(x[0] + x[1])
a0 = f_array[0] - a1*x[0] - a2*x[0]**2 - a3*x[0]**3
a = [a1, 2*a2, 3*a3] #coefficients of f'
# find the zeros of the f' polynomial (using the quadratic formula)
b = a2**2 - 3*a1*a3
X1 = (-a2 + sqrt(b))/(3*a3)
X2 = (-a2 - sqrt(b))/(3*a3)
print('roots = ', X1, X2)
# plot the results
plt.rc('pdf',fonttype=3) # for proper subsetting of fonts
plt.rc('axes',linewidth=0.5) # thin axes; the default for lines is 1pt
x = np.linspace( 0.05, 0.15, 500)
plt.plot(x, a0 +a1*x + a2*x**2 + a3*x**3, 'k--', label='Poly. approx.')
plt.plot(x, (1.11 + 1.11*x**(-0.18))/(1 - x), 'k', label='Objective func.')
plt.axis([0.05, 0.15, 3.09,3.17])
legend = plt.legend(loc='upper center', shadow=False, fontsize='large')
plt.ylabel("f")
plt.xlabel("x/D")
plt.show()
poly_root = [0.097620387704, 0.0985634827486, 0.0969340736066, \
0.098775097463, 0.102426814371, 0.101638472077, \
0.0991941169039, 0.095873175811]
print('polynomial root std. deviation = ', np.std(poly_root))
Both the equal interval search and the polynomial approximation seem to yeild results that agree to 3 decimal places. Both also return values squarely reside within our expected bounds (being between 0.09 and 0.10), which we determined by a qualitative examination of the plot of the original objective function.
The roots of our polynomial are real and unique. On the bounds $0.05<x<0.15$ the minimum of the function by the polynomial approximation is the first root, at $x=0.99$. Though this agrees well with the result obtained via the equal interval search method for this probelm, the polynomial approximation method appears to be sensitive to the initial guess points of $x$, even using the relatively accurate 4-point cubic approximation. The standard deviation of ten sets of randomly sampled $x$ points on our bounds is 0.0021.
The equal interval search might also suffer from choice of an initial point x, on the objective funtion, especially if the objective function is multi-modal. Practically speaking the hole size $d$ must be on the bounds $0<d<1$ to make phsyical sense. If $d \geq D$ or $d \leq 0$ the search function will return divide by zero errors. However, within the bounds of 0.05 and 0.15 the equal interval search function returns an average value of 0.0979.
In [ ]: