Solid mechanics simulations are mainly based on the finite element method that adopts the Lagrangian description—in which computational mesh is carved in the material body and deforms with the body. In post-processing, analysts naturally form their interpretation from results at Lagrangian mesh, for example, displacements and reaction forces at nodes, stresses and strains at integration points, nodes, or centroids of elements. Those results are attached to specific material points through Lagrangian mesh and are called results in Lagrangian representation. Investigating Lagrangian results satisfies analysts’ needs in solid mechanics problems most of time.

But in some dynamic or quasi-static problems in which solid material behaves like a flow, we may occasionally be more curious about results at some specific spatial locations rather than material points. For example, in a roll-to-roll process where a web travels along a designed path line for some functional processing from a unwinding roll to a rewinding roll, engineers may specifically be interested in the tension, traveling speed, temperature, or other physical quantities of web at some spatial locations such as some rollers since they may serve as parameters of control system or relate to product quality. Results of material expressed at spatial points are called results in Eulerian representation.

Abaqus has the ‘probe’ function in the visualization module which enables users to approximately get results at spatial locations through Lagrangian mesh in individual solution frame. However, it does not provide a convenient way of extracting the time history of results at spatial locations throughout the whole simulation. Yet, the time history of Eulerian results in these simulation is more valuable as it reveals when and how the steady state is achieved, which some times exposes physical information about the characteristics of the problem.

This project develops two general-purpose Abaqus Python scripts based on two methods that can automatically extract the time history of field outputs at some spatial location of interest. The developer wishes the scripts can save some labor and willpower of analysts so their energy can be reserved for pure problem solving rather than tedious external processing.

Physical Quantities in Lagrangian and Eulerian Representations

Essentially, physical quantities in solid mechanics are attached to material bodies to describe the state or property of material. Quantities of a material body can vary among its material points and change in time, thus mathematically, a physical quantity \(Q\) can be expressed as a function of material point \(p\) and time \(t\):

$$Q=f(p, t).$$

To distinguish each material point since there are infinite number of material points in a deformable body, we need a clever labeling system, so that when we talk about a quantity value, we know which specific material point we are referring to and when we want to express a quantity value for a specific material point we know how to refer it. Depending on whether the labeling system stays the same or changes along the course of deformation, we have the Lagrangian representation and Eulerian representation. Those are concepts from continuum mechanics, and I’ll summarize them as a brief reminder here. The summary is from a mechanical engineer’s point of view, with few sacrifice of mathematical rigorousness.

Lagrangian Representation

In Lagrangian description, the labeling system of material points remains the same during the course of deformation.

To set up the labeling system, we select a time instant and use the configuration of material body at that instant as a reference configuration. We build up a coordinates system in space at reference instant and label each material point by the coordinates it occupies at that instant. So each material point is mapped to a set of coordinates in the reference configuration by a mapping \(R\):

$$\boldsymbol{X}=(X_1, X_2, X_3)=R(p).$$

The inverse mapping \(R^{-1}\) maps the coordinates in reference configuration back to a material point \(p\):

$$p=R^{-1}(\boldsymbol{X})=R^{-1}(X_1, X_2, X_3).$$

When the material deforms in space and time, each material point probably moves to a new place. When we describe physical quantities at later instants, we stay with the label we’ve set up in the reference configuration. The label—the coordinates of material in reference configuration is called Lagrangian or material coordinates. Since we stay with this one labeling system during the whole course of deformation, the label of each material point remains same, or we say the Lagrangian coordinates deforms with the material body. Then, we can replace the argument of material point in the physical quantity function by the Lagrangian coordinates due to their one-to-one correspondence:

$$Q=f(p, t)=f(R^{-1}(\boldsymbol{X}), t)=\widehat{f}(\boldsymbol{X}, t).$$

As a result, the physical quantity is expressed as a function of the Lagrangian (material) coordinates and time, namely \(\widehat{f}\). This is the physical quantity in Lagrangian representation.

Eulerian Representation

In Eulerian description, the material point labeling system is constantly changing with time.

In Eulerian representation, we label each material point by the spatial coordinates it occupies at each instant \(t\), rather than at one reference instant in Lagrangian description. Since at different time instants, the spatial coordinates the material point occupies is probably different due to the deformation, the label of one material point at different instants is different. Therefore, the mapping between material point and spatial coordinates is time-dependent, namely \(\chi_{t}\):

$$\boldsymbol{x}=(x_1, x_2, x_3)=\chi_{t}(p).$$

The subscript \(t\) indicates the time dependence. This above mapping also has a physical significance—it actually describes the motion of material since it gives the location that each material point moves to at each time instant \(t\).

Similarly, the inverse mapping maps a set of spatial coordinates back to the material point that occupies that spatial coordinates, but at the specific time instant \(t\):

$$p=\chi_{t}^{-1}(\boldsymbol{x})= \chi_{t}^{-1}(x_1, x_2, x_3).$$

Again, we can replace the argument of material point in the function of a physical quantity:

$$Q=f(p, t)=f(\chi_{t}^{-1}(\boldsymbol{x}), t)=\bar{f}(\boldsymbol{x}, t).$$

Then the physical quantity is expressed as a function of the Eulerian (spatial) coordinates and time, namely \(\bar{f}\). This is the physical quantity in Eulerian representation.

A visualization of results in Lagrangian and Eulerian representations are shown below.

Transformation between Lagrangian and Eulerian Representations

If we replace the argument \(p\) in the labeling system of Eulerian description, i.e. the description of motion (\(\ref{motion}\)) with the inverse mapping adopted in the Lagrangian description, we’ll have the tool to transform the quantity expression between two representations. The mapping \(\chi_R\) transforms Lagrangian coordinates to Eulerian coordinates:

$$\boldsymbol{x}=\chi_{t}(p)=\chi_{t}(R^{-1}(\boldsymbol{X}))=\chi_{_{R}}(\boldsymbol{X}, t).\label{motion} \tag{1}$$

Noted, the indication of dependence on time \(t\) is moved from the subscript to an argument inside the function.

The inverse mapping \(\chi_R^{-1}\) transforms Eulerian coordinates to Lagrangian coordinates:

$$\boldsymbol{X}=\chi_{_{R}}^{-1}(\boldsymbol{x}, t).$$

Then, the transformations between a quantity in Lagrangian and Eulerian representations are:

$$Q=\widehat{f}(\boldsymbol{X}, t)=\widehat{f}(\chi_{_{R}}^{-1}(\boldsymbol{x},t),t)=\bar{f}(\boldsymbol{x}, t),$$

$$Q= \bar{f}(\boldsymbol{x}, t)=\bar{f}(\chi_{_{R}}(\boldsymbol{X}, t), t)=\widehat{f}(\boldsymbol{X}, t).\label{L2E} \tag{2}$$

The equal sign in above two expressions means two sides of each sign represent the same physical quantity and have identical value rather than identical form.

Extracting Time History of Eulerian Results

In pre-processing solid mechanics simulations, we mesh material bodies. This is when Abaqus creates Lagrangian coordinates. Then Abaqus solves the system of equations obtained by the finite element procedure. The resultant quantities are in two categories: node-based quantities and element-based quantities. Node-based quantities, namely \(Q_N\), are scalar and vector results such as displacements, temperature, etc., which usually are the direct variables to be solved at nodes in the final system of equations. Element-based quantities are higher order (\(\geq2\)) tensor results such as strains, stresses, etc., which are then calculated at integration points from results of direct variables by deformation measures and constitutive law. We name this type of quantities \(Q_{IP}\) indicating they are calculated at integration points. In most solid mechanics simulations Abaqus calculates all those results in Lagrangian representation, which can be represented as

$$Q_{N}=\widehat{f}(\boldsymbol{X}_{N}, t),\label{qn} \tag{3}$$

$$Q_{IP}=\widehat{f}(\boldsymbol{X}_{IP}, t).\label{qip} \tag{4}$$

If we want to have the time history of quantities in Eulerian representation at a requested spatial coordinates \(\boldsymbol{x}_{0}\), in other words, we want

$$Q=\bar{f}(\boldsymbol{x}_0, t),$$

which can be transformed from Lagrangian results by equation (\(\ref{L2E}\)). Thus we need to find every material point labeled by \(\boldsymbol{X}\) that occupies \(\boldsymbol{x}_0\) at every instant \(t\) of simulation by

$$\boldsymbol{X}=\chi_{_{R}}^{-1}(\boldsymbol{x}_0, t),$$

and extract quantity value for every frame \(t\) by

$$Q=\widehat{f}(\boldsymbol{X}, t).$$

In the simplest situation, where the quantities we’re interested in are some node-based result \(Q_{N}\) or element-based result \(Q_{IP}\), and we got lucky that material point \(\boldsymbol{X}\) at some instant \(t\) turns out to coincide with a node \(\boldsymbol{X}_{N}\) or a integration point \(\boldsymbol{X} _{IP}\), then we can get quantity value at \(t\) by equation (\(\ref{qn}\)) and (\(\ref{qip}\)), respectively.

However, a much greater chance is the returned material point \(\boldsymbol{X}\) doesn’t coincide with either a node nor a integration point. In these cases, an ideal treatment would be to find the element that encloses the requested spatial location \(\boldsymbol{x}_{0}\), then, for node-based quantities, interpolate the result at spatial location from results at nodes, and for element-based results, first extrapolate results from integration points out to nodes, then interpolate the result from nodes to the spatial location. This is the best accuracy of Eulerian results we can get from a finite element analysis. To completely implement this process requires a meticulous implementation of the shape function of each different types of elements.

Here I implemented two workaround methods. The first method takes advantage of the tool of spatial-point-based path provided in the Abaqus post-processing. This method basically adopted the above method of obtaining values at spatial locations, which should have a better accuracy. The second is based on finding the closest node to the spatial location and approximately taking values from or around the node rather than extrapolation and interpolation. This sacrifices the accuracy to some extent depending on mesh density. Each method has their limitations and advantages so I ultimately implemented both.

Path Method

Abaqus, in post-processing, allows users to create a spatial points based path and extract results from the path in each single solution frame. Thus, we can create a super short path centered at our requested \(\boldsymbol{x}_0\) to take advantage of this built-in tool. The extracted results from path have multiple values since Abaqus at least extracts results at path’s start and end points. We take the average value of results obtained from path as the result for \(\boldsymbol{x}_0\). This extraction process is automated to loop over all solution frames to obtain the time history result, from which a XYData object and XYPlot will be generated by the script.

The method should gives the best accuracy of Eulerian result we can expect from a finished ODB file. Also, the method runs fast since the data extraction from path is implemented internally by Abaqus.

The main limitation is the path has to intersect with some material body, from which the field output is extracted. There exists instances where as the deformation progresses, the path that intersects with material initially no long later. If this happens the method fails to extract field value. Another minor limitation is the method indifferently extracts results from material that intersects with the path and does not distinguish instance of part in a model. Thereby, if we only want to extract results at a requested \(\boldsymbol{x}_0\) from a specific instance of part of the model, we cannot use this method. Lastly,

in cases the method works fine, the accuracy of extraction is controlled by length of the short path, which is made as an input in the script. Users should examine the intersection of path and material body to make sure an appropriate condition.

To bypass these limitations, I implemented a second method, which is based on finding the closest node and extract results based on the node with approximation.

Closest Node Method

This method is based on finding the node \(\boldsymbol{X}_{N}^{0}\) from the input instance of the model that is closest to the requested \(\boldsymbol{x}_0\), at each solution frame \(t\), then, extracting the value of requested field output with the following approximation:

• for node-based results, the method directly returns the result from the closest node \(\boldsymbol{X}_{N}^{0}\)

$$Q=\bar{f}(\boldsymbol{x}_0, t)\approx \widehat{f}(\boldsymbol{X}_N^0, t),$$

• for element-based results, the method further finds all elements that shares the closest node \(\boldsymbol{X}_{N}^{0}\), and average results at each integration point \(\boldsymbol{X}_{IP}^{i}\) of all those elements with a total number of integration points of \(\Sigma\boldsymbol{X}_{IP}^i\) , regardless of full integration or reduced integration, and for beam, membrane, shell elements only considers in-plane integration points

$$Q=\bar{f}(\boldsymbol{x}_0, t)\approx \frac{1}{\Sigma\boldsymbol{X}_{IP}^i} \Sigma\, \widehat{f}(\boldsymbol{X}_{IP}^i, t).$$

The advantage of this method is the there always exists a closest node \(\boldsymbol{X}_{N}^{0}\) to the spatial point \(\boldsymbol{x}_{0}\), even if the material moves out of it. The second example problem next shows a case where even the material has moved out of \(\boldsymbol{x}_{0}\), the result extracted from \(\boldsymbol{X}_{N}^{0}\) is still want we want.

The drawback of this method is its expensive running time. Finding the closest node requires the current coordinates of each node, however, this output labeled by ‘COORD’ in Abaqus is not a default output. In this case, the current coordinates of each node have to be calculated from initial coordinates and displacements at each solution frame, which takes time. Thus, a potential node label series is made as an input to reduce the number of nodes that’s in the searching range, of which the current coordinates calculation are required. If the Abaqus ODB file has ‘COORD’ output, the script will directly use this for computation, which would significantly reduce running time.

Example problems

Here I’ll demonstrate two examples in which Eulerian results are of special interest and the two Python scripts are used to extract the time history of the Eulerian results. Those particular solid mechanics problems usually involve a steady state during the deforming process, when achieved, the material state like stresses and strains remain constant in space but vary among material points. Eulerian results can directly reflect if and how a steady state is achieved. Another minor motive of using Eulerian results is measurements in those problems are usually taken at spatial locations by placing sensors there. Thus, in validating finite element models, it’s more convenient if we can directly compare Eulerian results from simulations to test values.

Dynamic Simulation of Nip Mechanics in Web Transiting

A couple of nip rollers is widely used in a roll-to-roll process to regulate web velocity, adjust and separate web line tension, feed ink, laminate web and etc. In studies of the effect of compliance of nip cover on web velocity, we run 2D implicit dynamic simulations and are particularly interested in the velocity and tension in the upstream, nip contact zone, and the downstream of the nip couple. In the following example, the web is applied with different tensions at upstream and downstream ends. The two rollers have same rubber cover and the top roller is applied with a velocity to transmit web ahead. We used the path method script to extracte the time history of \(V_1\) and \(S_{11}\) at three locations: \(\color{red}{(-1, 0, 0)}\), \(\color{green}{(0,0,0)}\), and \(\color{blue}{(1,0,0)}\) shown in the following animation.

The script extracted Eulerian results through the whole simulation steps and frames and generated XYData objects, which are plotted in the following. As the XYData objects are created, we can also export their data to other programs, such as Excel, to plot and post-process data.

The closest node method was also tried, although it runs much slower than the path method due to the calculation of current coordinates of each node in each frame, the results are identical with that of the path method.

Dynamic Simulation of Lateral Mechanics of Web Camber

In studies of lateral mechanics of web with camber, we are specifically interested in the lateral deformation of cambered web transporting on rollers. To measure the lateral deformation of web we evenly placed five digital micrometer sensors S1~S5 on the edge of web in the test span, shown in the below schematic. We conducted implicit dynamic simulations on the problem. More information about the problem can be found in this post.

In post-processing, we hope to have a direct comparison with test results, thus, we extract \(U_3\) at five corresponded locations, shown in the following animation. Although the path method would be more efficient, it fails in this problem because the web moves laterally and no longer intersects with path created initially. The closest node based method was used.

The time history of Eulerian results show that \(U_3\) achieves a steady value at the end of the simulation and varies between those five locations. The deformed camber web in the test span has a S-shape.

How to Use

Source codes of the two Python scripts and instructions can be found in my GitHub repo:

https://github.com/yaor42/EulerianExtractor