1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
/*
Incompressible Navier Stokes
with Taylor-Hood Finite element
No linearity : Newton methode
continuation on Reynols Number
Mesh adaptation
*/
real reymax = 1600; // ok < 125000
func BCu1=4*x*(1-x);
mesh Th=square(8,8);
fespace Xh(Th,P2);
fespace Mh(Th,P1);
fespace XXMh(Th,[P2,P2,P1]);
XXMh [u1,u2,p];
XXMh [v1,v2,q];
real epsr=1e-6;
macro div(u1,u2) (dx(u1)+dy(u2))//
macro grad(u1,u2) [dx(u1),dy(u2)]//
macro ugrad(u1,u2,v) (u1*dx(v)+u2*dy(v)) //
macro Ugrad(u1,u2,v1,v2) [ugrad(u1,u2,v1),ugrad(u1,u2,v2)]//
solve Stokes ([u1,u2,p],[v1,v2,q],solver=UMFPACK) =
int2d(Th)( ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
+ dx(u2)*dx(v2) + dy(u2)*dy(v2) )
- p*q*epsr
- p*div(v1,v2)-q*div(u1,u2)
)
+ on(3,u1=BCu1,u2=0)
+ on(1,2,4,u1=0,u2=0);
Xh uu1=u1,uu2=u2;
plot(coef=0.2,cmm=" [u1,u2] et p ",p,[uu1,uu2],wait=1);
Xh psi,phi;
solve streamlines(psi,phi) =
int2d(Th)( dx(psi)*dx(phi) + dy(psi)*dy(phi))
+ int2d(Th)( -phi*(dy(u1)-dx(u2)))
+ on(1,2,3,4,psi=0);
plot(psi,wait=1);
int i=0;
real nu=1./100.;
real dt=0.1;
real alpha=1/dt;
/* NL
varf vNS ([u1,u2,p],[v1,v2,q],solver=Crout,init=i) =
int2d(Th)(
alpha*( u1*v1 + u2*v2)
+ nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
+ dx(u2)*dx(v2) + dy(u2)*dy(v2) )
- p*q*eprs
- p*div(v1,v2)-q*div(u1,u2)
+ Ugrad(u1,u2,u1,u2)'*[v1,v2]
)
+ on(3,u1=1,u2=0)
+ on(1,2,4,u1=0,u2=0)
Newtow DF(u) u1 = DF(u)u - F(u)
*/
XXMh [up1,up2,pp];
varf vDNS ([u1,u2,p],[v1,v2,q]) =
int2d(Th)(
+ nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
+ dx(u2)*dx(v2) + dy(u2)*dy(v2) )
- p*q*epsr
- p*div(v1,v2)-q*div(u1,u2)
+ Ugrad(u1,u2,up1,up2)'*[v1,v2]
+ Ugrad(up1,up2,u1,u2)'*[v1,v2]
)
+ on(3,u1=BCu1,u2=0)
+ on(1,2,4,u1=0,u2=0);
;
varf vNS ([u1,u2,p],[v1,v2,q]) = // DF(u)u - F(u)
int2d(Th)(
Ugrad(up1,up2,up1,up2)'*[v1,v2]//'
)
+ on(3,u1=BCu1,u2=0)
+ on(1,2,4,u1=0,u2=0);
;
for(real re=100;re<=reymax;re *=2)
{
real lerr=0.1;
if(re>8000) lerr=0.05;
if(re>10000) lerr=0.01;
for(int step=0;step<2;step++)
{
Th=adaptmesh(Th,[u1,u2],p,err=lerr,nbvx=100000,abserror=0, cutoff=0.01);
//plot(Th,wait=0);
[u1,u2,p]=[u1,u2,p];
[up1,up2,pp]=[up1,up2,pp];
for (i=0;i<=20;i++)
{
nu =1./re;
up1[]=u1[];
real[int] b = vNS(0,XXMh);
matrix Ans=vDNS(XXMh,XXMh);
set(Ans,solver=UMFPACK);
u1[] = Ans^-1*b;
b=u1[]-up1[];
cout << " iter = "<< i << " " << b.l2 << " rey = " << re << endl;
if(b.l2<1e-6) break;
// uu1=u1;uu2=u2;
//plot(coef=0.2,cmm=" [u1,u2] et p ",p,[uu1,uu2]);
} ;
}
uu1=u1;uu2=u2;
streamlines;
plot(coef=0.2,cmm="rey="+re+" [u1,u2] et p ",psi,[uu1,uu2],wait=0,nbiso=20,ps="cavity-"+re+".ps");
}
|
