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 first and the second elements as shown in the above diagram. The method of finite 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 stiffness matrix was assembled from the two elements stiffness matrices, the system equations given by \(K D = F\) were solved using the direct linear system solver.
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 finite elements method for each of the different 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 effect of the element distortion on the accuracy of the finite element result using the element selected.
The idea is that a good finite element should produce the same displacement at its nodes regardless of how it was fitted to the physical region in place. This report was to determine if the element selected would still produce good results when deformed. The first 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 first 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 specific 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\).
| element # element node # | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| 1 | 9 | 11 | 3 | 1 | 10 | 7 | 2 | 6 |
| 2 | 11 | 13 | 5 | 3 | 12 | 8 | 4 | 7 |
There is also a need to lookup the global coordinates \(\{x_i,y_i\}\) given the global node number. This is needed when finding 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.
| 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\) |
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 stiffness matrix.
Using the method in the project handout, the following table was generated.
| 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\) |
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.
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*}
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*}
Each of the above derivatives is evaluated and used in the Matlab code.
\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 find the Jacobian and its determinant and also find \(\Gamma \) and the matrix \(B_2\). To find 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 stiffness 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 stiffness matrices \(k_{\text{elem}}\) are then combined to make the global stiffness 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 stiffness 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.
This is a discussion of how the stress at the elements nodes was found. Initially the direct method was used to find 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.
The following diagram shows this layout more clearly.
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 different points for illustration
Now that the mapping between \(\xi ,\eta \) and \((r,s)\) is determined, the next step was to find \((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 first, and then the following extrapolation formula is applied
\begin{equation*} \sigma _p = \sum _{i=1}^{4} N_i \sigma _i \end{equation*}
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 find the stress at point \((\xi ,\eta )=(0,-1)\), the first 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 find \(\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.
The results are listed by the angle \(\theta \) used to distort the elements with. For each angle, the deflection 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 first given that compares Matlab and ANSYS result. This short summary contains only the deflection at the bottom right corner, which is node \(13\). After this, the full result and stress plots are given.
| \(x\) (meter) | \(y\) (meter) | |
| Matlab | \(-0.15\) | \(-1.02895\) |
| ANSYS | \(-0.15\) | \(-1.0289\) |
| 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\) |
The following figure shows the deformation found
The following table shows Matlab result for the direct stress \(\sigma _x\) found at each node for each element.
| 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\) |
The following shows the direct stress contour generated in Matlab
The following figure shows the deformation found
The following table shows ANSYS result for the direct stress \(\sigma _x\) found at each node for each element.
The following shows the direct stress contour generated in ANSYS
| \(x\) (meter) | \(y\) (meter) | |
| Matlab | \(-0.150262\) | \(-1.02555\) |
| ANSYS | \(-0.15026\) | \(-1.0256\) |
| global node # | \(x\) (meter) | \(y\) (meter) |
| \(1\) | \(0.000000\) | \(-0.008808\) |
| \(2\) | \(0.066221\) | \(-0.096603\) |
| \(3\) | \(0.115141\) | \(-0.356575\) |
| \(4\) | \(0.138750\) | \(-0.653161\) |
| \(5\) | \(0.146026\) | \(-1.012233\) |
| \(6\) | \(0.000000\) | \(0.000000\) |
| \(7\) | \(0.000596\) | \(-0.326056\) |
| \(8\) | \(0.001059\) | \(-1.018939\) |
| \(9\) | \(0.000000\) | \(-0.000457\) |
| \(10\) | \(-0.064844\) | \(-0.096198\) |
| \(11\) | \(-0.108540\) | \(-0.300195\) |
| \(12\) | \(-0.139066\) | \(-0.648279\) |
| \(13\) | \(-0.150262\) | \(-1.025550\) |
The following figure shows the deformation found
The following table shows Matlab result for the direct stress \(\sigma _x\) found at each node for each element.
| global node # | \(x\) | \(y\) | \(\sigma _x\) N/m2 |
| \(9\) | \(0.0000\) | \(0.0000\) | \(-299.049\) |
| \(11\) | \(4.7321\) | \(0.0000\) | \(-148.195\) |
| \(3\) | \(5.2679\) | \(2.0000\) | \(144.369\) |
| \(1\) | \(0.0000\) | \(2.0000\) | \(300.912\) |
| \(10\) | \(2.5000\) | \(0.0000\) | \(-223.622\) |
| \(7\) | \(5.0000\) | \(1.0000\) | \(-1.913\) |
| \(2\) | \(2.5000\) | \(2.0000\) | \(222.641\) |
| \(6\) | \(0.0000\) | \(1.0000\) | \(0.932\) |
| center | \(2.5000\) | \(1.0000\) | \(-0.491\) |
| Gauss point 1 | \(1.0566\) | \(0.4226\) | \(-154.111\) |
| Gauss point 2 | \(3.9434\) | \(0.4226\) | \(-104.520\) |
| Gauss point 3 | \(3.9434\) | \(1.5774\) | \(101.897\) |
| Gauss point 4 | \(1.0566\) | \(1.5774\) | \(154.772\) |
| global node # | \(x\) | \(y\) | \(\sigma _x\) N/m2 |
| \(11\) | \(4.7321\) | \(0.0000\) | \(-146.183\) |
| \(13\) | \(10.0000\) | \(0.0000\) | \(-0.994\) |
| \(5\) | \(10.0000\) | \(2.0000\) | \(2.746\) |
| \(3\) | \(5.2679\) | \(2.0000\) | \(142.734\) |
| \(12\) | \(7.5000\) | \(0.0000\) | \(-73.588\) |
| \(8\) | \(10.0000\) | \(1.0000\) | \(0.876\) |
| \(4\) | \(7.5000\) | \(2.0000\) | \(72.740\) |
| \(7\) | \(5.0000\) | \(1.0000\) | \(-1.724\) |
| center | \(7.5000\) | \(1.0000\) | \(-0.424\) |
| Gauss point 1 | \(6.0566\) | \(0.4226\) | \(-67.181\) |
| Gauss point 2 | \(8.9434\) | \(0.4226\) | \(-18.150\) |
| Gauss point 3 | \(8.9434\) | \(1.5774\) | \(18.803\) |
| Gauss point 4 | \(6.0566\) | \(1.5774\) | \(64.831\) |
The following shows the direct stress contour generated in Matlab
The following figure shows the deformation found
The following table shows ANSYS result for the direct stress \(\sigma _x\) found at each node for each element.
The following shows the direct stress contour generated in ANSYS
| \(x\) (meter) | \(y\) (meter) | |
| Matlab | \(-0.146118\) | \(-0.998214\) |
| ANSYS | \(-0.14612\) | \(-0.99821\) |
| global node # | \(x\) (meter) | \(y\) (meter) |
| \(1\) | \(0.000000\) | \(-0.013054\) |
| \(2\) | \(0.065955\) | \(-0.096602\) |
| \(3\) | \(0.115155\) | \(-0.384450\) |
| \(4\) | \(0.132363\) | \(-0.638148\) |
| \(5\) | \(0.137945\) | \(-0.972694\) |
| \(6\) | \(0.000000\) | \(0.000000\) |
| \(7\) | \(0.001094\) | \(-0.322637\) |
| \(8\) | \(0.002043\) | \(-0.985167\) |
| \(9\) | \(0.000000\) | \(0.003878\) |
| \(10\) | \(-0.063350\) | \(-0.095375\) |
| \(11\) | \(-0.102487\) | \(-0.265899\) |
| \(12\) | \(-0.132998\) | \(-0.628986\) |
| \(13\) | \(-0.146118\) | \(-0.998214\) |
The following figure shows the deformation found
The following table shows Matlab result for the direct stress \(\sigma _x\) found at each node for each element.
| global node # | \(x\) | \(y\) | \(\sigma _x\) N/m2 |
| \(9\) | \(0.0000\) | \(0.0000\) | \(-297.609\) |
| \(11\) | \(4.4226\) | \(0.0000\) | \(-138.871\) |
| \(3\) | \(5.5774\) | \(2.0000\) | \(128.858\) |
| \(1\) | \(0.0000\) | \(2.0000\) | \(302.166\) |
| \(10\) | \(2.5000\) | \(0.0000\) | \(-218.240\) |
| \(7\) | \(5.0000\) | \(1.0000\) | \(-5.007\) |
| \(2\) | \(2.5000\) | \(2.0000\) | \(215.512\) |
| \(6\) | \(0.0000\) | \(1.0000\) | \(2.279\) |
| center | \(2.5000\) | \(1.0000\) | \(-1.364\) |
| Gauss point 1 | \(1.0566\) | \(0.4226\) | \(-152.145\) |
| Gauss point 2 | \(3.9434\) | \(0.4226\) | \(-101.010\) |
| Gauss point 3 | \(3.9434\) | \(1.5774\) | \(94.076\) |
| Gauss point 4 | \(1.0566\) | \(1.5774\) | \(153.623\) |
| global node # | \(x\) | \(y\) | \(\sigma _x\) N/m2 |
| \(11\) | \(4.4226\) | \(0.0000\) | \(-129.964\) |
| \(13\) | \(10.0000\) | \(0.0000\) | \(-6.206\) |
| \(5\) | \(10.0000\) | \(2.0000\) | \(9.722\) |
| \(3\) | \(5.5774\) | \(2.0000\) | \(123.499\) |
| \(12\) | \(7.5000\) | \(0.0000\) | \(-68.085\) |
| \(8\) | \(10.0000\) | \(1.0000\) | \(1.758\) |
| \(4\) | \(7.5000\) | \(2.0000\) | \(66.611\) |
| \(7\) | \(5.0000\) | \(1.0000\) | \(-3.232\) |
| center | \(7.5000\) | \(1.0000\) | \(-0.737\) |
| Gauss point 1 | \(6.0566\) | \(0.4226\) | \(-60.856\) |
| Gauss point 2 | \(8.9434\) | \(0.4226\) | \(-18.386\) |
| Gauss point 3 | \(8.9434\) | \(1.5774\) | \(19.792\) |
| Gauss point 4 | \(6.0566\) | \(1.5774\) | \(56.500\) |
The following shows the direct stress contour generated in Matlab
The following figure shows the deformation found
The following table shows ANSYS result for the direct stress \(\sigma _x\) found at each node for each element.
The following shows the direct stress contour generated in ANSYS
| \(x\) (meter) | \(y\) (meter) | |
| Matlab | \(-0.135417\) | \(-0.934503\) |
| ANSYS | \(-0.13542\) | \(-0.93450\) |
| global node # | \(x\) (meter) | \(y\) (meter) |
| \(1\) | \(0.000000\) | \(-0.017363\) |
| \(2\) | \(0.064453\) | \(-0.096794\) |
| \(3\) | \(0.110558\) | \(-0.415091\) |
| \(4\) | \(0.119874\) | \(-0.603430\) |
| \(5\) | \(0.124609\) | \(-0.901330\) |
| \(6\) | \(0.000000\) | \(0.000000\) |
| \(7\) | \(0.001251\) | \(-0.315104\) |
| \(8\) | \(0.002702\) | \(-0.916990\) |
| \(9\) | \(0.000000\) | \(0.008169\) |
| \(10\) | \(-0.061186\) | \(-0.093436\) |
| \(11\) | \(-0.093955\) | \(-0.220369\) |
| \(12\) | \(-0.120788\) | \(-0.592443\) |
| \(13\) | \(-0.135417\) | \(-0.934503\) |
The following figure shows the deformation found
The following table shows Matlab result for the direct stress \(\sigma _x\) found at each node for each element.
| global node # | \(x\) | \(y\) | \(\sigma _x\) N/m2 |
| \(9\) | \(0.0000\) | \(0.0000\) | \(-294.696\) |
| \(11\) | \(4.0000\) | \(0.0000\) | \(-120.055\) |
| \(3\) | \(6.0000\) | \(2.0000\) | \(96.964\) |
| \(1\) | \(0.0000\) | \(2.0000\) | \(304.372\) |
| \(10\) | \(2.5000\) | \(0.0000\) | \(-207.376\) |
| \(7\) | \(5.0000\) | \(1.0000\) | \(-11.546\) |
| \(2\) | \(2.5000\) | \(2.0000\) | \(200.668\) |
| \(6\) | \(0.0000\) | \(1.0000\) | \(4.838\) |
| center | \(2.5000\) | \(1.0000\) | \(-3.354\) |
| Gauss point 1 | \(1.0566\) | \(0.4226\) | \(-148.254\) |
| Gauss point 2 | \(3.9434\) | \(0.4226\) | \(-94.038\) |
| Gauss point 3 | \(3.9434\) | \(1.5774\) | \(77.871\) |
| Gauss point 4 | \(1.0566\) | \(1.5774\) | \(151.005\) |
| global node # | \(x\) | \(y\) | \(\sigma _x\) N/m2 |
| \(11\) | \(4.0000\) | \(0.0000\) | \(-98.498\) |
| \(13\) | \(10.0000\) | \(0.0000\) | \(-17.052\) |
| \(5\) | \(10.0000\) | \(2.0000\) | \(22.022\) |
| \(3\) | \(6.0000\) | \(2.0000\) | \(91.497\) |
| \(12\) | \(7.5000\) | \(0.0000\) | \(-57.775\) |
| \(8\) | \(10.0000\) | \(1.0000\) | \(2.485\) |
| \(4\) | \(7.5000\) | \(2.0000\) | \(56.759\) |
| \(7\) | \(5.0000\) | \(1.0000\) | \(-3.501\) |
| center | \(7.5000\) | \(1.0000\) | \(-0.508\) |
| Gauss point 1 | \(6.0566\) | \(0.4226\) | \(-47.876\) |
| Gauss point 2 | \(8.9434\) | \(0.4226\) | \(-19.267\) |
| Gauss point 3 | \(8.9434\) | \(1.5774\) | \(21.706\) |
| Gauss point 4 | \(6.0566\) | \(1.5774\) | \(43.404\) |
The following shows the direct stress contour generated in Matlab
The following figure shows the deformation found
The following table shows ANSYS result for the direct stress \(\sigma _x\) found at each node for each element.
The following shows the direct stress contour generated in ANSYS
| \(x\) (meter) | \(y\) (meter) | |
| Matlab | \(-0.129803\) | \(-0.901557\) |
| ANSYS | \(-0.12980\) | \(-0.90156\) |
| global node # | \(x\) (meter) | \(y\) (meter) |
| \(1\) | \(0.000000\) | \(-0.018726\) |
| \(2\) | \(0.063489\) | \(-0.097037\) |
| \(3\) | \(0.107149\) | \(-0.426310\) |
| \(4\) | \(0.113960\) | \(-0.585035\) |
| \(5\) | \(0.118879\) | \(-0.868383\) |
| \(6\) | \(0.000000\) | \(0.000000\) |
| \(7\) | \(0.001133\) | \(-0.311076\) |
| \(8\) | \(0.002731\) | \(-0.883753\) |
| \(9\) | \(0.000000\) | \(0.009386\) |
| \(10\) | \(-0.060301\) | \(-0.092294\) |
| \(11\) | \(-0.090400\) | \(-0.200709\) |
| \(12\) | \(-0.114927\) | \(-0.574883\) |
| \(13\) | \(-0.129803\) | \(-0.901557\) |
The following figure shows the deformation found
The following table shows Matlab result for the direct stress \(\sigma _x\) found at each node for each element.
| global node # | \(x\) | \(y\) | \(\sigma _x\) N/m2 |
| \(9\) | \(0.0000\) | \(0.0000\) | \(-292.998\) |
| \(11\) | \(3.8082\) | \(0.0000\) | \(-111.352\) |
| \(3\) | \(6.1918\) | \(2.0000\) | \(80.771\) |
| \(1\) | \(0.0000\) | \(2.0000\) | \(305.494\) |
| \(10\) | \(2.5000\) | \(0.0000\) | \(-202.175\) |
| \(7\) | \(5.0000\) | \(1.0000\) | \(-15.290\) |
| \(2\) | \(2.5000\) | \(2.0000\) | \(193.133\) |
| \(6\) | \(0.0000\) | \(1.0000\) | \(6.248\) |
| center | \(2.5000\) | \(1.0000\) | \(-4.521\) |
| Gauss point 1 | \(1.0566\) | \(0.4226\) | \(-146.283\) |
| Gauss point 2 | \(3.9434\) | \(0.4226\) | \(-90.990\) |
| Gauss point 3 | \(3.9434\) | \(1.5774\) | \(69.513\) |
| Gauss point 4 | \(1.0566\) | \(1.5774\) | \(149.676\) |
| global node # | \(x\) | \(y\) | \(\sigma _x\) N/m2 |
| \(11\) | \(3.8082\) | \(0.0000\) | \(-84.728\) |
| \(13\) | \(10.0000\) | \(0.0000\) | \(-22.141\) |
| \(5\) | \(10.0000\) | \(2.0000\) | \(27.249\) |
| \(3\) | \(6.1918\) | \(2.0000\) | \(79.463\) |
| \(12\) | \(7.5000\) | \(0.0000\) | \(-53.435\) |
| \(8\) | \(10.0000\) | \(1.0000\) | \(2.554\) |
| \(4\) | \(7.5000\) | \(2.0000\) | \(53.356\) |
| \(7\) | \(5.0000\) | \(1.0000\) | \(-2.633\) |
| center | \(7.5000\) | \(1.0000\) | \(-0.039\) |
| Gauss point 1 | \(6.0566\) | \(0.4226\) | \(-41.931\) |
| Gauss point 2 | \(8.9434\) | \(0.4226\) | \(-19.803\) |
| Gauss point 3 | \(8.9434\) | \(1.5774\) | \(22.719\) |
| Gauss point 4 | \(6.0566\) | \(1.5774\) | \(38.858\) |
The following shows the direct stress contour generated in Matlab
The following figure shows the deformation found
The following table shows ANSYS result for the direct stress \(\sigma _x\) found at each node for each element.
The following shows the direct stress contour generated in ANSYS
| \(x\) (meter) | \(y\) (meter) | |
| Matlab | \(-0.123001\) | \(-0.86157\) |
| ANSYS | \(-0.123\) | \(-0.86157\) |
| global node # | \(x\) (meter) | \(y\) (meter) |
| \(1\) | \(0.000000\) | \(-0.019943\) |
| \(2\) | \(0.062204\) | \(-0.097514\) |
| \(3\) | \(0.102304\) | \(-0.438312\) |
| \(4\) | \(0.107168\) | \(-0.562073\) |
| \(5\) | \(0.112704\) | \(-0.830774\) |
| \(6\) | \(0.000000\) | \(0.000000\) |
| \(7\) | \(0.000874\) | \(-0.305999\) |
| \(8\) | \(0.002574\) | \(-0.844631\) |
| \(9\) | \(0.000000\) | \(0.010255\) |
| \(10\) | \(-0.059355\) | \(-0.090697\) |
| \(11\) | \(-0.086450\) | \(-0.177473\) |
| \(12\) | \(-0.108133\) | \(-0.554224\) |
| \(13\) | \(-0.123001\) | \(-0.861570\) |
The following figure shows the deformation found
The following table shows Matlab result for the direct stress \(\sigma _x\) found at each node for each element.
| global node # | \(x\) | \(y\) | \(\sigma _x\) N/m2 |
| \(9\) | \(0.0000\) | \(0.0000\) | \(-290.606\) |
| \(11\) | \(3.5719\) | \(0.0000\) | \(-101.874\) |
| \(3\) | \(6.4281\) | \(2.0000\) | \(61.164\) |
| \(1\) | \(0.0000\) | \(2.0000\) | \(306.879\) |
| \(10\) | \(2.5000\) | \(0.0000\) | \(-196.240\) |
| \(7\) | \(5.0000\) | \(1.0000\) | \(-20.355\) |
| \(2\) | \(2.5000\) | \(2.0000\) | \(184.022\) |
| \(6\) | \(0.0000\) | \(1.0000\) | \(8.137\) |
| center | \(2.5000\) | \(1.0000\) | \(-6.109\) |
| Gauss point 1 | \(1.0566\) | \(0.4226\) | \(-143.860\) |
| Gauss point 2 | \(3.9434\) | \(0.4226\) | \(-87.902\) |
| Gauss point 3 | \(3.9434\) | \(1.5774\) | \(59.234\) |
| Gauss point 4 | \(1.0566\) | \(1.5774\) | \(148.092\) |
| global node # | \(x\) | \(y\) | \(\sigma _x\) N/m2 |
| \(11\) | \(3.5719\) | \(0.0000\) | \(-70.386\) |
| \(13\) | \(10.0000\) | \(0.0000\) | \(-27.809\) |
| \(5\) | \(10.0000\) | \(2.0000\) | \(32.546\) |
| \(3\) | \(6.4281\) | \(2.0000\) | \(69.339\) |
| \(12\) | \(7.5000\) | \(0.0000\) | \(-49.097\) |
| \(8\) | \(10.0000\) | \(1.0000\) | \(2.369\) |
| \(4\) | \(7.5000\) | \(2.0000\) | \(50.943\) |
| \(7\) | \(5.0000\) | \(1.0000\) | \(-0.523\) |
| center | \(7.5000\) | \(1.0000\) | \(0.923\) |
| Gauss point 1 | \(6.0566\) | \(0.4226\) | \(-35.405\) |
| Gauss point 2 | \(8.9434\) | \(0.4226\) | \(-20.508\) |
| Gauss point 3 | \(8.9434\) | \(1.5774\) | \(24.022\) |
| Gauss point 4 | \(6.0566\) | \(1.5774\) | \(35.581\) |
The following shows the direct stress contour generated in Matlab
The following figures show the deformation found
The following table shows ANSYS result for the direct stress \(\sigma _x\) found at each node for each element.
The following shows the direct stress contour generated in ANSYS
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 element1 . These elements are known not to give good results when they distorted2 . 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 deflection reported by the finite 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 signification error in accuracy. The following graph shows how the error in the vertical deflection 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 fitting 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
Once distortion is added, it produced bad result in stress values when compared to ANSYS results.
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.
The following is the APDL script used for ANSYS analysis. ANSYS version 16.02, student version was used.
The following is the listing of the m file for Matlab implementation. For making the contour plot, an external m file was used from Mathworks file exchange called tricontf and is included in the zip file. This function is not listed here.
To run the Matlab program, the command is nma_project_EMA_471