86
3 Interaction d’un fluide incompressible avec une membrane élastique
64
return ( lambda *(r -1) );
65
}
66
67
// Probleme variationnel
68
problem IBM ([ u1 , u2 , p ] ,[ v1 , v2 , q ]) =
69
// Navier - Stokes
70
int2d ( Th )( u1 * v1 / dt + u2 * v2 / dt
71
+1/ Re *( dx ( u1 )* dx ( v1 )+ dy ( u1 )* dy ( v1 ) + dx ( u2 )* dx ( v2 )
+ dy ( u2 )* dy ( v2 )))
72
+ int2d ( Th )( dx ( p )* v1 + dy ( p )* v2 + q *( dx ( u1 )+ dy ( u2 )
) - 1e -10* p * q )
73
- int2d ( Th )( convect ([ u1n , u2n ],-dt , u1n )/ dt * v1 +
convect ([ u1n , u2n ],-dt , u2n )/ dt * v2 )
74
// Tenseur elastique eulerien
75
+ int2d ( Th )( Ep ( nabla ) *( dy ( phi )* dy ( phi )* dx ( v1 ) -dx (
phi )* dy ( phi ) *( dx ( v2 )+ dy ( v1 ))+ dx ( phi )* dx ( phi )*
dy ( v2 ))/ nabla * zet )
76
// Conditions aux limites
77
+ on (1 ,2 ,3 ,4 , u1 =0 , u2 =0) ;
78
79
80
// Boucle principale
81
real t =0 , taff =0;
82
for ( t =0; t < tmax ; t += dt ) {
83
cout << "t = " << t << " Volume variation (%) = " << ( vol -
vol0 )/ vol0 *100 << endl ;
84
// Resolution de Navier - Stokes
85
IBM ;
86
u1n = u1 ; u2n = u2 ;
87
// Transport de phi
88
aux = convect ([ u1 , u2 ],-dt , phi );
89
phi = aux ; vol = int2d ( Th )( phi <0) ;
90
nabla = sqrt (( dx ( phi )) ^2+( dy ( phi )) ^2+ eps1 ^2) ;
91
zet =1./( nabla * epsil )* zeta ( phi /( nabla * epsil ));
92
93
if (t > taff ) {
94
pp =p - p []. min ;
95
string vtkout = " resultats / memb_t ="+ t + ". vtk ";
96
savevtk ( vtkout , Th , pp , phi ,[ u1 , u2 ,0] , dataname = " Pressure
LevelSet Vitesse ");
97
taff = taff + dtaff ;
98
Th = adaptmesh ( Th , p );
99
}
100
}
Les figures ci-dessous correspondent à l’exécution de ce programme FreeFEM++ jusqu’à t = 0.4 ainsi que le profil de pression y = 1 associé.
Un code différences finies a aussi été implémenté en Matlab, et disponible en
ligne 2 . Il est construit à partir d’un code Navier-Stokes de Benjamin Seibold,
qui implémente une méthode de projection de type Chorin. Le paquetage
Level-Set de Baris Sumengen est aussi utilisé, où nous avons actualisé un
2. http://level-set.imag.fr
3 Interaction d’un fluide incompressible avec une membrane élastique
64
return ( lambda *(r -1) );
65
}
66
67
// Probleme variationnel
68
problem IBM ([ u1 , u2 , p ] ,[ v1 , v2 , q ]) =
69
// Navier - Stokes
70
int2d ( Th )( u1 * v1 / dt + u2 * v2 / dt
71
+1/ Re *( dx ( u1 )* dx ( v1 )+ dy ( u1 )* dy ( v1 ) + dx ( u2 )* dx ( v2 )
+ dy ( u2 )* dy ( v2 )))
72
+ int2d ( Th )( dx ( p )* v1 + dy ( p )* v2 + q *( dx ( u1 )+ dy ( u2 )
) - 1e -10* p * q )
73
- int2d ( Th )( convect ([ u1n , u2n ],-dt , u1n )/ dt * v1 +
convect ([ u1n , u2n ],-dt , u2n )/ dt * v2 )
74
// Tenseur elastique eulerien
75
+ int2d ( Th )( Ep ( nabla ) *( dy ( phi )* dy ( phi )* dx ( v1 ) -dx (
phi )* dy ( phi ) *( dx ( v2 )+ dy ( v1 ))+ dx ( phi )* dx ( phi )*
dy ( v2 ))/ nabla * zet )
76
// Conditions aux limites
77
+ on (1 ,2 ,3 ,4 , u1 =0 , u2 =0) ;
78
79
80
// Boucle principale
81
real t =0 , taff =0;
82
for ( t =0; t < tmax ; t += dt ) {
83
cout << "t = " << t << " Volume variation (%) = " << ( vol -
vol0 )/ vol0 *100 << endl ;
84
// Resolution de Navier - Stokes
85
IBM ;
86
u1n = u1 ; u2n = u2 ;
87
// Transport de phi
88
aux = convect ([ u1 , u2 ],-dt , phi );
89
phi = aux ; vol = int2d ( Th )( phi <0) ;
90
nabla = sqrt (( dx ( phi )) ^2+( dy ( phi )) ^2+ eps1 ^2) ;
91
zet =1./( nabla * epsil )* zeta ( phi /( nabla * epsil ));
92
93
if (t > taff ) {
94
pp =p - p []. min ;
95
string vtkout = " resultats / memb_t ="+ t + ". vtk ";
96
savevtk ( vtkout , Th , pp , phi ,[ u1 , u2 ,0] , dataname = " Pressure
LevelSet Vitesse ");
97
taff = taff + dtaff ;
98
Th = adaptmesh ( Th , p );
99
}
100
}
Les figures ci-dessous correspondent à l’exécution de ce programme FreeFEM++ jusqu’à t = 0.4 ainsi que le profil de pression y = 1 associé.
Un code différences finies a aussi été implémenté en Matlab, et disponible en
ligne 2 . Il est construit à partir d’un code Navier-Stokes de Benjamin Seibold,
qui implémente une méthode de projection de type Chorin. Le paquetage
Level-Set de Baris Sumengen est aussi utilisé, où nous avons actualisé un
2. http://level-set.imag.fr