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
|
mesh Th=square(8,8);
fespace Xh(Th,P2);
fespace Mh(Th,P1);
Xh u2,v2;
Xh u1,v1;
Mh p,q;
real epsr = 1e-8;
solve Stokes ([u1,u2,p],[v1,v2,q]) =
int2d(Th)( ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
+ dx(u2)*dx(v2) + dy(u2)*dy(v2) )
- p*q*epsr
+ p*dx(v1)+ p*dy(v2)
+ dx(u1)*q+ dy(u2)*q
)
+ on(3,u1=1,u2=0)
+ on(1,2,4,u1=0,u2=0);
plot(coef=0.2,cmm=" [u1,u2] et p ",p,[u1,u2],ArrowSize=0.5,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;
Xh up1,up2;
problem NS ([u1,u2,p],[v1,v2,q],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*epsr
+ p*dx(v1)+ p*dy(v2)
+ dx(u1)*q+ dy(u2)*q
)
+ int2d(Th) ( -alpha*convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 )
+ on(3,u1=1,u2=0)
+ on(1,2,4,u1=0,u2=0)
;
for (i=0;i<=20;i++)
{
up1=u1;
up2=u2;
NS;
if ( !(i % 10))
plot(coef=0.2,cmm=" [u1,u2] et p ",p,[u1,u2]);
} ;
if ( (i % 10))
plot(coef=0.2,cmm=" [u1,u2] et p ",p,[u1,u2]);
streamlines;
plot(psi,wait=1);
|
