# FEA: Newton-Raphson algorithm for Dirichlet BC non-linear static analysis

I have a standard non-linear static problem to solve with FEA. I am applying displacement boundary conditions only, no external force. The equation to solve look like these:

$$K(u)u=R(u)$$ Eq(1)

The terms $$K, R$$ are the stiffness matrix and residual vector(in this case internal force vector), in any standard FEA formulation, $$u$$ is the field variable(displacement). The relationship between $$K, R$$ can be given as per standard FEM definitions, $$\partial F/ \partial U=R$$(obtained from Taylor series expansion). Due to the non-linear nature, the equations cannot be solved directly, an incremental iterative approach is needed.

I am confused on how to set up the NR algorithm. All examples I see are breaking up the external load into steps, and run iterations in each step to check for the convergence $$Residual=External Force-Internal Force= 0$$. But I don't have an external force vector, then how should I proceed with the incremental steps? Moreover, in this case, $$External Force=0$$ so how should the convergence be checked in each iteration step?

• If you don't have any external force and also no displacement boundary conditions, the internal force is zero and nothing happens. If you have a displacement BC, that creates an external force. You can check out my notes at en.wikiversity.org/wiki/Nonlinear_finite_elements Mar 18 '19 at 19:49
• I do have an applied displacement. For e.g., if we consider a simple pull of a bar, one edge is fixed, on the other edge, I apply displacements, instead of force. However, I am not sure how to check for convergence and apply the constraints in this case. Thanks for the link, let me see if I can find what I am looking for there. Mar 18 '19 at 20:11
• The following link was very helpful, it answered some of my questions. However, I have a few more. How do you define the residual? Is it the same definition as I mentioned, then external force vector should be equal to zero? Secondly, why is $u_{r+1}=u_r- \delta u$ Shouldn't it be $u_{r+1}=u_r+ \delta u$? en.wikiversity.org/wiki/Nonlinear_finite_elements/… Mar 18 '19 at 21:17
• The negative sign comes from Newton's method. See commons.wikimedia.org/wiki/File:NewtonIteration_Ani.gif Mar 18 '19 at 21:41
• Well, I understood that, but it got me confused a bit. I have seen many cases where its calculated as $u_{r+1}=u_r+\delta u$. Probably depends on how they calculate the residual, $f_{int}-f_{ext}$ or $f_{ext}-f_{int}$ ?? Mar 18 '19 at 22:03

The discretized problem that you are trying to solve is $$\text{find}\,\, \mathbf{u}^{n+1} \,\, \text{such that} \,\, \mathbf{r}(\mathbf{u}^{n+1}) = \mathbf{0} \,\, \text{subject to the constraints} \,\, \mathbf{g}(\mathbf{u}^{n+1}) = \mathbf{0}$$ where $$\mathbf{u}^{n+1}$$ is the displacement at the end of load step $$n$$ and the residual is $$\mathbf{r}(\mathbf{u}^{n+1}) := \mathbf{f}^{\text{int}}(\mathbf{u}^{n+1}) - \mathbf{f}^{\text{ext}}(\mathbf{u}^{n+1}) \,.$$ In the absence of external forces, we have $$\mathbf{r}(\mathbf{u}^{n+1}) = \mathbf{f}^{\text{int}}(\mathbf{u}^{n+1}) \,.$$ Therefore, at the end of the load step, the displacements should be such that the internal forces are zero.

Newton's method can be used to find the values of $$\mathbf{u}^{n+1}$$ at which $$\mathbf{r} = \mathbf{0}$$.

We start with the solution at the beginning of load step $$n$$: $$\mathbf{u}^{n+1}_0 = \mathbf{u}^n$$

(Caveat: For implicit dynamic computations with the Newmark-$$\beta$$ method, a better initial guess is $$\mathbf{u}^{n+1}$$.)

As the iterations proceed, we will get better and better estimates of the quantity $$\mathbf{u}^{n+1}_k$$ where $$k$$ is the iteration number.

A Taylor series expansion of the residual about the current value of the displacement leads to $$\mathbf{0} = \mathbf{r}(\mathbf{u}^{n+1}_k) = \mathbf{r}(\mathbf{u}^{n+1}_{k-1}) + \frac{\partial \mathbf{r}(\mathbf{u}^{n+1}_{k-1})}{\partial \mathbf{u}} \cdot (\mathbf{u}^{n+1}_{k} - \mathbf{u}^{n+1}_{k-1}) + \dots$$ Reorganizing, we get $$\mathbf{u}^{n+1}_{k} = \mathbf{u}^{n+1}_{k-1} - \left[\frac{\partial \mathbf{r}(\mathbf{u}^{n+1}_{k-1})}{\partial \mathbf{u}}\right]^{-1} \cdot \mathbf{r}(\mathbf{u}^{n+1}_{k-1})$$ In the absence of external forces, we have \begin{align} \mathbf{u}^{n+1}_{k} &= \mathbf{u}^{n+1}_{k-1} - \left[\frac{\partial \mathbf{f}^{\text{int}}(\mathbf{u}^{n+1}_{k-1})}{\partial \mathbf{u}}\right]^{-1} \cdot \mathbf{f}^{\text{int}}(\mathbf{u}^{n+1}_{k-1}) \\ & = \mathbf{u}^{n+1}_{k-1} - \left[\mathbf{K}^{\text{int}}\right]^{-1} \cdot \mathbf{f}^{\text{int}}(\mathbf{u}^{n+1}_{k-1}) \end{align} The quantity $$\mathbf{K}^{\text{int}}$$ is called the tangent stiffness matrix. We continue iterating until a convergence criterion is satisfied - typically some sort of the energy norm computed from force and displacement.

Update:

1. Constraints:

I haven't talked about the process of applying the displacement constraints in the above procedure. There are several ways of applying these constraints. The most easily understood approach is to use Lagrange multipliers. A more common approach is to use a penalty method.

1. Computing the tangent stiffness:

The tangent stiffness matrix can be derived in several ways. For example, in the Belytchko et al. book on nonlinear FE, they proceed using the convected rate of the Kirchhoff stress. Simo's book on computation inelasticity tends to focus on algorithmically consistent tangent stiffness matrices. Deriving the tangent stiffness for complex material models and implementing it without errors is one of the major difficulties experience by authors of FE codes. Very few commercial codes get everything right.

• Thanks for the detailed explanation. Is the tangent stiffness matrix same as the stiffness matrix calculated with $u_{k-1}^{n+1}$ ? Mar 20 '19 at 18:43
• @Sauradeep: The tangent stiffness matrix has two terms - the material tangent stiffness and the geometric tangent stiffness. A detailed explanation can be found in many books. However, I've added a brief update to my answer to explain the basic idea. Mar 20 '19 at 22:47
• Thanks. I accepted your answer. However regarding the enforcement of constraints, can't it be done by using the standard method? In each iteration, I enforce the two conditions, displacement fixed($0$) on some nodes and applied value($\delta u$) in some nodes? The applied value will change with each incremental step. I just replace the corresponding rows in the $f^{int}$ vector with the values and make all non-diagonal terms in $K^{int}=zero$ and the diagonal terms $=one$ for each iteration. I don't understand why not. Mar 21 '19 at 2:30
• The approach one uses depends upon how general the code needs to be. For homogeneous displacement BCs you can use your approach. The values of diagonal terms have to be chosen such that the conditioning of the matrix does not suffer. Mar 21 '19 at 2:38
• It should work for inhomogenous BCs also, right? I have definite non-zero values for some nodes, which is basically the non-zero prescribed displacement. If I modify the corresponding rows of the internal force vector and update the corresponding stiffness matrix terms in each iteration, that should give me correct results. Mar 21 '19 at 3:29