Evaluation and interpretation of integrated shell forces from stress results of FE models

When conducting a certain structural design such as for bridges, one might need to model the geometry twice. One time in CAD (computer-aided design) programs and another in CAE (computer-aided engineering) programs, which requires a high modeling effort. This project aims to utilize the CAD model directly for analysis by exporting it as a volume model to a CAE program such as SOFiSTiK to perform a structural analysis.

SOFiSTiK is a powerful tool for evaluating stresses of 3D finite element volume models. When working with concrete slabs, internal forces are needed to determine the amount of reinforcement. Up to now, volume models in SOFiSTiK are not suitable for determining internal forces, but shell models are. For this purpose, a python tool was developed specifically for integrating the stresses from volume models into internal forces. These may then be applied to an equivalent shell model. Until now the python tool was only applicable for horizontal plates.

This work's purpose is to analyze and interpret the results from the python tool for different loading cases with the focus on temperature loads and extend the code to curved and arbitrarily oriented surfaces.

Introduction to the Python Tool

In this project, the focus is on the implementation in Python and SOFiSTiK. Nevertheless, a brief overview of the overall workflow of this project is worth mentioning. This section explains the work-flow and the different parts of the python code.

Work Flow of the project 

Step 1: Prepare Evaluation Surface Model of the BIM Volume Model of a Structure

In order to calculate the forces of the structure modeled in CAD software, one needs to create additionally the evaluation surface onto which the integrated shell forces will later be projected.

Evaluation Surfaces

Since having the results of internal forces on shell surfaces or beams is better for concrete design, the stresses from the volume model are integrated into shell forces on surface models called evaluation surfaces. The function of evaluation surfaces is to evaluate the internal force values, which are obtained by the volume model through sectional integration. The sectional planes, on which the stress integration is performed, are defined at the position of each node in the ES mesh and are perpendicular to the ES.

It is important to mention, that the evaluation surface is slightly smaller than the volume model because of the cutting sections. To obtain good calculation results, the cutting sections should not exceed the volume model. This is shown in the figures below. The left figure shows the volume model and its evaluation surface. The right figure represents a corner, illustrating the reason why the evaluation surface should be slightly smaller than the volume model. Furthermore it shows how the volume model is being cut by the sectional planes.

Step 2: Export of VM and ES into CAE

In this step, the VM and ES need to be exported to SOFiSTiK. To export from Siemens NX to SOFiSTiK, an external C# plugin was developed by the initiators of this work. In general, one can also use other means of export. This depends on used CAD and CAE software.

Step 3: Meshing and calculation

In SOFiSTiK the workflow is no different from a typical calculating procedure of the structure in CAE. Loads are applied to the VM, then the VM is meshed into 3D finite elements and a calculation is performed. The results of the calculation of the VM in SOFiSTiK are the stresses at the nodes and not forces, which is the main current restriction of VM in SOFiSTiK. Separately, the ES is meshed in SOFiSTiK into 2D finite (shell) elements. The figures show the shell model (left) of the respective structural elements of the VM (right) that are meant to be calculated with the Python code.

Step 4: Integration of Stresses

The ES mesh is the basis for the evaluation of internal forces, as the nodes of the mesh define the points of interest. The Python code then reads the different data from the database files and performs integrations at the nodes of the mesh. The results are written to the initial shell model to evaluate the amount of reinforcement.

For the calculation of internal forces the following information is read from SOFiSTiK:

1. Evaluation surface data: containing coordinates of integration points (nodes of the mesh).
2. Volume model data: includes coordinates and stresses at the nodes of tetrahedron elements for each load case.

The data is imported into Python from the database (.cdb) files from SOFiSTiK for both the volume model and the shell model.

Basic Work-Flow (Python)

Reading Data from Database Files

For reading data from the database files, SOFiSTiK created a library that may be imported in Python. This library allows us to read data from the database. As SOFiSTiK has a special architecture for storing data, we define appropriate classes in Python. The data is then stored in JSON files (separately for the volume and shell models). This makes the data easily accessible from the main Python script.

SOFiSTiK uses numbers to distinguish different groups for structural members and load cases, as can be seen in the figure below. The different colours correspond to different groups. In python, we may choose which groups and load cases to consider for the stress integration. 

This way we may assign only a few group numbers or specific load cases that should be analyzed.

Cutting Planes

Internal Forces

Step 5: Export and structural design

In the final step, the evaluated internal forces are exported to the shell model in SOFiSTiK. SOFiSTiK can perform a structural design of concrete slabs for shell models, where it can define the amount of reinforcement for concrete slabs. This may then be implemented in SOFiSTiK.

Evaluation and Interpretation of Results 

To judge the performance of the developed method, a comparison between the internal forces obtained from the volume model and the shell model of the bridge had been conducted. The analyzed load cases were thermal loads and earth pressure. Temperature loads are an important load case when considering bridge design. In general,  thermal analysis allows for evaluating the temperature distribution within a solid due to steady or transient heating or cooling or a temperature gradient. 

In the previous work done by Lu’s Thesis [2], it had been realized that the biggest difference between the results obtained from the volume model and the shell model was in regard to the temperature loads. That was possibly due to few reasons such as not having exact matching geometries between the shell and volume model, also the different methods utilized to calculate the results. The process of obtaining the results in the volume model, consist of the following: 

First, the results are extrapolated from the gauss points to the nodes, which are then averaged at the nodes, and finally, they are integrated on the cutting plane. This is illustrated in the figures below, in which the figure on the right just shows the averaging of extrapolated gauss point values at a node, and the figure on the left shows the cutting plane through the element, which represents the area of integration for this particular element and cutting plane. 

The process described above has a high possibility of yielding computational errors since until reaching the final result, there is a solid amount of data processing  (e.g. extrapolation, averaging, integration). This is of course not only limited to the temperature load cases but to every calculation that is being done using the python-tool. However, since the biggest differences appeared in the thermal load cases, they are more thoroughly discussed.  

Types of temperature load cases applied

These 4 different temperature load cases that are listed below had been applied inside SOFISTIK. Previously, these load cases had been applied on all the surfaces of the volume model where the thermal load acts. This had been modified, to having the loads defined on the group of the corresponding part of the structure. As a result, the load is applied to all the nodes of the quadrilateral elements of that particular group. This change in the definition of the application of the thermal load resulted in better performance of the volume model.

• Uniform heating +29 degrees.
• Uniform cooling down -26 degrees.
• Difference between top and bottom of 12.3 degrees. Here, the top is being the warmer surface.
• Difference between top and bottom of 8 degrees. Here, the bottom is being the warmer surface

Load application

When making a structural model and applying the load, one must think about how the structure will behave under a certain action. In the case of this project, the geometry was precisely defined for the volume model and simplified for the shell model. This could result in many possible discrepancies. Therefore it is important to be aware of the limitations of both methods and either try to reduce the discrepancies by making adjustments to the modeling method or understand how to interpret the results.

The way of applying the temperature load differs between the shell and the volume model. In the first case, the load is applied as a temperature gradient to the surface, where the sign of the applied load shows the direction of the gradient. In the case of the volume model, the load is applied separately on the surrounding surfaces using the SOFISTIK command hydra. It means that it should be manually defined which surface is warmer and which is cooler and define the specific temperature values in a way that the difference results in the same value as the temperature gradient in the shell model. Therefore it is important to define the local axes of the element very precisely.

Element thickness and the change of it also play an important role in the result interpretation. Since the stresses from the temperature appear because of the restricted strains, the results will vary for the elements with different stiffness.

In this project example, the bridge slab has a thickness change on the free edges of the span. In the volume model, the geometry is defined as it is, very precise, while in the shell model the slab edges are not modeled. For the uniform temperature load, the restricted strains will result in tension or compression forces and mostly uniform bending moment in the main span direction. Therefore the thickness change for these load cases most likely won't have a big influence on the results. However, it is not the case for the varying temperature load, because this action generates the bending moment.

Results 

Thermal Loads

Testing the method

First, the temperature load case for the uniform heating of +29 degrees had been discussed. Three different figures which are shown below, correspond to the shell model (right), volume model figure (left), volume model with mesh refinement at the bridge deck (middle).

The bending moment in the x-direction (span direction) had been analyzed. As it can be seen for the volume model without the mesh refinement, there exist big fluctuations of the bending moments on the bridge deck. This is simply because the existing mesh was not enough to capture the behavior of the structure properly. For the volume model with refinement at the bridge deck, it can be seen that these fluctuations had been reduced greatly after the mesh refinement.
Comparing the results of the shell (right figure) and volume model (middle figure), it can be seen that the results of the volume model are quite comparable to the shell model and shows a reasonable correlation. The results of the evaluation face still show some inconsistencies, which could be solved by refining the mesh even more. However, it is important to remember that every refinement is resulting in increased calculation cost. That is why this has always been a crucial task of an engineer to find a compromise between the level of certainty and the calculation effort.

Some other differences can be seen in the results, for example, the higher bending moment in the corners of the slab in the shell model, which is not traced in the volume model. For evaluation, it is important to acknowledge that the evaluation face does not have the edge results, because it is smaller than the real structure, as it was described in the previous chapters. Therefore the local effect of the corners can be excluded from the resulting domain.

Summary of the results

In the following figures, only the results for the refined mesh in the volume model are shown. The figures below represent the comparison of the results for the bending moment in the X-direction between the volume model and the shell model. It can be seen that the values on the inner part of the surface have very similar values. As it was discussed before, the corner effects of the shell model shouldn't be taken into account, because they cannot be properly compared.

 Next, the bending moment in the Y-direction for the same load cases is analyzed. First, two load cases (LC70 and LC71) represent the uniform temperature load. Here it is clearly seen that the Y-direction is not the main because the bending moment is much lower than in the X-direction (main span direction). Comparing the results between the volume model and the shell model, the higher difference can be traced. In the middle of the span, the results on the evaluation surface are on average 40% higher (absolute value) than on the shell model. This could be due to the differences in the modeled geometry, like the larger width of the slab which results in higher absolute moment values. In that case, the results of the volume model show more conservative results.

The load cases LC 72 and LC 73 have a higher chance of the bending moment in the Y-direction for the volume model in comparison to the shell model. The results in the middle of the span show the absolute difference of 25% in relation to the shell model. Having a careful look at the edge results of the evaluation face one can see that the absolute values are slightly lower for the volume model comparing to the shell model. This also could be the reason for the difference in geometry. As it was discussed before the load application technique was adjusted in order to get a more similar deformation pattern in the structure. However, the local effects and strains in the area where the slab has a thickness change were not considered during the modeling. Therefore the results could be affected by that. 

As it was discussed in this chapter, the modeled geometry has a big impact on the stress results. By modeling the structure one should be aware of the assumptions made and the limitations of the modeling method.

Earth Pressure

The main earth pressures are active earth pressure and in-rest earth pressure, but the one from backfill under traffic actions is the most unfavorable.

The results of the internal force of the earth pressure show that the internal forces of the volume model exceed its counterpart in the face model. These discrepancies can be explained in a similar manner to what had been described above in relation to the calculation method. 

Structural Design

The previous analysis showed that the temperature load is not causing a big difference for the bending moment in X-direction (main span direction) and a small (rather conservative) difference for the bending moment in Y-direction. Since in the previous work, described by Lu Gup in his master thesis, the temperature load cases were mentioned as the ones having the higher discrepancy in the results, other load cases were not further analyzed.

The next step is to perform the design load combination and get the results for the amount of reinforcement necessary. The figures below show the results for the upper and lower reinforcement of the volume and the shell model.

Analyzing the results for the amount of the reinforcement, a bigger difference between volume and shell model can be seen. The values on the evaluation face show smoother results, while the shell model shows the higher peak value next to the supports. It could also be the reason for the discrepancy lacked information about local stresses next to the supports since the evaluation face cannot capture the results on the edges and along with the supports.

For the lower reinforcement the results also visibly differ. Interesting is that the in volume model shows the same reinforcement amount necessary for both directions, while in the shell model the main direction (X-direction) can be clearly tracked. The reason for such discrepancy in the results could be caused by the geometry, load application technique, or specific load case. Some deeper analysis is necessary to find what causes the differences in the reinforcement calculation.

Extension of the Python Code

Up to this point, the python tool was able to perform the computation of internal forces for horizontal plates only. In other words, the orthogonal direction of the mid-plane was assumed to be parallel to the global z-direction. Also, the plate was assumed flat (no curvature in the geometry).

This section will explain the extension of the code to treat curved surfaces with arbitrary orientation.

A detailed description of how to use and run the code is available on the GitLab repository under the following link:

https://gitlab.lrz.de/ge75yov/mssm-project.git

Extension to arbitrarily oriented plates

The algorithm is assuming the local z-direction to coincide with the global z-direction. To adapt the code for arbitrary local z-direction, we use simple continuum mechanics operations. 

The code runs a small unit where cutting planes are defined for every point in the mesh of the ES. In this section, the local directions of the cutting planes are defined. These are gathered in the so-called coefficient matrix which is passed to the central part computing the internal forces. For the evaluation of internal forces, we fall back on continuum mechanics calculus. Normal and shear stresses are evaluated by using the stress vector and local directions of the cutting plane (from the coefficient matrix). 

With $//<![CDATA[ \begin{array}{l}\sigma\end{array} //]]>$ being the physical stress in direction of the normal $//<![CDATA[ \begin{array}{l}n\end{array} //]]>$ , $//<![CDATA[ \begin{array}{l}t\end{array} //]]>$ the stress vector, $//<![CDATA[ \begin{array}{l}\sigma^{ij}\end{array} //]]>$ the stress tensor and  $//<![CDATA[ \begin{array}{l}\tau\end{array} //]]>$ the shear stress in direction of the in-plane vector $//<![CDATA[ \begin{array}{l}p\end{array} //]]>$.

After gathering the different stresses, we may get the forces by multiplying with the area of the stress triangle and moments by additionally multiplying by the lever arm (see section Internal Forces in Introduction to the Python Tool).

Previously the lever arm was assumed parallel to the global z-direction. We now extend it to be in the local z-direction of the cutting plane. This is performed by using the information from the coefficient matrix for the local z-direction. Now we may project the distance between the mid-plane and the center of gravity of the stress triangle in the local directions. 

Inclined rectangular plate example

To benchmark the results of python code, the plates were additionally modeled and calculated in SOFISTIK as 2D plates meshed with QUAD elements. These models are called Shell Models (SM). In the picture below one can compare the bending moments in local direction of three plates evaluated with the help of python code and compare them to the local moments of SM.

The results can also be compared to a reference value for the moments at the center of the plate taken from the original Master's thesis by Max Bronner, which are: m_xx = 5.75 [kN/m] and m_yy = 8.11 [kN/m].

Influence of refinement of the Volume Model

 Additionally, we compared the results of Python force integration on the Volume models calculated with the different mesh sizes. For the first model H_min = 0.25 [m], for the second one H_min = 0.125 [m], where H_min is the maximal length of the mesh volume elements. The results can be seen in the figure below. As one can see, the results converge to the reference values from the Master's thesis by Max Bronner.

Extension to curved surfaces

Previously the cutting planes have been determined by determining a local coordinate system based on the global geometry of the plate (i.e. the corner nodes of the plate were used). This approach will be omitted and replaced by a more general approach. 

To avoid determining the local coordinate system of every point of interest, we may use information stored in the SOFiSTiK database file of the ES. The database contains information about local coordinates of every node based on the averaged local coordinates of the adjoining elements. We added the collection of this information in the appropriate script and exported the information to the corresponding JSON file. Back in the main script we may now read the information from the JSON file and assign the local coordinate systems to the respective nodes to form the cutting planes. The rest of the process for calculating internal forces remains unchanged. This way we may define arbitrarily shaped surfaces and compute internal forces at the nodes of the mesh of the corresponding ES.

Curved surface example

In the figure below the comparison of the results is shown for the SM and the VM, calculated with the modified python code.

Arch 2

As the Evaluation Surface and Shell Model are meshed in SOFISTIK separately, it is difficult to compare resulting forces on two models. For better comparison, an additional arch plate was modeled in Rhino7 CAD software. The Finite Element meshes for ES and SM were created in Rhino and exported to SOFISTIK. It is important to mention that mesh size was chosen to be of the same size as half of the length of integrating sectional planes. This way, just by deleting one row of the Finite Elements on four edges of the arch we could create the mesh of ES. Meshing is shown in the left figure below. At the same time, VM was meshed automatically in SOFISTIK.

Plate arch parameters are:

• Span: 10 [m]
• Width: 5 [m]
• Height: 2 [m]
• Thickness: 0.4 [m]
• Loading: 20 [kN/m2]

In the figure on the middle deformed structure is shown. In the right figure below section planes that were used for integration of the shell forces are shown.

In the figure below the comparison of the results is shown for the SM and the VM, calculated with the modified python code.

References

[1] Max Bronner. Verfahren zur berechnung der schalenschnittgroeßen eines fevolumenmodells im brueckenbau. Master's thesis, TUM, 2018.

[2] Guo, Lu. Method For Reinforced Concrete Design In Bridge Construction Based On FE Volume Model. Master's thesis, TUM 2020.