Fracture Formation and FEM
We’ll be looking at the phenomena of fracture formation, and modeling it using the finite element method based on this thesis by P. Federl. The general idea is that we have a material which consists of a surface with two layers that are connected. The bottom layer grows, causing the top layer to deform. When the amount of deformation passes a threshold, a fracture is formed. You can assume most if not all of the material presented is from the thesis, while any mistakes are my own.
Tensors and Growth
The growth is modeled using a growth tensor of rank 2. Think of it as a matrix, or an extension of a vector. A vector has magnitude and one direction. A tensor of rank 2 has magnitude and 2 directions associated with it. If we have a vector and we want to obtain another vector, we can multiply it by a scalar or a matrix. The matrix can change the direction, while the scalar cannot (except for transforming it in the opposite direction). Here is more info on tensors if you’re curious, they’re pretty wild. We can solve for the velocity of growth at any point from this matrix by solving for the velocity components of the following set of linear equations, where T is the growth tensor. The velocity components can then be calculated by a given position.
So for any point on the surface, we can get a growth vector, and this gives us a field of velocities.

The above is from the thesis, and shows you different types of growth that can occur on a surface. Isogonic growth is growth where the rate of growth is the same in all directions at all locations. If I apply this to a circle, it will still be a circle. Uniform anisotropic growth is that the growth is the same at all locations, but changes given a direction. If I apply this to a circle, it will become deform into another shape. Note that the growth operates on the shape such that the shape is at the origin.
In this example we’ll simply use growth along the x axis, where R is the rate of expansion.
And the corresponding growth tensor is simply
Creating the Surface
We begin by creating a triangulated surface with semi-random triangles. Because of the nature of FEM, uniform triangulation and random triangles can have undesirable effects. With random triangles, some points can be too close. Place the nodes on the surface randomly, and repel them. I used the separation function we made in the Boids algorithm, adding a drag force proportional to velocity. Run this until all of the nodes are in equilibrium. You can implement this as running the loop until all of the velocities are zero, or you can choose to run the loop for a certain amount of time as the particles separate (but this does not guarantee the particles will fully separate evenly). Now triangulate the nodes using a Delauney triangulation.
This gives you a triangular mesh on a plane. Project this plane in the normal direction to get a second layer. Each triangle now becomes a wedge element. Each wedge element has six nodes.
We are going to focus on this element for now. Then we will see how to combine all of them to calculate the stress on the nodes.
Finite Element Method – Qualitative
The finite element method can provide approximate numerical solutions to a wide range of engineering problems. This is most useful when encountering problems where an analytic solution is difficult to find, or simplifying assumptions cannot be made. We begin first with something familiar: finite differences. These are also used to find approximate numerical solutions. If we have a shape, we can approximate it using a grid to cover the shape. More points (a finer grid) produces a more accurate covering. The finite element method is similar, only it uses triangles instead and can much better approximate the region.
The above figure is from here, which provides a good introduction to the ideas of FEM. When we created smooth particle hydrodynamic simulations, we approximated continuous functions and properties over a region using kernel functions and interpolation. Similarly with FEM we divide the region into nodal points and use interpolation functions to approximate the property of interest over the region.
The FEM allows one to divide a region into elements, calculate the solutions for individual elements, and then combine the element solutions to form a soltuion for the region. In our case we find the stiffness of each element, and then combine them to obtain the stiffness of the entire material.
There are a variety of ways to calculate the properties of the individual elements. We choose the direct approach. Other approaches include the variational approach, and the weighted residuals approach. We divide the region into elements, calculate the element properties with an interpolation function with the element nodal points, then combine each of these element matrices into a giant matrix to represent the entire system. We impose the boundary conditions (any values we know about the system and the nodes) and then solve the matrix equation for the system.
Finite Element Method – Quantitative
So from the perspective of an element, the amount of force generated is proportional to the displacement of the nodes. Sound familiar? In the matrix form is looks similar to F=kx, Hooke’s law to determine force for a spring.
Where F is the force of the element, K is its stiffness, and q is its nodal displacements. For each element we have six nodes. Each node has three degrees of freedom, in the x direction, y direction, and z direction. This gives us a total of 18 degrees of freedom. [F] is an 18 x 1 vector with F_x, F_y, and F_z for each node. Similarly [q] is an 18×1 vector with the corresponding q_x, q_y, q_z for each node.
If we want to determine the value of some function over the entire element, it is helpful to define a coordinate system within the element. These coordinates would be called isoparametric coordinates. We can define three variables to integrate over in the wedge.
We can now use interpolation with the six nodes to determine the value of any point within the element, as a summation of the value at the node times a weighting function N. The weight is called a shape function, and is some value that depends on the isoparametric coordinates, and is inbetween 0 and 1. Mathematically this looks like
And if we have a matrix of N_i * I and concatanate them for i = 1, …, 6, where ‘I’ is the identity matrix, we have an 3×18 matrix for N. Given V is the 1×18 vector of concatanated values for each of the nodes, we can represent the above as a matrix equation [v] = [N][V]. If Q is the 1×18 vector of position points of all the nodes, we have [x,y,z]=[N][Q], converting the coordinates in the element to cartesian coordinates. Given a displacement [u,v,w] in an element, we can interpolate it over the nodes in an element using the shape function N.
Where u, v, w, and N are dependent on the isoparametric coordinates.
Before contintuing onward to calculate the element stiffness matrix, it is helpful to have a better understanding of the shape function. It should be pretty familiar in concept, and fortunately it is much simpler than the analog we used in SPH.
For any of the shape functions, you can vary the coordinates and see how it changes within the element. We can integrate over these coordinates to find [K_e], the element stiffness matrix. We can approximate the integral using the Gaussian integration technique. Instead of computing the continuous integral, we compute the discrete sum over select points that are weighted.
Now that we know how to go about integrating, we need to determine what exactly to integrate over. To find [K] we integrate over two other values, the material stiffness matrix [D], and the matrix to relate displacement with strain, [B]. The material is uniform throughout, so the stiffness matrix is I * constant. We can derive strain in each direction as a summation over the shape fuctions, since we can relate the shape functions and the displacement. There is also the Jacobian matrix [J] that transforms between coordinate systems (cartesian and isoparametric). The final equation for [K_e] is
where we use n Gaussian integration points.
Note we can rewrite the matrix equation for stiffness as [F] = [K][P]+[R] if we notice [q] = [P] – [Q] , where [Q] is the undeformed coordinates. So to obtain the global stiffness equations we must find
For a system of n nodes, F_g is a vector of length 3n, since each node has three force components. Similarly the coordinate vector [P_g] has length 3n. The global stiffness matrix [K_g] is simply adding the appropriate values at the correct indices of the element stiffness matrices.
In our simulation the bottom layer is fixed and the top layer is free to move. We set some initial strain (the growth of the bottom layer), insert the boundary conditions (force on all nodes is zero), and solve the equations (calculate the nodal coordinates of the free nodes. We now solve the familiar equation
Our set of linear equations is pretty sparse, since we created the global stress matrix by simply combining the element matrices. Because of this, we can use the conjugate gradient method to solve the set of linear equations.
We’ll take a brief moment to discuss a few concepts.
And here is an algorithm for it, from Wikipedia.



- repeat



- if rk+1 is sufficiently small then exit loop end if



- end repeat
- The result is xk+1
We calculate the eigenvalues and eigenvectors of the stress tensor on the mesh to determine the principle stresses at given points. If the largest of these stresses exceeds the material’s stress, a fracture is formed.
There are a few ways to represent the fracture. One way is to delete the element. This can produce fractures, though it is not realistic. Another method is to subdivide the elements along the principle stress plane.
To optimize we can use dynamic subdivision. If a given element is a candidate for fracture, we subdivide and recalculate the stresses. Once an element is small enough, a fracture can form on element relaxation.
More info will be posted soon!






