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:
- 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.
- 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.
