Boundary condtions on nonlinear FEM time integration

I'm using the finite element method to obtain the time response of a structure harmonically excited. I'm using a linear displacement function to obtain the stiffness matrix and the consistent mass matrix, i.e., I have this set of equations:

$$[M]\{\ddot x\} + [K]\{x\} = \{F(t)\}$$

What is the best way to apply the boundary conditions on the $[M]$ and $[K]$ matrices?

Typically, I have found two ways:

1. Remove the rows/columns whose degree of freedoms are constrained
2. Modify the matrices setting to zero the rows and columns and placing ones on the diagonal to maintain non-singularity.

I prefer the second way, since we can adjust it to account for non zero displacement constrains.

The problem here is that I'm using the Newmark method to integrate the system and the biggest difficult is that we are dealing with geometrical nonlinearities, which means $[K]$ is the tangent stiffness matrix. Since this is a nonlinear problem, we modify the Newmark method to account for the nonlinearities and use the Newton-Raphson method to find the equilibrium position at time $t(i+1)$ when marching forward in time. Usually, the Newton-Raphson method is fast enough to find the equilibrium position. However, if the initial estimate is not good enough, it will not converge.

The biggest problem here is that, when modifying the matrices, an error residue on any degree of freedom found by the Newton-Raphson is propagated and the final answer diverges. Any ideas on how to go around this problem?

I'm trying to keep it simple, without too much details on the Newmark method or the FEM method. If necessary, however, we can discuss it.

• The behavior you describe is very unusual because the solution at t(i) is typically a relatively close approximation to that at t(i+1). Are you sure that your tangent K is equal to the gradient of your internal force vector? As a test of this, I usually compute the tangent K by applying finite differences to the internal forces and then compare that to the analytical tangent K. – Bill Greene Oct 26 '16 at 17:29
• Good idea. I checked the tangent K with a few examples I had, but did not tried to calculate the gradient from the internal forces. Apparently, it is working, but it may be just for the example I have. I will try what you suggest. – Thales Oct 30 '16 at 2:45
• If you using Newmark method or not should not make a problem. It does nor matters as well whether you do (1) or (2), the solution of the system of equations should be the same. You have to set initial boundary conditions consistently with applied boundary conditions when you calculating residual. I prefer (2), for its simplicity. More interesting is how you modify Newmark for geometry nonlinear problem? – likask May 25 '17 at 22:25
• @likask I saw it in a discipline on my university. I have the algorithm on my class notes, a little big to write here. email me and when I have some time, I will type it for you (my manuscript is in portuguese). However, you can check this paper. It is not exactly what I did, but it gives enough references to do dynamic analysis on nonlinear structures. – Thales May 26 '17 at 13:08
• @Thales If I understand correctly you going to use Newmark method as it is, but modify approximation base with some hierarchy of bubble functions, i.e. higher order finite elements. I did similar thing for 3d mofem.eng.gla.ac.uk/mofem/html/elastic_dynamics.html However, what I did is not perfect, it would be much better if base functions only takes positive values, like NURBS or better Bernstein-bézier, speeding up assemble and removing problems with negative values in mass matrix, look at the second reference in mofem.eng.gla.ac.uk/mofem/html/citelist.html – likask May 26 '17 at 14:33

Have you tried using SOR (successive Over Relaxation scheme). Newton Raphson will face sever issue if one of the following is occuring in your case. Practicle Considerations for Newton Raphson.

• I did not tried that. I will try to implement and see if it helps. Thanks for the great idea. – Thales Oct 30 '16 at 2:43

What I would do (have done) if I had to solve a nonlinear time stepping problem..

1. Do time discretization in the strong form of PDE (use Newmark or whatever you want to)

2. Find the weak form of the time-discretized strong form of the PDE to find $$u^{n+1}$$. You know $$u^n$$ and $$u^{n-1}$$ from initial conditions.

3. Linearize the nonlinear system for $$u^{n+1}$$.

4. For each Time step to find $$u^{n+1}$$ - Take a newton raphson guess (preferably $$u^n$$, so your guess is very close to your solution and don't diverge) and impose the Dirichlet boundary condition on it. Since you are only solving for free degrees of freedom, you will not touch the solution at the dirichlet nodes. If you have Neumann condition, it will occur naturally in the weak form.