I am trying to solve a multiphysic problem using finite elements and a Newton Raphson solution scheme. I have two non-linear subsystems that are coupled bi-directionally.
The first subsystem includes an equation that I solve for $p_f$:
$\nabla u - \nabla a^3 \nabla p_f = \nabla a^3 B$
This equation becomes very non-linear, due to the parameter $a$ with an exponent of 3! $B$ are body forces and $u$ is the solution of the second subsystem that I am trying to solve. This second subsystem includes the Stokes equation that depends on the gradient of $p_f$. Therefore, both system are coupled bi-directionally.
I am using a Newton-Raphson solution scheme for each subsystem. First, I solve for $p_f$. After convergence occurred, I pass this solution to the second subsystem, which I solve using a Newton Raphson scheme, too. Then, the solution $u$ is passed to the first subsystem ... until the solutions of both subsystems do not change anymore. See the flow chart below.
This procedure works very well so far. However, if the values of the parameter $a$ become very low (thus the gradient of $a$ becomes very, very low), there is no convergence anymore! The reason seems to be that the coupling between both subsystem becomes very strong. A small change in one subsystem creates a significant changes in the other one (especially where the gradient of $a$ is very low). I tried a variety of good initial conditions, however, since the first subsystem reacts very sensitiv, convergence has been impossible so far.
In the figure below, you can see how the solution of this approach ("Segregated Approach"-blue line) is different from the "true" solution ("Fully Coupled"-red line). The "true solution" is the result of a benchmark, which I calculated by solving a single system of equations that includes the equations of both subsystems.
The greatest differences between both solutions are located where the gradient of $a$ becomes very small. See the figure below (in this figure the values of $p_f$ and $a$ are normalized to range between 0 and 1):
So: Is a Newton-Raphson solution scheme the wrong method for solving a strongly coupled system, if there are low gradients in a subsystem? Could my method be fixed to solve the problem that I described? If yes, how? Are there alternative solution schemes for this problem?