Stress Singularities

In my last entry, I wrote about how to check whether the mesh was of sufficient quality to ensure that the mesh would not cause divergence in the solver. However there are situations which arise where regardless of how many elements you have in the model, or how perfect the mesh is, the stress results will be inaccurate at best and divergent at worst. These locations can arise due to what are known as stress singularities.

There are a couple reasons where stress singularities occur. The first place would be at the location of bad elements. If there are elements in the mesh which have very high aspect ratio or Jacobian ratios, the stress value at these elements can incorrectly report as much higher than they actually are. This will cause a spike in the stress values of the surrounding elements as well because the stress value is averaged over an area. However if you followed the steps outlined in my last post, this should not be an issue.

A second place where stress singularities can occur is in regions where the physical geometry has a sharp corner which carries a significant loading. These areas are referred to as “re-entrant” corners, and according to the theory of elasticity, the stress at these points approaches infinity. In simpler terms, stress is equal to force/area.

σ = F/A

When the area tends towards zero, the stress tends towards infinity. Area at points or edges is very small, so the stress value will be very high.

A third place where stress singularities can occur is at fixtures and boundary conditions. Fixtures have an infinite stiffness value whereas the material you are analyzing does not. This can also result in an infinite stress values near the fixture. When analyzing a model, it is very important to take into account how the model will be restrained in the real world and whether the stress in areas near your fixtures are of interest. This particular problem will be addressed in my next post.

We will begin with a study of stress at sharp corners:

Figure 1: A 90 degree bracket with a sufficiently refined mesh to provide reasonably accurate results. There is a fixture at the top restricting all 3 degrees of freedom for each node (solid tetrahedral elements in SolidWorks Simulation have only 3 translational degrees of freedom each), as well as a load that is normal to the fixed face.

Figure 2: Running the study, we find the reported stress value is a maximum at 91.2 MPa

Figure 3: The mesh superimposed on the results.

Figure 4: The mesh is further refined

Figure 5: The maximum reported stress at the corner is now 103.2 MPa.

Figure 6

The difference between the mesh in figure 1 and figure 4 is that the global mesh density has been increased. While in general this should increase the accuracy of the stress and strain results, this also comes at a significant cost to solver computational time. Reducing the element mesh size by a factor of two can increase the number of elements by a factor of 8, and each high quality element has 10 nodes, each with 3 degrees of freedom. While many of the nodes are shared, you can see how this can substantially increase the solver time as each degree of freedom requires its own equation to define its displacement.

Instead, there is an option to refine the mesh locally, and is available for points, edges, faces, and components. This is referred to as a “Mesh Control” in SolidWorks Simulation. To apply a mesh control, right click on the mesh icon > “Apply Mesh Control”.

Figure 7: Right click on the “Mesh” icon > “Apply Mesh Control

Figure 8: The re-entrant edge is selected

Figure 9: The first case has a mesh control of 0.75 mm applied to the edge. You can see that the mesh is refined substantially compared to the surrounding global mesh element size.

Figure 10: Running the analysis, the maximum stress at the edge reports as 140.4 MPa.

Figure 11

Figure 12: The mesh control is refined to 0.5mm

Figure 13: Stress is now 167.9 MPa

Figure 14

Figure 15: 0.25mm Mesh Control

Figure 16: 282.9 MPa

Figure 17

Figure 18: 0.1mm Mesh Control

Figure 19: 425.7MPa

Figure 20

At this point, you can see that the results are divergent. Every time the mesh is refined, the stress values increase significantly for the same loading value. The results at this edge are in-accurate and most likely divergent.

We can resolve this issue by ignoring the results in this corner, however in analyzing this right angle bracket, this may be an area of interest as this is where the most stress will occur given the fixture and loading conditions. The second option is to physically modify the geometry such that the area of the edge or corner is not very small or approaching zero at that corner.

The simplest way to modify the geometry such that the area is not approaching zero at any edge is to add a tangent fillet.

Figure 21: A small fillet is added to the corner

Figure 22: The global mesh is set to 2mm, and no mesh refinement (mesh control) is applied at the fillet

Figure 23: The stress value reports as 214.1 MPa

Figure 24

Figure 25: A mesh control of 0.75mm is added to the fillet face

Figure 26: Maximum reported stress on the face is 228.1 MPa

Figure 27

Figure 28: Mesh control reduced to 0.5mm

Figure 29: Stress reported as 259.9 MPa

Figure 30

Figure 31: 0.25mm Mesh Control

Figure 32: Stress reported as 270.8 MPa

Figure 33

Figure 34: 0.1mm Mesh Control

Figure 35: 283.4 MPa Stress

Figure 36

In the second example where the fillet was added to the corner, the stress value continued to climb slowly, however this is expected. Although the stress results increase as the element size decreases, the % rate of change in reported stress also begins to decrease from iteration to iteration. This is considered a converging result.

In contrast, the first example without the fillet exhibited behavior such that the % difference in reported stress values continued to climb with each successive mesh refinement, indicating that the results were divergent.

Figure 37: Reported stress values versus element size for the bracket without the fillet

Figure 38: Reported stress values versus element size for the bracket with the fillet. As the mesh is successively refined, the % increase in reported stress values is also decreasing with each iteration indicating convergence.

In my last entry, I wrote about how to obtain accurate stress results by ensuring that there were at least two elements across the thickness in areas of interest. However to obtain accurate results, the quality of the mesh itself must be taken into account.

In order to talk about mesh quality, there are a couple important metrics that need to be taken into consideration.

The aspect ratio of an ideal tetrahedral element is 1.0. This is a ratio of the longest edge to the shortest normal dropped from a vertex to the opposite face, normalized with respect to the shortest normal dropped from a vertex to the opposite face of a perfect tetrahedral element. A general rule of thumb is to not have more than 10% of the elements with an aspect ratio higher than 10. Extremely large values >> 40 should be closely examined to determine where they exist and whether the stress results in those areas are of interest or not.

Figure 1: This is an example of a first order tetrahedral solid element with an aspect ratio of 1.0. The aspect ratio is the ratio of the longest edge to the shortest normal dropped from a vertex to the opposite surface normalized with respect to a perfect tetrahedral element. In this case, the shortest normal has a value of 1.0 because it is normalized with respect to itself.

Figure 2: On the left is a picture of a tetrahedral element with an aspect ratio of 1.0. The element on the left has a much larger aspect ratio.

In SolidWorks Simulation, you can create a mesh plot showing the values of the aspect ratio within the mesh after it has been generated.

Figure 3: Right click on the Mesh icon in the Simulation tree> “ Create Mesh Plot”

Figure 4: You have the choice to create a mesh plot, an aspect ratio plot, or a jacobian plot.

Figure 5: We see that the maximum aspect ratio in this mesh is 3.79, well within acceptable limits

The second metric used to determine mesh quality is the Jacobian Ratio. This method is only available for second order (High quality) mesh elements. Parabolic (second order) elements are able to map curvilinear geometry more accurately than the first order linear elements. The midside nodes are placed on the actual geometry of the model, and in extremely sharp or curved boundaries, the edges can cross over each other. This can result in a negative jacobian ratio which will cause the solver to fail.

Figure 6: The Jacobian is a measure of the curvature of the edge and distortion at the mid-side node.

Figure 7: This is an example of a 2-D representation of an element with a Negative Jacobian ratio. The curvature of the geometry that the element was trying to map was too great for the size of the element, causing the edge to collapse in on itself creating a negative jacobain ratio. This will cause the solver to fail.

As with the first order elements and the aspect ratio, the jacobian ratio of a perfect tetrahedral element with linear edges is 1.0. The jacobian ratio of an element increases as the curvature of the edges increase and are calculated at the selected number of Gaussian Points for each tetrahedral element. In general, elements with a jacobian ratio less than 40 are acceptable.

You can also create a Mesh Check Plot similar to the Aspect ratio check in SolidWorks.

Figure 8: The Jacobian check plot shows that the only areas where there are elements with a non-unity Jacobian value are areas where there is actual curvature of the geometry. Here the maximum value of the Jacobian is 1.095, and there is no need to further refine the mesh.

Mesh Density

When performing an FEA analysis, the goal is generally to obtain results which are accurate while taking the least amount of time. The amount of time that a study takes to solve is directly related to the number of degrees of freedom in the study. With second order solid tetrahedral elements there are 10 nodes per element, each with three translational degrees of freedom. Each degree of freedom has an associated equation which must be solved for displacement. For second order triangular 2-D shell elements, there are 6 nodes, and each node has 6 degrees of freedom, three translational, and three rotational.

The first unknown in linear static FEA that is solved is the displacement of each node. Based on this, the stresses and strains are calculated. The stress values at nodes in an FEA study are calculated at Gauss, or Quadrature points in the element, and then averaged with the stress values from the surrounding elements (http://pathfinder.scar.utoronto.ca/~dyer/csca57/book_P/node44.html). While the displacements are solved explicitly at the nodes, the stresses are an averaged value and if there are insufficient stress values present in an area, the stress value averaged at the node can be inaccurate.

How many elements are necessary to obtain an “accurate” result? This depends on several things, such as what degree of accuracy you need as well as the geometry itself. (More on the geometry in my next post). However I will say that stresses near sharp corners or boundary conditions (fixtures, connectors, etc) can be inaccurate. The recommended starting point for “accurate” results according to SolidWorks is two second order tetrahedral elements across the thickness in all directions. However this is just a starting point, and you will have to decide for yourself whether the results obtained are sufficiently refined for your needs. This also means that in areas where you are not interested in the stress results that you do not necessarily need two elements across the thickness, reducing computation time.

Let’s examine the relation of number of elements across the thickness versus stress:

Figure 1: A thin plate with a hole in it. There is a “sensor” placed at the location where the highest stress in the plate will occur

Figure 2: Standard fixtures and a tensile load are applied to the plate. The global element size is set to a value which is larger than the thickness of the plate, resulting in just one element across the thickness and skewed elements (aspect ratio>1).

Figure 3: If you examine the cross section of the plate, you can see that there is only 1 “triangular edge” along the thickness.

Figure 4

Figure 5: The plate is re-meshed with a target element size of 1.50mm. Although the thickness is 2.5mm, because the elements are not forced to keep an aspect ratio of 1 (An aspect ratio of 1 results in a perfect equilateral triangular tetrahedron), there are two elements across the thickness.

Figure 6

Figure 7: checking the results at the sensor point, we find that the stress value has increased by almost 10%

Figure 8: The mesh is further refined to 0.75mm

Figure 9: Three elements across the thickness

Figure 10: Increase in stress value by halving the mesh size results in a 0.7% increase in reported stress values.

Thickness of Plate = 2.5mm | Element Size | Number of elements across thickness | Von Mises Stress (Mpa) | % Increase in Stress |

Case # | ||||

Case 1 | 4.372 mm | 1 | 395.8 | N/A |

Case 2 | 1.5mm | 2 | 432.5 | 9.27 |

Case 3 | 0.75mm | 3 | 435.6 | 0.72 |

As we can see, the difference between a single element across the thickness and two elements across the thickness is almost 10%. However further refining the mesh to include three elements across the thickness results in only a 0.72% increase in reported stress value.

In order to calculate displacements, stresses, strains, or any other results of interest (modal frequencies, temperatures, etc) the FEA process can be used to obtain desired results. SolidWorks linear static FEA can be conducted under the following assumptions:

1. The material is linear

2. Structural deformations are small, and there is no rigid body motion

3. The applied loads are static

The basic FEA process:

Preprocessing-

The model is split up into a series of discrete (or finite) elements for which a matrix of equations can be created. SolidWorks simulation takes the CAD geometry and will automatically create a mesh for you. In most cases it is necessary to manually adjust the mesh size and refine the mesh in areas of interest.

Solution:

Computation of the mathematical model. The solver runs after you have defined your material, fixtures, and loads. The solver constructs a system of equations from the elements based on these parameters and solves for them either directly or iteratively.

Postprocessing:

Analyzing the results. This is where you evaluate the results and decide whether they are sufficiently accurate for your needs.

Linear static FEA process begins with taking the geometry and discretizing it into a series of smaller elements. Currently only basic shapes have analytical solutions. CAD geometry is often complex and must be broken down, or discretized, into a series of continuous elements which can be solved for displacements and subsequently stresses and strains. There are many different types of elements that FEA programs use to handle different classes of problems. SolidWorks has three main types of elements, 3-D solid tetrahedral, 2-D triangular shell, and 2-D beam elements.

There are two types of solid tetrahedral elements in SolidWorks Simulation:

First order (Draft Quality) Tetrahedral: 4 nodes

Second order (High Quality) Tetrahedral: 10 nodes

This is a 3-D tetrahedral Second Order Solid Element. It features 10 nodes (4 vertice and 6 midside).

Image: http://web.mit.edu/calculix_v1.6/CalculiX/ccx_1.6/doc/ccx/img91.png

The tetrahedral high quality element has 3 translational degrees of freedom per node for a total of 30 degrees of freedom per individual element. Rotation is accounted for because the element is 3 dimensional so the 3 rotational degrees of freedom are omitted from each node to reduce computational time.

–

In order to simplify some classes of problems where the geometry has a constant thickness and a length to thickness ratio greater than 20:1, such as sheet metal, a second type of element which has fewer nodes is often used to reduce computational time. This type of element is called the Triangular Shell element.

There are two types of triangular shell elements available in SolidWorks Simulation:

First order (Draft Quality) Triangular Shell: 3 nodes

Second order (High Quality) Triangular Shell: 6 nodes

This is a 2-D triangular draft quality shell element. Unlike the high quality element, this has no mid-side nodes, for a total of 3 nodes.

The triangular shell elements have 6 degrees of freedom (3 translation and 3 rotation) for each node to account for rotation and moments applied to the element.

–

Lets compare the results from these various types of elements to see how accurately they represent geometry.

Figure 3: This geometry has been discretized using 3-D draft quality tetrahedral elements. Note that it does not capture curvature well due to the linear edges

Figure 4: The same geometry using the same global element size as in Figure 3. The second order elements allow for parabolic displacement of the edges and increase the ability to more accurately map the geometry

Figure 5: Using 2-D first order triangular shell elements to discretize the geometry. The global element size is the same as with the solid elements.

Figure 6: Using 2-D second order triangular shell elements to discretize the geometry. The global element size is held constant, however analogous to the difference between first order and second order tetrahedral elements, the second order shell elements are able to more accurately represent the geometry than first order shell elements.

In general, it can be helpful to run an analysis using the first order draft quality elements initially to ensure that your analysis has been setup properly. Due to the nature of many complex simulations, the solve time can often be substantial and running the analysis with a large global element size and draft quality elements can help determine whether your problem setup is correct prior to running the full refined study.

I hope to put up a series of instructional posts regarding the analysis tools available in SolidWorks. SolidWorks purchased the rights to COSMOSWorks as well as COSMOS FloWorks and has since integrated these analysis tools into their CAD software. There are some advantages to having an integrated CAD/analysis tool such as eliminating the necessity to convert your files to an intermediate file format when you want to perform analysis on your components. Changes are immediate and automatic, and you can also avoid problems with file conversion from intermediate formats which are unfortunately quite common between CAD packages.

I have been working as an applications engineer performing FEA and CFD analysis for the last two years at this point and wanted to provide some tips and tricks that can be used to help perform accurate analysis on your designs.

While it is not necessary to have a full understanding of the theory behind FEA and CFD to perform an analysis, it does help during problem setup to understand what you are telling the software to do. Remember, Garbage in, garbage out.

For some basic FEA theory, MIT has an excellent paper about the fundamentals here:

Tuning a PID system can be daunting, especially when you do not have a mathematical model to describe the system. What Ziegler-Nichols Tuning offers is a heuristic method that can be used to determine the PID values of a system for which you can perform tuning safely on.

The basic steps are:

Increase the P gain until the system becomes a constant amplitude stable oscillatory system, this value is referred to as Ku or the Ultimate Gain.

The oscillation period, Tu, is then found and the combination of these two values are used to set the P,I, and D values according to the following rules:

Ziegler–Nichols method[1]

Control Type Kp Ki Kd

classic PID[2] 0.65 * Ku 0.5 * Tu 0.12 * Tu

Using these values, the system should be close to critically damped. In my previous post, I posted the program that I have written to implement the PID loop on my inverted Pendulum. I foudn that the Ku, or ultimate gain value was around 57. I also included a couple lines of code which converted values to strings so that they could be exported via the serial connection to a datalogging application.

I have found a third party open source data logger called “Gobetweeno” which among other things, the capability to record the serial output from the Arduino and record it to a .csv file with a time stamp. I have used this data to find the oscillatory period of the system (Tu). Graphing this data in excel, I find that the Oscillatory period of the system is around .45 seconds.

Inputting this into the Ziegler-Nichols Tuning method with Ku=57, Tu = 0.45 seconds, the values are:

Control Type Kp Ki Kd

classic PID[2] 37.05 0.225 0.045

Since these values are non-integer, I will have to re-define the variables as double precision.

Now that I have a cart which is able to move under its own power, it’s time to implement PID control. After research, I

have found that PID control stands for Proportional, Integral, Derivative control.

As with simple control systems, the inverted pendulum has a Process Variable(the angle of the pendulum relative to the ground), a Set Point ( what you want your process variable to be, in this case 90 degrees relative to the horizontal), and a Manipulated Variable (the current applied to the motor which allows the process variable to be manipulated by).

When it comes to the PID portion, the proportional part means just that, the output (Manipulated Variable) will be directly proportional to the input (Process Variable). In the case of my inverted pendulum, this can be thought of as simply turning on the motor in the opposite direction when the pendulum starts to fall one way, thus forcing it back into the upright position.

However if you only use proportional control, the result can be at best an oscillatory system. At worst the system will be unstable divergent and will fly out of control.

If the rate of response from proportional output is not sufficient, the integral term can be added to increase the current applied to the motor over time to try and get to the Set Point more quickly.

However the combination of Proportional and Integral terms can result in overshoot as the integral term sums the errors and accelerates the response towards the Set Point. In order to prevent this, a damping term, the Derivative term can be added to the system. This will calculate the rate of change of the error and reduce the effect of the other two terms which is most noticeable as the system gets close to the Set Point.

By implementing all three, the system can achieve critical damping, which is generally a desired result. In this case, I wish to use all three terms to achieve this critical damping, however determining the values of these constants can be difficult especially if I have no mathematical model of the system.

To this end, I plan to implement a Heuristic tuning method known as the Ziegler-Nichols Tuning method.

Before I started writing the program, I again created a Flow Chart to help guide myself through how the program should flow.

Here is the first program I have written incorporation PID control in a function call so as to reduce the size of the program on chip.

#define INA 53 \\defining the variables

#define INB 51

#define PWM2 13

#define POT 0

int errorcw = 0; \\Initializing the variables

int errorccw = 0;

int correction = 0;

double kp = 57; \\Setting the Proportional Gain

double ki = 0; \\Setting the Intergral Gain

double kd = 0; \\Setting the Derivative Gain

int PID = 0;

int error = 0;

int last_error = 0;

int angle = 0;

int sp = 839; \\Defining the Set Point as 839 (90Degrees Horizontal)

void setup()

{

Serial.begin(9600); \\Start Serial output to record data

pinMode(INA, OUTPUT); \\Defining the Pins as output or input

pinMode(INB, OUTPUT);

pinMode(PWM2, OUTPUT);

pinMode(POT, INPUT);

}

void loop()

{

char buffer[5]; \\Defining a string output for gobetweeno data logging

angle = analogRead(POT);

Serial.print(“#S|PIDLOG|[“);

Serial.print(itoa((angle), buffer, 10));

Serial.println(“]#”); \\Defining the string output for Gobetweeno Data log

if(angle < sp) \\Clockwise Motor loop

{

errorcw = (sp-angle);

correction = calcPid(errorcw);

analogWrite(PWM2, correction);

digitalWrite(INA, LOW);

digitalWrite(INB, HIGH);

}

if(angle > sp) \\Counter Clock wise Motor Loop

{

errorccw = (angle-sp);

correction = calcPid(errorccw);

analogWrite(PWM2, correction);

digitalWrite(INA, HIGH);

digitalWrite(INB, LOW);

}

return;

}

int calcPid(int error) \\PID Fuction call to reduce size of program on disk

{

PID = (kp*(error)) + (ki*(error + last_error)) + (kd*(error – last_error));

last_error = error;

return constrain(PID, 0, 255); \\scale the output PID value between

\-255 for the motor output

}

Currently I have found the Proportional gain value that works best for my system to be around 57. Any more and the system can become divergent, and much less than this and the system does not respond fast enough to keep the pendulum upright.

The next step will be to proceed with Ziegler-Nichols Tuning so that I can incorporate Integral and Derivative contol into the pendulum so that I can achieve critical damping.