3.6 Illustrations numériques et exemples de codes de calcul
85
12
real tmax =10;
13
real dtaff =0.2;
14
real CFL =0.2;
15
real dt =5. e -3;
16
real eps1 =0.001;
// petit parametre ( calcul normale )
17
18
// Maillage
19
mesh Th = square (n ,m ,[ xm *x , ym * y ], flags =1) ;
20
fespace Vitesses ( Th , P2 );
21
fespace Pression ( Th , P1 );
22
fespace LevelSet ( Th , P2 );
23
24
Vitesses u1 , u2 , v1 , v2 , u1n , u2n ;
25
Pression p ,q , pp ;
26
LevelSet phi , M1 , M2 , S , phiinit , nabla , N1 , N2 , zet , aux ;
27
28
// Definition de l ’ ellipse initiale
29
func phi0 = sqrt ((x - xm /2) ^2/( A ^2) +(y - ym /2) ^2/( B ^2) ) -1;
30
phi = phi0 ;
31
u1n =0; u2n =0;
32
33
// Calcul de la fonction distance a cette ellipse
34
// NB : les versions recentes de FreeFEM ++ implementent un
35
// calcul de la distance a une interface decrite par une
36
// ligne de niveau qui peut remplacer ces quelques lignes .
37
Pression h1 = hTriangle ;
38
real h = h1 []. max ; // Taille maximale d ’un triangle du maillage
39
real epsil = LARGINT * xm / n ;
40
int iterinit ;
41
real TT = CFL * h ; // pas de temps pour la re - initialisation
42
for ( iterinit =1; iterinit < 50* LARGINT ; iterinit = iterinit +1)
43
{
44
nabla =( dx ( phi )) ^2+( dy ( phi )) ^2;
45
S = phi /( sqrt ( phi ^2+ h * h * nabla )); // approximation de
signe ( phi )
46
M1 = dx ( phi ) /( sqrt ( nabla + eps1 ^2) ); // approximation de la
normale
47
M2 = dy ( phi ) /( sqrt ( nabla + eps1 ^2) );
48
phi = convect ([ - S * M1 ,-S * M2 ],TT , phi )+ TT * S ;
49
};
50
51
// Initialion de l ’ etirement
52
phi = ETIRINI * phi ;
53
nabla = sqrt (( dx ( phi )) ^2+( dy ( phi )) ^2+ eps1 ^2) ;
54
real vol0 = int2d ( Th )( phi <0) ;
55
real vol = vol0 , pslice ;
56
57
// fonction zeta
58
func real zeta ( real r ) {
59
return ((( r > -1) &&( r <1) ) ?0.5*(1+ cos ( pi * r )) :0) ;
60
}
61
62
// r --> E ’(r) loi deformation / contrainte
63
func real Ep ( real r ) {