2D Elasticity using Finite Element Analysis: Vector Field Problems

In the previous cases, the field of interested, e.g., Temperature, Concentration, Membrane displacement were scalar fields. However, many real-life problems involve vector field, e.g., displacement field in solid mechanics or velocity field in fluid mechanics. In this course, we will look at only linear elastic solid mechanics problems. 

Assumptions of theory of linear elasticity:
1. small deformations
2. behavior of material is linear
3. dynamical effects are neglected (however, they could be included)
4. no gaps or overlaps can occur during deformation (displacements are single valued)

Regarding Assumption-1. Generally, the deformations are not visible to the naked eye -- something that would be very clear when you are TA for the solid mechanics lab. Such problems come under linear elasticity. When the deformations are large, but the material is locally linear elastic (e.g., large deformations of a beam) we call this geometric nonlinearities. For Assumption-2, when the material is non-linear to begin with (for example, rubber) then the problem is non-linear elastic. In some other cases, the material undergo permanent plastic deformations and such problems are also non-linear problems. There are many other constitutive relations in solids, complex fluids and geotechnical engineering. Regarding Assumption-3, no problem is truly static. We can model a problem to be a static problem if  the time in which the load is applied is large compared to the period associated with the lowest frequency, i.e., load is applied slowly and slow is compared to the period of oscillation.  Regarding Assumption-4, there are some other problems such as fracture when crack is formed, we can have a two-valued displacement function, and gaps can be formed. 

Below, we focus only on linear elastic problems, specifically, homogenous, isotropic, materials, and how to solve them with FEA.

Kinematic of small deformations

Consider that an elastic deformation of object in terms of displacement, 

However, displacement in itself does not quantify deformations, which is obtained in terms of strains. 
	(Subscript[ϵ, xx  ])= ((∂ (Subscript[u, x]))/(∂ x)) is to quantify elongation along the x direction
	(Subscript[ϵ, xy  ])= ((∂ (Subscript[u, y]))/(∂ y)) is to quantify elongation along the y direction
	(Subscript[ϵ, xy ])= (1/2)(((∂ (Subscript[u, x]))/(∂ y))+((∂ (Subscript[u, y]))/(∂ x))) is to quantify 1/2 the angle between x and y axis after deformation
Though strain is a tensor, in finite element analysis the strain is arranged in terms of a vector (in sense of a matrix)

These can be in terms of displacements be expressed as

where the (Subscript[∇, s]) is the gradient matrix operator can be expressed as

Description of stress and tractions

Stress in simplest terms is stress per unit area. However, a more correct definition will describe stress as a tensor. 

			Stress components. 
	Note that the over bar (Overscript[□, ->]) indicates that the force per unit area on that particular face. There is a more technical definition for this vector which is called as traction
	
	
	I will not do a detailed derivation, but it can be shown that on the face with outer normal:
		n = (Subscript[n, x]) i + (Subscript[n, y]) j
	the traction (vector force per unit area) can be expressed as:
		(Subscript[t, x])=(Subscript[σ, xx])(Subscript[n, x])+ (Subscript[σ, xy ])(Subscript[n, y])
(Subscript[t, y])=  (Subscript[σ, xy ])(Subscript[n, x])+ (Subscript[σ, yy])(Subscript[n, y])
	which in a more compact form be written as:
		t = τ n,
	where τ is what is the stress tensor.

Equilibrium equation

Problem definition. (a) Element from a particular location in an elastic body. (b) Forces acting on the element. 
It could be easily shown from force balance that

Here, (Overscript[b, ->]) is the body force (force per unit volume, e.g., gravity). In terms of components it could be written as:

which could be expressed in FEA matrix notation as:

and the equilibrium equations are:

Now, the last thing we require is constitute relation between stress and strain. 

Please note that (Overscript[(Subscript[σ, x]), ->]) = (Subscript[σ, x])i + (Subscript[σ, xy])j  and  (Overscript[(Subscript[σ, y]), ->]) = (Subscript[σ, xy])i + (Subscript[σ, y])j in a somewhat artificial vectorial notation. Consequently, the equilibrium equation can also be written in a compact vectorial form as:
	∇(Overscript[·(Subscript[σ, x]), ->]) + (Subscript[b, x])= 0  and 
	∇(Overscript[·(Subscript[σ, y]), ->]) + (Subscript[b, y])= 0 , 
where all the quantities are defined earlier.

Constitutive equation

For a linear, elastic, homogeneous, isotropic, element we get:

where for, 

In Plane stress: (Subscript[σ, zz])= (Subscript[σ, xz  ])= (Subscript[σ, yz])= 0  (for example, thin plate). 
In Plane strain: (Subscript[ϵ, zz ])= (Subscript[ϵ, zx])= (Subscript[ϵ, zy])= 0.
Here, z is the out of plane direction. 

We are finally ready to combine equilibrium relation and the constitutive relation to provide us with the strong form. 

Note that pretty much all concepts expressed here are similar to those discussed in the strong form for the scalar field. The boundary conditions are of two types here also:
	1. Essential boundary condition or the displacement boundary condition where displacement is specified on the boundary (Subscript[Γ, u])(d).
	2. Natural boundary condition or the traction boundary condition where the tractions (vector force per unit area) is specified on the boundary (Subscript[Γ, t])(c).

Weak form corresponding to the strong form

The weight functions in the current case is a vector function (can also be thought of as virtual displacement.) However, apart from the fact that many multiplications with the weight function in the current case are dot products, as opposed to scalar, multiplication, the basic idea of generating the weak form is the same in the case of elasticity. 


where, the weight functions are of the form:
	
As in the earlier cases, w vanishes on the essential boundaries (Subscript[Γ, u]). Now, using the Green’s theorem on this weak form, we get:

Adding these two equations, we get:

Here, we note that the total boundary Γ = (Subscript[Γ, t]) ⋃(Subscript[Γ, u]) and (Subscript[Γ, u])⋂ (Subscript[Γ, t]) = Null, and since w vanish on (Subscript[Γ, u]), we the more simplified version. We note, however, that:
	(Overscript[(Subscript[σ, x]), ->]) · (Overscript[n, ->]) = (Overscript[(Subscript[t, x]), ->]) and 
(Overscript[(Subscript[σ, y]), ->]) · (Overscript[n, ->]) = (Overscript[(Subscript[t, y]), ->]) 
due to which the weak form becomes:

This vectorial weak form is what we will now convert in the more relevant matrix form:
	
	Finally, converting the vector weak form in the matrix form we get the weak form to look as:
	
	Note that this form is a general form independent of if the material is linear elastic. However, in the case of the linear elastic material we note that:
		σ = D ϵ , 
where D  is as defined earlier for the plane stress and the plane strain case.  Hence, the final statement of weak for for the linear elastic systems is:

All the notation is pretty much that same as what we had discussed earlier. Note that (H^1) corresponds to all functions which are (C^0)continuous (no jumps but kinks possible) but which are (Subscript[L, 2]) integrable, i.e., ∫ (u^2)dΩ is finite.

FInite element implementation for the weak form

Most of the steps are the same as before. We first discretize the domain

Since there are two degrees of freedom per node in the x and y direction, they are expressed in matrix notation as:
.
As before, the weight (virtual displacements) and the actual displacements are expressed as:

Note, however, that the element shape function matrix is:

The first row is staggered with respect to the second row because we have every alternate element of d to be in the x and y directions, respectively, also for every element as:


The integral, as earlier, is now obtained as a sum of integrals over individual elements as:
	
We note also that:

where the strain-displacement matrix (B^e) is defined as:
	
The derivatives of the weight functions are:

After putting all the equations together and recalling that (d^e)=(L^e)d, and (w^eT)= (w^T)(L^eT) we get

Using the same ideas as before, we get the following:
Element stiffness matrix:

Element external force matrix:

The weak form can now be written as:

Then by using this compact notation we obtain:
This could also be re-written as:
 ,  the residue,  r  = Kd - f.
Doing the partitioning with respect to the essential (E) and natural (F) nodes, as in the previous cases, we now get:
	
Since (Subscript[w, e]) = 0  and (Subscript[w, f]) is arbitrary, it is clear that (Subscript[r, F])= 0. The above equations, hence, could also be written as:

where the partitioning of the stiffness matrix is done in the same manner as in the earlier case, by numbering the essential nodes (E) first followed with the number of natural (F) nodes.

FEA using three-node triangular element

The ideas are pretty much the same as before, only because of the vector nature of the displacement there are some technical modifications. For example, for the element e as described earlier, 
The displacement field in terms of shape functions is

Hence, the strains are:

where  and  Using the shape functions which we have derived earlier in the Notes (based on Fish and Belytschko), the strain vector becomes

where (Subscript[x, IJ]^e)= (Subscript[x, I]^e)- (Subscript[x, J]^e), which defined the (B^e) matrix of the element, which in the current case is constant for a given element and not a function of x or y. Going forward, the stiffness matrix in the current case becomes:
 
Since everything is constant, 
	
Element Body Force Matrix
The element body force matrix is given by

There are two ways of evaluating this matrix:
1. by direct numerical integration from the expressions for b and N
2. by interpolating b, usually with a linear function, and integrating the results in closed form. Note: even in the case above for direct numerical integration such kind of interpolation can also be done. Also, sometimes, the forces are experimentally only obtainable at a discrete set of points and so interpolation needs to be done to obtain a complete field. 

The most convenient way to do interpolation for b is by using the same interpolation function (N^e). Using the triangular coordinates since it is convenient to do that way, we obtain

where (Subscript[b, x I]) and (Subscript[b, y I]) are the x and y components of the body force at node I. Substituting, we get

using the integration techniques for triangular element described earlier in the Notes. 
Boundary Force Matrix
The boundary force matrix is given by

Here, too we will have to use some kind of interpolation for t on (Subscript[Γ, t]).  

Triangular three-node element showing nodal displacements and nodal forces 

For the purposes of explanation, let us assume that the edge 1-2 is the boundary edge. Along this edge, (Subscript[N, 3]^e)vanishes, since it is zero and 1 and 2 and the function is linear between 1-2. As result, the displacement can be expressed in terms of (Subscript[N, 1]^(2L))and (Subscript[N, 2]^(2L)). Here, the superscript 2L corresponds to two-node linear element. 

Hence, the integral becomes:

Note here that the arc length along the boundary dΓ = l dξ .For the triangular element ξ ∈ [0,1]. The integral could be now carried out easily to provide

As expected, the traction forces on the internal node 3 of the element are zero in both x and y direction. When the tractions are constant, i.e., (Subscript[t, x1])= (Subscript[t, x2])= (Overscript[(Subscript[t, x]), -]) and (Subscript[t, y1])= (Subscript[t, y2])= (Overscript[(Subscript[t, y]), -]), we obtain

which shows that the total forces (assume uniform thickness of unity, else multiply with the thickness b) are split equally between the two nodes.

Generalization of boundary conditions

Though we specified the strong form of the equations earlier, there is more richness in the manner in which the traction and boundary conditions can be specified on the boundary since there are two conditions (Subscript[u, x]) and (Subscript[u, y]) and (Subscript[t, x]) and (Subscript[t, y]). Thus we can have combinations of traction and displacement boundary conditions on the same boundary and not a clear separation. 

The basic rules about the natural and essential boundary conditions in the vector for are as follows. 
	1. At every point of the boundary Γ one need to have one condition (traction or displacement) each in two perpendicular directions. 
	2. One cannot have traction and displacement condition simultaneously specified in the same direction.
	
These two conditions can mathematically be specified in a simple manner as follows. 

On the total boundary Γ we can have specification of the following:

The weak form can be derived by an appropriate choice of (Subscript[w, x]) and (Subscript[w, y]) on the boundary. Note that the same component of traction and displacement cannot be prescribed on any part of the boundary and so:

Also, since for each component, either the traction or the displacement has to be specified

Note that (Subscript[t, x]) and (Subscript[u, x]) (and likewise (Subscript[t, y])and (Subscript[u, y])) are work conjugates and such you cannot specificy work-conjugates at the same point. 

Example to Illustrate boundary conditions

Consider the problem of plate with a hole. 

All the boundary conditions are natural boundary conditions with either zero or some specified tractions at all the boundaries. However, if we solve the problem purely as a natural (traction) boundary condition problem then, the resulting stiffness matrix will be singular, because translations and rotation of the plate are not prevented. To prevent that we can specify, for example:
	(Subscript[u, xA])= (Subscript[u, yA ])= (Subscript[u, yB])= 0,
which prevents all rigid body translations and rotations. A more elegant (and also computationally efficient) way it so recognize the symmetry of the problem as shown in (b) of the figure above. Along any line of symmetry (FG or HK) since the displacements normal to the axis will form mirror images, there will be either gap or overlaps as shown in (c) due to which the compatibility condition (specified at the beginning of the notes will be violated.) Hence, the displacements normal to the axis should be zero. Also due to symmetry the shear traction (traction component along the axis of symmetry) should also vanish. This will give us:

Thus we now have a modified boundary value problem where we model only a quarter of the domain with the following boundary conditions in addition:


Note: The relation between the traction and stress on the inner curve are just a little bit tricky. You have to note that if the stress at the boundary is:

σ = {{σxx, σxy},{σxy, σyy}};
MatrixForm[σ]

Then the traction vector is: t = σn, where the normal n  = cos θ i + sin θ j. So the traction vector is:

n = {Cos[θ], Sin[θ]};
MatrixForm[n]
t = σ.n;
MatrixForm[t]

So, the traction is zero at every point. Hence depending on how the normal to that surface is the relation  could be somewhat involved. We do not explicitly need to put this condition in the weak form. This is just for your information. 

It is very easy to implement all these concepts in FEniCS. We will solve problem using FEniCS and triangular element as a part of demo and also as an exercise problem. 

Example using a Quadrilateral element.




			
Figure for the problem and the FEA mesh that is  used. Note that the FEA mesh is the same as used for the heat conduction problem.
				
The problem will be solved using the following steps for a quadrilateral element. The ideas of shape functions and how to do the integration using the ξ and η coordinates are given below. 

Step-1: Gather all the material and geometric properties of the structure

The constitutive D matrix is:

The coordinate matrix is:

Step-2: Shape functions and the derivatives
The shape functions in the iso-parametric space are:

The Jacobian matrix is:

The determinant of the Jacobian and the inverse are:

The strain-displacement matrix is:

The strain-displacement matrix will be integrated using 2x2 Gauss quadrature. As noted in the previous lectures, this is not exact, but as seen there the comparison with the exact solution is quite reasonable. 

Step-3: Stiffness matrix 

The stiffness matrix for the element is now obtained as 

We calculate the stiffness matrix at each Gauss point and then sum up to get the final element stiffness matrix. 

At the Gauss point ξ = - (1/(Sqrt[3])) and η = - (1/(Sqrt[3])), the values are as follows:

Thus, the strain-displacement matrix at the Gauss point is given as:

So, the stiffness contribution coming from the Gauss point: ((Subscript[ξ, 1]), (Subscript[η, 1])) is

After the process is repeated for all the Gauss-Points you will obtain:

Step-4: Force vector

There are no body forces in the current case and as a result the only contribution to the forces come from the boundary tractions. However, only the top boundary 1-4 has non-zero traction and hence only that will contribute to it. 


From the above figure it is clear that edge 1-4 corresponds to ξ = -1 in the parent coordinates. Also note that in the current case, since the length of edge 1-4 in both parent and the material coordinates is 2: dΓ = dη. Hence, the integral for the boundary traction can be evaluated as

Note that the only shape functions that do not vanish along edge 1-4 are (Subscript[N, 1]^(4Q)) and (Subscript[N, 4]^(4Q )). These non-zero shape functions would be linear function of η. Also, interestingly, the integral of the shape function over the edge is equal to 1. Forming the boundary force matrix and also accounting for the unknown reactions at node 1 and 2 (nodes that are fixed) we obtain:

Step-5: Solve the system of equations
The system of equations is:

The reduced system of equations is

which yields

Step-6: Post-processing 
The stress and strains at the Gauss points are:


and are evaluated at individual points are

The stress that is obtained is now used to check various material failure criteria. One of the most well-known criteria for plastic criteria is by evaluating the so called von Mises stress which should be (Subscript[les, HY])s than the yield stress (Subscript[σ, Y])


Here, (Subscript[σ, 1]) and (Subscript[σ, 2])are the principal stresses and have the values as given in the second equation below. 

We will now solve this same problem in FEniCS in the Jupyter notebook and obtain the von-Mises stress.

Expression for relation between stress and strain for plane stress and plane strain

(* for plane stress *)
eq[1] = ϵxx - (((σxx)/Y)- ν ((σyy)/Y));
eq[2] = ϵyy - (- ν ((σxx)/Y)+ ((σyy)/Y));
eq[3] = ϵxy - ((σxy)/(Y/(2 (1 + ν))));

Note that shear modulus G = (Y/(2 (1 + ν))), where Y is the Young’s modulus and ν is the Poisson’s ratio for linear, isotropic, elastic material.  The compliance matrix is of the form
	ϵ = C σ
And the stiffness matrix D  = (C^(-1))

c = {{(1/Y), ((-ν)/Y), 0}, {((-ν)/Y), (1/Y), 0}, {0, 0, ((2 (1 + ν))/Y)}};
MatrixForm[c]

d = Inverse[c];
MatrixForm[%]//FullSimplify

This form is the same as given above.

(x^y) (a/b)= (1/y)
Integrate[

Integrate[Sin[x], x]

(Overscript[7, ->]) + ((((a/b)))^2)

ComplexContourPlot[AbsArg[(z^2-I)/(z^3+I)],{z,-3-3 I,3+3 I},Contours->30]

VectorPlot[{2 x^2-y^2,3 x y},{x,-5,5},{y,-5,5}]