Get the crack tip trajectory

This post focuses on the crack tip tracking method. While I simulate several fatigue models, the crack tip trajectories are hard to track in the model. Previously, I just pick some points and decide the crack length and rate in several images. This useful method is not very effective when I get several video clips in crack place. Several image process tools such as Fiji can deal with the image fast and easy, but I still feel it is hard to realize my goal. So I try to write a simple tool to achieve the object. Since tracking the crack tip uses the technique of image processing, I use the Mathematica for this purpose.

The main function used in Mathematica is the function called “ImageFeatureTrack”. This function will track some feature points in a sequence of images, and return the coordinates of these points in each image. Thus, the idea is not very difficult, we import the video and transfer it into the image sequence, then we specify the tracking points with the coordinates in the image and use “ImageFeatureTrack” function to get the coordinates in each frame. The last step is to turn the image coordinates into the real model coordinates.

The first step is to import the file and to turn into a sequence of images. Additionally, the frame numbers are added to each frame to help inspect the image sequence,


Snippet: An Abaqus inp example

An example to show the inp file structure. The example demonstrates the stretch of a square plate with the dimension $1 \times 1 \mathrm{m^2}$. The plate material is steel with elastic stiffness of $2.1 \times 10^{11} \mathrm{MPa}$. The material density is $7800 \mathrm{kg/m^3}$. The left boundary is fixed and $0.1 \mathrm{m}$ displacement is applied to the right boundary. The plane strain is assumed in the model. Here are the model and the simulation results.


Plastic strain in MD simulation

This post provides the approach to extract the plastic strain in MD simulation. We can compute the atom displacement as the trajectory from the MD simulation, but the strain is the concept in continuum mechanics, it is really hard to compute the tensor in a discrete system. Recently, some tools provide this function based on a theory that accounting the effects of the neighbor atoms. While considering the relative displacement with the neighbor atom and using the weight function, the deformation gradient tensor can be calculated. Therefore, the strain tensor is obtained. Plastic strain can be expressed as,

$$\boldsymbol{\varepsilon}^p  = \boldsymbol{\varepsilon} – \boldsymbol{\varepsilon}^e $$

so two strain tensor should be generated in the simulation data: total strain tensor $\boldsymbol{\varepsilon}$ and the elastic strain tensor $\boldsymbol{\varepsilon}^e$. Obviously, LAMMPS data file does not provide such information. Here I use the OVITO to get this data.

OVITO provide two modifiers, the first one is the “Atomic strain” and the second one is the “Elastic strain calculation”. Both of them provide the deformation gradient tensor and the strain tensor data. But I found there are some bugs for strain tensor output in OVITO, so I use the deformation gradient tensor data and computing the strain tensor with a Python script.


Snippet: Plot format in Wolfram

This snippet demonstrates the polt format that I often used in Wolfram Mathematica.

Plot[Sin[x], {x, -Pi, Pi}, 
     Mesh->None, PlotRange->{{-Pi, Pi}, {-2, 2}}, 
     AspectRatio->9/16, PlotTheme->"Scientific", FrameStyle->Black, 
     TicksStyle->Directive[FontSize->20], ImageSize->Large,
     AxesLabel->Automatic, LabelStyle->Directive[FontSize->20,
     FontFamily->"Times New Roman"], Epilog->{Text[Style[
     ToExpression[ "y = 12 \\frac{\\partial^2 f(x)}{\\partial x^2}", 
     TeXForm, HoldForm], FontFamily->"Times New Roman", FontSize->20, 
     Bold, Black], Scaled[{.25, .80}]]}]

Here is the diagram,

Here is the Wolfram notebook (Ua98n5.nb) used in this post.

Notes: Simple tensor algebra III

Previously, the vector space $\mathbb{V}$ and three vector product operations were considered in Simple tensor algebra I and Simple tensor algebra II. But we do not talk much about the tensor. Here, the tensor will be introduced.

So, what is a tensor? The following figure shows a tensor in a 2D Cartesian coordinate system.


The vector is a tensor, and it is a first-order (1st-order) tensor. Actually, the scalar is a tensor as well, we can consider the scalar as a zeroth-order (0th-order) tensor.

Let $ a, b$ be two scalars in $ \mathbb{R}$, we put these two scalars (0th-order tensor) together with $ \{ a, b \}$ to form a vector (1st-order tensor). So, if we can put two vectors $ \boldsymbol g, \boldsymbol f $ together, then we will get a 2nd-order tensor. The tensor product $ \otimes $ plays such role here ($ \boldsymbol g \otimes \boldsymbol f $).

This is not a tutorial, please refer to the books for tensor algebra.


Notes: Simple tensor algebra II

The vector space $ \mathbb{V}$ is introduced in Simple tensor algrbra I, several rules are provided in it. To reduce the abstract, we represent the vectors $ \boldsymbol{x}, \boldsymbol{y}$ with the scalar coefficients $ \{{x_1},{x_2},{x_3} \}$ and $ \{{y_1},{y_2},{y_3} \}$ in $ \mathbb{E}^3$, respectively. The composition of two vectors $ \boldsymbol{x}, \boldsymbol{y} \in\mathbb{E}^3$ is expanded by,

Dot product ($ \cdot $):

$${\boldsymbol{x}} \cdot {\boldsymbol{y}} = {x_i}{{\boldsymbol{g}}^i} \cdot {y_j}{{\boldsymbol{g}}^j} = {x_i}{y_j}{{\boldsymbol{g}}^i} \cdot {{\boldsymbol{g}}^j} = {x_i}{y_j}{g^{ij}} = {x^i}{y_i} = \left[ {\begin{array}{*{20}{c}}
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
\end{array}} \right] = \sum\limits_{i = 1}^3 {{x_i}} {y_i}$$

The dot product is also called scalar product or inner product under Cartesian coordinates and it provides the Euclidean magnitudes of these two vectors and the cosine of the angle between them.

Cross product ($ \times$):

$${\boldsymbol{x}} \times {\boldsymbol{y}} = {x_i}{{\boldsymbol{g}}^i} \times {y_j}{{\boldsymbol{g}}^j} = {x_i}{y_j}{{\boldsymbol{g}}^i} \times {{\boldsymbol{g}}^j} = {x_i}{y_j}{\varepsilon ^{ijk}}g{{\boldsymbol{g}}_k} = \left| {\begin{array}{*{20}{c}}
\end{array}} \right|$$

The cross product is also called vector product, and it represents the area of a parallelogram with the sides of these two vectors. Now, a new composition for two vectors $ \boldsymbol{x}, \boldsymbol{y} \in\mathbb{E}^3$ is introduced here.

This is not a tutorial, please refer to the books for tensor algebra.


Notes: Simple tensor algebra I

The vector space $ \mathbb{V} $ with the operation (+) over a field of real number $ \mathbb{R} $ is an abelian group with a scalar multiplication. So, suppose $  \boldsymbol{x} $, $ \boldsymbol{y} $ and $ \boldsymbol{z} \in \mathbb{V}$, they satisfy the following conditions:

Abelian group:

  • Closure:$ \boldsymbol{x} + \boldsymbol{y} \in \mathbb{V}$,
  • Associativity: $ ( \boldsymbol{x} + \boldsymbol{y} )+ \boldsymbol{z}=  \boldsymbol{y} + ( \boldsymbol{x}+ \boldsymbol{z} ) $,
  • Identity: $ \exists \boldsymbol{0} \in \mathbb{V}, \forall \boldsymbol{x} \in  \mathbb{V} $ such that $ \boldsymbol{0} + \boldsymbol{x} = \boldsymbol{x} $ and $ \boldsymbol{x} + \boldsymbol{0} = \boldsymbol{x} $,
  • Invertibility $ \exists ! ( -\boldsymbol{x} ) \in \mathbb{V}, \forall \boldsymbol{x} \in  \mathbb{V} $ such that $ ( -\boldsymbol{x} ) + \boldsymbol{x} = \boldsymbol{0} $ and $ \boldsymbol{x} + ( -\boldsymbol{x} ) = \boldsymbol{0} $,
  • Commutativity: $ \boldsymbol{x} + \boldsymbol{y} =  \boldsymbol{y} + \boldsymbol{x} $.

Scalar multiplication:

  • $ \forall \alpha, \beta \in \mathbb{R} $, $ \forall \boldsymbol{x} \in \mathbb{V} $ such that $ ( \alpha \beta ) \boldsymbol{x} = \alpha ( \beta \boldsymbol{x} ) $,
  • $\forall \alpha, \beta \in \mathbb{R} $, $ \forall \boldsymbol{x}, \boldsymbol{y} \in \mathbb{V}  $ such that $ \alpha ( \boldsymbol{x} + \boldsymbol{y} ) = \alpha  \boldsymbol{x} + \alpha \boldsymbol{y} $ and $ ( \alpha + \beta )  \boldsymbol{x}  = \alpha  \boldsymbol{x} + \beta \boldsymbol{x} $,
  • $ 1 \boldsymbol{x} = \boldsymbol{x} $.

This is not a tutorial, please refer to the books for tensor algebra.


Build VTK 8.1.0 with Qt in Mac

This is the second time I try to build the VTK in Mac. One month ago, I took a lot of time to find how to compile it with CMake and built it successfully, but it is much easier to forget the process. Here I record the entire compiling process.

First step is to download the VTK from the VTK official website. Then, unzip the package to a folder. Here I put it into the folder called “VTK-8.1.0”. By the way, you should also install CMake from the CMake website (Please ensure that the Qt and Xcode are installed in your Mac). Now, create another folder called “VTK_Build”. This is the target folder to save the files generated by CMake. Open CMake and set the “source code” and “binaries” folder:


Snippet: Matrix with Euler angles

This snippet posts an Excel to calculate the matrix with Euler angles. To describe the rotation of a crystal frame, commonly, we use three angles: $ \alpha$, $ \beta$, and $ \gamma$. They are also known as Bunge angles.

The matrix to rotate the crystal frame in reference frame is described by Euler angles as below:

$$ \displaystyle \mathbf{B}=\begin{bmatrix}
\cos(\alpha)\cos(\gamma) \\
– \sin(\alpha)\sin(\gamma)\cos(\beta)
\sin(\alpha)\cos(\gamma) \\
-\cos(\alpha)\sin(\gamma) \\
– \sin(\alpha)\cos(\gamma)\cos(\beta)
-\sin(\alpha)\sin(\gamma) \\
+ \cos(\alpha)\cos(\gamma)\cos(\beta)
\end{matrix}& \cos(\gamma)\sin(\beta)\\
\sin(\alpha)\sin(\beta) & -\cos(\alpha)\sin(\beta) & \cos(\beta)
\end{bmatrix} $$

In this Excel (Cf06b0.xlsx), you just need to input three Euler angles and the vector in reference frame. It will give the matrix and the inverse matrix with these three Euler angles. Also, it will give the vector in crystal frame.


This snippet is very useful for debugging the code with the rotation of Euler angles.

Add a 3D convex hull in ABAQUS

Convex hull is the smallest envelope that contains the points set. It is used to construct the grain in grain-based model. There are several methods to generate the convex hull data, but it may need some efforts if you want to put it into ABAQUS model. Here, I post a simple method which is suitable for programming.

Generate convex hull data

Several methods can be used to generate the convex hull data. I am used to generate the convex hull data with my own code. But if you do not like write code, you can use Mathematica or Matlab to generate the data.

For example, in Mathematica, you can use following code to generate a 3D convex hull:

pts = RandomReal[{-1, 1}, {20, 3}];