Engineering at Alberta Courses » Quadrilateral Elements

Two Dimensional Solid Elements: Quadrilateral Elements

Quadrilateral Elements

Linear Quadrilateral Elements

It is usually easier to discretize a domain by using triangular elements; however, when the domain to be discretized has a regular shape, quadrilateral elements can be used to offer a more regular displacement discretization, which, in some cases, might offer a better approximation to the displacement shape. The bilinear quadrilateral element offers a bilinear displacement approximation where the displacement within the element is assumed to vary bilinearly (linear in two directions within the element). This element, however, behaves in a stiff manner when it is used to model linear strains (bending strains) since, as will be shown in the derivation, it is not possible to model pure bending (pure linear normal strain components in one direction) without introducing additional shear strains (usually termed parasitic shear strains).

Consider the plane quadrilateral shown in Figure 7. The bilinear displacement function for this element is assumed to have the following form:

(1) $\begin{equation*} u_{2D} = \begin{pmatrix} u \\ v \\ \end{pmatrix} = \begin{pmatrix} a_0 + a_1X_1 + a_2X_2 + a_3X_1X_2 \\ b_0 + b_1X_1 + b_2X_2 + b_3X_1X_2 \\ \end{pmatrix} \end{equation*}$

where $a_0$ , $a_1$ , $a_2$ , $a_3$ , $b_0$ , $b_1$ , $b_2$ , and $b_3$ are eight generalized degrees of freedom of the element. The displacement function in terms of the nodal degrees of freedom has the form:

$u_{2D} = \begin{pmatrix} u \\ v \\ \end{pmatrix} = \begin{pmatrix} N_1u_1 + N_2u_2 + N_3u_3 + N_4u_4 \\ N_1v_1 + N_2v_2 + N_3v_3 + N_4v_4 \\ \end{pmatrix}$

There are two ways to find the shape functions $N_i$ . The first way is to replace the generalized degrees of freedom with the nodal degrees of freedom in the first equation. The multipliers of the nodal degrees of freedom would then be the shape function. The following mathematica code does that for you:

View Mathematica Code

Clear[a, b]
u = a0 + a1*X1 + a2*X2 + a3*X1*X2;
v = b0 + b1*X1 + b2*X2 + b3*X1*X2;
Coordinates = {{X1 -> -a, X2 -> -b}, {X1 -> a, X2 -> -b}, {X1 -> a, X2 -> b}, {X1 -> -a, X2 -> b}};
Eq1 = u /. Coordinates[[1]];
Eq2 = u /. Coordinates[[2]];
Eq3 = u /. Coordinates[[3]];
Eq4 = u /. Coordinates[[4]];
Eq5 = v /. Coordinates[[1]];
Eq6 = v /. Coordinates[[2]];
Eq7 = v /. Coordinates[[3]];
Eq8 = v /. Coordinates[[4]];
sol = Solve[{Eq1 == u1, Eq2 == u2, Eq3 == u3, Eq4 == u4, Eq5 == v1, Eq6 == v2, Eq7 == v3, Eq8 == v4}, {a0, a1, a2, a3, b0, b1, b2, b3}]
u = u /. sol[[1]];
v = v /. sol[[1]];
N1 = FullSimplify[Coefficient[u, u1]]
N1 = FullSimplify[Coefficient[v, v1]]
N2 = FullSimplify[Coefficient[u, u2]]
N2 = FullSimplify[Coefficient[v, v2]]
N3 = FullSimplify[Coefficient[u, u3]]
N3 = FullSimplify[Coefficient[v, v3]]
N4 = FullSimplify[Coefficient[u, u4]]
N4 = FullSimplify[Coefficient[v, v4]]

View Python Code

import sympy as sp
from sympy import *
a,a0,a1,a2,a3,b,b0,b1,b2,b3,X1,X2,u1,u2,u3,u4,v1,v2,v3,v4 = symbols("a a_0 a_1 a_2 a_3 b b_0 b_1 b_2 b_3 X_1 X_2 u_1 u_2 u_3 u_4 v_1 v_2 v_3 v_4")
u = a0+a1*X1+a2*X2+a3*X1*X2
v = b0+b1*X1+b2*X2+b3*X1*X2
Coordinates = [{X1:-a,X2:-b},{X1:a,X2:-b},{X1:a,X2:b}, {X1:-a, X2:b}]
display("Coordinates: ",Coordinates)
Eq1 = u.subs(Coordinates[0])
Eq2 = u.subs(Coordinates[1])
Eq3 = u.subs(Coordinates[2])
Eq4 = u.subs(Coordinates[3])
Eq5 = v.subs(Coordinates[0])
Eq6 = v.subs(Coordinates[1])
Eq7 = v.subs(Coordinates[2])
Eq8 = v.subs(Coordinates[3])
s = solve((Eq1-u1,Eq2-u2,Eq3-u3,Eq4-u4,Eq5-v1,Eq6-v2, Eq7-v3, Eq8-v4),
          (a0,a1,a2,a3,b0,b1,b2,b3))
display("Solve: ",s)
u = u.subs(s).expand()
v = v.subs(s).expand()
display("u =",u,"v =",v)
N1_u = u.coeff(u1)
N1_v = v.coeff(v1)
N2_u = u.coeff(u2)
N2_v = v.coeff(v2)
N3_u = u.coeff(u3)
N3_v = v.coeff(v3)
N4_u = u.coeff(u4)
N4_v = v.coeff(v4)
display("N1u1 =",N1_u,"N1v1 =",N1_v,"N2u2 =",N2_u,
        "N2v2 =",N2_v,"N3u3 =",N3_u,"N3v3 =",N3_v,
        "N4u4=",N4_u,"N4v4 =",N4_v)

An alternate way is directly by realizing that $N_1$ has to satisfy that the condition that at node 1, $N_1=1$ , and is equal to zero on the lines $a-X_1=0$ and $b-X_2=0$ . The same can be applied for N_2 and N_3. These result in the following expressions for the shape functions:

$N_1 = \frac{(b - X_2)(a - X_1)}{4ab}$

$N_2 = \frac{(b - X_2)(a + X_1)}{4ab}$

$N_3 = \frac{(b + X_2)(a + X_1)}{4ab}$

$N_4 = \frac{(b + X_2)(a - X_1)}{4ab}$

The distribution of the shape functions on the element are illustrated in Figure 8.

Figure 8. Shape functions distribution in the bilinear quadrilateral element

The strain associated with the assumed displacement field has the following form:

$\varepsilon = \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \\ \end{pmatrix} = \begin{pmatrix} \frac{\partial u}{\partial X_1} \\ \frac{\partial v}{\partial X_2} \\ \frac{\partial u}{\partial X_2} + \frac{\partial v}{\partial X_1} \\ \end{pmatrix}$

$= \begin{pmatrix} \frac{\partial N_1}{\partial X_1} & 0 & \frac{\partial N_2}{\partial X_1} & 0 & \frac{\partial N_3}{\partial X_1} & 0 & \frac{\partial N_4}{\partial X_1} & 0 \\ 0 & \frac{\partial N_1}{\partial X_2} & 0 & \frac{\partial N_2}{\partial X_2} & 0 & \frac{\partial N_3}{\partial X_2} & 0 & \frac{\partial N_4}{\partial X_2} \\ \frac{\partial N_1}{\partial X_2} & \frac{\partial N_1}{\partial X_1} & \frac{\partial N_2}{\partial X_2} & \frac{\partial N_2}{\partial X_1} & \frac{\partial N_3}{\partial X_2} & \frac{\partial N_3}{\partial X_1} & \frac{\partial N_4}{\partial X_2} & \frac{\partial N_4}{\partial X_1} \\ \end{pmatrix} \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \\ u_4 \\ v_4 \\ \end{pmatrix}$

$= Bu_e$

Where, the explicit representation of the matrix $B$ is:

$B = \begin{pmatrix} -\frac{b - X_2}{4ab} & 0 & \frac{b - X_2}{4ab} & 0 & \frac{b+ X_2}{4ab} & 0 & -\frac{b + X_2}{4ab} & 0 \\ 0 & -\frac{a - X_1}{4ab} & 0 & -\frac{a + X_1}{4ab} & 0 & \frac{a + X_1}{4ab} & 0 & \frac{a - X_1}{4ab} \\ -\frac{a - X_1}{4ab} & -\frac{b - X_2}{4ab} & -\frac{a + X_1}{4ab} & \frac{b + X_2}{4ab} & \frac{a + X_1}{4ab} & \frac{b + X_2}{4ab} & \frac{a - X_1}{4ab} & -\frac{b + X_2}{4ab} \\ \end{pmatrix}$

Shear Locking

A quick glance at the matrix $B$ shows that the bilinear quadrilateral can model bending (linear strains) since the strains $\varepsilon_{11}$ and $\varepsilon_{22}$ contain linear expressions in $X_1$ and $X_2$ . However, it will be shown here that to model linear strains, the element predicts an associated shear strain $\gamma_{12}$ as well (termed parasitic shear strains). For simplicity, we will calculate the strains based on Equation 1. The strains have the form:

$\varepsilon_{11} = a_1 + a_3X_2$

$\varepsilon_{22} = b_2 + b_3X_1$

$\gamma_{12} = a_2 + a_3X_1 + b_1 + b_3X_2$

If we now assume pure bending state of strain with $a_3\neq 0$ . Then, the shear strain cannot be equal to zero since $a_3$ appears in the expression for $\gamma_{12}$ ! This behaviour is termed shear locking and causes the element to be too stiff under pure bending.

Stiffness Matrix

The constitutive relationship of plane linear elastic materials is defined using a $3\times 3$ matrix $C$ that depends on whether the material is in a plane strain or a plane stress state. For the sake of illustration, the plane strain state is chosen here, the material constitutive relationship matrix in that case is:

$\begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \\ \end{pmatrix} = \frac{E}{(1+\nu)(1-2\nu)} \begin{pmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \frac{1-2\nu}{2} \\ \end{pmatrix} \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \\ \end{pmatrix}$

The stiffness matrix of the linear strain triangle can be evaluated using the following integral (assuming a constant thickness $t$ :

$K^e = t \int_{-b}^b \int_{-a}^a B^T CB dX_1 dX_2$

The stiffness matrix has the dimensions of $8 \times 8$ , and the following Mathematica code can be utilized to view its components:

View Mathematica Code

Shapefun=Table[0,{i,1,4}];
Shapefun[[1]]=(b-x2)(a-x1)/4/a/b;
Shapefun[[2]]=(b-x2)(a+x1)/4/a/b;
Shapefun[[3]]=(b+x2)(a+x1)/4/a/b;
Shapefun[[4]]=(b+x2)(a-x1)/4/a/b;
B=Table[0,{i,1,3},{j,1,8}];
Do[B[[1,2i-1]]=B[[3,2i]]=D[Shapefun[[i]],x1];B[[2,2i]]=B[[3,2i-1]]=D[Shapefun[[i]],x2],{i,1,4}];
Cc=Ee/(1+nu)*{{(1-nu)/(1-2nu),nu/(1-2nu),0},{(nu)/(1-2nu),(1-nu)/(1-2nu),0},{0,0,1/2}};
Kbeforeintegration=t*FullSimplify[Transpose[B].Cc.B];
K=Integrate[Kbeforeintegration,{x2,-b,b},{x1,-a,a}];
K//MatrixForm

View Python Code

import sympy as sp
from sympy import *
a, b, x1, x2, nu, E, t = symbols("a b x_1 x_2 nu E t")
Shapefun = Matrix([(b-x2)*(a-x1)/4/a/b,(b-x2)*(a+x1)/4/a/b,(b+x2)*(a+x1)/4/a/b,(b+x2)*(a-x1)/4/a/b])
B = zeros(3,8)
for i in range(4):
    B[0,2*i] = B[2,2*i+1] = Shapefun[i].diff(x1)
    B[1,2*i+1] = B[2,2*i] = Shapefun[i].diff(x2)
display("B:", B)
Cc = E/(1 + nu)*Matrix([[(1-nu)/(1-2*nu),nu/(1-2*nu),0], 
                         [nu/(1-2*nu),(1-nu)/(1-2*nu), 0], 
                         [0, 0, 1/2]])  
KBeforeintegration = t * B.transpose()*Cc*B
K = Matrix([[simplify(integrate(KBeforeintegration[i,j],(x1,-a,a),(x2,-b,b)))for i in range (8)] for j in range (8)])
display("K: ", K)

Nodal Loads

The nodal loads due to distributed body forces and traction forces on the boundaries are equally distributed among the nodes. The same procedure followed for the triangular elements can be repeated to obtain the distribution shown in Figure 9. The following Mathematica code can be utilized for the calculations producing the distributions in Figure 9.

View Mathematica Code

Shapefun=Table[0,{i,1,4}];
Shapefun[[1]]=(b-x2)(a-x1)/4/a/b;
Shapefun[[2]]=(b-x2)(a+x1)/4/a/b;
Shapefun[[3]]=(b+x2)(a+x1)/4/a/b;
Shapefun[[4]]=(b+x2)(a-x1)/4/a/b;
Nn=Table[0,{i,1,2},{j,1,8}];
Do[Nn[[1,2i-1]]=Nn[[2,2i]]=Shapefun[[i]],{i,1,4}];
tn={t1,t2};
rb={rb1,rb2};
fetraction=Integrate[(Transpose[Nn].tn/.x1->-a),{x2,-b,b}]//MatrixForm
febodyforces=Integrate[(Transpose[Nn].rb),{x2,-b,b},{x1,-a,a}]//MatrixForm

View Python Code

import sympy as sp
from sympy import *
a, b, x1, x2, t1, t2, rb1, rb2 = symbols("a b x_1 x_2 t_1 t_2 rb_1 rb_2")
Shapefun = Matrix([(b-x2)*(a-x1)/4/a/b,(b-x2)*(a+x1)/4/a/b,(b+x2)*(a+x1)/4/a/b,(b+x2)*(a-x1)/4/a/b])
Nn= zeros(2,8)
for i in range(4):
    Nn[0,2*i] = Nn[1,2*i+1] = Shapefun[i]
display("Nn:", Nn)
tn = Matrix([t1,t2])
rb = Matrix([rb1,rb2])
integrand1 = (Nn.transpose()*tn).subs(x1,-a)
integrand2 = (Nn.transpose()*rb)
fetraction = Matrix([simplify(integrate(integrand1[i],(x2,-b,b)))for i in range (8)])
febodyforces = Matrix([simplify(integrate(integrand2[i],(x1,-a,a),(x2,-b,b)))for i in range (8)])
display("fetraction: ", fetraction)
display("febodyforces: ", febodyforces)

Figure 9. Nodal forces in a bilinear quadrilateral element with a constant unit thickness due to (a) constant body forces vectors, (b) constant traction vector on one side.

Quadratic Quadrilateral Elements

By introducing four additional nodes in the mid-sides of a quadrilateral element, the displacement shape within the element can have a quadratic form, and the parasitic shear stiffness of the bilinear quadrilateral can be avoided. The quadratic quadrilateral offers such advantage with the price of higher computational time due to having additional degrees of freedom. Consider the plane quadrilateral shown in Figure 10. The trial displacement function has the following form:

$u_{2D} = \begin{pmatrix} u \\ v \\ \end{pmatrix} = \begin{pmatrix} a_0 + a_1X_1 + a_2X_2 + a_3X_1X_2 + a_4X_1^2 + a_5X_2^2 + a_6X_1X_2^2 + a_7X_1^2X_2 \\ b_0 + b_1X_1 + b_2X_2 + b_3X_1X_2 +b_4X_1^2 + b_5X_2^2 + b_6X_1X_2^2 +b_7X_1^2X_2 \\ \end{pmatrix}$

where $a_0$ , $a_1$ , $a_2$ , $a_3$ , $a_4$ , $a_5$ , $a_6$ , $a_7$ , $b_0$ , $b_1$ , $b_2$ , $b_3$ , $b_4$ , $b_5$ , $b_6$ and $b_7$ are 16 generalized degrees of freedom of the element. The displacement function in terms of the nodal degrees of freedom has the form:

$u_{2D} = \begin{pmatrix} u \\ v \\ \end{pmatrix} = \begin{pmatrix} N_1u_1 + N_2u_2 + N_3u_3 + N_4u_4 + N_5u_5 + N_6u_6 + N_7u_7 + N_8u_8 \\ N_1v_1 + N_2v_2 + N_3v_3 + N_4v_4 + N_5v_5 + N_6v_6 + N_7v_7 + N_8v_8 \\ \end{pmatrix}$

Figure 10. Quadratic quadrilateral element

Following the procedures shown in the previous sections, the shape functions can be shown to have the following forms:

$N_1 = -\frac{(b-X_2)(a-X_1)(1 + \frac{X_1}{a} + \frac{X_2}{b})}{4ab}$

$N_2 = -\frac{(b-X_2)(a+X_1)(1 - \frac{X_1}{a} + \frac{X_2}{b})}{4ab}$

$N_3 = -\frac{(b+X_2)(a+X_1)(1 - \frac{X_1}{a} - \frac{X_2}{b})}{4ab}$

$N_4 = -\frac{(b+X_2)(a-X_1)(1 + \frac{X_1}{a} - \frac{X_2}{b})}{4ab}$

$N_5 = \frac{(b-X_2)(a-X_1)(a + X_1)}{2a^2b}$

$N_6 = \frac{(a + X_1)(b - X_2)(b + X_2)}{2ab^2}$

$N_7 = \frac{(b+X_2)(a-X_1)(a + X_1)}{2a^2b}$

$N_8 = \frac{(a - X_1)(b - X_2)(b + X_2)}{2ab^2}$

The distribution of the shape functions on the element are illustrated in Figure 11.

Figure 11. Shape functions distribution on the quadratic quadrilateral

Stiffness Matrix

The stiffness matrix of the linear strain triangle can be evaluated using the following integral (assuming a constant thickness $t$ :

$K^e = t \int_{-b}^b \int_{-a}^a B^T CB dX_1 dX_2$

The stiffness matrix has the dimensions of $16 \times 16$ , and the following Mathematica code can be utilized to view its components:

View Mathematica Code

Shapefun=Table[0,{i,1,8}];
Shapefun[[1]]=(b-x2)(a-x1)/4/a/b*-(1+x1/a+x2/b);
Shapefun[[2]]=(b-x2)(a+x1)/4/a/b*-(1-x1/a+x2/b);
Shapefun[[3]]=(b+x2)(a+x1)/4/a/b*-(1-x1/a-x2/b);
Shapefun[[4]]=(b+x2)(a-x1)/4/a/b*-(1+x1/a-x2/b);
Shapefun[[5]]=(b-x2)(a-x1)(a+x1)/2/a^2/b;
Shapefun[[6]]=(a+x1)(b-x2)(b+x2)/2/a/b^2;
Shapefun[[7]]=(b+x2)(a-x1)(a+x1)/2/a^2/b;
Shapefun[[8]]=(a-x1)(b-x2)(b+x2)/2/a/b^2;
B=Table[0,{i,1,3},{j,1,16}];
Do[B[[1,2i-1]]=B[[3,2i]]=D[Shapefun[[i]],x1];B[[2,2i]]=B[[3,2i-1]]=D[Shapefun[[i]],x2],{i,1,8}];
B//MatrixForm
Cc=Ee/(1+nu)*{{(1-nu)/(1-2nu),nu/(1-2nu),0},{(nu)/(1-2nu),(1-nu)/(1-2nu),0},{0,0,1/2}};
Kbeforeintegration=t*FullSimplify[Transpose[B].Cc.B];
K=Integrate[Kbeforeintegration,{x2,-b,b},{x1,-a,a}];
K//MatrixForm

View Python Code

import sympy as sp
from sympy import *
a, b, x1, x2, E, t, nu = symbols("a b x_1 x_2 E t nu")
Shapefun = Matrix([(b-x2)*(a-x1)/4/a/b*-(1+x1/a+x2/b),
                   (b-x2)*(a+x1)/4/a/b*-(1-x1/a+x2/b),
                   (b+x2)*(a+x1)/4/a/b*-(1-x1/a-x2/b),
                   (b+x2)*(a-x1)/4/a/b*-(1+x1/a-x2/b),
                   (b-x2)*(a-x1)*(a+x1)/2/a**2/b,
                   (a+x1)*(b-x2)*(b+x2)/2/a/b**2,
                   (b+x2)*(a-x1)*(a+x1)/2/a**2/b,
                   (a-x1)*(b-x2)*(b+x2)/2/a/b**2])
B = zeros(3,16)
for i in range(8):
    B[0,2*i] = B[2,2*i+1] = Shapefun[i].diff(x1)
    B[1,2*i+1] = B[2,2*i] = Shapefun[i].diff(x2)
display("B:", B)
Cc = E/(1 + nu)*Matrix([[(1-nu)/(1-2*nu),nu/(1-2*nu),0], 
                         [nu/(1-2*nu),(1-nu)/(1-2*nu), 0], 
                         [0, 0, 1/2]])  
KBeforeintegration = t * B.transpose()*Cc*B
K = Matrix([[simplify(integrate(KBeforeintegration[i,j],(x1,-a,a),(x2,-b,b)))for i in range (16)] for j in range (16)])
display("K: ", K)

Nodal Forces

Assuming that the distributed body forces vector per unit mass is $b={b_1,b_2}$ , the mass density is $\rho$ , and the traction vector per unit area on the left side is $t_n={t_1,t_2}$ , then, the lumped nodal forces due to the distributed body forces can be obtained as follows:

$f^e = \int_e N^T \rho b dX = \frac{\rho t A}{3} \begin{pmatrix} -\frac{b_1}{4} \\ -\frac{b_2}{4} \\ -\frac{b_1}{4} \\ -\frac{b_2}{4} \\ -\frac{b_1}{4} \\ -\frac{b_2}{4} \\ -\frac{b_1}{4} \\ -\frac{b_2}{4} \\ b_1 \\ b_2 \\ b_1 \\ b_2 \\ b_1 \\ b_2 \\ b_1 \\ b_2 \\ \end{pmatrix}$

The lumped nodal forces due to the distributed traction vector has the following form:

$f^e = \int_{\partial e} N^T t_n dS = \int_{-b}^{b} N^T t_n |_{X_1=a} dX_2 = \frac{bt}{3} \begin{pmatrix} t_1 \\ t_2 \\ 0 \\ 0 \\ 0 \\ 0 \\ t_1 \\ t_2 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 4t_1 \\ 4t_2 \\ \end{pmatrix}$

The nodal loads due to a constant distributed body forces vector and the traction vector on the left side are shown in Figure 12. Notice the surprising result of having small forces applied on the corner nodes opposite in direction to the applied body forces.

View Mathematica Code

Shapefun=Table[0,{i,1,8}];
Shapefun[[1]]=(b-x2)(a-x1)/4/a/b*-(1+x1/a+x2/b);
Shapefun[[2]]=(b-x2)(a+x1)/4/a/b*-(1-x1/a+x2/b);
Shapefun[[3]]=(b+x2)(a+x1)/4/a/b*-(1-x1/a-x2/b);
Shapefun[[4]]=(b+x2)(a-x1)/4/a/b*-(1+x1/a-x2/b);
Shapefun[[5]]=(b-x2)(a-x1)(a+x1)/2/a^2/b;
Shapefun[[6]]=(a+x1)(b-x2)(b+x2)/2/a/b^2;
Shapefun[[7]]=(b+x2)(a-x1)(a+x1)/2/a^2/b;
Shapefun[[8]]=(a-x1)(b-x2)(b+x2)/2/a/b^2;
Nn=Table[0,{i,1,2},{j,1,16}];
Do[Nn[[1,2i-1]]=Nn[[2,2i]]=Shapefun[[i]],{i,1,8}];
tn={t1,t2};
rb={rb1,rb2};
fetraction=Integrate[(Transpose[Nn].tn/.x1->-a),{x2,-b,b}]//MatrixForm
febodyforces=Integrate[(Transpose[Nn].rb),{x2,-b,b},{x1,-a,a}]//MatrixForm

View Python Code

import sympy as sp
from sympy import *
a, b, x1, x2, t1, t2, rb1, rb2 = symbols("a b x_1 x_2 t_1 t_2 rb_1 rb_2")
Shapefun = Matrix([(b-x2)*(a-x1)/4/a/b*-(1+x1/a+x2/b),
                   (b-x2)*(a+x1)/4/a/b*-(1-x1/a+x2/b),
                   (b+x2)*(a+x1)/4/a/b*-(1-x1/a-x2/b),
                   (b+x2)*(a-x1)/4/a/b*-(1+x1/a-x2/b),
                   (b-x2)*(a-x1)*(a+x1)/2/a**2/b,
                   (a+x1)*(b-x2)*(b+x2)/2/a/b**2,
                   (b+x2)*(a-x1)*(a+x1)/2/a**2/b,
                   (a-x1)*(b-x2)*(b+x2)/2/a/b**2])
Nn = zeros(2,16)
for i in range(8):
    Nn[0,2*i] = Nn[1,2*i+1] = Shapefun[i]
display("Nn:", Nn)
tn = Matrix([t1,t2])
rb = Matrix([rb1,rb2])
integrand1 = (Nn.transpose()*tn).subs(x1,-a)
integrand2 = (Nn.transpose()*rb)
fetraction = Matrix([simplify(integrate(integrand1[i],(x2,-b,b)))for i in range (16)])
febodyforces = Matrix([simplify(integrate(integrand2[i],(x1,-a,a),(x2,-b,b)))for i in range (16)])
display("fetraction: ", fetraction)
display("febodyforces: ", febodyforces)