Navier solution of a simply supported rectangular plate

navier_plate_interpolate.py

# coding: utf-8
"""Navier solution of a simply supported rectangular plate.

Author: Juan Luis Cano Rodríguez <juanlu001@gmail.com>

References
----------
* Timoshenko, S., & Woinowsky-Krieger, S. (1959). "Theory of plates and
  shells" (Vol. 2, p. 120). New York: McGraw-hill.

"""
from __future__ import division

import numpy as np
from numpy import pi, sin, cos
from scipy import integrate, interpolate

# In case JIT optimization is desired
try:
    from numba import jit
except ImportError:
    jit = lambda f: f

## Initial data

# Plate geometry
a = 1  # m
b = 1  # m
h = 50e-3  # m

# Material properties
E = 69e9  # Pa
nu = 0.35

# Series terms
MAX_M = 16
MAX_N = 16

# Computation points
NUM_POINTS = 1000

# Load
def q(x, y, a, b):
    """Uniform distributed load.

    """
    return 10.0  # Pa


## Computation

# Flexural rigidity
D = h**3 * E / (12 * (1 - nu**2))

# Coefficient
def a_mn(q, a, b, mm, nn):
    """Navier series coefficient.

    Parameters
    ----------
    q : function f(x, y, a, b)
        Load function.
    a, b : float
        Dimensions of the plate.
    mm, nn : integer
        Indices of the coefficient

    """
    def _inner(x, y, a, b, mm, nn):
        return q(x, y, a, b) * sin(mm * pi * x / a) * sin(nn * pi * y / b)

    res, _ = integrate.dblquad(_inner, 0, a, lambda x: 0, lambda x: b,
                               args=(a, b, mm, nn))
    return 4 * res / (a * b)

# Displacement field
x = np.linspace(0, a, num=NUM_POINTS)
y = np.linspace(0, b, num=NUM_POINTS)
xx, yy = np.meshgrid(x, y)

ww = np.zeros_like(xx)
for mm in range(1, MAX_M):
    for nn in range(1, MAX_N):
        ww += (a_mn(q, a, b, mm, nn) / (mm**2 / a**2 + nn**2 / b**2)**2 *
               sin(mm * pi * xx / a) * sin(nn * pi * yy / b) /
               (pi**4 * D))

# Creation of an interpolating function
w = interpolate.RegularGridInterpolator((x, y), ww)

# Maximum displacement
print("Maximum displacement: %.3e m" % w((a / 2, b / 2)))

# Solution plotting
import matplotlib.pyplot as plt

plt.contourf(xx, yy, ww)
plt.colorbar()
plt.xlabel("$x$ (m)")
plt.ylabel("$y$ (m)")
plt.savefig("navier_plate.png")

Solution

How can I run this?