Télécharger la plaquette Cast3M


fr en

Télécharger ale_mecaflu.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  ale_mecaflu.dgibi
  2. 'OPTI' echo 0 ;
  3. ************************************************************************
  4. ************************************************************************
  5.  
  6. ***************************************************************
  7. *                                                             *
  8. * NOM         : ale_mecaflu.dgibi                             *
  9. * DESCRIPTION : fiche de validation CASTEM2000                *
  10. *               Mécanique des Fluides                         *
  11. *               Equations de Navier-Stokes en description ALE *
  12. *               (Arbitraire Lagrange-Euler)                   *
  13. *               géométrie initiale : cavité rectangulaire (2D)*
  14. * FONCTIONS                                                   *
  15. *   TESTEES   : 'NS' option 'CENTREE' & 'ALE' incompressible  *
  16. * (entre             formulation non conservative             *
  17. *   autres)          éléments QUA4                            *
  18. *               'FORME' pour les déplacements du maillage     *
  19. *               'LAPN' option 'EF' 'IMPL'                     *
  20. *                      pour le calcul de la vitesse du        *
  21. *                      maillage                               *
  22. *                                                             *
  23. * AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)      *
  24. *               e-mail : gounand@semt2.smts.cea.fr            *
  25. *                                                             *
  26. ***************************************************************
  27. *                                                             *
  28. * VERSION     : v1, 29/12/97, version initiale                *
  29. * HISTORIQUE  : v1, 29/12/97, création                        *
  30. * HISTORIQUE  :                                               *
  31. * HISTORIQUE  :                                               *
  32. * HISTORIQUE  :                                               *
  33. *                                                             *
  34. ***************************************************************
  35. *                                                             *
  36. * Priere de completer la partie HISTORIQUE apres modification *
  37. * du jeu de donnees afin de faciliter la maintenance          *
  38. *                                                             *
  39. ***************************************************************
  40. *
  41. * DESCRIPTION DETAILLEE :
  42. *
  43. * Géométrie :
  44. * ---------
  45. *   A l'instant initial, le fluide est au repos dans une
  46. * cavité rectangulaire LxH ouverte (comme une gouttière).
  47. * Le fluide est newtonien,visqueux et incompressible.
  48. * Il est soumis à l'accélération de la pesanteur.
  49. *      
  50. *    z
  51. *      ^
  52. *      |
  53. * 
  54. *      |                |
  55. *      |                |                
  56. *     R|----------------|Q
  57. *      |                |    |
  58. *      |                |    |  g
  59. *      |H               |    |
  60. *      |     L          |    \/
  61. *      ------------------ - - - - - -> x
  62. *     O                  P
  63. * On modélise son comportement par les équations de
  64. * Navier-Stokes en description ALE (les noeuds du maillage
  65. * sont mobiles). On désigne par v, la vitesse absolue du
  66. * fluide aux noeuds du maillage et par w, la vitesse absolue
  67. * de déplacement des noeuds du maillage.
  68. * Dans la description ALE, seul le terme de convection des
  69. * équations de Navier-Stokes change, il s'écrit :
  70. *   (v-w).grad(v)
  71. * Si w=0, on retrouve la description eulérienne classique.
  72. * Si w=v, le terme de convection disparait => description
  73. * lagrangienne.
  74. * Dans le cas général, il nous faut une équation qui precrit w.
  75. * 
  76. * Condition initiale :
  77. * ------------------
  78. * A l'instant t=0, on imprime un pulse de pression sur la
  79. * surface libre (RQ) de la forme :
  80. *           P(t) = A Dirac(t) cos(pi x / L)
  81. * Celle-ci se met à osciller librement.
  82. * 
  83. * Conditions aux limites :
  84. * ----------------------
  85. * - contrainte nulle sur la surface ;
  86. * - v.n = 0 sur les bords (il faut que Q et R puissent se
  87. * déplacer).
  88. * 
  89. * Prescription de w :
  90. * -----------------
  91. * Trois facons de prescrire w ont été implémentés
  92. * (voir aussi plus loin les procédures correspondantes) :
  93. * - LAGR : lagrangienne : w=v
  94. * - SLAG : surface lagrangienne : w|surf=v|surf   
  95. *                            et   Laplacien(w)=0 sur le
  96. *                            domaine.
  97. * - CVS2LF : on impose w // Oz ie les points du maillage
  98. * doivent se déplacer uniquement verticalement.
  99. * Ceci nous amène à résoudre une équation de convection 1D sur la
  100. * surface qui est discrétisée en DF centrées <=> EF centrés.
  101. * 
  102. * Tests effectués :
  103. * ---------------
  104. * - on prescrit w par CVS2LF (le plus compliqué) ;
  105. * - on vérifie que le volume se conserve (les schémas utilisés
  106. * ne conservent PAS, a priori, le volume) ;
  107. * - on compare les ordonnées des points Q et R en fonction du
  108. * temps (zQ(t) et zR(t)) avec les profils obtenus par une
  109. * autre simulation numérique (Ramaswamy).
  110. *
  111. *                         Voir rapport CEA/DMT 97/371
  112. *                               pour plus de précisions
  113. *
  114. *************************************************************
  115. * Références : - Gounand, S. 1997
  116. *                Simulation Numérique d'écoulements à surface
  117. *                libre.
  118. *                Rapport CEA/DMT 97/371.
  119. *              - Ramaswamy, B. 1990
  120. *                Numerical simulation of unsteady viscous 
  121. *                free surface flow
  122. *                J. Comp. Phys. Vol. 90 pp. 396-430. 
  123. *              - Ramaswamy, B. & Kawahara, M. 1987
  124. *                Lagrangian finite element analysis applied
  125. *                to viscous free surface fluid flow.
  126. *                Int. J. for Num. Meth. in Fluids 
  127. *                Vol. 7 pp. 953-984. 
  128. *                                                             
  129. ***************************************************************
  130. *
  131. * JEU DE DONNEE :
  132. *
  133. prg = ('CHAINE' 'Non Linear Oscillations (NS ALE)') ;
  134. **** Options
  135. * court = VRAI  : le cas-test s'effectue en une trentaine
  136. *                 de secondes.
  137. * court = FAUX  : le cas-test est plus long. On ne raffine 
  138. *                 pas le maillage mais on allonge le temps
  139. *                 de simulation, ce qui permet une "vrai"
  140. *                 comparaison sur zQ(t) et zR(t) avec les
  141. *                 résultat de notre ami Ramaswamy.
  142. * graph = VRAI  : on trace les graphiques (historiques, champs)
  143. *                 On met graph = faux pour la fiche
  144. *                 de validation.
  145. * inta  = VRAI si on est en interactif ;  0 sinon
  146. *
  147. * N.B. si inta et graph sont VRAIs, on a droit au film des
  148. *      champs de vitesses en fin d'exécution...
  149. *
  150. *
  151. * discconv : option de discrétisation des termes de
  152. *            convection pour les équations de Navier-Stokes
  153. *            info 'NS' ;  pour les choix possibles
  154. *
  155. court = VRAI ;
  156. graph = FAUX ;
  157. inta  = FAUX ;
  158. typelem = 'QUA4' ;
  159. discconv = 'CENTREE' ;
  160. *
  161. * Méthodes de déplacement du maillage (vmail) :
  162. * LAGR - déplacement lagrangien du maillage
  163. * SLAG - surface lagrangienne, vitesse intérieure interpolée 
  164. *        par un laplacien.
  165. *  TOTAL : le contour entier est lagrangien (w=v)
  166. *  INTBOR : la surface libre est lagrangienne (w|surf = v|surf)
  167. *           Sur le fond, w=0
  168. *           Sur les bords : interpolation linéaire de w
  169. *  BORLIB : meme chose que INTBOR sauf :
  170. *           Sur les bords : condition de Neumann pour le
  171. *                           laplacien
  172. * CVS2LF - w // axe z sur la surface => eqn conv resolue
  173. *          par un schéma Leap-Frog (disc centrée dérivée 1°)
  174. *          Options : idem SLAG sauf TOTAL.
  175. *
  176. *vmail = 'LAGR' ;  optlapn = 'TOTAL' ;
  177. *vmail = 'SLAG' ;
  178. *   optlapn = 'TOTAL' ;
  179. *   optlapn = 'INTBOR' ;
  180. *   optlapn = 'BORLIB' ;
  181. vmail = 'CVS2LF' ;
  182.    optlapn = 'INTBOR' ;
  183. *   optlapn = 'BORLIB' ;
  184. 'OPTION' 'DIME' 2 'ELEM' typelem ;
  185. 'OPTION' 'ECHO' 1 ;
  186. 'SI' ('NON' inta) ;
  187.   'OPTION' 'TRAC' 'PS' ;
  188. 'SINON' ;
  189. **od = 'TEXTE' 'OPTION DONN 3' ;
  190. 'FINSI' ;
  191.  
  192. *************************************************************
  193. **** Définition de quelques (!) procédures :
  194. *
  195.  
  196. ** **** D'intéret général :
  197.  
  198. * -> LOG10 :
  199. *     LOG base10
  200. * -> CALCERR :
  201. *     Erreur relative entre deux champoints (norme Linfini)
  202. * -> TRMAIL :
  203. *     Tracé d'un maillage
  204. * -> TRVIT :
  205. *     Tracé d'un chpoint de vitesse VECT SOMMET
  206. * -> TRCH :
  207. *     Tracé d'un chpoint SCAL
  208.  
  209. ** **** Spécifique à l'ALE :
  210.  
  211. * Elles sont en grand nombre car il faut gérer soi-meme les
  212. * dicrétisations temporelle et spatiale (domaine mobile).
  213. * On ne peut pas utiliser (pour l'instant ?) EXEC pour s'en
  214. * occuper.
  215. * En outre, on propose plusieurs méthodes pour déplacer le
  216. * maillage.
  217.  
  218. * -> CREMAIL :
  219. *     Création de la cavité et de tous les sous-maillages
  220. *     nécessaires.
  221. *     Tout ce qui concerne la géométrie doit se trouver ici
  222. * -> CDXHAUT :
  223. *     Création des deltax pour la surface libre (ici $haut)
  224. *     pour utilisation dans la procedure CVS2LF
  225. *     (Différences Finies)
  226. * -> INITVAR :
  227. *     Initialisation de la table des variables
  228. *     Variables aux instants n+1/2 et n-1/2 (si nécessaire)
  229. * -> INITPDT :
  230. *     Initialisation de la table des pas de temps
  231. *     Tout ce qui concerne le calcul du pas de temps
  232. * -> INITHIST :
  233. *     Initialisation de la table des historiques
  234. *     Toutes les listes de scalaires
  235. * -> INITSAUV :
  236. *     Initialisation de la table des sauvegardes
  237. *     Sauvegarde des champoints dans une table indicée 
  238. *     par un entier associé à un de pas de temps
  239. * -> MAJPDT :
  240. *     Mise à jour de la table des pas de temps
  241. *     i.e calcul du nouveau
  242. * -> MAJHIST :
  243. *     Mise à jour des différents historiques
  244. * -> MAJSAUV :
  245. *     On sauvegarde les différents champoints dans une table
  246.  
  247. * -> SUMMARY :
  248. *     Affiche un résumé des paramètres
  249. * -> TRTOUT :
  250. *     Trace tous les champs
  251. ** NB : pour simplifier, on ne transmet pas les paramètres à
  252. **      SUMMARY et TRTOUT. C'est de la mauvaise programmation...
  253. * -> MESI :
  254. *     Affichage du message pour la boucle d'itérations
  255. * -> MEST :
  256. *     Affichage du message pour l'avance d'un pas de temps
  257.  
  258. * -> CALCON :
  259. *     Calcul de la contrainte appliquée sur la surface au premier
  260. *     pas de temps
  261. * -> PREPAS :
  262. *     Faire le premier pas de temps
  263. * -> CVIT :
  264. *     Faire un pas de temps : termes sources      en n-1/2
  265. *                             matrice-masse       en n
  266. *                             divergence-pression en n+1/2
  267. * -> LAGR :
  268. *     vitesse du fluide (v) =>
  269. *     vitesse du maillage (w) en LAGRangien (soit w=v)
  270. * -> SLAG :
  271. *     on choisit w|surf = v|surf (Surface LAGrangienne)
  272. *     on interpole ensuite w sur le reste du maillage par
  273. *     une équation au laplacien
  274. * -> CVS2LF :
  275. *     on calcule w|surf avec v|surf en imposant w|surf // Oy
  276. *     => eqn convection sur la surface.
  277. *     Elle est discrétisée en différences finies centrées
  278. *     d'où le nom : ConVection Surface ordre 2 (schéma de
  279. *     type Leap-Frog)
  280. *     on interpole ensuite w sur le reste du maillage par
  281. *     une équation au laplacien
  282. * ---> IVECT :
  283. *     Interpolation d'un vecteur entre deux points sur
  284. *     un segment. Un chpoint VECT => 2 chpos SCAL
  285. * ---> LAPC :
  286. *     Résolution d'une équation scalaire en implicite ;
  287. *     Laplacien + conditions Dirichlet ou Neumann
  288. *
  289. * -> ETENDRE :
  290. *     Procédure qui étend un champoint VECT SOMMET a tous les
  291. *     points (FACE, CENTRE...) du maillage total
  292. * -> CINVPDTM :
  293. *     Calcul du pas de temps dit "de maillage"
  294.  
  295. ** **** Post-traitement :
  296.  
  297. * -> FILMVIT :
  298. *     Construit une somme de déformées pour faire un film
  299. *     avec les champs de vitesse sauvegardés
  300. * -> TRHIST
  301. *     Tracé des historiques (la non plus, les paramètres ne
  302. *     sont pas transmis).
  303.  
  304. *
  305. *     LOG base10
  306. *
  307. 'DEBPROC' LOG10 machin*'FLOTTANT' ;
  308. 'SI' (machin '&lt;EG' 0.D0) ; a = -365.25 ;
  309.   'MESSAGE' 'Procédure LOG10 : donnée <= 0' ;
  310. 'SINON' ; a = (('LOG' machin) / ('LOG' 10.)) ;
  311. 'FINSI' ;
  312. 'FINPROC' a ;
  313.  
  314. *
  315. *     Erreur relative entre deux champoints (norme Linfini)
  316. *
  317. 'DEBPROC' CALCERR vitp1*'CHPOINT  ' vit*'CHPOINT  ' ;
  318. 'ARGUMENT' cle*'MOT     ' ;
  319. 'SI' ('EGA' cle 'ECHELLE') ;
  320.   'ARGUMENT' vref*'FLOTTANT' ;
  321.   'SI' (vref '&lt;EG' 0.D0) ;
  322.      'ERREUR' 'L échelle doit etre positive' ;
  323.   'FINSI' ;
  324. 'SINON' ;
  325.   mavp1 = 'MAXIMUM' vitp1 'ABS' ;
  326.   mav   = 'MAXIMUM' vit 'ABS' ;
  327.   pmav  = 'PROG' mavp1 mav ;
  328.   mimav = 'MINIMUM' pmav ;
  329.   mamav = 'MAXIMUM' pmav ;
  330.   'SI' (mamav '&lt;EG' 0.D0) ;
  331.      'ERREUR' 'Les deux champs sont nuls' ;
  332.   'FINSI' ;
  333.   'SI' (mimav '&lt;EG' 0.D0) ;
  334.      'MESSAGE' 'Procédure CALCERR : un des champs est nul' ;
  335.      'MESSAGE' '                    => option MOYENNE' ;
  336.      cle = 'MOYENNE' ;
  337.   'FINSI' ;
  338.   'SI' ('EGA' cle 'MINMAX') ;
  339.     vref = mimav ;
  340.   'SINON' ;
  341.     'SI' ('EGA' cle 'MOYENNE') ;
  342.       vref = (mimav '+' mamav ) '/' 2.0D0 ;
  343.     'SINON' ;
  344.       'ERREUR' 'Mauvaise option' ;
  345.     'FINSI' ;
  346.   'FINSI' ;
  347. 'FINSI' ;
  348. 'FINPROC' (('MAXIMUM' (vitp1 '-' vit) 'ABS') '/' vref) ;
  349.  
  350. *
  351. *     Tracé d'un maillage
  352. *
  353. 'DEBPROC' TRMAIL mailp*'MAILLAGE' tphyp/'FLOTTANT' ;
  354. 'SI' ('EXISTE' tphyp) ;
  355.   'TRACER' mailp 'TITR' ('CHAINE' 'Maillage à t=' tphyp) ;
  356. 'SINON' ;
  357.   'TRACER' mailp 'TITR' 'Maillage' ;
  358. 'FINSI' ;
  359. 'FINPROC' ;
  360.  
  361. *
  362. *     Tracé d'un chpoint de vitesse VECT SOMMET
  363. *
  364. 'DEBPROC' TRVIT vitp*'CHPOINT ' vitp2/'CHPOINT ' ;
  365. 'ARGUMENT' mailp*'MAILLAGE' ;
  366. titvp = 'Vitesses' ;
  367. 'REPETER' OPT ;
  368.   'ARGUMENT' cle/'MOT     ' ;
  369.   'SI' ('NON' ('EXISTE' cle)) ;
  370.     'QUITTER' OPT ;
  371.   'FINSI' ;
  372.   'SI' ('EGA' cle 'AMPLIF') ;
  373.     'ARGUMENT' amp*'FLOTTANT' ;
  374.   'FINSI' ;
  375.   'SI' ('EGA' cle 'TEMPS') ;
  376.     'ARGUMENT' tphp*'FLOTTANT' ;
  377.     titvp = ('CHAINE' titvp ' à t=' tphp) ;
  378.   'FINSI' ;
  379. 'FIN' OPT ;
  380. 'SI' ('NEG' ('TYPE' amp) 'FLOTTANT') ;
  381.   mvp = ('MAXIMUM'  vitp 'ABS') ;
  382.   'SI' (mvp '&lt;EG' 0.D0) ;
  383. *    'ERREUR' 'Champ de vitesses nul' ;
  384.     'MESSAGE' 'TRVIT : Champ de vitesses nul' ;
  385.     'QUITTER' TRVIT ;
  386.   'SINON' ;
  387.     xp yp = 'COORDONNEE' mailp ;
  388.     rp = 'PROG' (('MAXIMUM' xp) '-' ('MINIMUM' xp))  (('MAXIMUM' yp) '-' ('MINIMUM' yp)) ;
  389.     amp = ('MAXIMUM' rp) '/' (15.D0 '*' mvp) ;
  390.   'FINSI' ;
  391. 'FINSI' ;
  392. ungp = 'VECTEUR' vitp amp 'UX' 'UY' 'JAUNE' ;
  393. 'SI' ('EXISTE' vitp2) ;
  394.   ungp2 = 'VECTEUR' vitp2 amp 'UX' 'UY' 'ROUGE' ;
  395.   ungp = (ungp2 'ET' ungp) ;
  396. 'FINSI' ;
  397. 'TRACER'  ungp mailp 'TITR' titvp ;
  398. 'FINPROC';
  399.  
  400. *
  401. *     Tracé d'un chpoint SCAL
  402. *
  403. 'DEBPROC' TRCH vvp*'CHPOINT ' ;
  404. 'SI' ('NEG' ('DIME' ('EXTRAIRE' vvp 'COMP')) 1) ;
  405.   'ERREUR' 'Il faut un champ à UNE composante' ;
  406. 'FINSI' ;
  407. tvvp = 'EXTRAIRE' vvp 'TYPE' ;
  408. 'SI' ('EGA' tvvp 'CENTRE  ') ;
  409.   'ARGUMENT' $tdefp*'TABLE   ' cle*'MOT     ' ;
  410.   'SI' ('EGA' cle 'TRICHE') ;
  411.     vvp = ('ELNO' $tdefp ('KCHT' $tdefp 'SCAL' tvvp vvp)) ;
  412.     mvvp = ($tdefp . 'MAILLAGE') ;
  413.     tm = mvvp ;
  414.   'SINON' ;
  415.     'SI' ('EGA' cle 'NORMAL') ;
  416.       tm = ($tdefp . 'MAILLAGE') ;
  417.       mvvp = 'MODELISER' tm 'THERMIQUE' ;
  418.       vvp  = ('KCHA' $tdefp 'CHAM' vvp) ;
  419.     'SINON' ;
  420.         'ERREUR' 'Mauvaise option' ;
  421.     'FINSI' ;
  422.   'FINSI' ;
  423. 'SINON' ;
  424.   'SI' ('EGA' tvvp 'FACE    ') ;
  425.     'ERREUR' 'Les CHPOINTs FACE ne sont pas traités...' ;
  426.   'SINON' ;
  427.     'SI' ('NON' ('EGA' tvvp 'SOMMET  ')) ;
  428.       'MESSAGE' 'CHPOINT non créé par KCHT...' ;
  429.     'FINSI' ;
  430.     'ARGUMENT' mvvp*'MAILLAGE' ;
  431.     tm = mvvp ;
  432.   'FINSI' ;
  433. 'FINSI' ;
  434. 'ARGUMENT' nomvp/'MOT     ' tphp/'FLOTTANT' cmvvp/'MAILLAGE' ;
  435. 'ARGUMENT' nbisov/'ENTIER  ' ;
  436. 'SI' ('EXISTE' nomvp) ;
  437.   titvp = ('CHAINE' nomvp ' ; chpoint SCAL ' tvvp ) ;
  438. 'SINON' ;
  439.   titvp = ('CHAINE' 'Chpoint SCAL ' tvvp) ;
  440. 'FINSI' ;
  441. 'SI' ('EXISTE' tphp) ;
  442.   titvp = ('CHAINE' titvp ' à t=' tphp) ;
  443. 'FINSI' ;
  444. 'SI' ('NON' ('EXISTE' cmvvp)) ;
  445.   cmvvp = 'CONTOUR' tm ;
  446. 'FINSI' ;
  447. maxvvp = ('MAXIMUM' vvp 'ABS') ;
  448. minvvp = ('MINIMUM' vvp 'ABS') ;
  449. 'SI' (maxvvp '&lt;EG' (minvvp 'ABS')) ;
  450. *    'ERREUR' 'Champ constant' ;
  451.     'MESSAGE' ('CHAINE' 'TRCH : Champ constant = ' maxvvp ) ;
  452.     'QUITTER' TRCH ;
  453. 'FINSI' ;
  454. 'SI' ('EXISTE' nbisov) ;
  455.   'TRACER' vvp mvvp cmvvp nbisov 'TITR' titvp ;
  456. 'SINON' ;
  457.   'TRACER' vvp mvvp cmvvp 'TITR' titvp ;
  458. 'FINSI' ;
  459. 'FINPROC' ;
  460.  
  461. * 
  462. *     Création de la cavité et de tous les sous-maillages
  463. *     nécessaires.
  464. *     Tout ce qui concerne la géométrie doit se trouver ici
  465. * 
  466. 'DEBPROC' CREMAIL $cav*'TABLE   ' ;
  467. 'ARGUMENT'  l*'FLOTTANT'  h*'FLOTTANT' ;
  468. 'ARGUMENT' nl*'ENTIER  ' nh*'ENTIER  ' ;
  469. 'ARGUMENT' egeom*'FLOTTANT' ;
  470. p00 = 0.D0 0.D0 ;  pl0 = l 0.D0 ;
  471. p0h = 0.D0 h  ;    plh = l h  ;
  472. dxg = 1.0D0 ;  dxd = 1.0D0 ;  dyb = 1.0D0 ;  dyh = 1.0D0 ;
  473. 'REPETER' OPT ;
  474.   'ARGUMENT' cle/'MOT     ' ;
  475.   'SI' ('NON' ('EXISTE' cle)) ;
  476.     'QUITTER' OPT ;
  477.   'FINSI' ;
  478.   'SI' ('EGA' cle 'DX') ;
  479.     'ARGUMENT' dgp*'FLOTTANT' ddp*'FLOTTANT' ;
  480.     dxmoy = l '/' nl ;
  481.     dxg = dgp '*' dxmoy ;  dxd = ddp '*' dxmoy ;
  482.   'FINSI' ;
  483.   'SI' ('EGA' cle 'DY') ;
  484.     'ARGUMENT' dbp*'FLOTTANT' dhp*'FLOTTANT' ;
  485.     dymoy = h '/' nh ;
  486.     dyb = dbp '*' dymoy ;  dyh = dhp '*' dymoy ;
  487.   'FINSI' ;
  488. 'FIN' OPT ;
  489. bas   = pl0 'DROIT' (-1 '*' nl) p00 'DINI' dxd 'DFIN' dxg ;
  490. haut  = p0h 'DROIT' (-1 '*' nl) plh 'DINI' dxg 'DFIN' dxd ;
  491. rwall = plh 'DROIT' (-1 '*' nh) pl0 'DINI' dyh 'DFIN' dyb ;
  492. lwall = p00 'DROIT' (-1 '*' nh) p0h 'DINI' dyb 'DFIN' dyh ;
  493. $cav . 'dxg' = dxg ;  $cav . 'dxd' = dxd ;
  494. $cav . 'dyb' = dyb ;  $cav . 'dyh' = dyh ;
  495. mt = 'DALLER' haut rwall bas lwall ;
  496. $mt = 'DOMA' mt 'IMPR' ;
  497. mt = ($mt . 'MAILLAGE') ;  'TASSER' mt ;
  498. cmt = ('CONTOUR' mt)  ;  lmt = chan ligne $mt.maillage ;
  499. *haut  = 'CHANGER' haut SEG2   ;  bas   = 'CHANGER' bas SEG2 ;
  500. *lwall = 'CHANGER' lwall SEG2 ;   rwall = 'CHANGER' rwall SEG2 ;
  501. $cmt   = 'DOMA' cmt   'INCL' $mt egeom ;
  502. $lmt   = 'DOMA' lmt   'INCL' $mt egeom ;
  503. $haut  = 'DOMA' haut  'INCL' $mt egeom ;
  504. $bas   = 'DOMA' bas   'INCL' $mt egeom ;
  505. $lwall = 'DOMA' lwall 'INCL' $mt egeom ;
  506. $rwall = 'DOMA' rwall 'INCL' $mt egeom ;
  507. 'ELIMINATION'  (mt 'ET' lmt 'ET' bas 'ET' rwall 'ET' haut  'ET' lwall) egeom ;
  508. contini = 'COULEUR' cmt 'ROUG' ;
  509.  
  510. $cav . 'l' = l ;  $cav . 'nl' = nl ;
  511. $cav . 'h' = h ;  $cav . 'nh' = nh ;
  512. $cav . 'volini' = 'MESURE' mt ;
  513. $cav . '$mt' = $mt ;
  514. $cav . '$cmt'    = $cmt ;    $cav . '$lmt'  = $lmt ;
  515. $cav . '$haut'   = $haut ;   $cav . '$bas'  = $bas ;
  516. $cav . '$lwall'  = $lwall ;  $cav .'$rwall' = $rwall ;
  517. $cav . 'contini' =  (contini 'PLUS' (egeom egeom)) ;
  518. $cav . 'cnfini'  = FORM ;
  519. *   Définition des conditions aux limites pour Navier-Stokes
  520. *   VN nulle sur lwall, rwall et bas
  521. clr = 'MANUEL' 'CHPO'  (lwall 'ET' rwall) 1 '1UN' 0.D0 'NATU' 'DISCRET';
  522. cb  = 'MANUEL' 'CHPO'  (bas) 1 '2UN' 0.D0 'NATU' 'DISCRET' ;
  523. $cav . 'CLIM' = (clr 'ET' cb) ;
  524. 'MESSAGE' 'Procédure CREMAIL : table MAILLAGE créée' ;
  525. 'FINPROC' $cav ;
  526.  
  527. *
  528. *     Création des deltax pour la surface libre (ici $haut)
  529. *     pour utilisation dans la procedure CVS2LF
  530. *     (Différences Finies)
  531. *
  532. 'DEBPROC' CDXHAUT $h*'TABLE   ' ;
  533. h = ($h . 'MAILLAGE') ;
  534. xc = 'KCHT' $h 'SCAL' 'SOMMET' ('COORDONNEES' 1 h) ;
  535. xh = 'PROG' ;
  536. 'REPETER' BBB ($h . 'NPTD') ;
  537.   ppp = h 'POIN' (&BBB) ;
  538.   xh = (xh 'ET' ('PROG' ('EXTRAIRE' xc 'SCAL' ppp))) ;
  539. 'FIN' BBB ;
  540. psg  = 'LECT' 1 'PAS' 1 (($h . 'NPTD') '-' 1) ;
  541. psd  = 'LECT' 2 'PAS' 1  ($h . 'NPTD') ;
  542. $h . 'psgg' = 'LECT' 1 'PAS' 1 (($h . 'NPTD') '-' 2) ;
  543. $h . 'psm'  = 'LECT' 2 'PAS' 1 (($h . 'NPTD') '-' 1) ;
  544. $h . 'psdd' = 'LECT' 3 'PAS' 1 ($h . 'NPTD') ;
  545. $h . 'dxh' = ('EXTRAIRE' psd xh) '-' ('EXTRAIRE' psg xh) ;
  546. 'FINPROC' $h ;
  547.  
  548. *
  549. *     Initialisation de la table des variables
  550. *     Variables aux instants n+1/2 et n-1/2 (si nécessaire)
  551. *
  552. 'DEBPROC' INITVAR $v*'TABLE   ' ;
  553. 'ARGUMENT' $c*'TABLE   ' ;
  554. $mt = $c . '$mt' ;
  555. * Champoints
  556. $v . 'unm05'    = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
  557. $v . 'wnm05'    = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
  558. $v . 'wnp05i'   = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
  559. $v . 'unp05i'   = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
  560. $v . 'unp05im1' = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
  561. $v . 'anp05i'   = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
  562. $v . 'pnp05i'   = 'KCHT' $mt 'SCAL' 'CENTRE' 0.D0 ;
  563. $v . 'dnp05i'   = 'KCHT' $mt 'SCAL' 'CENTRE' 0.D0 ;
  564. * Configuration
  565. $v . 'deltaxn' = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
  566. $v . 'cnfn'    = 'FORME' ;
  567. $v . 'cnfnm05' = 'FORME' ;
  568. $v . 'cnfnp05i' = 'FORME' ;
  569. 'MESSAGE' 'Procédure INITVAR : table VARIABLE créée' ;
  570. 'FINPROC' $v ;
  571.  
  572. *
  573. *     Initialisation de la table des pas de temps
  574. *     Tout ce qui concerne le calcul du pas de temps
  575. *
  576. 'DEBPROC' INITPDT $t*'TABLE   ' ;
  577. 'ARGUMENT' alfa*'FLOTTANT' gamma*'FLOTTANT' ;
  578. 'ARGUMENT' dtini*'FLOTTANT' poini*'ENTIER  ' ;
  579. * Entiers
  580. $t . 'nupdt' = 0 ;
  581. $t . 'iter' = 0 ;
  582. * Réels
  583. $t . 'tphyn' = 0.D0 ;   $t . 'tpscpu' = 0.D0 ;
  584. $t . 'alfa'  = alfa ;
  585. $t . 'gamma' = gamma ;  $t . 'poini'  = poini ;
  586. $t . 'conv'  = dtini ;  $t . 'diff'   = dtini ;
  587. $t . 'mail'  = dtini ;  $t . 'moy'    = 0.D0 ;
  588. $t . 'min'   = dtini ;
  589. 'MESSAGE' 'Procédure INITPDT : table PASDETPS créée' ;
  590. 'FINPROC' $t ;
  591.  
  592. *
  593. *     Initialisation de la table des historiques
  594. *     Toutes les listes de scalaires
  595. *
  596. 'DEBPROC' INITHIST $h*'TABLE   ' ;
  597. * Listes
  598. $h . 'ltps'  = 'PROG' ;
  599. $h . 'liter' = 'PROG' ;
  600. $h . 'lzq'   = 'PROG' ;  $h . 'lzr'  = 'PROG' ;
  601. $h . 'ldv'   = 'PROG' ;  $h . 'lvol' = 'PROG' ;
  602. $h . 'lpdt'  = 'TABLE' 'LPDT' ; hl = $h . 'lpdt' ;
  603. hl . 'conv'  = 'PROG' ;  hl . 'diff' = 'PROG' ;
  604. hl . 'mail'  = 'PROG' ;  hl . 'moy'  = 'PROG' ;
  605. hl . 'min'   = 'PROG' ;
  606. 'MESSAGE' 'Procédure INITHIST : table HIST créée' ;
  607. 'FINPROC' $h ;
  608.  
  609. *
  610. *     Initialisation de la table des sauvegardes
  611. *     Sauvegarde des champoints dans une table indicée 
  612. *     par un entier associé à un pas de temps
  613. *
  614. 'DEBPROC' INITSAUV $s*'TABLE   ' ;
  615. * Entiers
  616. $s . 'ntrc' = 0 ;
  617. * Listes
  618. $s . 'lttrc' = 'PROG' ;
  619. * Tables
  620. $s . 'DELTAX'   = 'TABLE' 'DELTAX' ;
  621. $s . 'VITESSE'  = 'TABLE' 'VITESSE' ;
  622. *$s . 'PRESSION' = 'TABLE' 'PRESSION' ;
  623. *$s . 'DIVV'    = 'TABLE' 'DIVV' ;
  624. 'MESSAGE' 'Procédure INITSAUV : table SAUV créée' ;
  625. 'FINPROC' $s ;
  626.  
  627. *
  628. *     Mise à jour de la table des pas de temps
  629. *     i.e calcul du nouveau
  630. *
  631. 'DEBPROC' MAJPDT $p*'TABLE   ' ;
  632. dtcr = ($p . 'alfa')  '*' ($p . 'conv') ;
  633. dtdr = ($p . 'alfa')  '*' ($p . 'diff') ;
  634. dtmr = ($p . 'gamma') '*' ($p . 'mail') ;
  635. dti = ('MINIMUM' ('PROG' dtcr dtdr dtmr)) ;
  636. dtmoy = ($p . 'moy') ;
  637. npdt = ($p . 'nupdt') ;  poi = ($p . 'poini') ;
  638. $p . 'moy' = (((npdt '+' poi) '*' dtmoy) '+' dti) / (npdt '+' poi '+' 1) ;
  639. $p.'min' = ('MINIMUM' ('PROG' dti ($p . 'moy'))) ;
  640. 'FINPROC' $p ;
  641.  
  642. *
  643. *     Mise à jour des différents historiques
  644. *
  645. 'DEBPROC' MAJHIST $h*'TABLE   ' $c*'TABLE   ' ;
  646. 'ARGUMENT' $p*'TABLE   ' $v*'TABLE   ' ;
  647. $h . 'ltps'  = (($h . 'ltps')  'ET' ('PROG' ($p . 'tphyn'))) ;
  648. $h . 'liter' = (($h . 'liter') 'ET' ('PROG' ($p . 'iter'))) ;
  649. r = (($c . '$haut' . 'MAILLAGE') 'POINT' 'INITIAL') ;
  650. q = (($c . '$haut' . 'MAILLAGE') 'POINT' 'FINAL') ;
  651. hini = ($c . 'h') ;  volini = ($c . 'volini') ;
  652. $h . 'lzr' = (($h . 'lzr') 'ET'  ('PROG' (('COORDONNEE' 2 r) '-' hini))) ;
  653. $h . 'lzq' = (($h . 'lzq') 'ET'  ('PROG' (('COORDONNEE' 2 q) '-' hini))) ;
  654. $h . 'ldv' = (($h . 'ldv') 'ET'  ('PROG' ('MAXIMUM' ($v . 'dnp05i') 'ABS'))) ;
  655. $h . 'lvol' = (($h . 'lvol') 'ET'  ('PROG' ((('MESURE' ($c . '$mt' . 'MAILLAGE')) '-' volini)  '/' volini))) ;
  656. $hl = ($h . 'lpdt') ;
  657. $hl . 'conv' = (($hl . 'conv') 'ET' ('PROG' ($p . 'conv'))) ;
  658. $hl . 'diff' = (($hl . 'diff') 'ET' ('PROG' ($p . 'diff'))) ;
  659. $hl . 'moy'  = (($hl . 'moy')  'ET' ('PROG' ($p . 'moy'))) ;
  660. $hl . 'mail' = (($hl . 'mail') 'ET' ('PROG' ($p . 'mail'))) ;
  661. $hl . 'min'  = (($hl . 'min')  'ET' ('PROG' ($p . 'min'))) ;
  662. 'FINPROC' $h ;
  663.  
  664. *
  665. *     On sauvegarde les différents champoints dans une table
  666. *
  667. 'DEBPROC' MAJSAUV $s*'TABLE   ' ;
  668. 'ARGUMENT' $v*'TABLE   ' $p*'TABLE   ' ;
  669. $s . 'ntrc' = (($s . 'ntrc') '+' 1) ;
  670. $s . 'lttrc' = (($s . 'lttrc') 'ET' ('PROG' ($p . 'tphyn'))) ;
  671. $s . 'DELTAX' . ($s . 'ntrc')   = 'COPIE' ($v . 'deltaxn') ;
  672. $s . 'VITESSE' . ($s . 'ntrc')  = 'COPIE' ($v . 'unp05i') ;
  673. *$s . 'PRESSION' . ($s . 'ntrc') = 'COPIE' ($v . 'pnp05i') ;
  674. *$s . 'DIVV' . ($s . 'ntrc') = 'COPIE' ($v . 'dnp05i') ;
  675. 'FINPROC' $s ;
  676.  
  677. *
  678. *     Affiche un résumé des paramètres
  679. *
  680. 'DEBPROC' SUMMARY ;
  681. 'SAUTER' 'LIGNE' ;
  682. 'MESSAGE' ('CHAINE' 'Programme : ' prg) ;
  683. 'SAUTER' 'LIGNE' ;
  684. 'MESSAGE' ('CHAINE' '***  GEOMETRIE :') ;
  685. 'MESSAGE' ('CHAINE' 'Boite carrée : ' l ' x' h) ;
  686. 'MESSAGE' ('CHAINE' '      soient :  ' nl 'x' nh '='  (nl '*' nh) ' éléments ' typelem) ;
  687. *'MESSAGE' ('CHAINE' '      soient :  ' (nl '*' nh '*' 4)
  688. * ' éléments linéaires') ;
  689. 'SAUTER' 'LIGNE' ;
  690. 'MESSAGE' ('CHAINE' '***  PHYSIQUE :') ;
  691. 'MESSAGE' ('CHAINE' 'Gravité : ' ('COORDONNEES' 2 g) ' uz ; viscosité :' nu) ;
  692. 'MESSAGE' ('CHAINE' 'Intensité Dirac Pression : ' pt0) ;
  693. 'SAUTER' 'LIGNE' ;
  694. 'MESSAGE' ('CHAINE' '***  NUMERIQUE :') ;
  695. 'MESSAGE' ('CHAINE' 'Mode de déplacement du maillage : ' vmail ' ' optlapn) ;
  696. 'MESSAGE' ('CHAINE' 'Discrétisation ' discconv ' des termes de convection.');
  697. 'MESSAGE' ('CHAINE' ' Alfa  = ' alfa) ;
  698. 'MESSAGE' ('CHAINE' ' Gamma = ' gamma) ;
  699. 'MESSAGE' ('CHAINE' ' Cvg (bcl iter) = ' cvg) ;
  700. 'MESSAGE' ('CHAINE' ' Poids (pdt1) = ' poini) ;
  701. 'SAUTER' 'LIGNE' ;
  702. 'MESSAGE' ('CHAINE' '***  TEMPS :') ;
  703. 'MESSAGE' ('CHAINE' 'Temps physique : ' ($pdt . 'tphyn')) ;
  704. 'MESSAGE' ('CHAINE' 'Temps CPU      : ' ($pdt . 'tpscpu')) ;
  705. 'FINPROC' ;
  706.  
  707. *
  708. *     Trace tous les champs
  709. *
  710. 'DEBPROC' TRTOUT ;
  711. $mt = ($cavite . '$mt') ; mtm = ($mt . 'MAILLAGE') ;
  712. cini = ($cavite . 'contini') ;
  713. TRMAIL (mtm 'ET' cini) ($pdt . 'tphyn') ;
  714. 'SI' ('NON' inta) ;
  715.   TRVIT  ($var . 'unm05') (('CONTOUR' mtm) 'ET' cini) 'AMPLIF' ampvit 'TEMPS' ($pdt . 'tphyn') ;
  716. 'SINON' ;
  717.   TRVIT  ($var . 'unm05') ($var . 'wnm05') (('CONTOUR' mtm) 'ET' cini) 'AMPLIF' ampvit 'TEMPS' ($pdt . 'tphyn') ;
  718. 'FINSI' ;
  719. TRCH   ($var . 'pnp05i') $mt 'NORMAL' 'Pression' ($pdt . 'tphyn') ;
  720. 'FINPROC' ;
  721.  
  722. *
  723. *     Affichage du message pour la boucle d'itérations
  724. *
  725. 'DEBPROC' MESI $v*'TABLE   ' $p*'TABLE   ' ;
  726. 'SI' ('EGA' ($p . 'iter') 1) ;
  727.   'MESSAGE' ('CHAINE' 'Iter'/1 'Dtconv'/7 'Dtdiff'/21  'Dtmaill'/35 'Dt'/49 'Erreur (log)'/61) ;
  728. 'FINSI' ;
  729. 'MESSAGE' ('CHAINE' ($p . 'iter')*3 ($p . 'conv')/5  ($p . 'diff')/19 ($p . 'mail')/33 ($p . 'min')/47  (LOG10 ($v . 'erri'))/61 ) ;
  730. 'FINPROC' ;
  731.  
  732. *
  733. *     Affichage du message pour l'avance d'un pas de temps
  734. *
  735. 'DEBPROC' MEST $v*'TABLE   ' $p*'TABLE   ' ;
  736. 'MESSAGE' ('CHAINE' '----- Nupdt'/1 'Tphy'/15 'Dt'/29  'Max(div)'/43 'Min(pres)'/57) ;
  737. 'MESSAGE' ('CHAINE' ($p . 'nupdt')*9 ($p . 'tphyn')/13 ($p . 'min')/27 ('MAXIMUM' ($v . 'dnp05i') 'ABS')/41 ('MINIMUM' ($v . 'pnp05i') 'ABS')/55) ;
  738. 'FINPROC' ;
  739.  
  740. *
  741. *     Calcul de la contrainte appliquée sur la surface au premier
  742. *     pas de temps
  743. *
  744. 'DEBPROC' CALCON $h*'TABLE   ' ac*'FLOTTANT' ;
  745. lini = 'ABS' ( ('COORDONNEE' 1  (($h . 'MAILLAGE') 'POINT' 'INITIAL')) '-' ('COORDONNEE' 1  (($h . 'MAILLAGE') 'POINT' 'FINAL'))) ;
  746. xhaut = 'COORDONNEE' 1 ($h . 'CENTRE') ;
  747. xts = 'NOMC' 'UX' ('KCHT' $h 'SCAL' 'CENTRE' 0.D0)  'NATU' 'DISCRET' ;
  748. yts = 'NOMC' 'UY'  (ac '*' ('COS' ((180.D0 '/' lini) '*' xhaut)))  'NATU' 'DISCRET' ;
  749. 'FINPROC' ('KCHT' $h 'VECT' 'CENTRE' (xts 'ET' yts)) ;
  750.  
  751. *
  752. *     Faire le premier pas de temps
  753. *
  754. 'DEBPROC' PREPAS $v*'TABLE   ' $p*'TABLE   ' ;
  755. 'ARGUMENT' $c*'TABLE   ' ;
  756. 'ARGUMENT' tsu*'CHPOINT ' nu*'FLOTTANT' g*'POINT   ' ;
  757. $mt = $c . '$mt' ;  $haut = $c . '$haut' ;
  758. rv = 'EQEX' $mt 'OPTI' 'ALE' 'CENTREE' 'ZONE' $mt   'OPER' 'NS' nu g 'CN' 'INCO' 'UN' 'ZONE' $haut 'OPER' 'TOIM' 'PSURF' 'INCO' 'UN' ;
  759. rv . 'PRESSION' = 'EQPR' $mt  'ZONE' $mt   'OPER' 'PRESSION' 0.D0 ;
  760. rvp = rv . 'PRESSION' ;
  761. rv . 'INCO' = 'TABLE' 'INCO' ;
  762. rv . 'INCO' . 'UN' = 'COPIE' ($v . 'unm05') ;
  763. rv . 'INCO' . 'CN' = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
  764. rv . 'INCO' . 'PSURF' = 'KOPS' tsu '/' ($p . 'min') ;
  765. rv . 'CLIM' = 'COPIE' ($c . 'CLIM') ;
  766. ******
  767. * Calcul des matrices masse diagonales
  768. 'KDIA' rv ;
  769. rvp . 'CLIM' = rv . 'CLIM' ;
  770. rv . 'KIZD' . 'UN' = ('KOPS' (rv . 'KIZD' . 'UN') 'CLIM' (rv . 'CLIM') 1) ;
  771. rvp . 'DIAGV'    = (rv . 'KIZD' . 'UN') ;
  772.  
  773. * Calcul des matrices de la divergence 
  774. rvp . 'MATC'     = 'KMAB' rvp ;
  775. rvp . 'PRESSION' = 'KCHT' (rvp . 'DOMAINE')  'SCAL' 'CENTRE' 0.0D0 ;
  776. * Calcul des termes source 
  777. nbop = 'DIMENSION' (rv . 'LISTOPER') ;
  778. 'REPETER' BLOP nbop ;
  779.   nomper  = 'EXTRAIRE' &BLOP (rv . 'LISTOPER') ;
  780.   notable = 'MOT' ('TEXTE' ('CHAINE' &BLOP nomper)) ;
  781.   ('TEXTE' nomper) (rv . notable) ;
  782. 'FIN' BLOP ;
  783. 'OUBLIER' $mt 'MATESI' ;  'OUBLIER' $mt 'XXPSOML' ;
  784. 'OUBLIER' $mt 'XXDXDY' ;  'OUBLIER' $mt 'XXVOLUM' ;
  785. dt = ($p . 'min') ;
  786. rv . 'PASDETPS' . 'DELTAT' = dt ;
  787.  
  788. * Calcul de la pression
  789. rvp . 'DELTAT' = dt ;
  790. dm1f = 'KOPS' ('KOPS' (rv . 'KIZG' . 'UN') '/'  (rv . 'KIZD' . 'UN'))  '-'  ('KOPS' (rv . 'INCO' . 'UN')  '/' dt) ;
  791. P = 'KMF' (rvp . 'MATC') dm1f ;
  792. rvp . 'FORCE' = 'KCHT' (rvp . 'DOMAINE') 'SCAL' 'CENTRE' 0.D0 ;
  793. 'KRES' rvp P 'BETA' (rvp . 'KBETA') (rvp . 'BETA') 'PIMP' (rvp . 'KPIMP') (rvp . 'PIMP') ;
  794. rvp . 'PRESSION' = P ;
  795. rvp . 'GRADP'    = 'KMTP' (rvp . 'MATC')   (rvp . 'PRESSION') ;
  796. rv . 'KIZG' . 'UN' = ('KOPS' (rv . 'KIZG' . 'UN')  '-'  (rvp . 'GRADP')) ;
  797. * Avance en temps 
  798. $v . 'unp05im1' = ('COPIE' ($v . 'unp05i')) ;
  799. $v . 'unp05i'   = ('KOPS' (rv . 'INCO' . 'UN')  '-'  ('KOPS' ('KOPS' (rv . 'KIZG' . 'UN') '/' (rv . 'KIZD' . 'UN')) '*' dt)) ;
  800. $v . 'pnp05i'  = rvp . 'PRESSION' ;
  801. $v . 'dnp05i'  = ('KMF' (rvp . 'MATC') ($v . 'unp05i')) ;
  802. $p . 'conv' = rv . 'PASDETPS' . 'DTCONV' ;
  803. $p . 'diff' = rv . 'PASDETPS' . 'DTDIFU' ;
  804. $p . 'iter' = ($p . 'iter' '+' 1) ;
  805. FINPROC $v $p ;
  806.  
  807. *
  808. *     Faire un pas de temps : termes sources      en n-1/2
  809. *                             matrice-masse       en n
  810. *                             divergence-pression en n+1/2
  811. *
  812. 'DEBPROC' CVIT $c*'TABLE   ' $v*'TABLE   ' $p*'TABLE   ' ;
  813. 'ARGUMENT' nu*'FLOTTANT' g*'POINT   ' dc*'MOT     ' ;
  814. $mt = $c . '$mt' ;
  815. rv = 'EQEX' $mt 'OPTI' 'ALE' dc 'ZONE' $mt   'OPER' 'NS' nu g 'CN' 'INCO' 'UN' ;
  816. rv . 'PRESSION' = 'EQPR' $mt  'ZONE' $mt 'OPER' 'PRESSION' 0.D0 ;
  817. rvp = rv . 'PRESSION' ;
  818. rv . 'INCO' = 'TABLE' 'INCO' ;
  819. rv . 'INCO' . 'UN' = ('COPIE' ($v . 'unm05')) ;
  820. rv . 'INCO' . 'CN' = ('KOPS' (rv . 'INCO' . 'UN')  '-'  ($v . 'wnm05')) ;
  821. rv . 'PASDETPS' . 'PDTIMP' = ($p . 'min') ;
  822. rv . 'CLIM' = ($c . 'CLIM') ;
  823. *****
  824. * Calcul des matrices masse diagonales
  825. 'FORME' ($v . 'cnfn') ;
  826. 'KDIA' rv ;
  827. rvp . 'CLIM' = rv . 'CLIM' ;
  828. rv . 'KIZD' . 'UN' = ('KOPS' (rv . 'KIZD' . 'UN') 'CLIM' (rv . 'CLIM') 1) ;
  829. rvp . 'DIAGV'    = (rv . 'KIZD' . 'UN') ;
  830.  
  831. * Calcul des matrices de la divergence 
  832. 'FORME' ($v . 'cnfnp05i') ;
  833. rvp . 'MATC'     = 'KMAB' rvp ;
  834. rvp . 'PRESSION' = 'KCHT' (rvp . 'DOMAINE')  'SCAL' 'CENTRE' 0.0D0 ;
  835. * Calcul des termes source 
  836. 'FORME' ($v . 'cnfnm05') ;
  837. 'NS' (rv . '1NS') ;
  838. 'OUBLIER' $mt 'MATESI' ;  'OUBLIER' $mt 'XXPSOML' ;
  839. 'OUBLIER' $mt 'XXDXDY' ;  'OUBLIER' $mt 'XXVOLUM' ;
  840. dt = ($p . 'min') ;
  841. rv . 'PASDETPS' . 'DELTAT' = dt ;
  842.  
  843. * Calcul de la pression
  844. 'FORME' ($v . 'cnfnp05i') ;
  845. rvp . 'DELTAT' = dt ;
  846. dm1f = 'KOPS' ('KOPS' (rv . 'KIZG' . 'UN') '/'  (rv . 'KIZD' . 'UN'))  '-'  ('KOPS' (rv . 'INCO' . 'UN')  '/' dt) ;
  847. P = 'KMF' (rvp . 'MATC') dm1f ;
  848. rvp . 'FORCE' = 'KCHT' (rvp . 'DOMAINE') 'SCAL' 'CENTRE' 0.D0 ;
  849. 'KRES' rvp P 'BETA' (rvp . 'KBETA') (rvp . 'BETA') 'PIMP' (rvp . 'KPIMP') (rvp . 'PIMP') ;
  850. rvp . 'PRESSION' = P ;
  851. rvp . 'GRADP'    = 'KMTP' (rvp . 'MATC')   (rvp . 'PRESSION') ;
  852. rv . 'KIZG' . 'UN' = ('KOPS' (rv . 'KIZG' . 'UN')  '-'  (rvp . 'GRADP')) ;
  853. * Avance en temps 
  854. $v . 'unp05im1' = ('COPIE' ($v . 'unp05i')) ;
  855. $v . 'unp05i'   = ('KOPS' (rv . 'INCO' . 'UN')  '-'  ('KOPS' ('KOPS' (rv . 'KIZG' . 'UN') '/'  (rv . 'KIZD' . 'UN')) '*' dt)) ;
  856. $v . 'pnp05i'  = rvp . 'PRESSION' ;
  857. $v . 'dnp05i'  = ('KMF' (rvp . 'MATC') ($v . 'unp05i')) ;
  858. $p . 'conv' = rv . 'PASDETPS' . 'DTCONV' ;
  859. $p . 'diff' = rv . 'PASDETPS' . 'DTDIFU' ;
  860. 'FINPROC' $v $p ;
  861.  
  862. *
  863. *     vitesse du fluide (v) =>
  864. *     vitesse du maillage (w) en Lagrangien
  865. *
  866. 'DEBPROC' LAGR veuler*'CHPOINT '  'ARGUMENT' $c/'TABLE   ' dumby1/'MOT     ' dumby2/'CONFIGUR' ;
  867. 'FINPROC' ('COPIE' veuler) ;
  868.  
  869. *
  870. *     on choisit w|surf = v|surf (Surface LAGrangienne)
  871. *     on interpole ensuite w sur le reste du maillage par
  872. *     une équation au laplacien
  873. *
  874. 'DEBPROC' SLAG veuler*'CHPOINT ' $c*'TABLE   ' ;
  875. 'ARGUMENT' optlap*'MOT     ' conflap*'CONFIGUR' ;
  876. $haut  = ($c . '$haut')  ;  $bas   = ($c . '$bas') ;
  877. $lwall = ($c . '$lwall') ;  $rwall = ($c . '$rwall') ;
  878. $cmt = ($c . '$cmt') ;
  879. $mt  = ($c . '$mt') ;
  880. haut  = ($haut . 'MAILLAGE') ;  bas   = ($bas . 'MAILLAGE') ;
  881. lwall = ($lwall . 'MAILLAGE') ; rwall = ($rwall . 'MAILLAGE') ;
  882. sauvform = 'FORME' ;
  883. 'FORME' conflap ;
  884. 'SI' ('EGA' optlap 'INTBOR') ;
  885.   vithx = 'KCHT' $haut 'SCAL' 'SOMMET' ('EXCO' 'UX' veuler) ;
  886.   vithy = 'KCHT' $haut 'SCAL' 'SOMMET' ('EXCO' 'UY' veuler) ;
  887.   vlwx vlwy = IVECT $lwall veuler ;
  888.   vrwx vrwy = IVECT $rwall veuler ;
  889.   vcx = 'KCHT' $cmt  'SCAL' 'SOMMET' 0.D0 vithx vlwx vrwx ;
  890.   vmx = LAPC $mt vcx ;
  891.   vcy = 'KCHT' $cmt  'SCAL' 'SOMMET' 0.D0 vithy vlwy vrwy ;
  892.   vmy = LAPC $mt vcy ;
  893. 'SINON' ;
  894.   'SI' ('EGA' optlap 'BORLIB') ;
  895.     vithx = 'KCHT' $haut 'SCAL' 'SOMMET' ('EXCO' 'UX' veuler) ;
  896.     vithy = 'KCHT' $haut 'SCAL' 'SOMMET' ('EXCO' 'UY' veuler) ;
  897.     vcx = 'KCHT' $cmt  'SCAL' 'SOMMET' 0.D0 vithx ;
  898.     vmx = LAPC $mt vcx ;
  899.     vcy = vithy 'ET' ('KCHT' $bas 'SCAL' 'SOMMET' 0.D0) ;
  900.     vmy = LAPC $mt vcy ;
  901.   'SINON' ;
  902.     'SI' ('EGA' optlap 'TOTAL') ;
  903.       vcx = 'KCHT' $cmt 'SCAL' 'SOMMET' ('EXCO' 'UX' veuler) ;
  904.       vmx = LAPC $mt vcx ;
  905.       vcy = 'KCHT' $cmt 'SCAL' 'SOMMET' ('EXCO' 'UY' veuler) ;
  906.       vmy = LAPC $mt vcy ;
  907.     'SINON' ;
  908.       'ERREUR' 'Pas la bonne option pour le laplacien...' ;
  909.     'FINSI' ;
  910.   'FINSI' ;
  911. 'FINSI' ;
  912. vitx = 'NOMC' 'UX' vmx 'NATU' 'DISCRET' ;
  913. vity = 'NOMC' 'UY' vmy 'NATU' 'DISCRET' ;
  914. 'FORME' sauvform ;
  915. 'FINPROC' ('KCHT' $mt 'VECT' 'SOMMET' (vitx 'ET' vity)) ;
  916.  
  917. *
  918. *     on calcule w|surf avec v|surf en imposant w|surf // Oy
  919. *     => eqn convection sur la surface.
  920. *     Elle est discrétisée en différences finies centrées
  921. *     d'où le nom : ConVection Surface ordre 2 (schéma de
  922. *     type Leap-Frog)
  923. *     on interpole ensuite w sur le reste du maillage par
  924. *     une équation au laplacien
  925. *
  926. 'DEBPROC' CVS2LF veuler*'CHPOINT ' $c*'TABLE   ' ;
  927. 'ARGUMENT' optlap*'MOT     ' conflap*'CONFIGUR' ;
  928. $haut  = ($c . '$haut') ;   $bas   = ($c . '$bas') ;
  929. $lwall = ($c . '$lwall') ;  $rwall = ($c . '$rwall') ;
  930. $cmt = ($c . '$cmt') ;
  931. $mt = ($c . '$mt') ;
  932. haut  = ($haut . 'MAILLAGE') ;  bas   = ($bas . 'MAILLAGE') ;
  933. lwall = ($lwall . 'MAILLAGE') ; rwall = ($rwall . 'MAILLAGE') ;
  934. ** Calcul de withy, withx est nul
  935. * Calcul des paramètres pour le schéma aux différences finies
  936. nblh = ($haut . 'NELD') ;  nbph = ($haut . 'NPTD') ;
  937. dxh = ($haut . 'dxh')  ;  psm = ($haut . 'psm') ;
  938. psg = ($haut . 'psgg') ;  psd = ($haut . 'psdd') ;
  939. hc = 'KCHT' $haut 'SCAL' 'SOMMET' ('COORDONNEE' 2 haut) ;
  940. uc = 'KCHT' $haut 'SCAL' 'SOMMET' ('EXCO' 'UX' veuler) ;
  941. vc = 'KCHT' $haut 'SCAL' 'SOMMET' ('EXCO' 'UY' veuler) ;
  942. h = 'PROG' ;  u = 'PROG' ;  v = 'PROG' ;  sdtl = 'PROG' ;
  943. 'REPETER' BBB nbph ;
  944.   ppp = haut 'POINT' (&BBB) ;
  945.   h = h 'ET' ('PROG' ('EXTRAIRE' hc 'SCAL' ppp)) ;
  946.   u = u 'ET' ('PROG' ('EXTRAIRE' uc 'SCAL' ppp)) ;
  947.   v = v 'ET' ('PROG' ('EXTRAIRE' vc 'SCAL' ppp)) ;
  948. 'FIN' BBB ;
  949. vi   = 'EXTRAIRE' v   psm ;  ui = 'EXTRAIRE' u psm ;
  950. *vip1   = 'EXTRAIRE' v   psd ;  uip1 = 'EXTRAIRE' u psd ;
  951. *vim1   = 'EXTRAIRE' v   psg ;  uim1 = 'EXTRAIRE' u psg ;
  952. hi   = 'EXTRAIRE' h   psm ;
  953. hip1 = 'EXTRAIRE' h   psd ;
  954. him1 = 'EXTRAIRE' h   psg ;
  955. dxi  = 'EXTRAIRE' dxh psm ;  dxim1 = 'EXTRAIRE' dxh psg ;
  956. * points intermédaires
  957. dxi2 = dxi '*' dxi   ;  dxim12 = dxim1 '*' dxim1 ;
  958. dxp  = dxi '+' dxim1 ;  dxf    = dxim1 '*' dxi   ;
  959. dxm  = dxi '-' dxim1 ;
  960. hxi = ((((dxim12 '*' hip1) '-' (dxi2 '*' him1) '/' dxp) '+' ( dxm    '*' hi  )) '/' dxf) ;
  961. lwithy = vi '-' (ui '*' hxi) ;
  962. * 1er point
  963. lwithy = 'INSE' lwithy 1 ('EXTRAIRE' v 1) ;
  964. * dernier point
  965. lwithy = (lwithy 'ET' ('PROG' ('EXTRAIRE' v nbph))) ;
  966.   sdtl  = (ui '/' dxim1) '+' (ui '/' dxi) ;
  967.   dtmax = 1.D0 '/' ('MAXIMUM' sdtl) ;
  968.   'MESSAGE' 'Procédure CVS2LF : dtmax =' dtmax ;
  969. *withyg = 'MANUEL' 'CHPO' ($haut . 'SOMMET') 1 'SCAL' lwithy ;
  970. withyg = 'MANUEL' 'CHPO' ('CHAN' 'POI1' haut) 1 'SCAL' lwithy ;
  971. withy = 'KCHT' $haut 'SCAL' 'SOMMET' withyg ;
  972. * Interpolation de wity a partir de withy sur le reste du 
  973. * maillage
  974. sauvform = 'FORME' ;
  975. 'FORME' conflap ;
  976. 'SI' ('EGA' optlap 'INTBOR') ;
  977.   vlwx vlwy = IVECT $lwall veuler ;
  978.   vrwx vrwy = IVECT $rwall veuler ;
  979.   vcy = 'KCHT' $cmt  'SCAL' 'SOMMET' 0.D0 withy vlwy vrwy ;
  980.   vmy = LAPC $mt vcy ;
  981. 'SINON' ;
  982.   'SI' ('EGA' optlap 'BORLIB') ;
  983.     vcy = withy 'ET' ('KCHT' $bas 'SCAL' 'SOMMET' 0.D0) ;
  984.     vmy = LAPC $mt vcy ;
  985.   'SINON' ;
  986.     'ERREUR' 'Pas la bonne option pour le laplacien...' ;
  987.   'FINSI' ;
  988. 'FINSI' ;
  989. witx = 'NOMC' 'UX' ('KCHT' $mt 'SCAL' 'SOMMET' 0.D0)  'NATU' 'DISCRET' ;
  990. wity = 'NOMC' 'UY' vmy 'NATU' 'DISCRET' ;
  991. 'FORME' sauvform ;
  992. 'FINPROC' ('KCHT' $mt 'VECT' 'SOMMET' (witx 'ET' wity)) ;
  993.  
  994. *
  995. *     Interpolation d'un vecteur entre deux points sur
  996. *     un segment. Un chpoint VECT => 2 chpos SCAL
  997. *
  998. 'DEBPROC' IVECT $dom*'TABLE   ' chv*'CHPOINT' ;
  999. dom = ($dom . 'MAILLAGE') ;  nbd = ($dom . 'NPTD') ;
  1000. cyd = 'COORDONNEES' 2 dom ;
  1001. yd  = 'PROG' ;
  1002. 'REPETER' BBB nbd ;
  1003.    ppp = dom 'POINT' (&BBB) ;
  1004.    yd = yd 'ET' ('PROG' ('EXTRAIRE' cyd 'SCAL' ppp )) ;
  1005. 'FIN' BBB ;
  1006. di = dom 'POINT' 'INITIAL' ;  df = dom 'POINT' 'FINAL' ;
  1007. yi   = 'EXTRAIRE' yd 1 ;
  1008. yf   = 'EXTRAIRE' yd nbd ;
  1009. vdxi = 'EXTRAIRE' chv 'UX' di ;
  1010. vdxf = 'EXTRAIRE' chv 'UX' df ;
  1011. vdyi = 'EXTRAIRE' chv 'UY' di ;
  1012. vdyf = 'EXTRAIRE' chv 'UY' df ;
  1013. ax = (vdxf '-' vdxi) '/' (yf '-' yi) ;
  1014. bx = vdxi '-' (ax '*' yi) ;
  1015. lvdx = 'PROG' 'LINE' 'A' ax 'B' bx yd ;
  1016. ay = (vdyf '-' vdyi) '/' (yf '-' yi) ;
  1017. by = vdyi '-' (ay '*' yi) ;
  1018. lvdy = 'PROG' 'LINE' 'A' ay 'B' by yd ;
  1019. dom1 = chan 'POI1' dom ;
  1020. vdx = 'MANUEL' 'CHPO' dom1 1 'SCAL' lvdx 'NATURE' 'DISCRET' ;
  1021. vdy = 'MANUEL' 'CHPO' dom1 1 'SCAL' lvdy 'NATURE' 'DISCRET' ;
  1022. 'FINPROC' vdx vdy ;
  1023.  
  1024. *
  1025. *     Résolution d'une équation scalaire en implicite ;
  1026. *     Laplacien + conditions Dirichlet ou Neumann
  1027. *
  1028. 'DEBPROC' LAPC $dom*'TABLE   ' condlim*'CHPOINT ' ;
  1029. rt = 'EQEX' $dom 'OPTI' 'EF' 'IMPL' 'ZONE' $dom 'OPER' 'LAPN' 1.D0 'INCO' 'SCAL' ;
  1030. rt . 'INCO' = 'TABLE' 'INCO' ;
  1031. rt . 'INCO' . 'SCAL' = 'KCHT' $mt 'SCAL' 'SOMMET' 0.D0 ;
  1032. rt . 'CLIM' = condlim ;
  1033. EXEC rt ;
  1034. *  'OUBLIER' $mt 'MATESI' ;
  1035. *  'OUBLIER' $mt 'XXDXDY' ; 'OUBLIER' $mt 'XXVOLUM' ;
  1036. *  'OUBLIER' $mt 'XXDIAME' ; 'OUBLIER' $mt 'XXDIEMIN' ;
  1037. 'FINPROC' (rt . 'INCO' . 'SCAL') ;
  1038.  
  1039. *
  1040. *     Procédure qui étend un champoint VECT SOMMET a tous les
  1041. *     points (FACE, CENTRE...) du maillage total
  1042. *
  1043. 'DEBPROC' ETENDRE $c*'TABLE   ' depli*'CHPOINT ' ;
  1044. $mt = ($c . '$mt') ;  $lmt = ($c . '$lmt') ;
  1045. deplic = 'NOEL' $mt depli ;
  1046. deplil = 'KCHT' $lmt 'VECT' 'SOMMET' depli ;
  1047. deplif = 'NOEL' $lmt deplil ;
  1048. deplif = 'KCHT' $mt  'VECT' 'FACE'   deplif ;
  1049. deplic = 'KCHT' $mt  'VECT' 'CENTRE' deplic ;
  1050. 'FINPROC' (depli 'ET' deplic 'ET' deplif) ;
  1051.  
  1052. *
  1053. *     Calcul du pas de temps dit "de maillage"
  1054. *
  1055. 'DEBPROC' CINVPDTM dvit*'CHPOINT ' $mt*'TABLE   ' ;
  1056. gdvitx = ('EXCO' 'UX' ('KOPS' ('KCHT' $mt 'SCAL' 'SOMMET' ('EXCO' 'UX' dvit)) 'GRAD' $mt)) ;
  1057. gdvity = ('EXCO' 'UY' ('KOPS' ('KCHT' $mt 'SCAL' 'SOMMET' ('EXCO' 'UY' dvit)) 'GRAD' $mt)) ;
  1058. lmg = 'PROG' ('MINIMUM' gdvitx) ('MINIMUM' gdvity) ;
  1059. 'FINPROC' ('MINIMUM' lmg) ;
  1060.  
  1061. *
  1062. *     Construit une somme de déformées pour faire un film
  1063. *     avec les champs de vitesse sauvegardés
  1064. *
  1065. 'DEBPROC' FILMVIT $s*'TABLE   ' $v*'TABLE   ' ;
  1066. indmax = ($s . 'ntrc') ;
  1067. $s . 'anim' = 'TABLE' 'ANIM' ;
  1068. 'REPETER' BCL1 indmax ;
  1069.   ind = &BCL1 ;
  1070.   dx = ($s . 'DELTAX' . ind) ;
  1071.   vv = 'VECTEUR' ($s . 'VITESSE' . ind)  ampvit 'UX' 'UY' 'JAUNE' ;
  1072.   $s . 'anim' . ind =  ('DEFORME' ($cavite . '$mt' . 'MAILLAGE') dx 1.D0 vv) ;
  1073. 'FIN' BCL1 ;
  1074. 'FINPROC' (@STBL ($s . 'anim')) ;
  1075. *
  1076. *     Tracé d'historiques
  1077. *
  1078. 'DEBPROC' TRHIST ;
  1079. ltps  = ($hist . 'ltps') ;
  1080. liter = ($hist . 'liter') ; ldv  = ($hist . 'ldv') ;
  1081. lvol  = ($hist . 'lvol') ;  ltps = ($hist . 'ltps') ;
  1082. lzr   = ($hist . 'lzr') ;   lzq  = ($hist . 'lzq') ;
  1083. $hl = ($hist . 'lpdt') ;
  1084. lpdt  = ($hl . 'min') ;
  1085. *evpdt = 'EVOL' 'MANU' 'TEMPS' ltps 'PDT' lpdt ;
  1086. *'DESSIN' evpdt 'MIMA' 'TITR' 'Pas de temps' ;
  1087.  
  1088. ltco = ('LOG' (($hl . 'conv') '*' ($pdt . 'alfa')))  '/' ('LOG' 10.) ;
  1089. ltdi = ('LOG' (($hl . 'diff') '*' ($pdt . 'alfa')))  '/' ('LOG' 10.) ;
  1090. ltma = ('LOG' (($hl . 'mail') '*' ($pdt . 'gamma')))  '/' ('LOG' 10.) ;
  1091. ltmo = ('LOG' ($hl . 'moy')) '/' ('LOG' 10.) ;
  1092. ltmi = ('LOG' ($hl . 'min')) '/' ('LOG' 10.) ;
  1093. evco = 'COULEUR' ('EVOL' 'MANU' 'TEMPS' ltps 'dtconv' ltco) 'BLEU' ;
  1094. evdi = 'COULEUR' ('EVOL' 'MANU' 'TEMPS' ltps 'dtdiff' ltdi) 'ROUG' ;
  1095. evma = 'COULEUR' ('EVOL' 'MANU' 'TEMPS' ltps 'dtdiff' ltma) 'JAUN' ;
  1096. evmo = 'COULEUR' ('EVOL' 'MANU' 'TEMPS' ltps 'dtdiff' ltmo) 'VERT' ;
  1097. evmi = 'COULEUR' ('EVOL' 'MANU' 'TEMPS' ltps 'dtdiff' ltmi) 'BLAN' ;
  1098. tabt = 'TABLE' ;  tabt . 'TITRE' = 'TABLE' ;
  1099. tabt . 1 = 'TIRL' ;  tabt . 'TITRE' . 1 = 'Pdt diffusion' ;
  1100. tabt . 2 = 'TIRR' ;  tabt . 'TITRE' . 2 = 'Pdt maillage' ;
  1101. tabt . 3 = 'TIRM' ;  tabt . 'TITRE' . 3 = 'Pdt moyen' ;
  1102. tabt . 4 = 'NOLI' ;  tabt . 'TITRE' . 4 = 'Pdt minimum' ;
  1103. tabt . 5 = 'TIRC' ;  tabt . 'TITRE' . 5 = 'Pdt convection' ;
  1104.  
  1105. * On ne dessine pas le pdt de conv en lagrangien (infini)
  1106. 'SI' ('NON' ('EGA' vmail 'LAGR')) ;
  1107.   evtot = (evdi 'ET' evma 'ET' evmo 'ET' evmi 'ET' evco) ;
  1108. 'SINON' ;
  1109.   evtot = (evdi 'ET' evma 'ET' evmo 'ET' evmi) ;
  1110. 'FINSI' ;
  1111. 'DESSIN' evtot 'MIMA' 'TITR' 'Pas de temps' 'LEGE' tabt ;
  1112.  
  1113. maxy = (('ENTIER' ('MAXIMUM' liter)) '+' 1.) ;
  1114. tity = 'CHAINE' 'Nb Iter (Total = ' ('ENTIER'  ('SOMM' liter))  ' ;  Cpu = ' ($pdt . 'tpscpu') ')' ;
  1115. evit = 'EVOL' 'MANU' 'TEMPS' ltps 'NITER' liter ;
  1116. 'DESSIN' evit 'YBOR' 1. maxy 'MIMA' 'TITR' tity ;
  1117.  
  1118. ldv  = (('LOG' ldv) '/' ('LOG' 10.)) ;
  1119. evdv = 'EVOL' 'MANU' 'TEMPS' ltps 'DIVV (log)' ldv ;
  1120. 'DESSIN' evdv 'MIMA' 'TITR' 'Div v' ;
  1121.  
  1122. 'FINPROC' ;
  1123.  
  1124. *                 
  1125. **  Fin de la définition des procédures
  1126. **************************************************************
  1127.  
  1128. **************************************************************
  1129. **  Début du 'main'
  1130. *
  1131. **** Définition des constantes
  1132. * l, h, nl, nh : données pour la boite carrée
  1133. * dh, db, dg, dd : densités (facultatives) pour la boite
  1134. * egeom : erreur géométrique pour 'ELIMINATION' 
  1135. * cvg : critère de convergence pour la boucle d'itérations
  1136. * nitmax : nombre maximum d'itérations
  1137. * minnz : zéro pour ne pas avoir de pbs avec / , LOG ...
  1138. * nu, g, rho : paramètres pour l'eau 
  1139. * minnu, ming : squizzage de certains termes dans NS
  1140. * alfa : paramètre pour NS
  1141. * gamma : coeff. gamma = 1 => le maillage dégénère
  1142. * psauv : fréquence de sauvegarde des chpoints 
  1143. * tmax : temps d'arret
  1144. * pt0 : coeff intensité du cos de pression sur la surface
  1145. * dtini : 1er pas de temps
  1146. * poids du dt initial dans dtmoy
  1147. l = 4.8D0 ;  nl = 14 ;
  1148. h = 4.0D0 ;  nh = 22 ;
  1149. db = 3.0D0 ;  dh = 0.3D0 ;
  1150. dg = 2.0D0 ;  dd = 0.5D0 ;
  1151. egeom = 1.D-4 ;
  1152. cvg  = 1.D-6 ;
  1153. npdtmax = -1 ;
  1154. nitmax  = -1 ;
  1155. minnz = 1.D-13 ;
  1156. nu   = 1.D-2 ;      minnu = nu '*' minnz ;
  1157. g = (0.D0 -1.D0) ;  ming = (0.D0 0.D0) ;
  1158. rho  = 1.D3 ;
  1159. alfa = 0.8D0 ;
  1160. gamma = 0.01D0 ;
  1161. psauv = 1 ;
  1162. *pt0 = -1.D0 ;   dtini = 0.25D0 ;
  1163. pt0 = -1.D0 ;   dtini = 2.5D-6 ;
  1164. poini = 0 ;
  1165. 'SI' ('NON' court) ;
  1166.   tmax = 4.0D0 ;  ptrc = 1.0D0 ;  ttrc = ptrc ;
  1167. *  tmax = 60.0D0 ;  ptrc = 10.0D0 ;  ttrc = ptrc ;
  1168. *  fsauv = 'CHAINE'  'nlt' ('ENTIER' tmax) vmail optlapn
  1169. *   'a' ('ENTIER' (alfa '*' 10)) 
  1170. *   'c' ('ENTIER' ((LOG10 cvg) '-' 0.5D0)) 
  1171. *   'g' ('ENTIER' ((LOG10 gamma) '-' 0.5D0)) 'po' poini '.sauv' ;
  1172. 'SINON' ;
  1173.   tmax = 0.25D0 ;  ptrc = 0.05D0 ;  ttrc = ptrc ;
  1174. *  tmax = 0.25D0 ;  ptrc = 0.25D0 ;  ttrc = ptrc ;
  1175. *  fsauv = '/test4/gounand/ale/dumb' ;
  1176. 'FINSI' ;
  1177.  
  1178. **** Création des maillages
  1179. $cavite = 'TABLE' 'MAILLAGE' ;
  1180. $cavite = CREMAIL $cavite l h nl nh egeom ;
  1181. *  'DX' dg dd 'DY' db dh ;
  1182.  
  1183. 'SI' ('EGA' vmail 'CVS2LF') ;
  1184.   CDXHAUT ($cavite . '$haut') ;
  1185. 'FINSI' ;
  1186.  
  1187. **** Initialisation des tables diverses
  1188. $var  = 'TABLE' 'VARIABLE' ;
  1189. $pdt  = 'TABLE' 'PASDETPS' ;
  1190. $hist = 'TABLE' 'HIST' ;
  1191. $sauv = 'TABLE' 'SAUV' ;
  1192. INITVAR $var $cavite ;
  1193. INITPDT $pdt alfa gamma dtini poini ;
  1194. INITHIST $hist ;
  1195. 'SI' (graph 'ET' inta) ;
  1196.     INITSAUV $sauv ;
  1197. 'FINSI' ;
  1198. 'TEMPS' 'ZERO' ;  SUMMARY ;  'TEMPS' 'ZERO' ;
  1199.  
  1200.  
  1201. **** Premier pas de temps
  1202. * Calcul de la contrainte à appliquer sur la surface
  1203. rits = CALCON ($cavite . '$haut') pt0 ;
  1204. * Calcul des champs à l'instant dtini
  1205. *           Peu d'avance en temps => pas de diffusion
  1206. *           Dirac => champ de pression du à g négligeable
  1207. $var $pdt = PREPAS $var $pdt $cavite rits minnu ming ;
  1208. * Pour ne pas avoir une estimation de pas de temps avec 0
  1209. $var . 'wnp05i' = ('TEXTE' vmail)  ($var . 'unp05i') $cavite optlapn ($cavite . 'cnfini') ;
  1210. $var . 'erri' =  CALCERR ($var . 'unp05i') ($var . 'unp05im1') 'MINMAX' ;
  1211. $pdt . 'nupdt' = (($pdt . 'nupdt') '+' 1) ;
  1212. $var . 'unm05' = ('COPIE' ($var . 'unp05i')) ;
  1213. $var . 'wnm05' = ('COPIE' ($var . 'wnp05i')) ;
  1214. ampvit = l '/' (10.D0 '*' ('MAXIMUM' ($var . 'unm05') 'ABS')) ;
  1215. MESI $var $pdt ;  MEST $var $pdt ;
  1216. 'SI' graph ;
  1217.   TRTOUT ;
  1218.   'SI' inta ;
  1219.     MAJSAUV $sauv $var $pdt ;
  1220.   'FINSI' ;
  1221. 'FINSI' ;
  1222.  
  1223. ***** Boucle sur les pas de temps
  1224. 'REPETER' TITER npdtmax ;
  1225. **** Boucle d'itérations
  1226. * Initialisation
  1227.   $pdt . 'iter' = 0 ;
  1228.   $var . 'unp05i' = ('COPIER'  ($var . 'unm05')) ;
  1229.   $var . 'wnp05i' = ('COPIER'  ($var . 'wnm05')) ;
  1230.   'FORME' ($var . 'cnfnm05') ;
  1231.   $var . 'cnfnp05i' = 'FORME' ;
  1232. * Estimation du pas de temps pour le maillage
  1233.   pm = CINVPDTM ($var . 'wnm05') ($cavite . '$mt') ;
  1234.   $pdt . 'mail' = (-1.0D0 '/' pm) ;
  1235.   MAJPDT $pdt ;
  1236. * Boucle
  1237.   'REPETER' BITER nitmax ;
  1238. * On bouge le maillage -> xn+1/2(i)
  1239.     'FORME' ($var . 'cnfnp05i') ;
  1240.     $var . 'wnp05i' = ('TEXTE' vmail)  ($var . 'unp05i')  $cavite optlapn ($cavite . 'cnfini') ;
  1241.     dxns2 =  (ETENDRE $cavite  ('KOPS' (0.5D0 '*' ($pdt . 'min'))  '*'   ($var . 'wnp05i'))) ;
  1242.     'FORME' ($var . 'cnfn') ;
  1243.     $var . 'cnfnp05i' = 'FORME' dxns2 ;
  1244. * Calcul de un+1/2(i)
  1245.     $var $pdt = CVIT $cavite $var $pdt nu g discconv ;
  1246.     $var . 'erri' = (CALCERR ($var . 'unp05i')  ($var . 'unp05im1') 'MINMAX') ;
  1247.     $pdt . 'iter'   = (($pdt . 'iter') '+' 1) ;
  1248.     MESI $var $pdt ;
  1249.     'SI' (($var . 'erri') '<' cvg) ;
  1250.       'QUITTER' BITER ;
  1251.     'FINSI' ;
  1252.   'FIN' BITER ;
  1253. * On avance d'un pas de temps :
  1254. * estimation des maillages xn+1/2 
  1255. *                       et xn+1 avec le dernier u calculé
  1256.   'FORME' ($var . 'cnfnp05i') ;
  1257.   $var . 'wnp05i' = ('TEXTE' vmail)  ($var . 'unp05i')  $cavite optlapn ($cavite . 'cnfini') ;
  1258.   dxns2 = (ETENDRE $cavite  ('KOPS' (0.5D0 '*' ($pdt . 'min'))  '*'  ($var . 'wnp05i'))) ;
  1259.   $var . 'deltaxn' = (($var . 'deltaxn')  '+' (2.0D0 '*' dxns2)) ;
  1260.   $pdt . 'nupdt'   = (($pdt . 'nupdt') '+' 1) ;
  1261.   $pdt . 'tphyn'   = (($pdt . 'tphyn') '+' ($pdt . 'min')) ;
  1262.   $var . 'unm05' = ('COPIE' ($var . 'unp05i')) ;
  1263.   $var . 'wnm05' = ('COPIE' ($var . 'wnp05i')) ;
  1264.   'FORME' ($var . 'cnfn') ;
  1265.   $var . 'cnfnm05' = 'FORME' dxns2 ;
  1266. * Sauvegarde éventuelle des champs à n-1/2
  1267.   'SI' (('MULT' ($pdt . 'nupdt') psauv) 'ET' graph 'ET' inta) ;
  1268.     MAJSAUV $sauv $var $pdt ;
  1269.   'FINSI' ;
  1270.   $var . 'cnfn'    = 'FORME' dxns2 ;
  1271. * Sauvegarde du reste (volume, pos des points) à n
  1272.   MAJHIST $hist $cavite $pdt $var ;
  1273.   'SI' ((($pdt . 'tphyn') '>EG' ttrc) 'ET' graph) ;
  1274.     TRTOUT ; ttrc = ttrc '+' ptrc ;
  1275.   'FINSI' ;
  1276.   MEST $var $pdt ;
  1277.   'SI' (($pdt . 'tphyn') '>' tmax) ;
  1278.     'QUITTER' TITER ;
  1279.   'FINSI' ;
  1280.   'MENAGE' ;
  1281. 'FIN' TITER ;
  1282.  
  1283. TABTPS = 'TEMPS' 'NOEC' ;
  1284. $pdt . 'tpscpu' = TABTPS.'TEMPS_CPU'.'INITIAL' ;
  1285.  
  1286. ** Sauvegarde
  1287. 'SI' ('EXISTE' ($cavite . '$mt') 'MATESI') ;
  1288.   'OUBLIER' ($cavite . '$mt') 'MATESI' ;
  1289. 'FINSI' ;
  1290. 'MENAGE' ;
  1291. *'OPTION' 'SAUV' fsauv ;
  1292. *'SAUVER' ;
  1293.  
  1294. *
  1295. **  Fin du 'main'
  1296. **************************************************************
  1297.  
  1298. **************************************************************
  1299. **  Début du post-traitement
  1300. *
  1301. 'SI' graph ;
  1302.   'SI' inta ;
  1303.     'FORME' ($cavite . 'cnfini') ;
  1304.     fv = FILMVIT $sauv $var ;
  1305.     'TRACER' fv 'ANIME' ;
  1306.     'OUBLIER' fv ;  MENAGE ;
  1307.   'FINSI' ;
  1308.   TRHIST;
  1309. 'FINSI' ;
  1310.  
  1311. *** Les tests...
  1312. ltps = ($hist . 'ltps') ;  lvol = ($hist . 'lvol') ;
  1313. lzr   = ($hist . 'lzr') ;   lzq  = ($hist . 'lzq') ;
  1314. ** Ici, on teste la conservation du volume du domaine
  1315. critere1 = 'MAXIMUM' lvol 'ABS' ;
  1316. * Résultat obtenu le 29/12/97 :
  1317. * critere1 = .9806970050855623D-14
  1318. limite1 = 1.D-10 ;
  1319. 'SI' graph ;
  1320.   evvol = 'EVOL' 'MANU' 'TEMPS'  ltps 'VOLUME' lvol ;
  1321.   'DESSIN' evvol 'AXES' 'MIMA'  'TITR' 'Ecart au volume initial' ;
  1322. 'FINSI' ;
  1323.  
  1324. ltrram = ('PROG' 0.0 0.4 0.9 1.55 2.35 3.1 3.9 4.7 5.65 6.75 7.75 8.6 9.25 9.95 10.65 11.35 11.9 12.45)  '*' (4. '/' 12.) ;
  1325. ltqram = ('PROG' 0.0 0.5 1. 1.6 2.25 2.9 3.7 4.4 5.1 6.1 7.1 7.9 8.6 9.2 9.9 10.5 11.1 11.6 12. 12.45)  '*' (4. '/' 12.) ;
  1326. lzrram = ('PROG' 0.0 0.45 1.05 1.6 2.4 2.95 3.5 3.8 3.95 3.95  3.75 3.25 2.8 2.15 1.6 0.9 0.3 -0.25)  '*' ((-2.) '/' 11.8) ;
  1327. lzqram = ('PROG' 0.0 0.75 1.35 2.3 3.25 4.15 5.15 5.8 6.3 6.5 6.4 5.85 5.25 4.5 3.8 3. 2.25 1.45 0.75 0.25)  '*' (2. '/' 11.8) ;
  1328. * Comparaison des profils zQ(t) et zR(t) de Ramaswami aux notres
  1329. * on les remet sur ltps par interpolation linéaire
  1330. lzrint = 'IPOL' ltps ltrram lzrram ;
  1331. lzqint = 'IPOL' ltps ltqram lzqram ;
  1332. * échelles pour r et q
  1333. echr = 'MINIMUM' ('PROG' ('MAXIMUM' lzr 'ABS') ('MAXIMUM' lzrram 'ABS') ('MAXIMUM' lzrint 'ABS')) ;
  1334. echq = 'MINIMUM' ('PROG' ('MAXIMUM' lzq 'ABS') ('MAXIMUM' lzqram 'ABS') ('MAXIMUM' lzqint 'ABS')) ;
  1335. echrq = 'MINIMUM' ('PROG' echr echq) ;
  1336. critere2 = ('MAXIMUM' (lzrint '-' lzr) 'ABS') '/' echrq ;
  1337. critere3 = ('MAXIMUM' (lzqint '-' lzq) 'ABS') '/' echrq ;
  1338. * Résultats obtenus le 29/12/97 :
  1339. * critere2 = .4174969242483980D-01
  1340. * critere3 = .9776923547482334D-01
  1341. limite2 = 0.05D0 ;
  1342. limite3 = 0.12D0 ;
  1343.  
  1344. 'SI' graph ;
  1345.   evzq = 'COULEUR'  ('EVOL' 'MANU' 'TEMPS' ltps 'zQ' lzq) 'ROSE' ;
  1346.   evzr = 'COULEUR'  ('EVOL' 'MANU' 'TEMPS' ltps 'zR' lzr) 'TURQ' ;
  1347.   evzqint = 'COULEUR'  ('EVOL' 'MANU' 'TEMPS' ltps 'zQint' lzqint) 'ROUG' ;
  1348.   evzrint = 'COULEUR'  ('EVOL' 'MANU' 'TEMPS' ltps 'zRint' lzrint) 'BLEU' ;
  1349.   evzqram = 'COULEUR'  ('EVOL' 'MANU' 'TEMPS' ltqram 'zQram' lzqram) 'JAUN' ;
  1350.   evzrram = 'COULEUR'  ('EVOL' 'MANU' 'TEMPS' ltrram 'zRram' lzrram) 'VERT' ;
  1351.   tabevzqr = 'TABLE' ;  tabevzqr . 'TITRE' = 'TABLE' ;
  1352.   tabevzqr . 1 = 'NOLI' ;
  1353.   tabevzqr . 'TITRE' . 1 = 'zQ (Castem 2000)' ;
  1354.   tabevzqr . 2 = 'NOLI' ;
  1355.   tabevzqr . 'TITRE' . 2 = 'zR (Castem 2000)' ;
  1356.   tabevzqr . 3 = 'TIRR' ;
  1357.   tabevzqr . 'TITRE' . 3 = 'zQ (Rama. interpolé)' ;
  1358.   tabevzqr . 4 = 'TIRR' ;
  1359.   tabevzqr . 'TITRE' . 4 = 'zR (Rama. interpolé)' ;
  1360.   tabevzqr . 5 = 'NOLI MARQ CARR' ;
  1361.   tabevzqr . 'TITRE' . 5 = 'zQ (Ramaswami)' ;
  1362.   tabevzqr . 6 = 'NOLI MARQ CARR' ;
  1363.   tabevzqr . 'TITRE' . 6 = 'zR (Ramaswami)' ;
  1364.   'DESSIN' (evzq 'ET' evzr 'ET' evzqint 'ET' evzrint 'ET' evzqram 'ET' evzrram)  'MIMA' 'TITR' 'Amplitude de Q et R'  'XBOR' 0. tmax 'LEGE' 'AXES' tabevzqr ;
  1365. 'FINSI' ;
  1366. *
  1367. ** Affichage des messages et des erreurs éventuelles
  1368. *
  1369.   'OPTION' 'ECHO' 0 ;
  1370.   SUMMARY ;
  1371.   'SAUTER' 2 'LIGNE' ;
  1372.   'MESSAGE' 'Debriefing :' ;
  1373.   'MESSAGE' '------------' ;
  1374.   'MESSAGE' 'Conservation relative du volume' ;
  1375.   'MESSAGE' 'Valeur / Limite :' ;
  1376.   'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere1 ' / ' limite1) ;
  1377.   'SAUTER' 'LIGNE' ;
  1378.   'MESSAGE' 'Comparaison d ordonnées (points extremes de la surface)' ;
  1379.   'MESSAGE' 'avec une autre simulation numérique (Ramaswamy)' ;
  1380.   'MESSAGE' 'zR(t) (Erreur / Limite) :' ;
  1381.   'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere2 ' / ' limite2) ;
  1382.   'MESSAGE' 'zQ(t) (Erreur / Limite) :' ;
  1383.   'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere3 ' / ' limite3) ;
  1384.   'SAUTER' 'LIGNE' ;
  1385.   'MESSAGE' ('CHAINE' 'Temps CPU / Temps CPU (29/12/97) (interactif) : ') ;
  1386.   'MESSAGE' ('CHAINE' ('ENTIER' (($pdt . 'tpscpu') '/' 100.0D0)) ' / ~28' ' secondes') ;
  1387. 'SI' ('OU' (critere1 '>' limite1) ('OU' (critere2 '>' limite2) (critere3 '>' limite3))) ;
  1388.   'ERREUR' 5 ;
  1389. 'FINSI' ;
  1390. 'MESSAGE' 'Tout s est bien passé' ;
  1391. 'OPTION' 'ECHO' 1 ;
  1392. 'SI' inta ;  'OPTION' 'DONN' 5 ;  'FINSI' ;
  1393. *
  1394. **  Fin du post-traitement
  1395. **************************************************************
  1396. 'FIN' ;
  1397.  
  1398.  
  1399.  
  1400.  
  1401.  
  1402.  
  1403.  
  1404.  
  1405.  
  1406.  
  1407.  
  1408.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales