5 my project report

For completion of the report, the problem statement is given below taken from the project handout.

The problem described above was solved for the following values of the angle

$\alpha$ (in degrees)

$\fbox{ \{0,15,30,45,50,55\} }$ Where

$\alpha$ is the distortion angle between the ﬁrst and the second elements as shown in the above diagram. The method of ﬁnite elements was implemented in Matlab to solve for the displacements at the nodes.

$4$ Gaussian points were used for the integration step (this is also called the

$2\times 2$ integration rule). After the global stiﬀness matrix was assembled from the two elements stiﬀness matrices, the system equations given by

$K D = F$ were solved using the direct linear system solver.

5.1.2 Geometry, nodes and element description

Two elements were used each having

$8$ nodes.

$4$ of these nodes are at the corners, and the other

$4$ are at the mid point of the element edges. The following diagram shows the global coordinates of the elements nodes. The origin of the global coordinate system is at the lower left corner as shown in the diagram below.

$L$ is the overall length of the beam, which is

$10$ meters, and

$h$ is the height of the beam (which is

$2$ meters).

The vertical and horizontal displacement of each node was solved for using the ﬁnite elements method for each of the diﬀerent values of the angle

$\alpha$ and the vertical displacement at the right bottom edge of the beam was compared to the expected theoretical value in order to see the eﬀect of the element distortion on the accuracy of the ﬁnite element result using the element selected.

The idea is that a good ﬁnite element should produce the same displacement at its nodes regardless of how it was ﬁtted to the physical region in place. This report was to determine if the element selected would still produce good results when deformed. The ﬁrst step was to map each element local node number to a global node number. Local element numbers go from

$1$ to

$8$ since there are only

$8$ nodes per element, but the global node numbers enumerates over all the nodes in all the elements.

Local element node numbering is made in the standard anti clock wise direction by numbering the corner nodes ﬁrst from

$1$ to

$4$ , followed by numbering the middle nodes also in the anti clock wise sense from

$5$ to

$8$ .

The following diagram shows the mapping between the local element node numbers and the global node numbers. The top diagram shows the global node numbers and the lower diagram shows each element local node numbers.

Based on the the above, the table elem_map_nodes was constructed. This table gives the global node number (the entries inside the table) for an element number (the row of the table) and the node number within that speciﬁc element (the column in the table). For example, for element

$1$ with local node number

$1$ the table above shows that the global node number is

$9$ .

There is also a need to lookup the global coordinates

$\{x_i,y_i\}$ given the global node number. This is needed when ﬁnding the Jacobian.

The following table called global_coordinates_tbl was constructed for this purpose. In this table,

$H=2,L=10$ and

$\Delta = \tan \theta \frac{H}{2}$ is the amount of shift in meters of the global node

$3$ and

$11$ . These are the only two nodes which shift location when changing the angle

$\alpha$ . When

$\alpha$ is zero, there will be no distortion of the elements.

There are two degrees of freedom at each node. These are

$u,v$ , representing the horizontal and vertical displacement of a node. Hence there is a need for a lookup table called elem_map_dofs which gives the degree of freedom number of each element’s local node. This table is used for assembling the global stiﬀness matrix.

The following diagram, generated in the Matlab program, gives the degree of freedom number corresponding to each element node (

$u,v$ ). These DOF numbers represent the position of the unknowns in the solution of the

$K D = F$ . This means that there are a total of

$26$ degrees of freedom initially. However, due to boundary conditions constraints, the total number of degrees of freedom reduces to

$22$ .

The above was a general description of the problem and the geometry and data structures used in the Matlab implementation. Next is a discussion of the analytical derivation and the post processing stage which starts after solving for the displacements, followed by discussion of how stress was calculated at the element nodes from the stress value at the four Gaussian points.

5.2 Theory and analytical derivation

5.2.1 Shape functions

Since an 8 node element is used, then there will be 8 shape functions. In ANSYS, the element used is called PLANE183. This element is a serendipity element, which means it has nodes only on the edges and no node in the middle of the element.

The following are the 8 shape functions used

$\begin{align*} f_1&=-\frac{1}{4} \left (1-\eta ^2\right ) (1-\xi )-\frac{1}{4} (1-\eta ) \left (1-\xi ^2\right )+\frac{1}{4} (1-\eta ) (1-\xi )\\ f_2&=-\frac{1}{4} \left (1-\eta ^2\right ) (\xi +1)-\frac{1}{4} (1-\eta ) \left (1-\xi ^2\right )+\frac{1}{4} (1-\eta ) (\xi +1)\\ f_3&=-\frac{1}{4} \left (1-\eta ^2\right ) (\xi +1)-\frac{1}{4} (\eta +1) \left (1-\xi ^2\right )+\frac{1}{4} (\eta +1) (\xi +1)\\ f_4&=-\frac{1}{4} \left (1-\eta ^2\right ) (1-\xi )-\frac{1}{4} (\eta +1) \left (1-\xi ^2\right )+\frac{1}{4} (\eta +1) (1-\xi )\\ f_5&=\frac{1}{2} (1-\eta ) \left (1-\xi ^2\right )\\ f_6&=\frac{1}{2} \left (1-\eta ^2\right ) (\xi +1)\\ f_7&=\frac{1}{2} (\eta +1) \left (1-\xi ^2\right )\\ f_8&=\frac{1}{2} \left (1-\eta ^2\right ) (1-\xi )\\ \end{align*}$

The global coordinates

$x,y$ are now expressed as functions of the natural coordinates

$\xi ,\eta$ using

$\begin{align*} x(\xi ,\eta ) &= \sum _{i=1}^M x_i f_i(\xi ,\eta ) \\ y(\xi ,\eta ) &= \sum _{i=1}^M y_i f_i(\xi ,\eta ) \end{align*}$

Where

$M$ is the number of nodes of each element (which is

$8$ ) and

$x_i,y_i$ are the global coordinates of these nodes. Expanding the above gives the result below.

The result below is shown for some element number

$k$ . In this expansion,

$x_i^k$ means the global

$x$ coordinate of the

$i^{\text{th}}$ node in the

$k^{\text{th}}$ element.

Similarly,

$y_i^k$ means the global

$y$ coordinate of the

$i^{\text{th}}$ node in the

$k^{\text{th}}$ element. These are read in the Matlab code using the global_coordinates_tbl table, which gives the global

$x,y$ coordinates of each global node and by using elem_map_nodes table to map the element node number to the global node number.

$\begin{align*} x(\xi ,\eta ) &=x_1^k\left (-\frac{1}{4} \left (1-\eta ^2\right ) (1-\xi )-\frac{1}{4} (1-\eta ) \left (1-\xi ^2\right )+\frac{1}{4} (1-\eta ) (1-\xi )\right )+x_2^k\left (-\frac{1}{4} \left (1-\eta ^2\right ) (\xi +1)-\frac{1}{4} (1-\eta ) \left (1-\xi ^2\right )+\frac{1}{4} (1-\eta ) (\xi +1)\right )+x_3^k\left (-\frac{1}{4} \left (1-\eta ^2\right ) (\xi +1)-\frac{1}{4} (\eta +1) \left (1-\xi ^2\right )+\frac{1}{4} (\eta +1) (\xi +1)\right )+x_4^k\left (-\frac{1}{4} \left (1-\eta ^2\right ) (1-\xi )-\frac{1}{4} (\eta +1) \left (1-\xi ^2\right )+\frac{1}{4} (\eta +1) (1-\xi )\right )+x_5^k\left (\frac{1}{2} (1-\eta ) \left (1-\xi ^2\right )\right )+x_6^k\left (\frac{1}{2} \left (1-\eta ^2\right ) (\xi +1)\right )+x_7^k\left (\frac{1}{2} (\eta +1) \left (1-\xi ^2\right )\right )+x_8^k\left (\frac{1}{2} \left (1-\eta ^2\right ) (1-\xi )\right )\\ y(\xi ,\eta ) &=y_1^k\left (-\frac{1}{4} \left (1-\eta ^2\right ) (1-\xi )-\frac{1}{4} (1-\eta ) \left (1-\xi ^2\right )+\frac{1}{4} (1-\eta ) (1-\xi )\right )+y_2^k\left (-\frac{1}{4} \left (1-\eta ^2\right ) (\xi +1)-\frac{1}{4} (1-\eta ) \left (1-\xi ^2\right )+\frac{1}{4} (1-\eta ) (\xi +1)\right )+y_3^k\left (-\frac{1}{4} \left (1-\eta ^2\right ) (\xi +1)-\frac{1}{4} (\eta +1) \left (1-\xi ^2\right )+\frac{1}{4} (\eta +1) (\xi +1)\right )+y_4^k\left (-\frac{1}{4} \left (1-\eta ^2\right ) (1-\xi )-\frac{1}{4} (\eta +1) \left (1-\xi ^2\right )+\frac{1}{4} (\eta +1) (1-\xi )\right )+y_5^k\left (\frac{1}{2} (1-\eta ) \left (1-\xi ^2\right )\right )+y_6^k\left (\frac{1}{2} \left (1-\eta ^2\right ) (\xi +1)\right )+y_7^k\left (\frac{1}{2} (\eta +1) \left (1-\xi ^2\right )\right )+y_8^k\left (\frac{1}{2} \left (1-\eta ^2\right ) (1-\xi )\right )\\ \end{align*}$

5.2.2 Finding the Jacobian

The Jacobian is evaluated at each Gaussian integration point during the integration step. It has the form

$\begin{align*} J &= \begin{pmatrix} \frac{\partial x}{\partial \xi } & \frac{\partial y}{\partial \xi }\\ \frac{\partial x}{\partial \eta } & \frac{\partial y}{\partial \eta } \end{pmatrix} \end{align*}$

$\begin{align*} \frac{\partial x}{\partial \xi } &= x_1 \left (\frac{1}{4} \left (1-\eta ^2\right )+\frac{1}{2} (1-\eta ) \xi +\frac{\eta -1}{4}\right )+x_2 \left (\frac{1}{4} \left (\eta ^2-1\right )+\frac{1}{2} (1-\eta ) \xi +\frac{1-\eta }{4}\right )+x_3 \left (\frac{1}{4} \left (\eta ^2-1\right )+\frac{1}{2} (\eta +1) \xi +\frac{\eta +1}{4}\right )+x_4 \left (\frac{1}{4} \left (1-\eta ^2\right )+\frac{1}{2} (\eta +1) \xi +\frac{1}{4} (-\eta -1)\right )+\frac{1}{2} \left (1-\eta ^2\right ) x_6-\frac{1}{2} \left (1-\eta ^2\right ) x_8-(1-\eta ) \xi x_5-(\eta +1) \xi x_7\\ \frac{\partial y}{\partial \xi } &= y_1 \left (\frac{1}{4} \left (1-\eta ^2\right )+\frac{1}{2} (1-\eta ) \xi +\frac{\eta -1}{4}\right )+y_2 \left (\frac{1}{4} \left (\eta ^2-1\right )+\frac{1}{2} (1-\eta ) \xi +\frac{1-\eta }{4}\right )+y_3 \left (\frac{1}{4} \left (\eta ^2-1\right )+\frac{1}{2} (\eta +1) \xi +\frac{\eta +1}{4}\right )+y_4 \left (\frac{1}{4} \left (1-\eta ^2\right )+\frac{1}{2} (\eta +1) \xi +\frac{1}{4} (-\eta -1)\right )+\frac{1}{2} \left (1-\eta ^2\right ) y_6-\frac{1}{2} \left (1-\eta ^2\right ) y_8-(1-\eta ) \xi y_5-(\eta +1) \xi y_7\\ \frac{\partial x}{\partial \eta } &= x_1 \left (\frac{1}{2} \eta (1-\xi )+\frac{1}{4} \left (1-\xi ^2\right )+\frac{\xi -1}{4}\right )+x_2 \left (\frac{1}{2} \eta (\xi +1)+\frac{1}{4} \left (1-\xi ^2\right )+\frac{1}{4} (-\xi -1)\right )+x_3 \left (\frac{1}{2} \eta (\xi +1)+\frac{1}{4} \left (\xi ^2-1\right )+\frac{\xi +1}{4}\right )+x_4 \left (\frac{1}{2} \eta (1-\xi )+\frac{1}{4} \left (\xi ^2-1\right )+\frac{1-\xi }{4}\right )-\eta (\xi +1) x_6-\eta (1-\xi ) x_8-\frac{1}{2} \left (1-\xi ^2\right ) x_5+\frac{1}{2} \left (1-\xi ^2\right ) x_7\\ \frac{\partial y}{\partial \eta } &= y_1 \left (\frac{1}{2} \eta (1-\xi )+\frac{1}{4} \left (1-\xi ^2\right )+\frac{\xi -1}{4}\right )+y_2 \left (\frac{1}{2} \eta (\xi +1)+\frac{1}{4} \left (1-\xi ^2\right )+\frac{1}{4} (-\xi -1)\right )+y_3 \left (\frac{1}{2} \eta (\xi +1)+\frac{1}{4} \left (\xi ^2-1\right )+\frac{\xi +1}{4}\right )+y_4 \left (\frac{1}{2} \eta (1-\xi )+\frac{1}{4} \left (\xi ^2-1\right )+\frac{1-\xi }{4}\right )-\eta (\xi +1) y_6-\eta (1-\xi ) y_8-\frac{1}{2} \left (1-\xi ^2\right ) y_5+\frac{1}{2} \left (1-\xi ^2\right ) y_7\\ \end{align*}$

The above is now used to ﬁnd the Jacobian and its determinant and also ﬁnd

$\Gamma$ and the matrix

$B_2$ . To ﬁnd the matrix

$B_3$ the derivatives of each shape function is taken w.r.t.

$\xi$ and

$\eta$ .

This below gives the result of this computation

$\begin{align*} \frac{\partial f_1}{\partial \xi } &= \frac{1}{4} \left (1-\eta ^2\right )+\frac{1}{2} (1-\eta ) \xi +\frac{\eta -1}{4}\\ \frac{\partial f_1}{\partial \eta } &= \frac{1}{2} \eta (1-\xi )+\frac{1}{4} \left (1-\xi ^2\right )+\frac{\xi -1}{4}\\ \frac{\partial f_2}{\partial \xi } &= \frac{1}{4} \left (\eta ^2-1\right )+\frac{1}{2} (1-\eta ) \xi +\frac{1-\eta }{4}\\ \frac{\partial f_2}{\partial \eta } &= \frac{1}{2} \eta (\xi +1)+\frac{1}{4} \left (1-\xi ^2\right )+\frac{1}{4} (-\xi -1)\\ \frac{\partial f_3}{\partial \xi } &= \frac{1}{4} \left (\eta ^2-1\right )+\frac{1}{2} (\eta +1) \xi +\frac{\eta +1}{4}\\ \frac{\partial f_3}{\partial \eta } &= \frac{1}{2} \eta (\xi +1)+\frac{1}{4} \left (\xi ^2-1\right )+\frac{\xi +1}{4}\\ \frac{\partial f_4}{\partial \xi } &= \frac{1}{4} \left (1-\eta ^2\right )+\frac{1}{2} (\eta +1) \xi +\frac{1}{4} (-\eta -1)\\ \frac{\partial f_4}{\partial \eta } &= \frac{1}{2} \eta (1-\xi )+\frac{1}{4} \left (\xi ^2-1\right )+\frac{1-\xi }{4}\\ \frac{\partial f_5}{\partial \xi } &= -(1-\eta ) \xi \\ \frac{\partial f_5}{\partial \eta } &= \frac{1}{2} \left (\xi ^2-1\right )\\ \frac{\partial f_6}{\partial \xi } &= \frac{1}{2} \left (1-\eta ^2\right )\\ \frac{\partial f_6}{\partial \eta } &= -\eta (\xi +1)\\ \frac{\partial f_7}{\partial \xi } &= -(\eta +1) \xi \\ \frac{\partial f_7}{\partial \eta } &= \frac{1}{2} \left (1-\xi ^2\right )\\ \frac{\partial f_8}{\partial \xi } &= \frac{1}{2} \left (\eta ^2-1\right )\\ \frac{\partial f_8}{\partial \eta } &= -\eta (1-\xi ) \end{align*}$

With the above, the matrix

$B_3$ matrix was calculated giving

$\begin{align*} B &= B_1 B_2 B3\\ \end{align*}$

And calculate the element stiﬀness matrix given by

$\begin{align*} k_{\text{elem}}&=\int B^{T}EB\, dV\\ &=\int _{-1}^{+1}\int _{-1}^{+1} B^{T}EB |J| \, dV \end{align*}$

The elements stiﬀness matrices

$k_{\text{elem}}$ are then combined to make the global stiﬀness matrix

$K$ and then the system

$K D = F$ was solved for the unknowns

$D$ which are the nodal displacements in the

$x$ and

$y$ direction.

In the above

$F$ is the load vector, which in this problem contains only one non-zero entry, which is the vertical load of

$-20 N$ at the middle of the right edge of the beam.

In the above, the matrix

$B_1$ is

$\begin{align*} B_1 &= \begin{bmatrix} 1&0&0&0\\ 0&0&0&1\\ 0&1&1&0 \end{bmatrix} \end{align*}$

And the matrix

$B_2$ is

$\begin{align*} B_2 &= \begin{bmatrix} \Gamma _{11} & \Gamma _{12}& 0&0\\ \Gamma _{21} & \Gamma _{22}& 0&0\\ 0&0& \Gamma _{11}&\Gamma _{12}\\ 0&0& \Gamma _{21}&\Gamma _{22} \end{bmatrix} \end{align*}$

Where

$\Gamma$ is the inverse of the Jacobian matrix

$\Gamma = J^{-1}$ . And the matrix

$B_3$ is

$\begin{align*} B_3 &= \begin{bmatrix} \frac{\partial f_1}{\partial \xi }&0 &\frac{\partial f_2}{\partial \xi }&0& \dots &\frac{\partial f_8}{\partial \xi }&0\\ \frac{\partial f_1}{\partial \eta }&0 &\frac{\partial f_2}{\partial \eta }&0& \dots &\frac{\partial f_8}{\partial \eta }&0\\ 0 &\frac{\partial f_1}{\partial \xi }&0 &\frac{\partial f_2}{\partial \xi }&\dots &0&\frac{\partial f_8}{\partial \xi }\\ 0 &\frac{\partial f_1}{\partial \eta }&0 &\frac{\partial f_2}{\partial \eta }&\dots &0&\frac{\partial f_8}{\partial \eta } \end{bmatrix} \end{align*}$

The following diagram shows the internal structure of the global stiﬀness matrix, found using the spy() command in Matlab. It shows the bands along the diagonals and illustrates how sparse the matrix is.

The above concludes the solve stage. The next stage is the post processing, where stress calculations are performed.

5.2.3 Stress recovery

This is a discussion of how the stress at the elements nodes was found. Initially the direct method was used to ﬁnd the stress at any point in the element.

This method was found to be accurate as long as there was no distortion. Once the angle was increased, this method did not produce stress results which agreed with ANSYS. Therefore, this method was not used, and instead a new implementation was made based on the extrapolation method.

This method is described in reference [1], pages 230-232. This method is more complicated that the direct method, but it is much more accurate. It uses two coordinates systems. The original natural coordinates system of the element

$(\xi ,\eta )$ , which extends from

$-1\dots 1$ across the length and height of the element, and a new coordinates systems called

$(r,s)$ which extends across what is called the Gaussian element.

Therefore, when

$\xi =\frac{1}{\sqrt 3}$ the value of

$r$ is one. And when

$\eta =\frac{1}{\sqrt 3}$ then

$s=1$ also.

Therefore, the relation between

$(r,s)$ coordinates system and

$(\xi ,\eta )$ coordinates system is as follows

$\begin{align*} r &= \xi \sqrt{3} \\ s &= \eta \sqrt{3} \end{align*}$

The following diagram shows the mapping at two diﬀerent points for illustration

Now that the mapping between

$\xi ,\eta$ and

$(r,s)$ is determined, the next step was to ﬁnd

$(r,s)$ at each point where the stress needs to be found at by extrapolating the stress value from the

$4$ Gaussian points to that point of interest. The reason for doing all of the above, is because

$(r,s)$ are used when evaluating the

$N_i$ shape functions described below, and not the original

$(\xi ,\eta )$ values.

Therefore, for each point where the stress needs to be found, say point

$p$ , its coordinates in the

$(r,s)$ system are found ﬁrst, and then the following extrapolation formula is applied

Where

$\sigma _i$ is the stress at the Gaussian point (which was found using the direct method based on the full

$8$ shape functions of the main element).

In the above,

$N_i$ are the shape functions used for extrapolation. These are not the same shape functions used in the original element. These shape functions are based on

$4$ node Gaussian element and are given by

$\begin{align*} N_1 &= \frac{1}{4}(1-r)(1-s) \\ N_2 &= \frac{1}{4}(1+r)(1-s) \\ N_3 &= \frac{1}{4}(1+r)(1+s) \\ N_4 &= \frac{1}{4}(1-r)(1+s) \\ \end{align*}$

For example, to ﬁnd the stress at point

$(\xi ,\eta )=(0,-1)$ , the ﬁrst step is to determine this point’s

$(r_p,s_p)$ coordinates. Since

$r_p=\sqrt{3}\xi$ then

$r_p=0$ and since

$s_p=\sqrt{3}\eta$ then

$s_p=-\sqrt{3}=$ . Applying the extrapolation formula above gives

$\begin{align*} \sigma _p &= \left ( \frac{1}{4}(1-r_p)(1-s_p) \right )\sigma _1 + \left ( \frac{1}{4}(1+r_p)(1-s_p) \right )\sigma _2 + \left ( \frac{1}{4}(1+r_p)(1+s_p) \right )\sigma _3 + \left ( \frac{1}{4}(1-r_p)(1+s_p) \right )\sigma _4 \\ &= \left ( \frac{1}{4}(1+\sqrt{3}) \right )\sigma _1 + \left ( \frac{1}{4}(1+\sqrt{3}) \right )\sigma _2 + \left ( \frac{1}{4}(1-\sqrt{3}) \right )\sigma _3 + \left ( \frac{1}{4}(1-\sqrt{3}) \right )\sigma _4 \\ &= 0.6830127 \sigma _1 + 0.6830127 \sigma _2 -0.1830127 \sigma _3 -0.1830127 \sigma _4 \end{align*}$

Since the stresses at four Gaussian points

$\sigma _i$ are known, the stress at the point

$p$ is now found from the above. The above method was found to produce more accurate result for

$\sigma _p$ than using the direct method to ﬁnd

$\sigma _p$ and the result found for the stress at the nodes agreed with those found by ANSYS.

The following diagram shows the

$(r,s)$ coordinates of all the element points used to calculate the stresses at using this method. A total of

$13$ points was used per element. These are the

$8$ nodes of the element, and also the center of the element and the Gaussian points themselves giving a total of

$13$ points. These are used to generate the stress contour. This was done for both elements. The generated stress contour agreed with ANSYS results.

5.3 Results

The results are listed by the angle

$\theta$ used to distort the elements with. For each angle, the deﬂection and

$\sigma _x$ stress contour and tables are generated using Matlab and also using ANYS to compare with side by side.

The following angles are used (in degrees)

$\{0,15,30,45,50,55\}$ . When trying to use

$60$ degrees distortion, ANSYS complained and gave number of computational error messages relating to the element shape. It is not clear why ANSYS did not accept such large angle, since the Matlab implementation worked. But since ANSYS did not produce result for this case, the angle

$55$ degrees was the maximum distortion used for both Matlab and ANSYS.

For each angle, a short summary of the result in the form of a table is ﬁrst given that compares Matlab and ANSYS result. This short summary contains only the deﬂection at the bottom right corner, which is node

$13$ . After this, the full result and stress plots are given.

5.3.1 No element distortion. zero angle

summary of result

Matlab result

The following table shows Matlab result for the direct stress

$\sigma _x$ found at each node for each element.

ANSYS result

The following table shows ANSYS result for the direct stress

$\sigma _x$ found at each node for each element.

5.3.2 15 degrees distortion

summary of result

Matlab result

The following table shows Matlab result for the direct stress

$\sigma _x$ found at each node for each element.

ANSYS result

The following table shows ANSYS result for the direct stress

$\sigma _x$ found at each node for each element.

5.3.3 30 degrees distortion

summary of result

Matlab result

The following table shows Matlab result for the direct stress

$\sigma _x$ found at each node for each element.

ANSYS result

The following table shows ANSYS result for the direct stress

$\sigma _x$ found at each node for each element.

5.3.4 45 degrees distortion

summary of result

Matlab result

The following table shows Matlab result for the direct stress

$\sigma _x$ found at each node for each element.

ANSYS result

The following table shows ANSYS result for the direct stress

$\sigma _x$ found at each node for each element.

5.3.5 50 degrees distortion

summary of result

Matlab result

The following table shows Matlab result for the direct stress

$\sigma _x$ found at each node for each element.

ANSYS result

The following table shows ANSYS result for the direct stress

$\sigma _x$ found at each node for each element.

5.3.6 55 degrees distortion

summary of result

Matlab result

The following table shows Matlab result for the direct stress

$\sigma _x$ found at each node for each element.

ANSYS result

The following table shows ANSYS result for the direct stress

$\sigma _x$ found at each node for each element.

5.4 Observations, discussion and conclusions

It is clear from the above result that the

$8$ node element used did not behave well at all when it became distorted.

The element used is a serendipity element¹ . These elements are known not to give good results when they distorted² . There are a number of distortions that an element could have, such as an edge distortion or angular distortion and others. In this report we only looked at angular distortion.

As the angular distortion increased, the result became less accurate. The inaccuracy also accelerated when the angle was above

$35^0$ based on the diagram generated below. At angle

$55^0$ , the vertical deﬂection reported by the ﬁnite element Matlab program (and ANSYS as well) was

$-0.86157$ meters where the theoretical result should be close to

$-1.03$ meters. This is over

$16\%$ error. A signiﬁcation error in accuracy. The following graph shows how the error in the vertical deﬂection changed as a function of the distortion angle.

From the above one can see that to keep the error below

$5\%$ , the angular distortion should not be more than about

$35^0$ as the element behaves very badly after that. ANSYS will not even accept the analysis when trying angle of

$60^0$ and gave number of errors relating to element shape. This means this element is not suitable for modeling physical regions which are highly irregular. This element is suitable for ﬁtting in physical region which has straight edges and mostly rectangular regions. Curved regions and other such regions would need to be modeled using other elements.

When it comes to post processing and stress recovery, which the term used for calculating stresses in the post processing stage, there are three methods that can be used. These are

When making the contour plots for

$\sigma _x$ , one can see that the stress along the same line between the two elements is no longer smooth as the angle increases. ANSYS shows this to be smooth transition in the stress contour, but this must be because ANSYS did averaging across element boundaries.

In the Matlab implementation, No stress averaging was made between nodes across elements, hence the distortion (contour lines) is more clear in the stress contour along the line between the two elements as the following plot shows when the angle is

$55^0$ . This shows clearly that stress across elements is not smooth and changes abruptly now.

The more angular distortion there is, the sharper this distortion in the stress contour between the two element became. This is due to the element becoming less accurate as it distorts. Comparing the above plot to the one when the angle was zero (no distortion) one can see that in the no distortion case the stress across the elements is smooth and has same values at the nodes connecting the two elements.

In conclusion, it is recommended that this element be used only when there is no distortion in the geometry.

5.4.1 References

5.5 Appendix

5.5.1 APDL used for ANSYS

The following is the APDL script used for ANSYS analysis. ANSYS version 16.02, student version was used.

5.5.2 Matlab source code

The following is the listing of the m ﬁle for Matlab implementation. For making the contour plot, an external m ﬁle was used from Mathworks ﬁle exchange called tricontf and is included in the zip ﬁle. This function is not listed here.

function nma_project_EMA_471()
%Final project Finite Elements option, EMA option.
%By Nasser M. Abbasi
%EMA 471, spring 2016, Univ. Of Wisconsin, Madison
%This files uses a Mathworks file exchange function
%which is included in the zip file and is in the same folder
%as this m file. Needed for making the contour plot.
 
close all; clc;
 
data      = PRE_PROCESSOR();    %set up data structure
solution  = SOLVE(data);        %assemble and solve KD=F 
 
%stress calculations and print results
POST_PROCESSOR(solution,data);  
 
end
%=============================================
function data = PRE_PROCESSOR()
%In this function we allocate the mapping tables and
%set all the problem parameters. All is saved in a struct data
 
%get folder name we are running from. Needed to save results
if(~isdeployed)
    data.baseFolder = fileparts(which(mfilename));
    cd(data.baseFolder);
end
 
%we are using 2 by 2 Gaussian rule for integration.
wt(1) = 1;           wt(2) = 1;
gs(1) = -0.57735027; gs(2)  = 0.57735027;
 
%allocate data parameters
data.wt = wt;
data.gs = gs;
data.num_elements = 2;
data.angle_degree = 0;  %change as needed
data.angle = data.angle_degree*pi/180;
 
%  global node numbering:
%
%  1     2      3       4         5
%  o-----o------o-------o---------o
%  |            |                 |
%  o 6          o7                o8
%  |            |                 |
%  o-----o------o-------o---------o
%  9     10     11     12         13
%
%
data.elem_map_nodes = [9,11,3,1,10,7,2,6;   %element (1)
                       11,13,5,3,12,8,4,7];  %element (2)
 
L          = 10; 
H          = 2;
data.L     = L;  %meter
data.H     = H;   %meter
data.shift = (H/2)*tan(data.angle); %meter
 
data.global_coord_tbl = [...  %global coordinates of the 13 nodes
    0,          H;         %node 1, top left corner
    L/4,        H;         %node 2, etc...
    L/2+data.shift,  H;    %if angle is 90 degrees, then shift=0
    (3/4)*L,    H;
    L,          H;         %last node in top of beam, node 5
    0,          H/2;       %node 6, middle of left edge of beam.
    L/2,        H/2;
    L,          H/2;
    0,          0;         %node 9, at origin, bottom left corner
    L/4,        0;
    L/2-data.shift,  0;
    (3/4)*L,    0;
    L,          0];
 
data.young_module  = 10^4;
data.mu            = 0.3;
data.E0            = data.young_module/(1-data.mu^2)*...
                     [1,data.mu,0;data.mu,1,0;0,0,(1-data.mu)/2];
 
data.elem_map_dofs = zeros(data.num_elements ,16);
display_diagram_of_dof(data.L,data.H);
 
%elem_map_dofs is used to merge the element k to the global K
for i=1:data.num_elements
    for j=1:8
        data.elem_map_dofs(i,2*j-1)= 2*data.elem_map_nodes(i,j)-1;
        data.elem_map_dofs(i,2*j)  = 2*data.elem_map_nodes(i,j);
    end
end
 
end
%====================================================
function solution = SOLVE(data)
%all the problem data description is now in struct data. This
%function assembles the stiffness matrix and solves KD=F
%and returns the solution
 
N        = size(data.global_coord_tbl,1);
k_global = zeros(2*N,2*N);
 
%do initial verification that shape functions adds to one.
 
%throws error if not verified
verify_shape_functions_sum_to_one(data.gs);  
%fprintf('verified shape functions ok....\n');
 
for k = 1:data.num_elements
    %obtain global coordinates for the node of this element
    [x_coord,y_coord] = find_XY_coordinates(k,...
                      data.elem_map_nodes,data.global_coord_tbl);
    k_elem = zeros(16,16); %allocate local element k
     
    for i = 1:2  %we are using 2 by 2 Gaussian integration rule
        xi = data.gs(i);
        for j = 1:2
            eta   = data.gs(j);
            J     = get_J(xi,eta,x_coord,y_coord);
            detJ  = det(J);
            if detJ<0
                error('Internal error. negative |J| detected.');
            end
            B = get_B(xi,eta,J);  %find the strain rate matrix
            v = B'*data.E0*B;
             
            for ii = 1:16  %numerical integration
                for jj = 1:16
                    if(ii<=jj) %only do upper diagonal, symmetry
                        k_elem(ii,jj) = k_elem(ii,jj)+...
                            data.wt(i)*data.wt(j)*v(ii,jj)*detJ;
                    end
                end
            end
             
        end
    end
     
    %First copy the upper part of k to the lower part,  symmetric
    k_elem  = k_elem + triu(k_elem,1)';
     
    %Merge it to global stiffness matrix
    k_global = merge_element_to_global(k, k_elem, ...
                                   k_global,data.elem_map_dofs);
end
 
%we need now to zero out rows/cols {1,11,12,17} from the 
%above, since these correspond to the fixed boundary conditions. 
%Make sure to keep the diagonal element. Since the boundary 
%conditions are zero at these, we do not have to patch the F 
%vector on the RHS as normally we would. put a 1 in the diagonal
 
fix   = [1,11,12,17];  %these are the fixed DOF to zero out
 
%----- COMMENTED OUT ------------------------
%below is one method to fix. It gives same as the next method. 
%But the next method below this is simpler since it keeps the 
%same sizes of data kept for reference
 
%k_global(fix,:)=[];
%k_global(:,fix)=[];
%r_global       = zeros(22,1);
%r_global(end)  = -20;  % Newton
%d              = k_global\r_global;  %solve for displacement
 
%----- END COMMENTED --------------------
 
%second method to fix K
 
for i = 1:length(fix)
    tmp                     =  k_global(fix(i),fix(i));
    k_global(fix(i),:)      = 0;
    k_global(:,fix(i))      = 0;
    k_global(fix(i),fix(i)) = tmp; %same effect as setting to 1.
end
 
r_global      = zeros(26,1);
r_global(16)  = -20;  % Newton
solution      = k_global\r_global;  %solve for displacement
 
%finished. Now write the solution to file to include it in report
write_nodal_solution_to_file(solution,data.baseFolder,...
                                              data.angle_degree);
 
figure(); %show the global stiffness matrix structure using spy
spy(k_global);
title('global stiffness matrix after fixing for B.C.',...
    'interpreter','Latex','Fontsize',11);
%print(gcf, '-dpdf', '-r600','../images/spy.pdf');
 
end
%=====================================================
function POST_PROCESSOR(solution,data)
%This is final stage. We find stress and make contour plots
%and draw the deflection shape of the beam
 
draw_deflection(solution,data);
generate_stress_diagram(solution,data);
end
%========================================
function write_nodal_solution_to_file(d,baseFolder,angle)
%d is the nodal solution vector. 26 by 1.
%write the X,Y solution to text file to use in report
 
d              = reshape(d,2,13)';
fileName       = [baseFolder sprintf(...
                       '/data/deformation_matlab_%d.tex',angle)];
fileName       = strrep(fileName, '/', filesep);
[fileID0,errMsg] = fopen(fileName,'w');
 
if fileID0<0
    fprintf('Error opening %s\n, the message is [%s\n]',...
                                                fileName,errMsg);
    error(errMsg);
end
 
fprintf(fileID0,'\\begin{table}[!htbp]\n');
fprintf(fileID0,'\\centering\n');
fprintf(fileID0,'\\captionsetup{width=.8\\textwidth}\n');
fprintf(fileID0,'\\begin{tabular}{|l|l|l|}\\hline\n');
fprintf(fileID0,'global node \\#& $x$ (meter)& $y$ (meter)\\\\\\hline\n');
for i=1:size(d,1)
    fprintf(fileID0,'$%d$ & $%7.6f$& $%7.6f$\\\\ \n',i,d(i,1),d(i,2));
end
fprintf(fileID0,'\\hline\n\\end{tabular}\n');
fprintf(fileID0,...
  '\\caption{Matlab result. nodal solutions, angle [$%d$] degree}\n',angle);
fprintf(fileID0,'\\end{table}\n');
fprintf(fileID0,'\\FloatBarrier\n');
 
fclose(fileID0);
 
end
%==========================================
function generate_stress_diagram(solution,data)
 
%find the stress at node of each element
element_stress_1 = stress_calculation(solution,data,1);
element_stress_2 = stress_calculation(solution,data,2);
 
%make one diagram of the overall stress contour across the beam
stress_diagram(data,element_stress_1,element_stress_2);
end
%=============================================
function element_stress = stress_calculation(solution,data,...
                                                element_number)
 
%first calculate stress at the 4 Gaussian points for
%element in order to use for extrapolation. These are ordered
%anticlock wise.
 
%used to store stress at the 4 Gaussian points
gauss_stress = zeros(4,1);  
 
%these are the r,s coordinates of the Gaussian element inside
%the element itself used for extrapolation. See report for
%more details. These are ordered anticlock wise, with one in
%the center. So there are 9 of them. Also, the Gaussian stress
%is added as well. So we end up with 9+4=13 total stress points
%for each element. This should be enough to make nice contuor with
 
z  = sqrt(3);
r  = [-z,z,z,-z,0,z,0,-z];
s  = [-z,-z,z,z,-z,0,z,0];
 
%these are the natural coordinates in xi,eta space used 
%to calculate the stress at Gaussian points, using the full 8 
%shape functions
z   =  1/sqrt(3);
xi  = [-z,z,z,-z];
eta = [-z,-z,z,z];
 
L  = data.L;
H  = data.H;
 
%this is used to find global coordinates of all stress points, for 
%contour plot this is center of element in global space
if element_number == 1
    X0=L/4;
    Y0=H/2;
else
    X0=(3/4)*L;
    Y0=H/2;
end
 
%find element nodal coordinates in global space
%this only has nodes. We will add the center and the gaussian 
%points also later to make contour plot
[x_coord,y_coord] = find_XY_coordinates(element_number,...
                     data.elem_map_nodes,data.global_coord_tbl);
 
%Now make matrix to store all stress result in for element 1. 
%We are finding stress at 13 points. 8 for nodes, one for center, 
%and the 4 Gaussian points. we need 4 columns. First two are 
%the x,y in global space of the point, and the stress. The 
%first column is just the point ID, for tracking. Not used for 
%plotting.
 
element_stress = zeros(13,4);
 
global_node_numbers = data.elem_map_nodes(element_number,:);
U   = solution(data.elem_map_dofs(element_number,:));
for i = 1:length(xi)
    J    = get_J( xi(i), eta(i), x_coord, y_coord); %jacobian
    B    = get_B( xi(i), eta(i), J); %strain rate matrix
    %find actual strain from displacements at nodes
    strain  = B  * U;  
    stress  = data.E0 * strain;
    gauss_stress(i) = stress(1); %use direct stress only
end
 
%Now do the extrapolation. See report
for i=1:length(r)
    element_stress(i,4)=0;
    element_stress(i,1)=global_node_numbers(i);
    element_stress(i,2)=x_coord(i);
    element_stress(i,3)=y_coord(i);
    for j=1:4 %extrapolation
        switch j
            case 1,
                f = (1/4)*(1-r(i))*(1-s(i));
            case 2,
                f = (1/4)*(1+r(i))*(1-s(i));
            case 3,
                f = (1/4)*(1+r(i))*(1+s(i));
            case 4,
                f = (1/4)*(1-r(i))*(1+s(i));
        end
        element_stress(i,4) = element_stress(i,4)+ f * ...
                                             gauss_stress(j);
    end
end
 
%now add the center and the 4 gaussian points we found before. 
%These come after the nodes. This is in order to improve the 
%contour plot by having more points.
element_stress(9,4) = 0;
element_stress(9,1) = -1;  %we do not have a global node number 
                           %for this. this is just a place holder
element_stress(9,2)=X0;
element_stress(9,3)=Y0;
for j=1:4 %extrapolation
  element_stress(9,4)= element_stress(9,4)+ 1/4 * gauss_stress(j);
end
 
%now add the acutal Gaussian stress found. This make it up to 
%13 points first find the global coordinates of the element 
%gaussian points
z = 1/sqrt(3);
gauss_global_coordinates=[X0-z*(1/4)*L,Y0-z*(1/2)*H;...
    X0+z*(1/4)*L,Y0-z*(1/2)*H;...
    X0+z*(1/4)*L,Y0+z*(1/2)*H;...
    X0-z*(1/4)*L,Y0+z*(1/2)*H];
 
for i=1:length(gauss_global_coordinates)
    element_stress(9+i,4)=gauss_stress(i);
    element_stress(9+i,1)=-1; %place holder
    element_stress(9+i,2)=gauss_global_coordinates(i,1);
    element_stress(9+i,3)=gauss_global_coordinates(i,2);
end
 
end
%==========================================
function stress_diagram(data,element_stress_1,element_stress_2)
 
%we are done!  Now we can make contour of first element stress
%This uses tricontf, which is a mathworks file exchange file
%since Matlab does not have such a function build in
figure;
 
% I commented out the stress average last minute. I think it is
%better NOT to do stress averging across elements, in order 
%to more clearly see the difference. I left the code here for 
%reference in case need to use it later
 
%----------- COMMENTED OUT ---------------------
%now do stress avergaing on the nodes that are between
%element one and two. These nodes have global node
%numbers of 2,7,11 which correspond to local nodes
% 2,6,3 for first element and nodes 1,8,4 for second element.
 
%element_stress_1(2,4) = (element_stress_1(2,4)+element_stress_2(1,4))/2;
%element_stress_1(6,4) = (element_stress_1(6,4)+element_stress_2(8,4))/2;
%element_stress_1(3,4) = (element_stress_1(3,4)+element_stress_2(4,4))/2;
 
%now that we averaged the stress, remove these entry from the second
%element before merging, since it is duplicate
%element_stress_2 = element_stress_2([2:3,5:7,9:end],:);
 
%----------- END COMMENTED OUT ----------------
 
x = [element_stress_1(:,2);element_stress_2(:,2)];
y = [element_stress_1(:,3);element_stress_2(:,3)];
z = [element_stress_1(:,4);element_stress_2(:,4)];
M = delaunay(x,y);
max_stress = max(z);
min_stress = min(z);
range_of_stress = linspace(min_stress,max_stress,15);
 
%this below uses mathworks file exchange function. 
%It is in the same folder
[~,h]=tricontf(x,y,M,z,range_of_stress,'-k');
%set(h,'edgecolor','none');
axis equal tight;
%hold on;
%[~,h]=tricont(x,y,M,z,range_of_stress,'-k');
colorbar;
 
title(sprintf('$\\sigma_x$ contour for angle $%d^o$',...
                                           data.angle_degree),...
                            'Fontsize',12,'interpreter','Latex');
 
xlabel('length of beam in meters','Fontsize',10,'interpreter','Latex');
ylabel('height in meters','interpreter','Latex','Fontsize',10);
 
%uncomment to write the plot
%print(gcf, '-dpdf', '-r600',...
%    sprintf('../images/stress_matlab_%d.pdf',data.angle_degree));
 
write_the_stress_table(data,element_stress_1,element_stress_2);
 
end
%==========================================
function write_the_stress_table(data,element_stress_1,...
                                                element_stress_2)
 
fileName         = [data.baseFolder sprintf(...
                                 '/data/stress_matlab_%d.tex',...
                    data.angle_degree)];
fileName         = strrep(fileName, '/', filesep);
[fileID0,errMsg] = fopen(fileName,'w');
 
if fileID0<0
    fprintf('Error opening %s\n, the message is [%s\n]',...
                                                fileName,errMsg);
    error(errMsg);
end
 
fprintf(fileID0,'\\begin{table}[!htbp]\n');
fprintf(fileID0,'\\centering\n');
fprintf(fileID0,'\\begin{minipage}{0.49\\textwidth}\n');
fprintf(fileID0,'\\centering\n');
fprintf(fileID0,'\\captionsetup{width=.95\\textwidth}\n');
fprintf(fileID0,'\\begin{tabular}{|l|l|l|l|}\\hline\n');
fprintf(fileID0,['global node \\# & $x$ & $y$ & $\\sigma_x$',...
                 '{\\footnotesize N/m^2} \\\\\\hline\n']);
for i=1:size(element_stress_1,1)
    if i==9
        fprintf(fileID0,...
               'center & $%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n',...
               element_stress_1(i,2),element_stress_1(i,3),...
               element_stress_1(i,4));
    elseif i==10
        fprintf(fileID0,['{\\footnotesize Gauss point 1}&',...
                       '$%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n'],...
            element_stress_1(i,2),element_stress_1(i,3),...
                                           element_stress_1(i,4));
    elseif i==11
        fprintf(fileID0,['{\\footnotesize Gauss point 2}&',...
                      '$%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n'],...
            element_stress_1(i,2),element_stress_1(i,3),...
                                           element_stress_1(i,4));
    elseif i==12
        fprintf(fileID0,['{\\footnotesize Gauss point 3}&',...
                        '$%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n'],...
            element_stress_1(i,2),element_stress_1(i,3),...
            element_stress_1(i,4));
    elseif i==13
        fprintf(fileID0,['{\\footnotesize Gauss point 4}&',...
                       '$%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n'],...
            element_stress_1(i,2),element_stress_1(i,3),...
            element_stress_1(i,4));
    else
        fprintf(fileID0,'$%d$ & $%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n',...
            element_stress_1(i,1),element_stress_1(i,2),...
            element_stress_1(i,3),...
            element_stress_1(i,4));
    end
end
fprintf(fileID0,'\\hline\n\\end{tabular}\n');
fprintf(fileID0,...
    ['\\caption{Matlab result. direct stress $\\sigma_x$ at',...
     ' each node, First element, angle [$%d$] degree}\n'],...
    data.angle_degree);
fprintf(fileID0,'\\end{minipage}\n');
fprintf(fileID0,'\\hfill\n');
 
fprintf(fileID0,'\\begin{minipage}{0.49\\textwidth}\n');
fprintf(fileID0,'\\centering\n');
fprintf(fileID0,'\\captionsetup{width=.95\\textwidth}\n');
 
fprintf(fileID0,'\\begin{tabular}{|l|l|l|l|}\\hline\n');
fprintf(fileID0,['global node \\# & $x$ & $y$ & $\\sigma_x$',...
                       ' {\\footnotesize N/m^2} \\\\\\hline\n']);
for i=1:size(element_stress_2,1)
    if i==9
        fprintf(fileID0,'center & $%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n',...
           element_stress_2(i,2),element_stress_2(i,3),...
           element_stress_2(i,4));
    elseif i==10
        fprintf(fileID0,['{\\footnotesize Gauss point 1}&',...
                      ' $%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n'],...
            element_stress_2(i,2),element_stress_2(i,3),...
            element_stress_2(i,4));
    elseif i==11
        fprintf(fileID0,['{\\footnotesize Gauss point 2}&',...
                      ' $%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n'],...
            element_stress_2(i,2),element_stress_2(i,3),...
            element_stress_2(i,4));
    elseif i==12
        fprintf(fileID0,['{\\footnotesize Gauss point 3}&',...
                       '$%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n'],...
            element_stress_2(i,2),element_stress_2(i,3),...
            element_stress_2(i,4));
    elseif i==13
        fprintf(fileID0,['{\\footnotesize Gauss point 4}&',...
                      ' $%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n'],...
            element_stress_2(i,2),element_stress_2(i,3),...
            element_stress_2(i,4));
    else
        fprintf(fileID0,'$%d$ & $%4.4f$ & $%4.4f$ & $%5.3f$ \\\\ \n',...
            element_stress_2(i,1),element_stress_2(i,2),...
            element_stress_2(i,3),...
            element_stress_2(i,4));
    end
end
fprintf(fileID0,'\\hline\n\\end{tabular}\n');
 
fprintf(fileID0,...
    ['\\caption{Matlab result. direct stress at each node,',...
     ' Second element, angle [$%d$] degree}\n'],...
      data.angle_degree);
fprintf(fileID0,'\\end{minipage}\n');
fprintf(fileID0,'\\end{table}\n');
fprintf(fileID0,'\\FloatBarrier\n');
 
fclose(fileID0);
 
end
%=======================================
function B = get_B(xi,eta,J)
%calculate the B matrix
 
B1 = [1,0,0,0;
      0,0,0,1;
      0,1,1,0];
 
gamma = 1/det(J) * [J(2,2) ,-J(1,2);
                    -J(2,1) ,J(1,1)];
 
B2    = [gamma, zeros(2,2);
         zeros(2,2),gamma];
 
Z  = zeros(2,1);
 
N1 = [dfdxi(1,xi,eta);
      dfdeta(1,xi,eta)];
 
N2 = [dfdxi(2,xi,eta);
      dfdeta(2,xi,eta)];
 
N3 = [dfdxi(3,xi,eta);
      dfdeta(3,xi,eta)];
 
N4 = [dfdxi(4,xi,eta);
      dfdeta(4,xi,eta)];
 
N5 = [dfdxi(5,xi,eta);
      dfdeta(5,xi,eta)];
 
N6 = [dfdxi(6,xi,eta);
      dfdeta(6,xi,eta)];
 
N7 = [dfdxi(7,xi,eta);
      dfdeta(7,xi,eta)];
 
N8 = [dfdxi(8,xi,eta);
      dfdeta(8,xi,eta)];
 
B3 = [N1,Z,  N2,Z,  N3,Z, N4,Z, N5,Z, N6,Z, N7,Z, N8,Z;
      Z,N1,  Z,N2,  Z,N3, Z,N4, Z,N5, Z,N6, Z,N7, Z,N8];
 
B  = B1*B2*B3;
end
%==========================================
function [x_coord,y_coord] = find_XY_coordinates(k,...
                                               elem_map_node,...
                                               global_coord_tbl)
%This function returns the x,y global coordinates of 
%specific element nodes
 
N       = size(elem_map_node,2);  %number of nodes in element
x_coord = zeros(N,1);  %x for this element
y_coord = zeros(N,1);  %y for this element
 
%collect this element node coordinates, go over each node 
%of this element and find its global x,y coordinates
for i = 1:N
  global_node_of_this_element_node = elem_map_node(k,i);
  x_coord(i) = global_coord_tbl(global_node_of_this_element_node,1);
  y_coord(i) = global_coord_tbl(global_node_of_this_element_node,2);
end
end
%======================================
function J = get_J(xi,eta,x_coord,y_coord)
 
J = [ ddxi(xi,eta,x_coord),  ddxi(xi,eta,y_coord);
      ddeta(xi,eta,x_coord), ddeta(xi,eta,y_coord)];
 
end
%==================================
function v = ddxi(xi,eta,c)
%find dx/d(xi) or dy/d(xi)
v = c(1)*(1/4*(1-eta^2)+(1/2)*(1-eta)*xi+(eta-1)/4)...
    +c(2)*(1/4*(eta^2-1)+(1/2)*(1-eta)*xi+(1-eta)/4)...
    +c(3)*(1/4*(eta^2-1)+(1/2)*(eta+1)*xi+(eta+1)/4)...
    +c(4)*(1/4*(1-eta^2)+(1/2)*(eta+1)*xi+(-eta-1)/4)...
    -c(5)*(1-eta)*xi...
    +c(6)*(1/2)*(1-eta^2)...
    -c(7)*xi*(eta+1)...
    -c(8)*(1/2)*(1-eta^2);
end
%==================================
function v = ddeta(xi,eta,c)
%find dx/d(eta) or dy/d(eta)
v = c(1)*(1/2*eta*(1-xi)+(1/4)*(1-xi^2)+(xi-1)/4)...
    +c(2)*((1/2)*(1+xi)*eta+1/4*(1-xi^2)+(-xi-1)/4)...
    +c(3)*((1/2)*eta*(1+xi)+1/4*(xi^2-1)+(xi+1)/4)...
    +c(4)*((1/2)*(1-xi)*eta+1/4*(xi^2-1)+(1-xi)/4)...
    -c(5)*(1/2)*(1-xi^2)...
    -c(6)*eta*(xi+1)...
    +c(7)*(1/2)*(1-xi^2)...
    -c(8)*eta*(1-xi);
end
%==================================
function v = dfdxi(shape_function_number,xi,eta)
%evaluate shape function at some x,y point
switch shape_function_number
    case 1
        v=(1/4)*(1-eta^2)+(1/2)*(1-eta)*xi+(eta-1)/4;
    case 2
        v=(1/4)*(eta^2-1)+(1/2)*(1-eta)*xi+(1-eta)/4;
    case 3
        v=(1/4)*(eta^2-1)+(1/2)*(eta+1)*xi+(eta+1)/4;
    case 4
        v=(1/4)*(1-eta^2)+(1/2)*(eta+1)*xi+(-eta-1)/4;
    case 5
        v=-(1-eta)*xi;
    case 6
        v=(1/2)*(1-eta^2);
    case 7
        v=-(eta+1)*xi;
    case 8
        v=(1/2)*(eta^2-1);
end
 
end
%==================================
function v = dfdeta(shape_function_number,xi,eta)
switch shape_function_number
    case 1
        v=(1/2)*eta*(1-xi)+(1/4)*(1-xi^2)+(xi-1)/4;
    case 2
        v=(1/2)*eta*(xi+1)+(1/4)*(1-xi^2)+(-xi-1)/4;
    case 3
        v=(1/2)*eta*(xi+1)+(1/4)*(xi^2-1)+(xi+1)/4;
    case 4
        v=(1/2)*eta*(1-xi)+(1/4)*(xi^2-1)+(-xi+1)/4;
    case 5
        v=(1/2)*(xi^2-1);
    case 6
        v=-eta*(xi+1);
    case 7
        v=(1/2)*(1-xi^2);
    case 8
        v=-eta*(1-xi);
end
 
end
%=============================
function k_global = merge_element_to_global(k,k_elem,...
                                         k_global,elem_map_dofs)
 
%assemble local element k to global K
 
for i=1:16
    for j =1:16
       global_i =  elem_map_dofs(k,i);
       global_j =  elem_map_dofs(k,j);
       k_global(global_i,global_j) = k_global(global_i,global_j)...
                                                    + k_elem(i,j);
    end
end
 
end
%======================
function draw_deflection(solution,data)
%This function is called at the end, after we have solved 
%the problem using finite elements and have nodal (x,y) 
%deformations. It plots the deformed shape against the undeformed
%original shape.
 
u = reshape(solution,2,13)';
L = data.L;
H = data.H;
 
figure();
%This is before demformation shape
 
top_line    = [0,L/4,L/2+data.shift,(3/4)*L,L;H,H,H,H,H];
left_line   = [0,0,0;H,H/2,0];
right_line  = [L,L,L;H,H/2,0];
bottom_line = [0,L/4,L/2-data.shift,(3/4)*L,L;0,0,0,0,0];
 
line(top_line(1,:),top_line(2,:),'LineStyle','--'); hold on;
line(left_line(1,:),left_line(2,:),'LineStyle','--');
line(right_line(1,:),right_line(2,:),'LineStyle','--');
line(bottom_line(1,:),bottom_line(2,:),'LineStyle','--');
 
%this is after adding deformation
line(top_line(1,:)+u(1:5,1)',top_line(2,:)+u(1:5,2)',...
                                               'Color','red');
line(left_line(1,:)+u([1 6 9],1)',left_line(2,:)+u([1 6 9],2)',...
                                                   'Color','red');
line(right_line(1,:)+u([5 8 13],1)',right_line(2,:)+...
                                   u([5 8 13],2)','Color','red');
line(bottom_line(1,:)+u(9:13,1)',bottom_line(2,:)+...
                                       u(9:13,2)','Color','red');
 
%There are the nodes. Draw nodes on top line
for i=1:size(top_line,2)
    plot(top_line(1,i)+u(i,1),top_line(2,i)+u(i,2),'bo');
end
 
%draw nodes on left
idx=5;
plot(left_line(1,2)+u(1+idx,1),left_line(2,2)+u(1+idx,2),'bo');
 
%draw nodes on right
idx=7;
plot(right_line(1,2)+u(1+idx,1),right_line(2,2)+u(1+idx,2),'bo');
 
%draw nodes in middle
idx=6;
plot(L/2+u(1+idx,1),H/2+u(1+idx,2),'bo');
 
%Draw nodes on bottom line
idx=8;
for i=1:size(bottom_line,2)
    plot(bottom_line(1,i)+u(i+idx,1),bottom_line(2,i)+...
                                                u(i+idx,2),'bo');
end
 
%draw dashed line between elements
line([bottom_line(1,3)+u(11,1),...
      L/2+u(7,1),...
      L/2+data.shift+u(3,1)],...
      [bottom_line(2,3)+u(11,2),...
      H/2+u(7,2),...
      H+u(3,2)],...
      'LineStyle','-.');
 
%put title, x,y arrows at (0,0) and save the image to include in
%document/report at end
 
title(sprintf('deflection: $y=%3.4f$ m, $x=%3.4f$ m, angle $%d^o$',...
    u(end,2),u(end,1),data.angle_degree),...
    'interpreter','Latex','Fontsize',11);
 
xlabel('Length in meters ($x$ direction)','interpreter',...
                                          'Latex','Fontsize',11);
ylabel('Height in meters ($y$ direction)','interpreter',...
                                          'Latex','Fontsize',11);
 
quiver(0,0,1,0,1,'MaxHeadSize',0.5,'Color','black');
text(1.1,.2,'$x$','interpreter','Latex');
quiver(0,0,0,1,1,'MaxHeadSize',0.5,'Color','black');
text(0.1,1.1,'$y$','interpreter','Latex');
axis equal;
xlim([-0.5,L+0.5]);
ylim([-1.5,H+1]);
grid;
 
%print(gcf, '-dpdf', '-r600',...
%sprintf('../images/deflection_matlab_%d.pdf',data.angle_degree));
 
end
 
%===============================
function v = get_shape_function(shape_function_number,xi,eta)
 
switch shape_function_number
    case 1
        v=-(1/4)*(1-eta^2)*(1-xi)-(1/4)*(1-eta)*(1-xi^2)+... 
                                            (1/4)*(1-eta)*(1-xi);
    case 2
        v=-(1/4)*(1-eta^2)*(1+xi)-(1/4)*(1-eta)*...
                                   (1-xi^2)+(1/4)*(1-eta)*(1+xi);
    case 3
        v=-(1/4)*(1-eta^2)*(1+xi)-(1/4)*(1+eta)*...
                                   (1-xi^2)+(1/4)*(1+eta)*(1+xi);
    case 4
        v=-(1/4)*(1-eta^2)*(1-xi)-(1/4)*(1+eta)*...
                                   (1-xi^2)+(1/4)*(1+eta)*(1-xi);
    case 5
        v=(1/2)*(1-eta)*(1-xi^2);
    case 6
        v=(1/2)*(1-eta^2)*(1+xi);
    case 7
        v=(1/2)*(1+eta)*(1-xi^2);
    case 8
        v=(1/2)*(1-eta^2)*(1-xi);
end
 
end
%=======================================
function verify_shape_functions_sum_to_one(gs)
for i=1:2
    xi = gs(i);
    for j=1:2
        eta   = gs(j);
        chk_1 = 0;
        for k=1:8 %sum all shape functions at this Gaussian point
            chk_1 = chk_1 + get_shape_function(k,xi,eta);
        end
        if chk_1 ~= 1
            error(['Internal error. sum of shape functions',...
                   ' not 1 at $\xi=%3.3f,\eta=%3.3f'],...
                xi,eta);
        end
    end
end
end
%=================================
function display_diagram_of_dof(L,H)
 
figure();
top_line    = [0,L/4,L/2,(3/4)*L,L;
    H,H,H,H,H];
left_line   = [0,0,0;
    H,H/2,0];
right_line  = [L,L,L;
    H,H/2,0];
bottom_line = [0,L/4,L/2,(3/4)*L,L;
    0,0,0,0,0];
middle_line = [L/2;H/2];
 
line(top_line(1,:),top_line(2,:)); hold on;
%axis equal;
xlim([-2,L+2]);
ylim([-.75,H+1]);
line(left_line(1,:),left_line(2,:));
line(right_line(1,:),right_line(2,:));
line(bottom_line(1,:),bottom_line(2,:));
 
k=0;
node_number=0;
for i=1:size(top_line,2)
    x=top_line(1,i); y=top_line(2,i);
    plot(x,y,'ro');
    quiver(x,y,0.5,0,1,'MaxHeadSize',2,'Color','black');
    k=k+1;
    node_number=node_number+1;
    text(x+.1,y-.1,sprintf('$%d$',node_number),...
                                       'interpreter','Latex',...
        'Fontsize',11,'Color','red');
    text(x+.3,y+.1,sprintf('$%d$',k),'interpreter',...
                                         'Latex','Fontsize',11);
    quiver(x,y,0,0.5,1,'MaxHeadSize',2,'Color','black');
    k=k+1;
    text(x,y+.6,sprintf('$%d$',k),'interpreter',...
                                          'Latex','Fontsize',11);
end
for i=1:size(left_line,2)
    x=left_line(1,i); y=left_line(2,i);
    plot(x,y,'ro');
    quiver(x,y,0.5,0,1,'MaxHeadSize',2,'Color','black');
    quiver(x,y,0,0.5,1,'MaxHeadSize',2,'Color','black');
end
k=k+1;
node_number=node_number+1;
text(.1,H/2-.1,sprintf('$%d$',node_number),'interpreter','Latex',...
    'Fontsize',11,'Color','red');
text(.5,H/2,sprintf('%d',k),'interpreter','Latex','Fontsize',11);
k=k+1;
text(x,H/2+.5,sprintf('%d',k),'interpreter','Latex','Fontsize',11);
 
for i=1:size(middle_line,2)
    x=middle_line(1,i); y=middle_line(2,i);
    plot(x,y,'ro');
    quiver(x,y,0.5,0,1,'MaxHeadSize',2,'Color','black');
    k=k+1;
    node_number=node_number+1;
    text(x+.1,H/2-.15,sprintf('$%d$',node_number),...
                                        'interpreter','Latex',...
        'Fontsize',11,'Color','red');
    text(x+.5,H/2,sprintf('%d',k),'interpreter',...
                                          'Latex','Fontsize',11);
    quiver(x,y,0,0.5,1,'MaxHeadSize',2,'Color','black');
    k=k+1;
    text(x+.1,H/2+.4,sprintf('%d',k),'interpreter',...
                                          'Latex','Fontsize',11);
end
 
for i=1:size(right_line,2)
    x=right_line(1,i); y=right_line(2,i);
    plot(x,y,'ro');
    quiver(x,y,0.5,0,1,'MaxHeadSize',2,'Color','black');
    quiver(x,y,0,0.5,1,'MaxHeadSize',2,'Color','black');
end
k=k+1;
node_number=node_number+1;
text(L+.1,H/2-.1,sprintf('$%d$',node_number),...
          'interpreter','Latex','Fontsize',11,'Color','red');
text(L+.5,H/2,sprintf('%d',k),'interpreter','Latex','Fontsize',11);
k=k+1;
text(L,H/2+.5,sprintf('%d',k),'interpreter','Latex','Fontsize',11);
 
 
for i=1:size(bottom_line,2)
    x=bottom_line(1,i); y=bottom_line(2,i);
    plot(x,y,'ro');
    quiver(x,y,0.5,0,1,'MaxHeadSize',2,'Color','black');
    k=k+1;
    node_number=node_number+1;
    text(x-.1,-.2,sprintf('$%d$',node_number),...
                                       'interpreter','Latex',...
        'Fontsize',11,'Color','red');
    text(x+.4,.1,sprintf('%d',k));
    quiver(x,y,0,0.5,1,'MaxHeadSize',2,'Color','black');
    k=k+1;
    text(x+.1,.4,sprintf('%d',k));
end
 
title({'global D.O.F. numbering used (black letters).',...
    'with associated global node numering (in red letters)'},...
    'interpreter','Latex','Fontsize',11);
 
%
%print(gcf, '-dpdf', '-r600', sprintf('../images/dof.pdf'));
%
end


global node # global node coordinates	$x_i$	$y_i$

1	$0$	$H$

2	$\frac{L}{4}$	$H$

3	$\frac{L}{2}+\Delta$	$H$

4	$\frac{3}{4}L$	$H$

5	$L$	$H$

6	$0$	$\frac{H}{2}$

7	$\frac{L}{2}$	$\frac{H}{2}$

8	$L$	$\frac{H}{2}$

9	$0$	$0$

10	$\frac{L}{4}$	$0$

11	$\frac{L}{2}-\Delta$	$0$

12	$\frac{3}{4}L$	$0$

13	$L$	$0$


element # element node #	node 1		node 2		node 3		node 4		node 5		node 6		node 7		node 8

$1$	$17$	$18$	$21$	$22$	$5$	$6$	$1$	$2$	$19$	$20$	$13$	$14$	$3$	$4$	$11$	$12$

$2$	$21$	$22$	$25$	$26$	$9$	$10$	$5$	$6$	$23$	$24$	$15$	$16$	$7$	$8$	$13$	$14$


global node #	$x$ (meter)	$y$ (meter)

$1$	$0.000000$	$-0.004650$
$2$	$0.065794$	$-0.096522$
$3$	$0.112725$	$-0.329300$
$4$	$0.140794$	$-0.655828$
$5$	$0.150000$	$-1.028950$
$6$	$0.000000$	$0.000000$
$7$	$-0.000000$	$-0.327050$
$8$	$-0.000000$	$-1.029100$
$9$	$0.000000$	$-0.004650$
$10$	$-0.065794$	$-0.096522$
$11$	$-0.112725$	$-0.329300$
$12$	$-0.140794$	$-0.655828$
$13$	$-0.150000$	$-1.028950$


global node #	$x$	$y$	$\sigma _x$ N/m2

$9$	$0.0000$	$0.0000$	$-300.000$
$11$	$5.0000$	$0.0000$	$-150.000$
$3$	$5.0000$	$2.0000$	$150.000$
$1$	$0.0000$	$2.0000$	$300.000$
$10$	$2.5000$	$0.0000$	$-225.000$
$7$	$5.0000$	$1.0000$	$-0.000$
$2$	$2.5000$	$2.0000$	$225.000$
$6$	$0.0000$	$1.0000$	$-0.000$
center	$2.5000$	$1.0000$	$-0.000$
Gauss point 1	$1.0566$	$0.4226$	$-154.904$
Gauss point 2	$3.9434$	$0.4226$	$-104.904$
Gauss point 3	$3.9434$	$1.5774$	$104.904$
Gauss point 4	$1.0566$	$1.5774$	$154.904$


global node #	$x$	$y$	$\sigma _x$ N/m2

$11$	$5.0000$	$0.0000$	$-150.000$
$13$	$10.0000$	$0.0000$	$-0.000$
$5$	$10.0000$	$2.0000$	$0.000$
$3$	$5.0000$	$2.0000$	$150.000$
$12$	$7.5000$	$0.0000$	$-75.000$
$8$	$10.0000$	$1.0000$	$0.000$
$4$	$7.5000$	$2.0000$	$75.000$
$7$	$5.0000$	$1.0000$	$0.000$
center	$7.5000$	$1.0000$	$0.000$
Gauss point 1	$6.0566$	$0.4226$	$-68.301$
Gauss point 2	$8.9434$	$0.4226$	$-18.301$
Gauss point 3	$8.9434$	$1.5774$	$18.301$
Gauss point 4	$6.0566$	$1.5774$	$68.301$

Chapter 5
my project report

5.1 Description of the problem, geometry and element description

5.1.1 Problem statement

5.1.2 Geometry, nodes and element description

5.2 Theory and analytical derivation

5.2.1 Shape functions

5.2.2 Finding the Jacobian

5.2.3 Stress recovery

5.3 Results

5.3.1 No element distortion. zero angle

summary of result

Matlab result

ANSYS result

5.3.2 15 degrees distortion

summary of result

Matlab result

ANSYS result

5.3.3 30 degrees distortion

summary of result

Matlab result

ANSYS result

5.3.4 45 degrees distortion

summary of result

Matlab result

ANSYS result

5.3.5 50 degrees distortion

summary of result

Matlab result

ANSYS result

5.3.6 55 degrees distortion

summary of result

Matlab result

ANSYS result

5.4 Observations, discussion and conclusions

5.4.1 References

5.5 Appendix

5.5.1 APDL used for ANSYS

5.5.2 Matlab source code


	$x$ (meter)	$y$ (meter)
Matlab	$-0.15$	$-1.02895$
ANSYS	$-0.15$	$-1.0289$

Chapter 5my project report

5.1 Description of the problem, geometry and element description

5.1.1 Problem statement

5.1.2 Geometry, nodes and element description

5.2 Theory and analytical derivation

5.2.1 Shape functions

5.2.2 Finding the Jacobian

5.2.3 Stress recovery

5.3 Results

5.3.1 No element distortion. zero angle

summary of result

Matlab result

ANSYS result

5.3.2 15 degrees distortion

summary of result

Matlab result

ANSYS result

5.3.3 30 degrees distortion

summary of result

Matlab result

ANSYS result

5.3.4 45 degrees distortion

summary of result

Matlab result

ANSYS result

5.3.5 50 degrees distortion

summary of result

Matlab result

ANSYS result

5.3.6 55 degrees distortion

summary of result

Matlab result

ANSYS result

5.4 Observations, discussion and conclusions

5.4.1 References

5.5 Appendix

5.5.1 APDL used for ANSYS

5.5.2 Matlab source code

Chapter 5
my project report