# Rectangular Plate Bending Element

This derivation follows the procedure given in Chapter 12 of "A First Course in the Finite Element Method, 4th Edition" by Daryl L. Logan.

We'll start by importing a few Python libraries that are useful for symbolic math, and initializing "pretty" printing.


```python
from sympy import symbols, Matrix, diff, integrate, simplify, factor, latex, init_printing
from IPython.display import display, Math
init_printing()
```

The plate width will be defined as $2b$, and the height will be $2c$ to be consistent with Figure 12-1. We'll set up some Sympy symbols to represent $b$ and $c$.


```python
b, c = symbols('b, c')
```

The plate is defined by four nodes specified in counter-clockwise order: i, j, m, and n. The local x-axis runs from node i toward node j, and the local y-axis runs from node i toward node n. Next we'll define the element's local displacement vector, $[d]$, at each node. There are 3 degrees of freedom at each node: $w$, $\theta_x$, and $\theta_y$.


```python
wi, theta_xi, theta_yi = symbols('w_i, theta_x_i, theta_yi')
wj, theta_xj, theta_yj = symbols('w_j, theta_xj, theta_yj')
wm, theta_xm, theta_ym = symbols('w_m, theta_xm, theta_ym')
wn, theta_xn, theta_yn = symbols('w_n, theta_xn, theta_yn')
d = Matrix([wi, theta_xi, theta_yi, wj, theta_xj, theta_yj, wm, theta_xm, theta_ym, wn, theta_xn, theta_yn])
display(Math('[d] = ' + latex(d)))
```

A 12-term polynomial displacement function will be assumed to define the out-of-plane displacement, w, at any point (x, y) in the plate's local coordinate system. The rotations about each axis are derivatives of this displacement:

$w = a_1 + a_2x + a_3y + a_4x^2 + a_5xy + a_6y^2 + a_7x^3 + a_8x^2y + a_9xy^2 + a_{10}y^3 + a_{11}x^3y + a_{12}xy^3$

$\theta_x = \frac{dw}{dy} = a_3 + a_5x + 2a_6y + a_8x^2 + 2a_9xy + 2a_{10}y^2 + a_{11}x^3 + 3a_{12}xy^2$

$\theta_y = -\frac{dw}{dx} = -a_2 - 2a_4x - a_5y - 3a_7x^2 - 2a_8xy - a_9y^2 - 3a_{11}x^2y - a_{12}y^3$

The negative sign on $\frac{dw}{dx}$ is required to be consistent with the right hand rule. These equations can be rewritten in matrix form as follows:

$[\psi] = [P][a]$

where $[\psi]$ is shorthand for  $\begin{bmatrix} w \\ \theta_x \\ \theta_y \end{bmatrix}$ and $[P]$ is defined as follows:


```python
x, y = symbols('x, y')
P = Matrix([[1, x, y, x**2, x*y, y**2, x**3, x**2*y, x*y**2, y**3, x**3*y, x*y**3],
            [0, 0, 1, 0, x, 2*y, 0, x**2, 2*x*y, 3*y**2, x**3, 3*x*y**2],
            [0, -1, 0, -2*x, -y, 0, -3*x**2, -2*x*y, -y**2, 0, -3*x**2*y, -y**3]])
display(Math('P = ' + latex(P)))
```

This equation for $[w]$ is valid for a single node. Evaluating $[P]$ at each node gives us a larger set of equations:

$[d] = [C][a]$

where $[C]$ is merely $[P]$ evaluated at each node, and $[d]$ is correpsondingly $[\psi]$ at each node. Knowing that the plate width is $2b$ and the plate height is $2c$, we can obtain the matrix $[C]$.


```python
C = Matrix([P, P, P, P])
C[0:3, 0:12] = C[0:3, 0:12].subs(x, 0).subs(y, 0)      # i-node @ x = 0, y = 0
C[3:6, 0:12] = C[3:6, 0:12].subs(x, 2*b).subs(y, 0)    # j-node @ x = 2b, y = 0
C[6:9, 0:12] = C[6:9, 0:12].subs(x, 2*b).subs(y, 2*c)  # m-node @ x = 2b, y = 2c
C[9:12, 0:12] = C[9:12, 0:12].subs(x, 0).subs(y, 2*c)  # n-node @ x = 0, y = 2c
display(Math('[C] = ' + latex(C)))
```

An important matrix that we will come back to later is the shape function matrix $[N]$, defined as:

$[N] = [P][C]^{-1}$

The closed form solution of $[N]$ for a rectangular plate is:


```python
N = P*C.inv()
display(Math('[N] = ' + latex(simplify(N))))
```

We can now solve for the $[a]$ matrix in terms of the nodal displacements:


```python
a = simplify(C.inv()*d)
display(Math('[a] = ' + latex(a)))
```

The next step is to define the curvature matrix:

$[\kappa] = \begin{bmatrix} -\frac{d^2w}{dx^2} \\ -\frac{d^2w}{dy^2} \\ -\frac{2d^2w}{dxdy} \end{bmatrix} = [Q][a]$

It should be recognized that $w/[a]$ is simply the first row of our $[P]$ matrix. Evaluating the derivatives in this expression gives $[Q]$ as follows:


```python
Q = Matrix([-diff(diff(P[0, :], x), x),
            -diff(diff(P[0, :], y), y),
            -2*diff(diff(P[0, :], x), y)])
display(Math('[Q] = ' + latex(Q)))
```

With $[Q]$ in hand we can now solve for the $[B]$ matrix which is essential for formulating the stiffness matrix $[k]$


```python
B = simplify(Q*C.inv())
display(Math('[B] = ' + latex(B)))
```

Now we form the constitutive matrix for isotropic materials, [D]. This matrix is analagous to the flexural stiffness of a beam EI.


```python
E, t, nu = symbols('E, t, nu')
Coef = E*t**3/(12*(1-nu**2))
D = Coef*Matrix([[1, nu, 0],
                 [nu, 1, 0],
                 [0, 0, (1-nu)/2]])
display(Math('[D] = ' + latex(D)))
```

Now we can calculate the stiffness matrix:

$[k] = \int_0^{2c} \int_0^{2b} [B]^T[D][B] dx dy$


```python
k = integrate(integrate(B.T*D*B, (x, 0, 2*b)), (y, 0, 2*c))
display(Math('[k] = {Et^3}/{12(1-\nu^2)}' + latex(simplify(k/Coef))))
```

The surface force matrix $[F_s]$ can be obtained from the shape function matrix. Since we're interested in the surface force matrix for uniform pressures in the direction of w,


```python
q = symbols('q')
Fs = integrate(integrate(N[0, :].T*q, (x, 0, 2*b)), (y, 0, 2*c))
display(Math('[F_s] = 4qcb' + latex(Fs/(4*q*c*b))))
print(Fs/(4*q*c*b))
```


```python

```

986c99e1495952c0dba475739b3da3b4ec45df55
ipynb
Jupyter Notebook
Derivations/Rectangular Plate Element.ipynb
geosharma/PyNite
b7419ac352b48d7cd92b1a34a30d0083e19cfb78
MIT
2019-04-12T05:30:43.000Z
2022-03-29T23:42:01.000Z
Derivations/Rectangular Plate Element.ipynb
mohashrafy/PyNite
efffccdbff6727d3b271ba2937e35892d9df8c00
MIT
2019-07-22T19:41:26.000Z
2022-03-31T22:18:32.000Z
Derivations/Rectangular Plate Element.ipynb
mohashrafy/PyNite
efffccdbff6727d3b271ba2937e35892d9df8c00
MIT
2019-02-07T12:02:01.000Z
2022-03-24T13:24:38.000Z
Qwen/Qwen-72B
1. YES
2. YES