In this blog post we will have a look at how the use of the package Lasso can help out the post processing of an LS-Dyna analysis. For this blog we will use the free version of Lasso and this can be downloaded from here .
The model
Some simple yet illustrative examples are often useful when learning new things. In this blog post we will use an example from DNV-RP-C208, which is a document describing recommended practises from DNV (Det Norske Veritas). There is no intention of recreating the results that can be found in DNV-RP-C208 for this example. If you are interested in that topic, a great paper by Marcus Lilja can be found here
The example is related to determining the buckling resistance of a “crane”. The geometry is shown below.
The finite element model can be seen below and contains approximately 3200 shell elements. The model is constrained along certain boundaries and this is also visible in the figure.
Buckle up for safety…
OK, let's start with the finite element techniques that will be utilized to determine the buckling resistance. This contains three steps:
- Set up and run a linear buckling analysis
- Extract the first buckling mode
- Use the first buckling mode as a geometric imperfection to a “full” nonlinear analysis (material nonlinear and geometric nonlinear)
OK, that looks nice. Of course, since this state of deformation represents an eigenvector, it is more the deformation shape that is of interest not the magnitude. After extracting the first mode it will be scaled before imported to the nonlinear analysis. To test out the model, let’s scale the shape of the first mode so that the maximum magnitude of the imperfection is 19.5 mm. The material properties for the nonlinear analysis can be seen below
| Parameter | Value |
|---|---|
| Density | 7850 kg/m^3 |
| Young's modulus | 210 GPa |
| Poisson's ratio | 0.3 |
| Yield strenght | 355 MPa |
| Tangent modulus | 1000 MPa |
These are the same properties that were used for the linear buckling analysis, but then only the elastic properties were used. The resulting von Mises stress can be seen below
Nice, then the “pure” LS-Dyna stuff can be considered as done and now for the part that this blog post is about
Post processing with Python and Lasso
The first thing that needs to be done is to install Lasso if you haven't already done that. This can be done using pip:
python -m pip install lasso-python
If this works out fine you are good to go. Let's test the installation by doing some imports. We will here also import the other packages that will be useful.
import os from lasso.dyna import D3plot, ArrayType, FilterType, Binout import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl mpl.style.use('classic') %matplotlib inline
Once that is done, the next step is to create some data to play around with. We will here use the nonlinear analysis described above, but change some of the input parameters. To do this, we simply create some files containing the parameters that we are interested in. First let's create some folders where we can store the results:
# Create folders for analysis wai = os.getcwd() Name = ['Run0'+str(i+1) for i in range(6)] for na in Name: os.mkdir(wai+os.sep+na)
Sweet, now we have six folders to store stuff in. The next step is to create the parameter-files:
# Paremeter 1: Initial pertubation of 19.5 mm and 30.8 mm Par1 = [0.0687, 0.1087] # Parameter 2: Tangent modulus of 1000, 3000, and 5000 MPa Par2 = [1000, 3000, 5000] i = 0 for pa1 in Par1: for pa2 in Par2: # Print parameter file f = open(wai+os.sep+Name[i]+os.sep+'Parameter.key','w') f.write('*KEYWORD\n') f.write('*parameter\n') f.write('rpert, %5.4f, ryield, %5.3f\n' % (pa1,pa2)) f.write('*END') f.close() i += 1
When opening one of the Parameter.key files in folder Run01, this is what it looks like:
*KEYWORD *parameter rpert, 0.0687, ryield, 1000.000 *END
The two parameters that we change in the six analysis are the initial perturbation (rpert) and the hardening parameter (ryield) for the material. A summary of the settings can be seen below:
| Analysis | Parameter 1 (rpert) | Parameter 2 (ryield) |
|---|---|---|
| Run01 | 19.5 mm | 1000 MPa |
| Run02 | 19.5 mm | 3000 MPa |
| Run03 | 19.5 mm | 5000 MPa |
| Run04 | 30.8 mm | 1000 MPa |
| Run05 | 30.8 mm | 3000 MPa |
| Run06 | 30.8 mm | 5000 MPa |
OK, let us start with comparing the reaction force in y-direction at the possition where the creane in attached to the ground. To do this, all we have to do is to load the binout file and find the reaction forces and plot them in a graph. This can be done with just a few lines of code.
# Plot Run01 to Run03 Name = ['Run0'+str(i) for i in [1,2,3]] for na in Name: binout = Binout(na+os.sep+"binout") y_force = np.sum(binout.read('spcforc','y_force'),axis=1) plt.plot(-y_force,lw=2,label=na) # Plot Run04 to Run06 Name = ['Run0'+str(i) for i in [4,5,6]] for na in Name: binout = Binout(na+os.sep+"binout") y_force = np.sum(binout.read('spcforc','y_force'),axis=1) plt.plot(-y_force,'--',lw=2,label=na) plt.legend(loc=4) plt.ylabel('Y-Force [N]') plt.xlabel('Analysis step') plt.grid()
and the result can be seen below:
That was really easy. Note that the solid lines represent the “small” perturbation and the dashed lines represent the “larger” perturbation. The maximum load does not seem to be affected to any greater extent by the change of perturbation and hardening parameter. However, what happened after the peak load looks to be influenced by the hardening parameter, as would be expected.
To further investigate the models, let's load some d3plot data and have a look at the distribution of plastic strain. Here we will only investigate two analyses: Run04 and Run06. Since Run06 has a higher value on the hardening parameter, it would be expected that the maximum value of plastic strain will be smaller for Run06 than for Run04. It is also expected that a higher number of elements will undergo plastic deformation. Let’s try to create a plot to show if our expectations are correct.
# Load d3plot. Only load the last state d3plot1 = D3plot("Run04/d3plot", state_filter={-1}) d3plot2 = D3plot("Run06/d3plot", state_filter={-1}) #Extract the plastic strain Plastic_strain1 = d3plot1.arrays[ArrayType.element_shell_effective_plastic_strain] Plastic_strain2 = d3plot2.arrays[ArrayType.element_shell_effective_plastic_strain] # Compute the element average over the different integration points pstrain1 = Plastic_strain1.mean(axis=2) pstrain2 = Plastic_strain2.mean(axis=2) plt.xlabel('Analysis step') plt.grid()
Sweet! Now we need to visualise this in some way that makes sense. One way to do this is to sort the plastic strain (from smallest value to largest value). By doing so, it is possible to find out how many of the elements that have not reached the plastic area. Also it will be visible if our expectation on the plastic strain mentioned above is correct or not. One again, not much code is needed to do this:
plt.plot(np.arange(len(pstrain1[0]))/len(pstrain1[0])*100 ,np.sort(pstrain1[0])*100,'b',lw=2,label='Run04') plt.plot(np.arange(len(pstrain2[0]))/len(pstrain2[0])*100, np.sort(pstrain2[0])*100,'g',lw=2,label='Run06') plt.grid() plt.legend(loc=2) plt.ylabel('Effective plastic strain [%]') plt.xlabel('Total fraction of element [%]') plt.xlim(80,102)
And here is the resulting graph.
OK. From this graph it is clear that the maximum value of plastic strain is higher for the analysis with lowest hardening parameter (Run04). Since the curves cross over at approximately 96% it is also clear that more elements has reached the plastic area. Another way to put this is that the plastic zone will be more localized for the model with lower value of the hardeing parameter, wish also was expected
Time to conclude
In this blog post we had our first look at Lasso. It will not be the last. Here focus was mostly on getting the data into the system, but in upcoming blog posts we will utilize Lasso to do more statistical analysis on data coming from LS-Dyna analysis.
That is all for this time. As always, please feel free to reach out to us on info@gemello.se. Happy coding and post processing!