Télécharger la plaquette Cast3M


fr en

Télécharger cl_E_4.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  cl_E_4.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *C***********************************************************************
  5. *C*                                                                     *
  6. *C* PROJET      : Opérateur CLMI                                        *
  7. *C* NOM         : cl_E_4  .dgibi                                        *
  8. *C* DESCRIPTION : Jeu de données pour le calcul de la couche limite     *
  9. *C*               turbulente en utilisant les relations de fermeture de *
  10. *C*               Cousteix - Cas de la plaque plane avec gradient de    *
  11. *C*               pression positif constant -Méthode de Head            *
  12. *C* LANGAGE     : Gibiane                                               *
  13. *C* AUTEUR      : Guillaume VENCO - DRN/DMT/SEMT/LTMF                   *
  14. *C*                                                                     *
  15. *C***********************************************************************
  16. *C*                                                                     *
  17. *C* APPELES     : Opérateur CLMI                                        *
  18. *C*                                                                     *
  19. *C***********************************************************************
  20. *C*                                                                     *
  21. *C*                                                                     *
  22. *C***********************************************************************
  23. *C*                                                                     *
  24. *C* VERSION     : 19/05/2000                                            *
  25. *C* CREATION    : 19/04/2000                                            *
  26. *C*                                                                     *
  27. *C***********************************************************************
  28. * 19/04/2000 Source: TTMF2 Nom: Cl_E_3.dgibi
  29. *Couche limite plaque plane turbulente 
  30. * avec gradient de pression négatif constant
  31. * Fermetures de HEAD
  32.  
  33. 'OPTION' dime 2 elem seg2 ;
  34. GRAPH=FAUX ;
  35. COMPLET=FAUX;
  36. *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  37. *>>>>>>>>>>>>>Maillage et domaine>>>>>>>>>>>>>>>>>>>>>>>>>>>
  38. *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  39. P1 = 0.15 0.;
  40. p2 =  1. 0. ;
  41. nb1 = 10;
  42. p1p2 = 'DROIT' p1 p2 'DINI' 0.0001 'DFIN' 0.01 ;
  43. *p1p2 = 'DROIT' p1 p2 20;
  44. mt = p1p2 ;
  45. * 'TRACER'  mt ;
  46.  
  47. *Définition du modèle
  48. Mmt = 'CHANGER'  mt 'QUAF' ;
  49. Mail = 'CHANGER' Mmt 'POI1';
  50. Pt0 = Mail 'ELEM' 1;
  51. $mt = 'MODELISER'  Mmt 'NAVIER_STOKES' 'LINE' ;
  52.  
  53. *Données 
  54. Ue0 = 1.;
  55. NU = 1.d-6 ;
  56. *Valeurs initiales données par la loi de puissance
  57. D20=4.33d-4;
  58. Q10=6.37*D20;
  59.  
  60. *Pas de temps et nombre de pas
  61. DT = 0.1;
  62.  
  63. NBIT =20;
  64. Si COMPLET;
  65. NBIT =200;
  66. Finsi;
  67.  
  68. *Champ de vitesse à l'extérieur
  69.  
  70. *vitesse en puissance
  71. K=1.;
  72. m=(1.);
  73. m1= 'DOMA' $mt 'SOMMET';
  74. X = 'COORDONNEE' 1 m1;
  75. *X1 = (X**m);
  76. *UEX = K*X1;
  77. X1 = ((X-0.15)**m);
  78. UEX = 2.-(K*X1);
  79. UE = 'KCHT' $mt scal 'SOMMET' UEX;
  80. X2 = ((X-0.15)**(m-1));
  81. DUEX = (-1.)*K*m*X2 ;
  82. DUEX=-1.;
  83. DUE = 'KCHT' $mt scal 'SOMMET' DUEX;
  84.  
  85. *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  86. *>>>>>>>>>>>>Création de la table RV et de la table RK>>>>>>
  87. *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  88.  
  89. *Dans CLMI le 1er coef correspond au type de fermeture
  90. *le 2ème coef correspond au n°de l'équation à traiter
  91. *1=QDM 2=Energie cinétique
  92.  
  93. rv = 'EQEX'  'OMEGA' 1.  'NITER' 1  'ITMA' 0   'FIDT' 1
  94.      'OPTI' 'EF' 'IMPL' 'CENTREE' 'BDF2'
  95.      'ZONE' $mt
  96.      'OPER' 'DFDT' 1. 'D2NM' 'D2N2' DT  'INCO' 'D2' ;
  97. *     'OPER' 'DFDT' 1. 'D2NM' DT  'INCO' 'D2' ;
  98.  
  99. rv = 'EQEX' rv
  100.      'OPTI' 'EF' 'IMPL' 'SUPG'
  101.       'ZONE'  $mt
  102.      'OPER'  'CLMI' 5 1 UE DUE 'D2NM' 'Q1NM'
  103.      'INCO' 'D2' ;
  104.  
  105. rv = 'EQEX' rv
  106.      'OPTI' 'EF' 'IMPL' 'CENTREE' 'BDF2'
  107.      'ZONE' $mt
  108.      'OPER' 'DFDT' 1. 'Q1NM' 'Q1N2' DT  'INCO' 'Q1' ;
  109. *     'OPER' 'DFDT' 1. 'D1NM' DT  'INCO' 'D1' ;
  110.  
  111. rv = 'EQEX' rv
  112.      'OPTI' 'EF' 'IMPL' 'SUPG'
  113.      'ZONE'  $mt
  114.      'OPER'  'CLMI' 5 3 UE DUE 'Q1NM' 'D2NM'
  115.      'INCO' 'Q1' ;
  116.  
  117. rv = 'EQEX' rv
  118.       'CLIM' 'D2' 'TIMP' pt0 D20
  119.       'CLIM' 'Q1' 'TIMP' pt0 Q10;
  120.  
  121.  
  122. rv.inco = 'TABLE'  inco ;
  123. rv.inco.'D2'=  'KCHT'  $mt scal sommet D20;
  124. rv.inco.'D2NM'=  'KCHT'  $mt scal sommet D20;
  125. rv.inco.'D2N2'=  'KCHT'  $mt scal sommet D20;
  126. rv.inco.'Q1'=  'KCHT'  $mt scal sommet Q10;
  127. rv.inco.'Q1NM'=  'KCHT'  $mt scal sommet Q10;
  128. rv.inco.'Q1N2'=  'KCHT'  $mt scal sommet Q10;
  129.  
  130. ltps = 'PROG' 0. ;
  131. lteta = 'PROG' 0. ;
  132. TPS = 0. ;
  133. nupt = 0 ;
  134. mts1 = 'DOMA'  $mt 'SOMMET' ;
  135. pt1 = 'POIN'  mts1 2 ;
  136. 'LISTE'  pt1 ;
  137.  
  138. *'OPTION'  donn 5 ;
  139.  
  140. 'REPETER'  BCLT nbit ;
  141. nupt=nupt+1;
  142. tps=tps+dt ;
  143. mess ' Pas de temps DT=' dt ' Temps = ' tps ;
  144. ltps=ltps et (prog tps) ;
  145. EXEC  rv ;
  146. teta='EXTRAIRE'  (rv.inco.'D2') 'SCAL' pt1 ;
  147. mess ' Teta=' teta ;
  148. lteta= lteta et ('PROG' teta);
  149.  
  150.  
  151. YNM = rv.'INCO'.'D2NM';
  152. ZNM = rv.'INCO'.'Q1NM';
  153.  
  154. rv.'INCO'.'D2N2' = rv.'INCO'.'D2NM';
  155. rv.'INCO'.'D2NM' = rv.'INCO'.'D2';
  156. rv.'INCO'.'Q1N2' = rv.'INCO'.'Q1NM';
  157. rv.'INCO'.'Q1NM' = rv.'INCO'.'Q1';
  158.  
  159. 'FIN'  bclt ;
  160.  
  161. li4 = 'PROG' 'EXPO' 'A' -1. 'B' 0. ltps;
  162. evt2 = 'EVOL' 'ROUGE' manu ltps li4 ;
  163. evt= evol manu ltps lteta;
  164.  
  165.  
  166.  
  167. *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  168. *>>>>>>>>>>>>>>>>>>>>>>> résultats>>>>>>>>>>>>>>>>>>>>>>>>>>
  169. *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  170. *évolution de la vitesse extérieure et du gradient
  171. vit = 'EVOL' 'CHPO' UE mt;
  172. gvit = 'EVOL' 'ROUG' 'CHPO' DUE mt;
  173. *quantité de mouvement
  174. D2S = (rv.inco.'D2');
  175. RD2S = UE*D2S/NU;
  176. ev1s = 'EVOL' 'ROUGE' 'CHPO' RD2S mt;
  177. *différence entre (d-d1)
  178. Q1S = (rv.inco.'Q1');
  179. RQ1S = UE*Q1S/NU;
  180. ev2s = 'EVOL' 'ROUGE' 'CHPO' RQ1S mt;
  181. *Epaisseur de déplacement
  182. D1S = (rv.'INCO'.'D1');
  183. RD1S = UE*D1S/NU;
  184. ev3s = 'EVOL' 'ROUGE' 'CHPO' RD1S mt;
  185. *évolution du frottement à la paroi
  186. CfS = rv.'INCO'.'CF';
  187. ev4s = 'EVOL' 'ROUGE' 'CHPO' CfS mt;
  188. *évolution du facteur de forme H final
  189. HS = rv.'INCO'.'H';
  190. ev5s = 'EVOL' 'ROUGE' 'CHPO' HS mt;
  191. *évolution du facteur H*
  192. HHS = Q1S*(D2S**(-1));
  193. ev6s = 'EVOL' 'ROUGE' 'CHPO' HHS mt;
  194.  
  195.  
  196. *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  197. *solution théorique cf:SCHLICHTING >>>>>>>>>>>>>>>>>>>>>>>>
  198. *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  199. RX = X*UE*(1./NU);
  200. *Quantité de mvt
  201. *----Rteta donné par Nikuradze
  202. RD2N = (0.01738/1.304)*(RX**(0.861)); 
  203. ev1n =  'EVOL' 'BLEU' 'CHPO' RD2N mt;
  204. *----Rteta par la loi de puissance
  205. RD2P = (1./1.e-6)*(0.036*X*(RX**(-1./5)))*UE;
  206. ev1p = 'EVOL' 'VERT' 'CHPO' RD2P mt;
  207. *différence (d-d1) théorique
  208. Q1P=(0.37-(1.2857*0.036))*(X*(RX**(-1./5)));
  209. RQ1P=UE*Q1P/NU;
  210. ev2p='EVOL' 'VERT' 'CHPO' RQ1P mt;
  211. *Epaisseur de déplacement
  212. *----Rdelta1 par la Nikuradze
  213. RD1N = 1.304*RD2N;
  214. ev3n = 'EVOL' 'BLEU' 'CHPO' RD1N mt;
  215. *----Rdelta1 par la loi de puissance
  216. RD1P = (9./7.)*RD2P;
  217. ev3p = 'EVOL' 'VERT' 'CHPO' RD1P mt;
  218. *Coefficient de frottement
  219. *----Formule par Nikuradze
  220. CfN = 0.02296*(RX**(-0.139));
  221. ev4n = 'EVOL' 'BLEU' 'CHPO' CfN mt;
  222. *----Formule par loi de puissance
  223. CfP = 0.0576*(RX**(-1./5.));
  224. ev4p = 'EVOL' 'VERT' 'CHPO' CfP mt;
  225. *facteur de forme
  226. *----H par Nikuradze
  227. HN = 'KCHT' $mt scal sommet 1.304;
  228. ev5n = 'EVOL' 'BLEU' 'CHPO' HN mt;
  229. *----H par loi de puissance
  230. HP = 'KCHT' $mt scal sommet 1.2857;
  231. ev5p = 'EVOL' 'VERT' 'CHPO' HP mt;
  232. *coef H* théorique
  233. HHP='KCHT' $mt scal sommet 7.89;
  234. ev6p= 'EVOL' 'BLEU' 'CHPO' HHP mt;
  235.  
  236. *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  237. *Courbes>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  238. *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  239. Si GRAPH ;
  240.  
  241. 'DESSIN'  (evt) 'TITR' 'Convergence de la valeur au dernier noeud'
  242.            'TITY' 'Valeurs de Y' 'TITX' 'Temps';
  243.  
  244. 'DESSIN' (vit 'ET' gvit)
  245.      'TITRE' 'Vitesse et gradient de vitesse extérieure';
  246.  
  247. TAB1 = 'TABLE' ;
  248. TAB1.'TITRE' = 'TABLE' ;
  249. TAB1.'TITRE' . 1 = MOT ' Simulation';
  250. TAB1.'TITRE' . 2 = MOT ' Nikuradze';
  251. TAB1.'TITRE' . 3 = MOT ' Loi de puissance';
  252.  
  253. TAB2 = 'TABLE' ;
  254. TAB2.'TITRE' = 'TABLE' ;
  255. TAB2.'TITRE' . 1 = MOT ' Simulation';
  256. TAB2.'TITRE' . 2 = MOT ' Loi de puissance';
  257.  
  258. 'DESSIN' (ev1s 'ET' ev1n 'ET' ev1p)
  259.        'TITRE' 'Epaisseur de quantité de mouvement'
  260.        'TITY' 'Rteta'
  261.        'TITX' 'Abscisses'
  262.        'GRIL'
  263.        'LEGE' TAB1;       
  264.  
  265. 'DESSIN' (ev2s 'ET' ev2p) 'TITRE' 'différence entre delta et delta1' 
  266.        'TITY' 'Rdelta-Rdelta1'
  267.        'TITX' 'Abscisses'
  268.        'GRIL'
  269.        'LEGE' TAB2;              
  270.  
  271. 'DESSIN' (ev3s 'ET' ev3n 'ET' ev3p)
  272.        'TITRE' 'Epaisseur de déplacement'
  273.        'TITY' 'Rdelta1'
  274.        'TITX' 'Abscisses'
  275.        'GRIL'
  276.        'LEGE' TAB1;
  277.  
  278. 'DESSIN' (ev4s 'ET' ev4n 'ET' ev4p)
  279.            'TITRE' 'Coefficient de frottement à la paroi'
  280.            'TITX' 'Abscisses'
  281.            'TITY' 'Cf'
  282.            'GRIL'
  283.            'LEGE' TAB1;
  284.  
  285. 'DESSIN' (ev5s 'ET' ev5n 'ET' ev5p)
  286.         'TITRE' 'Evolution du facteur de forme'
  287.         'TITX' 'Abscisses'
  288.         'TITY' 'H=delta1/theta'
  289.         'LEGE' TAB1;
  290.  
  291. 'DESSIN' (ev6s 'ET' ev6p) 'TITRE' 'Evolution du facteur H*'
  292.         'TITX' 'Abscisses'
  293.         'TITY' 'H'
  294.         'LEGE' TAB1;
  295. Finsi;
  296.  
  297. *calcul de l'erreur au sens de la norme L2
  298. diag = 'DOMA' $mt 'XXDIAGSI';
  299. r = (RD2S-RD2N)**2;
  300. r2= (RD2S-RD2P)**2;
  301. EL2N = 'SOMT'(diag*r);
  302. EL2P = 'SOMT'(diag*r2);
  303. 'MESSAGE' 'Erreur sur teta au sens de la norme L2=';
  304. 'LISTE'  EL2N;
  305. 'LISTE'  EL2P;
  306. r3=(RQ1S-RQ1P)**2;
  307. EL3P= 'SOMT'(diag*r3);
  308. 'MESSAGE' 'Erreur sur (delta-delta1) au sens de la norme L2=';
  309. 'LISTE'  EL3P;
  310. *SAUVEGARDE
  311. *'OPTION' 'SAUVER' '/u2/ttmf2/guillaume/c20/turb4H_sauv';
  312. *SAUVER;
  313.  
  314. **************************************************************************
  315. **************************************************************************
  316. *Facteur de forme
  317. x = 'EXTRAIRE' ev5s 'ABSC';
  318. hs = 'EXTRAIRE' ev5s 'ORDO';
  319. hn = 'EXTRAIRE' ev5n 'ORDO';
  320. hp = 'EXTRAIRE' ev5p 'ORDO';
  321. nlist = 'DIME' x;
  322. *'OPTION' echo 0;
  323. *'OPTION' impr 'turb4H_H.data';
  324. 'REPETER' BLIST nlist;
  325.    xi = 'EXTRAIRE' x &BLIST;
  326.    hsi = 'EXTRAIRE' hs &BLIST;
  327.    hni = 'EXTRAIRE' hn &BLIST;
  328.    hpi = 'EXTRAIRE' hp &BLIST;
  329. *  'MESSAGE' ('CHAINE' xi ' ' hsi ' ' hni ' ' hpi);
  330. 'FIN' BLIST;
  331. *'OPTION' impr 'Bidon';
  332. *Epaisseur de quantité de mouvement
  333. x = 'EXTRAIRE' ev1s 'ABSC';
  334. d2s = 'EXTRAIRE' ev1s 'ORDO';
  335. d2n = 'EXTRAIRE' ev1n 'ORDO';
  336. d2p = 'EXTRAIRE' ev1p 'ORDO';
  337. nlist = 'DIME' x;
  338. *'OPTION' echo 0;
  339. *'OPTION' impr 'turb4H_D2.data';
  340. 'REPETER' BLIST nlist;
  341.    xi = 'EXTRAIRE' x &BLIST;
  342.    d2si = 'EXTRAIRE' d2s &BLIST;
  343.    d2ni = 'EXTRAIRE' d2n &BLIST;
  344.    d2pi = 'EXTRAIRE' d2p &BLIST;
  345. *  'MESSAGE' ('CHAINE' xi ' ' d2si ' ' d2ni ' ' d2pi);
  346. 'FIN' BLIST;
  347. *'OPTION' impr 'Bidon';
  348. *Epaisseur de déplacement
  349. x = 'EXTRAIRE' ev3s 'ABSC';
  350. d1s = 'EXTRAIRE' ev3s 'ORDO';
  351. d1n = 'EXTRAIRE' ev3n 'ORDO';
  352. d1p = 'EXTRAIRE' ev3p 'ORDO';
  353. nlist = 'DIME' x;
  354. *'OPTION' echo 0;
  355. *'OPTION' impr 'turb4H_D1.data';
  356. 'REPETER' BLIST nlist;
  357.    xi = 'EXTRAIRE' x &BLIST;
  358.    d1si = 'EXTRAIRE' d1s &BLIST;
  359.    d1ni = 'EXTRAIRE' d1n &BLIST;
  360.    d1pi = 'EXTRAIRE' d1p &BLIST;
  361. *  'MESSAGE' ('CHAINE' xi ' ' d1si ' ' d1ni ' ' d1pi);
  362. 'FIN' BLIST;
  363. *'OPTION' impr 'Bidon';
  364.  
  365. *Coeff de frottement
  366. x = 'EXTRAIRE' ev4s 'ABSC';
  367. cfs = 'EXTRAIRE' ev4s 'ORDO';
  368. cfn = 'EXTRAIRE' ev4n 'ORDO';
  369. cfp = 'EXTRAIRE' ev4p 'ORDO';
  370. nlist = 'DIME' x;
  371. *'OPTION' echo 0;
  372. *'OPTION' impr 'turb4H_CF.data';
  373. 'REPETER' BLIST nlist;
  374.    xi = 'EXTRAIRE' x &BLIST;
  375.    cfsi = 'EXTRAIRE' cfs &BLIST;
  376.    cfni = 'EXTRAIRE' cfn &BLIST;
  377.    cfpi = 'EXTRAIRE' cfp &BLIST;
  378. *  'MESSAGE' ('CHAINE' xi ' ' cfsi ' ' cfni ' ' cfpi);
  379. 'FIN' BLIST;
  380. *'OPTION' impr 'Bidon';
  381.  
  382. 'FIN' ;
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  

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