4.6 Exercises

Exercise 7: Circular clamped plate with concentrated single load

Figure 62: Schematic representation of a circular clamped plate with concentrated single load. The plate is clamped at $ R=R_0 $ and the load $ P $ is applied at $ R=0 $.

Figure 62 show a rigidly clamped plate with radius $ R_0 $. The plate is loaded with a single load $ P $ in the plate centre. The differential equation for the deflection $ W(R) $ is given by: $$ \begin{align} \frac{d^3 W}{d R^3} + \frac{1}{R}\frac{d^2 W}{d R^2} - \frac{1}{R^2}\frac{d W}{d R} = \frac{P}{2 \pi D \cdot R} \tag{4.138} \end{align} $$

The plate stifness $ D $ is given by: $ D = \frac{E t^3}{12(1-\nu)} $, where $ E $ is the modulus of elasticity, $ \nu $ is the Poissons ratio and $ t $ is the plate thickness. The boundary conditions are given by: $$ \begin{align} W(R_0) = \frac{d W}{d R}(R_0) = 0, \quad \frac{d W}{d R}(0) = 0 \tag{4.139} \end{align} $$

In which the two first are because it is rigidly clamped, and the last due to the symmetry. We introduce dimmensionless variables: $ r=\frac{R}{R_0} $, $ \omega (R) = \frac{16 \pi D W}{P R_0^2} $, so that equation (4.138) can be written: $$ \begin{align} \frac{d^3 \omega}{d r^3} + \frac{1}{r}\frac{d^2 \omega}{d r^2} - \frac{1}{r^2}\frac{d \omega}{d r} = \frac{8}{r}, \quad 0 < r < 1 \tag{4.140} \end{align} $$

and equation (4.139): $$ \begin{align} \omega(1) = \frac{d \omega}{d r}(1) = 0, \quad \frac{d \omega}{d r}(0) = 0 \tag{4.141} \end{align} $$

The analytical solution is given by: $$ \begin{align} \omega(r) = r^2 \left[2 \, ln(r) - 1\right] +1, \quad \frac{d \omega}{d r}(r) = 4 \, r \, ln(r) \tag{4.142} \end{align} $$

Pen and paper:
The following problems should be done using pen and paper:
a) We introduce the inclination $ \phi (r) = -\frac{d \omega}{d r} $ and insert into equation (4.140): $$ \begin{align} \frac{d^2 \phi}{d r^2} + \frac{1}{r}\frac{d \phi}{d r} - \frac{\phi}{r^2} = - \frac{8}{r}, \tag{4.143} \end{align} $$

with boundary conditions: $$ \begin{align} \phi(0) = \phi(1) = 0 \tag{4.144} \end{align} $$

discretize equation (4.143) with central differences. Partion the interval $ [0, 1] $ into $ N $ segments, so that $ h = \Delta r = \frac{1}{N} $, which gives $ r = h \cdot i $, $ i = 0, 1, \dots , N $. The coefficient matrix should be tridiagonal. Write out the discretized equation for $ i=1 $, $ i $ and $ i = N-1 $. Write the expressions for the diagonals, and the right hand side.

b) Choose $ N=4 $ ($ h=\frac{1}{4} $) in a) and solve the corresponding system of equations.

c) We will now find $ \omega (r) $ by integrating the equation $ \frac{d \omega}{d r} = - \phi (r) $ using Heuns method. Since the right hand side is independenpt of $ \omega $ , Heuns method reduces to the trapes method. (The predictor is not necessary).

Find the $ \omega- $ values in the points as in b).

d) equation (4.143) can also be written: $$ \begin{align} \frac{d}{d r} \left[\frac{1}{r}\frac{d}{d r}\left[r \cdot \phi(r) \right]\right] \tag{4.145} \end{align} $$

Discretize equation (4.145) using equation (2.45). The coefficient matrix should be tridiagonal. Write out the discretized equation for $ i=1 $, $ i $ and $ i = N-1 $. Write the expressions for the diagonals, and the right hand side.

e) Solve equation (4.140) and (4.141) directly by introducing a new independant variable $ z = ln(r) $ and show that equation (4.140) can be written: $ \omega'''(z)-2\omega''(z)=8r^2 \equiv 8e^{2z} $. Next guess a particular solution on the form $ \omega_p (z)= k \cdot z \cdot e^{2z} $, where $ k $ is a constant. Lastly decide the constants using (4.141).

Programing:
a) Write a python program that solves b) and c) from the pen and paper exercise, numerically. Experiment with finer segmentation. If you want you can download the python scelleton clampedPlate.py and fill in where applicable.

Hint 1.

pen and paper:

b) Solutions: $ \phi_1=\frac{142}{105}=1.3524 $, $ \phi_2 =\frac{48}{35}=1.3714 $, $ \phi_3=\frac{6}{7}=0.8571 $

c) Since $ \omega_0 $ is unknown, but $ \omega_4 = \phi_4 = 0 $ is known, first find $ \omega_3 $, then $ \omega_2 $, $ \omega_1 $ and $ \omega_0 $. Use the $ \phi $-values obtained in a).

Solutions: $ \omega_0=\frac{94}{105}=0.8953 $, $ \omega_1 =\frac{61}{84}=0.7262 $, $ \omega_2=\frac{27}{70}=0.3857 $, $ \omega_2=\frac{3}{28}=0.1071 $

e) Solution given in equation (4.142)

Hint 2.

programming:

# src-ch2/clampedPlate.py

import scipy 
import scipy.linalg
import scipy.sparse
import scipy.sparse.linalg

#import matplotlib; matplotlib.use('Qt4Agg')
import matplotlib.pylab as plt
#plt.get_current_fig_manager().window.raise_()
import numpy as np

#### set default plot values: ####
LNWDT=2; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT


def solve_phi(N, h):
    
    i_vector = np.linspace(0, N, N + 1) # vector containing the indice number of all nodes
    i_mid = i_vector[1:-1] # vector containing the indice number of all interior nodes nodes
    
    # mainDiag = ?
    # subDiag = ?
    # superDiag = ? 
    
    # RHS = ?

    # A = scipy.sparse.diags([?, ?, ?], [?, ?, ?], format='csc')
    
    # Phi = scipy.sparse.linalg.spsolve(A, RHS)
    # Phi = np.append(0, Phi)
    # Phi = np.append(Phi, 0)
    
    return Phi

def solve_omega(Phi, N, h):
    
    Omega = np.zeros_like(Phi)
    Omega[-1] = 0
    #for i in range(N - 1, -1, -1):
        # i takes the values N-2, N-3, .... , 1, 0
        #Omega[i] = ?

    return Omega
    

N = 3
r = np.linspace(0, 1, N + 1)
h = 1./N

Phi = solve_phi(N, h)

Omega = solve_omega(Phi, N, h)

Omega_analytic = r[1:]**2*(2*np.log(r[1:]) - 1) + 1
Omega_analytic = np.append(1, Omega_analytic) # log(0) not possible

Phi_analytic = - 4*r[1:]*np.log(r[1:])
Phi_analytic = np.append(0, Phi_analytic) # log(0) not possible


fig, ax = plt.subplots(2, 1, sharex=True, squeeze=False)

ax[0][0].plot(r, Phi, 'r')
ax[0][0].plot(r, Phi_analytic, 'k--')

ax[0][0].set_ylabel(r'$\phi$')

ax[1][0].plot(r, Omega, 'r')
ax[1][0].plot(r, Omega_analytic, 'k--')
ax[1][0].legend(['numerical', 'analytical'], frameon=False)
ax[1][0].set_ylabel(r'$\omega$')
ax[1][0].set_xlabel('r')

plt.show()