Équation non stationnaire de la dynamique

Considérons l’équation de la dynamique, sous forme non stationnaire (i.e. dépendant du temps) : Mu.t / R C Cu.t / P C Ku.t / D F .t / (17.1) qui sous forme discrétisée par éléments finis sera réécrite : ŒM fqRg .t / C ŒC fqPg .t / C ŒK fqg .t / D fF g .t / (17.2) avec les conditions initiales fqg .t D 0/ D fqg0 et fqPg .t D 0/ D fqPg0 . On cherche à obtenir la discrétisation fqg du champ u et ses dérivées temporelles fqPg et fqRg à différents instants t et vérifiant l’équation de la dynamique. Pour résoudre un tel problème, on se propose de réaliser une discrétisation temporelle qui pourra être explicite ou implicite. Le schéma généralement utilisé pour la discrétisation temporelle est celui de Newmark, qui est présenté dans le cas de la résolution des équations différentielles au paragraphe C.2.5. 

Comparaison des méthodes explicite et implicite 

Concernant la méthode explicite, on retiendra d’abord que cette méthode nécessite moins de mémoire, est donc plus rapide et mieux adaptée aux problèmes de grandes tailles. Comme cette méthode nécessite moins de mémoire et des petits pas de temps (pour la stabilité), elle est bien adaptée au cas des chocs. Comme elle n’est que conditionnellement stable, elle est plutôt adaptée à la résolution élément par élément, donc au traitement local. Enfin, cette méthode est très robuste numériquement et permet de traiter le cas de non linéarités couplées. Concernant la méthode implicite, on retiendra que puisque cette méthode est inconditionnellement stable, elle est bien adaptée à la résolution de problèmes globaux (qui nécessitent qu’il y ait convergence). De plus, cette méthode est bien moins robuste numériquement que la précédente (pivots nuls, divergence)… Enfin et surtout, on se souviendra que c’est une méthode coûteuse en mémoire et en temps… La figure 17.1 donne, de manière imagée, les domaines d’utilisation des méthodes implicite et explicite

Exemple : un calcul de propagation avec FREEFEM++

Quelque soit le logiciel utilisé 1 , il y a toujours un grand nombre d’exemples traités fournis avec, sans compter le nombre d’ouvrages disponibles. Nous avons souhaité néanmoins présenter un listing (qui sera discuté oralement) correspondant à un calcul réalisé sous FREEFEM++. Ce logiciel gratuit vous permettra de faire vos calculs par éléments finis… mais c’est également un outil pédagogique sans équivalent : sa manipulation se fait en des termes proches de la formulation mathématique, et en cela, il est une passerelle parfaite entre la théorie et la pratique de codes plus « industriels ». Le fait de pouvoir discrétiser à loisir une formulation variationnelle que l’on rentre soi-même est également un argument qui nous semble un vrai plus, même si cela commence à exister dans des codes commerciaux (mais reste confidentiel). 2 Voici donc un petit listing sur lequel nous pourrons revenir de vive voix. Commençons donc, comme il se doit, par définir la géométrie du problème, paramétrée par a et b, que nous maillons :
real mu =1.; 2 real rog =0. , ro =1.; 3 / / t a i l l e du m a i l l a g e 4 int ncuts = 20; / / c o n s t r u c t i o n du d o m ai n e 7 real a =pi *0.8 , b =1.2* pi ;  border gam1 (t =0. , a ){x= t;y =0; label =1;}; 10 border gam2 (t =0. , b ){x= a;y= t; label =2;}; 11 border gam3 (t=a ,0.) {x= t;y= b; label =3;}; 12 border gam4 (t=b ,0.) {x =0.; y=t; label =4;}; 13 14 / / M a i l l a g e 15 mesh Th = buildmesh ( gam1 ( ncuts )+ gam2 ( ncuts )+ gam3 ( ncuts )+ gam4 ( ncuts )) ; 16 / / p l o t ( Th , p s = » me sh . e p s  » ) ; 17 / / p l o t ( Th , w a i t = 1 ) ; Nous définissons quelques valeurs afin de pouvoir définir le chargement (le chargement défini par la fonction w1 ne sera pas utilisée dans le calcul effectué). Notons comme il est aisé de définir des fonctions : 
/ / c o n d i t i o n s a u x l i m i t e s e s s e n t i e l l e s 
real r0 =a /10. , r , xc = a /3. , yc =3.* b /4.; 
/ / c h a r g e m e n t i n i t i a l s i n u s o i d a l 
 func real pertb ( real r) 23 { if(r < r0 ) { return (1.+ cos (pi* r/ r0 )) ;} 
 else return 0.;} 25 26 / / f u n c r e a l w1 ( r e a l x , r e a l y ) 
 func real w11 ( real x, real y) 
 { 
 return pertb ( sqrt ((x- xc ) ^2+(y- yc ) ^2) ); 
 } 
/ / c h a r g e m e n t c o n s t a n t s u r un c a r r e 
 / / p a s u t i l i s e s u r c e t e x e m p l e

func real w1 ( real c , real d )
{
if(c >(0.8* xc ) ) 37 {if (c <(1.2* xc )) {if (d >(0.8* yc )) 39 {if (d <(1.2* yc )) { return 0.8;}
else return 0.;}
else return 0.;}
else return 0.;}
else return 0.;
}
func real w0 ( real x, real y)
{ 4
return 0.* sin(x)* sin (y) ;
}

Ensuite, nous définissons l’espace des éléments finis utilisés, et les variables appartenant à ces espaces. L’analogie avec la « formulation mathématique » saute au yeux… et c’est pourquoi cet outil nous semble particulièrement intéressant sur le plan pédagogique. 

 / / E s p a c e s e l e m e n t s f i n i s 

 fespace Vh (Th , P2 ); 

 fespace Wh (Th , P1dc ); 

 Vh w , wa , wd , wda , wdd , wdda , fi ; 

 Wh sxz , syz ; 

Il s’agit d’un exemple non stationnaire, nous allons donc faire une boucle sur le temps. Dans FREEFEM++, nous entrons directement la formulation variationnelle du problème sous forme « mathématique »… le lien entre pratique des éléments finis et théorie est beaucoup plus clair. La boucle de temps est basée sur un schéma de Newmark implicite. Le temps « courant » est le temps t C t et le temps précédent est le temps t pour rester cohérent avec les notations du chapitre 17. Il s’agit d’un problème de déplacement hors plan, dont le champ inconnu est traditionnellement noté w. Ainsi, dans le listing, w représente qtCt , wa qt , wd PqtCt , wda Pqt , wdd RqtCt et wdda Rqt . fi ‘ est la fonction test.