The rectangular element formulation is simple and effective, but it has a major drawback: it is defined for a perfect rectangle whose sides are aligned with the global Cartesian axes (x, y). Real-world geometries are complex, and it is often inefficient or impossible to mesh a complex shape using only perfect rectangles.

If we try to use the shape functions derived above for a general, distorted quadrilateral shape, they will fail to satisfy the required properties (e.g., a shape function may not be zero along an opposite edge). This is because the functions are explicitly tied to the x and y coordinates and the element’s specific dimensions, a and b.

To analyze arbitrarily shaped quadrilateral elements, we must introduce a new coordinate system that is independent of the element’s shape in the global system. This leads us to the concept of Natural Coordinates.

1. The Natural Coordinate System

The core idea of the natural coordinate system is to map any arbitrarily shaped quadrilateral in the global (x, y) coordinate system to a single, perfectly square “parent” element in a local (r, s) coordinate system.

• Global System (Real Element): The physical element with coordinates x and y. It can be distorted, and its sides may not be parallel.
• Natural System (Parent Element): A perfect square defined by coordinates r and s, where both r and s range from -1 to +1.

This mapping transforms a complex geometry into a simple, standardized one, which greatly simplifies the formulation, especially the integration.

2. Isoparametric Formulation: Mapping and Shape Functions

The genius of the isoparametric formulation is to use the very same shape functions to define the geometric mapping as are used to interpolate the displacement field.

2.1 Shape Functions in Natural Coordinates

On the (r, s) parent square, the shape functions for the 4-node element have a simple and universal form, analogous to the product method used for the rectangle:

• \(N_1(r,s) = \frac{1}{4}(1-r)(1-s)\)
• \(N_2(r,s) = \frac{1}{4}(1+r)(1-s)\)
• \(N_3(r,s) = \frac{1}{4}(1+r)(1+s)\)
• \(N_4(r,s) = \frac{1}{4}(1-r)(1+s)\)

These functions are always the same, regardless of the shape of the real element in the (x, y) plane.

2.2 The Geometric Mapping

The connection between the two coordinate systems is established by using these shape functions to interpolate the global coordinates from the nodal coordinates:

\[ \begin{aligned} x(r,s) = \sum_{i=1}^{4} N_i(r,s) x_i,\\ y(r,s)= \sum_{i=1}^{4} N_i(r,s) y_i \end{aligned} \]

This mapping allows us to find the global (x, y) coordinates corresponding to any point given by its natural coordinates (r, s).

3. The Jacobian Matrix: Relating Derivatives

To form the B matrix, we need the derivatives of the shape functions with respect to the global coordinates, x and y. However, our shape functions are defined in terms of the natural coordinates, r and s. We can relate these derivatives using the chain rule of partial differentiation:

\[ \frac{\partial N_i}{\partial r} = \frac{\partial N_i}{\partial x} \frac{\partial x}{\partial r} + \frac{\partial N_i}{\partial y} \frac{\partial y}{\partial r} \] \[ \frac{\partial N_i}{\partial s} = \frac{\partial N_i}{\partial x} \frac{\partial x}{\partial s} + \frac{\partial N_i}{\partial y} \frac{\partial y}{\partial s} \]

This system can be written in matrix form:

\[ \begin{Bmatrix} \frac{\partial N_i}{\partial r} \\ \frac{\partial N_i}{\partial s} \end{Bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\ \frac{\partial x}{\partial s} & \frac{\partial y}{\partial s} \end{bmatrix} \begin{Bmatrix} \frac{\partial N_i}{\partial x} \\ \frac{\partial N_i}{\partial y} \end{Bmatrix} \]

The matrix in this equation is the Jacobian Matrix, J.

\[ \mathbf{J} = \begin{bmatrix} J_{11} & J_{12} \\ J_{21} & J_{22} \end{bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\ \frac{\partial x}{\partial s} & \frac{\partial y}{\partial s} \end{bmatrix} \]

Since we need the derivatives with respect to x and y to form the B matrix, we must invert this relationship:

\[ \begin{Bmatrix} \frac{\partial N_i}{\partial x} \\ \frac{\partial N_i}{\partial y} \end{Bmatrix} = \mathbf{J}^{-1} \begin{Bmatrix} \frac{\partial N_i}{\partial r} \\ \frac{\partial N_i}{\partial s} \end{Bmatrix} \]

The terms of the Jacobian can be calculated by differentiating the mapping equations (e.g., \(\displaystyle \frac{\partial x}{\partial r} = \sum \frac{\partial N_i}{\partial r} x_i\)).

6. Integration in Natural Coordinates

The final step is to transform the stiffness matrix integral into the natural coordinate system. The differential area element also transforms:

\[ dA = dx dy = \det(\mathbf{J}) dr ds \]

The determinant of the Jacobian, det(J), acts as a scaling factor between the differential areas in the two systems. The element stiffness matrix integral becomes:

\[ \mathbf{K} = t \int_{-1}^{1} \int_{-1}^{1} \mathbf{B}(r,s)^T \mathbf{E} \mathbf{B}(r,s) \det(\mathbf{J}) \, dr ds \]

This integral has constant limits from -1 to 1, regardless of the real element’s shape. This standardized form is ideal for numerical integration techniques, such as Gaussian Quadrature, which is the standard method for evaluating these integrals in FEM software.