Chapter 5
my project report

 5.1 Description of the problem, geometry and element description
 5.2 Theory and analytical derivation
 5.3 Results
 5.4 Observations, discussion and conclusions
 5.5 Appendix
  1. This report in one PDF (letter size)
  2. This report in one PDF (legal size)

5.1 Description of the problem, geometry and element description

5.1.1 Problem statement

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


pict

Figure 5.1: EMA project problem description

The problem described above was solved for the following values of the angle α  (in degrees)

|---------------------|
|                     |
--{0,15,30,45,50,55-}--|

Where α  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 × 2  integration rule). After the global stiffness matrix was assembled from the two elements stiffness matrices, the system equations given by KD  =  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).


pict

Figure 5.2: Global node coordinates using general coordinates, showing the distortion angle

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 α  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.


pict

Figure 5.3: Global and element node numbering used for project EMA

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










Table 5.1: elem_map_nodes table. Mapping element node to global nodes

There is also a need to lookup the global coordinates {x ,y }
   i  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 Δ  = tan θH2   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 α  . When α  is zero, there will be no distortion of the elements.





global node # global node coordinates xi  yi



1 0  H



2 L-
4   H



3 L-
2 + Δ  H



4 3
4L  H



5 L  H



6 0  H-
2



7 L
2-   H
2-



8 L  H2-



9 0  0



10 L
4-   0



11 L2 − Δ  0



12 3L
4  0



13 L  0




Table 5.2: global_coordinates_tbl. Mapping global node number to global coordinates

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


















Table 5.3: elem_map_dof table. Mapping local element DOF to global stiffness matrix locations

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 KD  =  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  .


pict

Figure 5.4: global DOF numbering used

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.


pict

Figure 5.5: 8-node element used for the finite elements

The following are the 8 shape functions used

pict

The global coordinates x,y  are now expressed as functions of the natural coordinates ξ,η  using

pict

Where M  is the number of nodes of each element (which is 8  ) and xi,yi  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, xki  means the global x  coordinate of the ith   node in the kth   element.

Similarly,  k
yi  means the global y  coordinate of the  th
i   node in the  th
k   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.

pict

5.2.2 Finding the Jacobian

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

pict

Each of the above derivatives is evaluated and used in the Matlab code.

pict

The above is now used to find the Jacobian and its determinant and also find Γ  and the matrix B2   . To find the matrix B3   the derivatives of each shape function is taken w.r.t. ξ  and η  .

This below gives the result of this computation

pict

With the above, the matrix B3   matrix was calculated giving

pict

And calculate the element stiffness matrix given by

pict

The elements stiffness matrices kelem   are then combined to make the global stiffness matrix K  and then the system KD   = 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 −  20N  at the middle of the right edge of the beam.

In the above, the matrix B1   is

pict

And the matrix B2   is

pict

Where Γ  is the inverse of the Jacobian matrix Γ =  J−1   . And the matrix B3   is

pict

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.


pict

Figure 5.6: Global stiffness matrix spy() output showing the bands

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 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 (ξ,η)  , which extends from −  1...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 ξ =  1√--
      3   the value of r  is one. And when η = √1-
      3   then s = 1  also.

The following diagram shows this layout more clearly.


pict

Figure 5.7: Stress recovery using r,s  coordinates system inside ξ,η

Therefore, the relation between (r,s)  coordinates system and (ξ,η )  coordinates system is as follows

pict

The following diagram shows the mapping at two different points for illustration


pict

Figure 5.8: Stress recovery using r,s  coordinates system inside ξ,η

Now that the mapping between ξ,η  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 Ni  shape functions described below, and not the original (ξ,η)  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

      ∑4
σp =     Ni σi
      i=1

Where σ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, Ni  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

pict

For example, to find the stress at point (ξ, η) = (0, − 1 )  , the first step is to determine this point's (rp,sp)  coordinates. Since      √ --
rp =   3ξ  then rp = 0  and since      √ --
sp =   3η  then        √ --
sp = −   3 =  . Applying the extrapolation formula above gives

pict

Since the stresses at four Gaussian points σ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 σp  than using the direct method to find σ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.


pict

Figure 5.9: Location of points used for stress recovery in the r,s  coordinate system

5.3 Results

The results are listed by the angle θ  used to distort the elements with. For each angle, the deflection and σ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.

5.3.1 No element distortion. zero angle

summary of result





x  (meter) y  (meter)
Matlab −  0.15  − 1.02895
ANSYS −  0.15  − 1.0289




Table 5.4: Short summary of test case zero degree distortion

Matlab result





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




Table 5.5: Matlab result. nodal solutions, angle [0  ] degree

The following figure shows the deformation found


pict

Figure 5.10: deflection found using 2 elements using Matlab, zero degree

The following table shows Matlab result for the direct stress σx  found at each node for each element.






global node # x  y  σx  N/m^2




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





Table 5.6: Matlab result. direct stress σx  at each node, First element, angle [0  ] degree




global node # x  y  σx  N/m^2




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





Table 5.7: Matlab result. direct stress at each node, Second element, angle [0  ] degree

The following shows the direct stress contour generated in Matlab


pict

Figure 5.11: Contour of direct stress found using 2 elements using Matlab, zero degree

ANSYS result

The following figure shows the deformation found


pict

Figure 5.12: deflection found using 2 elements using ANSYS, zero degree

The following table shows ANSYS result for the direct stress σx  found at each node for each element.

The following shows the direct stress contour generated in ANSYS


pict

Figure 5.13: Contour of direct stress found using 2 elements using ANSYS, zero degree

5.3.2 15 degrees distortion

summary of result





x  (meter) y  (meter)
Matlab −  0.150262  − 1.02555
ANSYS −  0.15026  − 1.0256




Table 5.8: Short summary of test case 15 degrees distortion

Matlab result





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




Table 5.9: Matlab result. nodal solutions, angle [15  ] degree

The following figure shows the deformation found


pict

Figure 5.14: deflection found using 2 elements using Matlab, 15 degrees

The following table shows Matlab result for the direct stress σx  found at each node for each element.






global node # x  y  σx  N/m^2




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





Table 5.10: Matlab result. direct stress σx  at each node, First element, angle [15  ] degree




global node # x  y  σx  N/m^2




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





Table 5.11: Matlab result. direct stress at each node, Second element, angle [15  ] degree

The following shows the direct stress contour generated in Matlab


pict

Figure 5.15: Contour of direct stress found using 2 elements using Matlab, 15 degrees

ANSYS result

The following figure shows the deformation found


pict

Figure 5.16: deflection found using 2 elements using ANSYS, 15 degrees

The following table shows ANSYS result for the direct stress σx  found at each node for each element.

The following shows the direct stress contour generated in ANSYS


pict

Figure 5.17: Contour of direct stress found using 2 elements using ANSYS, 15 degrees

5.3.3 30 degrees distortion

summary of result





x  (meter) y  (meter)
Matlab −  0.146118  − 0.998214
ANSYS −  0.14612  − 0.99821




Table 5.12: Short summary of test case 30 degrees distortion

Matlab result





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




Table 5.13: Matlab result. nodal solutions, angle [30  ] degree

The following figure shows the deformation found


pict

Figure 5.18: deflection found using 2 elements using Matlab, 30 degrees

The following table shows Matlab result for the direct stress σx  found at each node for each element.






global node # x  y  σx  N/m^2




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





Table 5.14: Matlab result. direct stress σx  at each node, First element, angle [30  ] degree




global node # x  y  σx  N/m^2




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





Table 5.15: Matlab result. direct stress at each node, Second element, angle [30  ] degree

The following shows the direct stress contour generated in Matlab


pict

Figure 5.19: Contour of direct stress found using 2 elements using Matlab, 30 degrees

ANSYS result

The following figure shows the deformation found


pict

Figure 5.20: deflection found using 2 elements using ANSYS, 30 degrees

The following table shows ANSYS result for the direct stress σx  found at each node for each element.

The following shows the direct stress contour generated in ANSYS


pict

Figure 5.21: Contour of direct stress found using 2 elements using ANSYS, 30 degrees

5.3.4 45 degrees distortion

summary of result





x  (meter) y  (meter)
Matlab −  0.135417  − 0.934503
ANSYS −  0.13542  − 0.93450




Table 5.16: Short summary of test case 45 degrees distortion

Matlab result





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




Table 5.17: Matlab result. nodal solutions, angle [45  ] degree

The following figure shows the deformation found


pict

Figure 5.22: deflection found using 2 elements using Matlab, 45 degrees

The following table shows Matlab result for the direct stress σx  found at each node for each element.






global node # x  y  σx  N/m^2




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





Table 5.18: Matlab result. direct stress σx  at each node, First element, angle [45  ] degree




global node # x  y  σx  N/m^2




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





Table 5.19: Matlab result. direct stress at each node, Second element, angle [45  ] degree

The following shows the direct stress contour generated in Matlab


pict

Figure 5.23: Contour of direct stress found using 2 elements using Matlab, 45 degrees

ANSYS result

The following figure shows the deformation found


pict

Figure 5.24: deflection found using 2 elements using ANSYS, 45 degrees

The following table shows ANSYS result for the direct stress σx  found at each node for each element.

The following shows the direct stress contour generated in ANSYS


pict

Figure 5.25: Contour of direct stress found using 2 elements using ANSYS, 45 degrees

5.3.5 50 degrees distortion

summary of result





x  (meter) y  (meter)
Matlab −  0.129803  − 0.901557
ANSYS −  0.12980  − 0.90156




Table 5.20: Short summary of test case 50 degrees distortion

Matlab result





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




Table 5.21: Matlab result. nodal solutions, angle [50  ] degree

The following figure shows the deformation found


pict

Figure 5.26: deflection found using 2 elements using Matlab, 50 degrees

The following table shows Matlab result for the direct stress σx  found at each node for each element.






global node # x  y  σx  N/m^2




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





Table 5.22: Matlab result. direct stress σx  at each node, First element, angle [50  ] degree




global node # x  y  σx  N/m^2




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





Table 5.23: Matlab result. direct stress at each node, Second element, angle [50  ] degree

The following shows the direct stress contour generated in Matlab


pict

Figure 5.27: Contour of direct stress found using 2 elements using Matlab, 50 degrees

ANSYS result

The following figure shows the deformation found


pict

Figure 5.28: deflection found using 2 elements using ANSYS, 50 degrees

The following table shows ANSYS result for the direct stress σx  found at each node for each element.

The following shows the direct stress contour generated in ANSYS


pict

Figure 5.29: Contour of direct stress found using 2 elements using ANSYS, 50 degrees

5.3.6 55 degrees distortion

summary of result





x  (meter) y  (meter)
Matlab −  0.123001  − 0.86157
ANSYS −  0.123  − 0.86157




Table 5.24: Short summary of test case 55 degrees distortion

Matlab result





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




Table 5.25: Matlab result. nodal solutions, angle [55  ] degree

The following figure shows the deformation found


pict

Figure 5.30: deflection found using 2 elements using Matlab, 55 degrees

The following table shows Matlab result for the direct stress σx  found at each node for each element.






global node # x  y  σx  N/m^2




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





Table 5.26: Matlab result. direct stress σx  at each node, First element, angle [55  ] degree




global node # x  y  σx  N/m^2




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





Table 5.27: Matlab result. direct stress at each node, Second element, angle [55  ] degree

The following shows the direct stress contour generated in Matlab


pict

Figure 5.31: Contour of direct stress found using 2 elements using Matlab, 55 degrees

ANSYS result

The following figures show the deformation found


pict

Figure 5.32: deflection found using 2 elements using ANSYS, 55 degrees

The following table shows ANSYS result for the direct stress σx  found at each node for each element.

The following shows the direct stress contour generated in ANSYS


pict

Figure 5.33: Contour of direct stress found using 2 elements using ANSYS, 55 degrees

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 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 350   based on the diagram generated below. At angle   0
55   , 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.


pict

Figure 5.34: Error as function of amount of of angular element distortion

From the above one can see that to keep the error below 5%  , the angular distortion should not be more than about 350   as the element behaves very badly after that. ANSYS will not even accept the analysis when trying angle of    0
60   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

  1. direct method. In this method the stress is found at arbitrary points by directly evaluating σ = EBu  at the point. Where in this, the u  is the nodal solution. This requires finding the Jacobian at point in question. This method is simple and works very well as long as there is no distortion.

    Once distortion is added, it produced bad result in stress values when compared to ANSYS results.

  2. Method of extrapolation. This method is described in reference [1], pages 230-232. In this method, the stress at the nodes of the element (or at any other arbitrary point) is found by extrapolating the stress values found at the four Gaussian points. This works well because the stress at the Gaussian points is the most accurate since these are also the integration points used. This was the method used in this project and worked very well. Results from the Matlab implementation all agreed with ANSYS result for the direct stress at the nodes.
  3. Patch Recovery. This is based on using a polynomial of same order as the shape functions and then using such polynomial on a small patch around the point on interest to find the stress. It would be interesting to compare this method in the future with the extrapolation method to see which is more accurate or easier to implement.

When making the contour plots for σ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 550   . This shows clearly that stress across elements is not smooth and changes abruptly now.


pict

pict

Figure 5.35: Stress contour, at 55 degrees

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.


pict

pict

Figure 5.36: Stress contour, at zero degrees

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

5.4.1 References

1
Concepts And Applications Of Finite Element Analysis. 4th edition. Robert D. Cook, David S. Malkus, Michael E. Plesha, Robert J. Witt. John Wiley & Sons. Inc.
2
Effect Of Element Distortions On The Performance Of Isoparametric Elements. Nam-Sua Lee, Klaus-Jurgen Bathe. Dept of Mechanical Engineering. MIT. International Journal For Numerical Methods In Engineering. Vol 36, 3553-3676 (1993)
3
ANSYS help manuals for APDL and general ANSYS use.

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 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