Télécharger ale_mecaflu.dgibi

Retour à la liste

Numérotation des lignes :

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

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