Hele shaw flow problem

Dear developer:
I’m solving for a hele shaw flow over a rectangular region (2×1). The boundary conditions on all four edges are no-flux boundary conditions. The code I developed fails to converge.Note that If the nonlinear term vr terms in the following equation is removed, the result agree well with the comsol result. My question is that how to handle vr terms so that the code converge. I would appreciate that if your could kindly help me to solve this problem.

I would appreciate that if your could kindly help me

The specific equation is shown below:
image

Where Tr is 0.95.
The initial conditions are as follows:

Some key codes are as follows:

varf form

image

Solving procedure

image

See attachment for the complete code.
The results of my simulation using comsol are shown below:

image
image
image

Best wish,
Guangyuan Wei
modify41.edp (2.0 KB)

I have tried to resolve the nonlinearity v_r(\rho) implicitly by iteration:
\int_\Omega\omega\frac{\rho^{(k+1)}-\rho^n}{dt} +\int_\Omega\nabla\omega\cdot\left(\frac{\rho^{(k)2}}{\mu_r(\rho^{(k)})}\nabla\Pi^{(k+1)}\right) +\int_\Omega\chi\bigl(v_r(\rho^{(k)})+r(\rho^{(k+1)}-\rho^{(k)})-\Pi^{(k+1)}\bigr)
+\int_\Omega\Lambda\nabla\chi\cdot\nabla\rho^{(k+1)}=0
where k is an iteration index. Then \rho^{(k)} is supposed to converge to \rho^{n+1} solving the implicit problem. The constant r has to be taken large enough, in particular larger than \sup v_r'(\rho), where v_r' stands for the derivative of v_r with respect to \rho. Moreover r must be greater than 6.*(1.-Tr) which is minus the minimum of v_r'.
Warning: when you write ( 8 / 3 ), the quotient is taken in the sense of integers, giving the result of 2. You must write ( 8. / 3. ).
Here is the code, that seems to work
heleshaw_new.edp (2.8 KB)

Dear François Bouchut,
Thank you very much for the great help you have given me in my difficulties. I will be greatly appreciated. I will consider your suggestions and solutions in depth.

Best wish,
Guangyuan Wei