The difference to the common model used is the f,g terms. The 3D FEM formulation has been derived and a new fortran program is written. The derivative matrix is symmetric. But the calculation result does not converge. I have tried every means to get a converging result. But it doesn't work, the residual norm always increases.