I have developed a 2D flow to solve the steady NS equations for a more complicated mesh. Once I got the code working, I have tried to test it against some analytic problem. I have used the developing region of plane Poiseuille flow (x \in[0,L] and y \in [0 H]). The inlet boundary condition is u1 = 1, u2=0, the upper and lower walls are non slip. I solve the steady problem using Newton-Raphson and once I got the results, I have realised that the inlet boundary condition is slightly modified (it is not really imposing Ux=1). I send you here a plot
I upload a test problem. The results are written in datosPlanePoiseuille.dat, where the first column is the coordinate y, the second column u1(x=0) and the third one u1(x=L). You can see that there is an overshoot near the walls at the inlet boundary condition whereas I would expect to have u1 = 1.
Your problem is discretized with 20 intervals in the y direction,
but you plot 100 values of y. It follows that the points that you plot on the boundary are just interpolation of the values 0 and 1.
If in your code you take nl = 21 then the .dat file is correct:
0 -2.683e-32 4.91667e-34
0.05 1 0.28025
0.1 1 0.531
0.15 1 0.75225
0.2 1 0.944
0.25 1 1.10625
0.3 1 1.239
0.35 1 1.34225
0.4 1 1.416
0.45 1 1.46025
0.5 1 1.475
0.55 1 1.46025
0.6 1 1.416
0.65 1 1.34225
0.7 1 1.239
0.75 1 1.10625
0.8 1 0.944
0.85 1 0.75225
0.9 1 0.531
0.95 1 0.28025
1 -2.65432e-32 4.91667e-34
The left corners are set to 0 because you impose 0 on upper and lower boundaries.
Note that the order in which you set the BC matters:
+ on(inlet, up1=1, up2=0)
+ on(upperWall, up1=0,up2=0)
+ on(lowerWall, up1=0, up2=0)
+ on(outlet,up2=0)
On the left corners you have two definitions of up1 : inlet and upper/lower wall.
The one that is applied is the last one in the list: up1=0.
If you change the order, as
+ on(upperWall, up1=0,up2=0)
+ on(lowerWall, up1=0, up2=0)
+ on(inlet, up1=1, up2=0)
+ on(outlet,up2=0)
Then the left corners take the value up1=1.
