* fichier : cyltest6.dgibi ************************************************************************ * Section : Fluides Non Stationnaire ************************************************************************ 'OPTI' 'ECHO' 0 ; ************************************************************************ * * * Ecoulement autour d'un cylindre circulaire * * Résolution des équations de Navier Stokes laminaires instationnaires * * * * Pompé sur cyltest.dgibi * * Ajouts Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF) * * mél : gounand@semt2.smts.cea.fr * * en date du 20/12/2007 * * teste une méthode couplée et une méthode de projection * * algébrique incrémentale : vérification du bilan des forces * * locales et de la force s'exercant sur le cylindre * * SG 08/10/2024 : passage aux quadratiques pour la discrétisation * * du cercle en TRIA RAFT * * Ajout double projection * * SG 31/01/2025 : Modifs tolerance tests car MacOS/Arm genere un * * maillage un peu different (TRIA/RAFT) * * Test simple projection aussi * ************************************************************************ 'SAUTER' 2 'LIGNE' ; 'MESSAGE' ' Execution de cyltest6.dgibi' ; 'SAUTER' 2 'LIGNE' ; * interact = FAUX ; graph = FAUX ; * ***************************************** * * * Procedure tracvit * * * ***************************************** 'DEBPROC' TRACVIT; 'ARGUMENT' RVX*'TABLE'; rv = rvx . 'EQEX' ; umax = 2. ; CHUN = 'VECTEUR' (rv.inco.un) (1. '/' (1. * umax)) UX UY JAUN; 'TRACER' CHUN mt cmt 'TITRE' 'Champ de vitesse' 'NCLK' ; as2 ama1 = 'KOPS' 'MATRIK'; 'RESPRO' as2 ama1; 'FINP'; ***************************************** * * * Proc Bloc de pression nul * * * ***************************************** 'DEBPROC' BLOCPNUL ; 'ARGUMENT' rvx*'TABLE' ; RV = RVX . 'EQEX'; * 'SI' ('NON' ('EXISTE' rv 'MATPNULL')) ; * * Creation d'un bloc de pression nulle * 'MESSAGE' ('CHAINE' 'Calcul du bloc de pression') ; $mt = rvx . 'DOMZ' ; rt = 'EQEX' 'OPTI' 'EF' 'IMPL' 'CENTREE' 'CENTREP1' 'ZONE' $mt 'OPER' 'KMAB' 0. 'INCO' 'UN' 'PN' ; Rt.'INCO' = RV . 'INCO' ; smbb matb = kmab rt . '1KMAB' ; chpv = 'KCHT' $mt 'VECT' 'SOMMET' 'COMP' '1UN' '2UN' (0. 0.) ; matpres = 'KOPS' 'CMCT' matb chpv matb ; rv . 'MATPNULL' = matpres ; 'SINON' ; matpres =rv . 'MATPNULL' ; 'FINSI' ; smbvide matvide = 'KOPS' 'MATRIK' ; 'FINPROC' smbvide matpres ; * * Options globales * 'OPTION' 'DIME' 2 'ELEM' 'QUA8' 'ECHO' 0 ; 'SI' interact ; 'OPTION' 'TRAC' 'X' ; 'SINON' ; 'OPTION' 'TRAC' 'PSC'; 'FINSI' ; 'OPTION' isov suli; disv = 'QUAF' ; DISP = 'CENTREP1'; * iraff : raffinement du maillage * mproj : 0 méthode couplée * 1 méthode simple projection * 2 méthode double projection * miter : méthode itérative ou directe (FAUX ou VRAI) iraff = 1 ; impkres = 0 ; * lok = VRAI ; * 'REPETER' imet 10 ; 'SAUTER' 1 'LIGN' ; iimet = &imet ; lprec = faux ; 'SI' ('EGA' iimet 1) ; mproj = 0 ; miter = FAUX ; 'FINSI' ; 'SI' ('EGA' iimet 2) ; mproj = 0 ; miter = VRAI ; 'FINSI' ; 'SI' ('EGA' iimet 3) ; mproj = 1; lprec = FAUX ; miter = FAUX ; 'FINSI' ; 'SI' ('EGA' iimet 4) ; mproj = 1 ; lprec = FAUX ; miter = VRAI ; 'FINSI' ; 'SI' ('EGA' iimet 5) ; mproj = 1 ; lprec = VRAI ; miter = FAUX ; 'FINSI' ; 'SI' ('EGA' iimet 6) ; mproj = 1 ; lprec = VRAI ; miter = VRAI ; 'FINSI' ; 'SI' ('EGA' iimet 7) ; mproj = 2 ; lprec = FAUX ; miter = FAUX ; 'FINSI' ; 'SI' ('EGA' iimet 8) ; mproj = 2 ; lprec = FAUX ; miter = VRAI ; 'FINSI' ; 'SI' ('EGA' iimet 9) ; mproj = 2 ; lprec = VRAI ; miter = FAUX ; 'FINSI' ; 'SI' ('EGA' iimet 10) ; mproj = 2 ; lprec = VRAI ; miter = VRAI ; 'FINSI' ; 'SI' ('EGA' mproj 0) ; chmet = 'CHAI' 'Methode directe' ; 'FINS' ; 'SI' ('EGA' mproj 1) ; chmet = 'CHAI' 'Methode simple projection' ; 'FINS' ; 'SI' ('EGA' mproj 2) ; chmet = 'CHAI' 'Methode double projection' ; 'FINS' ; 'SI' lprec ; chmet = 'CHAI' chmet ' preconditionnee' ; 'FINSI' ; 'MESS' chmet ; 'SI' miter ; 'MESSAGE' 'Solveur iteratif' ; 'SINON' ; 'MESSAGE' 'Solveur direct' ; 'FINSI' ; 'SAUTER' 1 'LIGN' ; *============================ * Construction du maillage: *============================ * Construction des points: nc = 4 ; nx = 11 ; ny = 6 ; nc = '*' nc iraff ; nx = '*' nx iraff ; ny = '*' ny iraff ; * p0 = 0. 0. ; p1 = 0.5 0. ; p2 = 0. 0.5 ; p3 = -0.5 0. ; p4 = 0. -0.5 ; * c1 = 'CERCLE' nc p1 p0 p2 ; c2 = 'CERCLE' nc p2 p0 p3 ; c3 = 'CERCLE' nc p3 p0 p4 ; c4 = 'CERCLE' nc p4 p0 p1 ; * ligc = c1 'ET' c2 'ET' c3 'ET' c4 ; ligc = 'INVERSE' ligc ; * * pA = -6. -6. ; pB = 15. -6. ; pC = 15. 6. ; pD = -6. 6. ; bas = pA 'DROIT' nx pB ; dro = pB 'DROIT' ny pC ; hau = pC 'DROIT' nx pD ; gau = pD 'DROIT' ny pA ; * lige = bas 'ET' dro 'ET' hau 'ET' gau ; * cntt = ligc 'ET' lige ; * mt = tria cntt ; xmt ymt = 'COORDONNEE' mt ; rmt = '**' ('+' ('**' xmt 2) ('**' ymt 2)) 0.5D0 ; *rmt = '*' rmt 3. ; rmt = '/' rmt 3. ; rmt = '/' rmt iraff ; mtr = 'RAFT' rmt mt ; 'SI' graph ; 'TRACER' mtr 'TITR' ('CHAINE' 'nbl=' ('NBEL' mtr)) ; 'FINSI' ; *'OPTION' 'DONN' 5 ; cercle = ligc ; * Création du domaine total: mt = mtr; cmt = 'CONTOUR' mt; *============================================== * Changement des éléments du maillage en QUAF: *============================================== _mt = 'CHANGER' mt 'QUAF'; _cmt = 'CHANGER' cmt 'QUAF'; _cercle = 'CHANGER' cercle 'QUAF'; _entree = 'CHANGER' gau 'QUAF'; _sortie = 'CHANGER' dro 'QUAF'; _hau = 'CHANGER' hau 'QUAF'; _bas = 'CHANGER' bas 'QUAF'; 'ELIM' 1.E-5 (_mt et _cercle et _entree 'ET' _sortie 'ET' _cmt 'ET' _hau 'ET' _bas); NN = 'NBEL' mt; 'MESSAGE' 'Nombre d elements du maillage :' ' ' NN; *======================================= * Formulation du domaine Navier Stokes: *======================================= $mt = 'MODE' _mt 'NAVIER_STOKES' disv ; $cmt = 'MODE' _cmt 'NAVIER_STOKES' disv ; $cercle = 'MODE' _cercle 'NAVIER_STOKES' disv ; $entree = 'MODE' _entree 'NAVIER_STOKES' disv ; $sortie = 'MODE' _sortie 'NAVIER_STOKES' disv ; $hau = 'MODE' _hau 'NAVIER_STOKES' disv ; $bas = 'MODE' _bas 'NAVIER_STOKES' disv ; * mt = 'DOMA' $mt 'MAILLAGE' ; cmt = 'DOMA' $cmt 'MAILLAGE' ; cercle = 'DOMA' $cercle 'MAILLAGE' ; entree = 'DOMA' $entree 'MAILLAGE' ; sortie = 'DOMA' $sortie 'MAILLAGE' ; hau = 'DOMA' $hau 'MAILLAGE' ; bas = 'DOMA' $bas 'MAILLAGE' ; *==================================== * Création des tables de résolution: *==================================== dt = .5D0 ; *rey = 30. ; rey = 100.; nu = 1./rey; arelax = 1.; * niter = 6 ; nitime = 3 ; * RV = 'EQEX' 'OMEGA' 1.0 'NITER' niter 'ITMA' 1 'FIDT' 1 ; * RV = 'EQEX' RV 'OPTI' 'EF' 'CENTREE' 'IMPL' * 'ZONE' $mt 'OPER' 'NS' 'INV_RE' 'INCO' 'UN' 'ZONE' $mt 'OPER' 'KONV' 1. 'UN' 'INV_RE' 'INCO' 'UN' 'OPTI' 'EF' 'CENTREE' 'IMPL' 'ZONE' $mt 'OPER' 'LAPN' 'INV_RE' 'INCO' 'UN' 'OPTI' 'EF' 'CENTREE' 'IMPL' 'ZONE' $mt 'OPER' 'DFDT' 1. 'UNM' 'DT' 'INCO' 'UN' 'OPTI' 'EF' 'IMPL' 'CENTREE' DISP 'ZONE' $mt 'OPER' 'KBBT' 1. 'INCO' 'UN' 'PN' ; 'SI' graph ; rv = 'EQEX' rv 'ZONE' $mt 'OPER' TRACVIT ; 'FINSI' ; * Implantation des conditions aux limites: RV = 'EQEX' RV 'CLIM' 'UN' 'UIMP' cercle 0. 'UN' 'VIMP' cercle 0. 'UN' 'UIMP' entree 1. 'UN' 'VIMP' entree 0. 'UN' 'VIMP' hau 0. 'UN' 'VIMP' bas 0. ; * * Choix de la méthode de résolution (couplée, projection) * et des solveurs (direct, itératif) * 'SI' ('EGA' mproj 0) ; * Solveur pour le système vitesse-pression rvm = rv . 'METHINV' ; rvm . 'SCALING' = 1 ; rvm . 'IMPINV' = impkres ; 'SI' ('NON' miter) ; rvm . 'TYPINV' = 1 ; 'SINON' ; * On rajoute le bloc de pression nul rv = 'EQEX' rv 'ZONE' $mt 'OPER' BLOCPNUL ; * rvm . 'TYPINV' = 4 ; rvm . 'LBCG' = 2 ; rvm . 'PRECOND' = 5 ; rvm . 'ILUTLFIL' = 3.D0 ; 'FINSI' ; 'SINON' ; rv . 'GPROJ' = 'TABLE' ; rv . 'GPROJ' . 'NOMVIT' = 'UN' ; rv . 'GPROJ' . 'NOMPRES' = 'PN' ; 'SI' ('NON' lprec) ; rv . 'GPROJ' . 'NOPREC' = VRAI ; 'FINSI' ; rv . 'GPROJ' . 'dblproj' = ('EGA' mproj 2) ; * * Solveur pour les composantes de la vitesse rvm = rv . 'METHINV' ; rvm . 'SCALING' = 1 ; rvm . 'IMPINV' = impkres ; 'SI' ('NON' miter) ; rvm . 'TYPINV' = 1 ; 'SINON' ; rvm . 'TYPINV' = 4 ; rvm . 'LBCG' = 2 ; rvm . 'PRECOND' = 3 ; 'FINSI' ; * Solveur pour l'équation de pression rv . 'GPROJ'. 'METHINV' = 'TABLE' 'METHINV' ; rvgm = rv . 'GPROJ' . 'METHINV' ; rvgm . 'SCALING' = 1 ; rvgm . 'IMPINV' = impkres ; 'SI' ('NON' miter) ; rvgm . 'TYPINV' = 1 ; 'SINON' ; rvgm . 'TYPINV' = 4 ; rvgm . 'LBCG' = 2 ; rvgm . 'PRECOND' = 3 ; * rvgm . 'PRECOND' = 5 ; * rvgm . 'ILUTLFIL' = 1.D0 ; 'FINSI' ; 'FINSI' ; * Création de la table des inconnues et initialisation: RV.INCO = 'TABLE' INCO; * RV.INCO.'UN' = 'KCHT' $mt 'VECT' 'SOMMET' (1.E-3 1.E-3); RV.INCO.'UNM' = 'KCHT' $mt 'VECT' 'SOMMET' (1.E-3 1.E-3); rv . 'INCO' . 'G' = (0. -9.8) ; RV.INCO.'DT' = dt ; RV.INCO.'INV_RE'= nu; RV.INCO.'PN' = 'KCHT' $mt 'SCAL' DISP 0.; * * Table avec les résultats en vitesse * res = 'TABLE' ; itres = 0 ; res . itres = 'COPIER' (RV.INCO.'UN') ; *************************************************** * * Boucle en temps faite main * 'TEMPS' 'ZERO' ; 'REPETER' iitime nitime ; rv . 'INCO' . 'UNM' = 'COPIER' (res . itres) ; EXEC RV; itres = '+' itres 1 ; res . itres = 'COPIER' (RV.INCO.'UN') ; 'FIN' iitime ; TABTPS = TEMP 'NOEC'; tcpu = TABTPS.'TEMPS_CPU'.'INITIAL' ; * *************************************************** lastime = '-' ('DIME' res) 1 ; vit = res . lastime ; vitm1 = res . ('-' lastime 1) ; umax = 2. ; * * * Créer les matrices * bidon matkon = 'KONV' (rv . '1KONV') ; bidon matlap = 'LAPN' (rv . '2LAPN') ; chdfdt matdfdt = 'DFDT' (rv . '3DFDT') ; * rt = 'EQEX' 'ZONE' $mt 'OPTI' 'EF' 'IMPL' 'CENTREE' DISP 'OPER' 'KMBT' -1. 'INCO' 'PN' 'UN' ; rt.'INCO' = rv . 'INCO' ; bidon matgrp = 'KMBT' (rt . '1KMBT') ; * pre = rv . 'INCO' . 'PN' ; nomavv = 'MOTS' 'UX' 'UY' ; nomapv = 'MOTS' '1UN' '2UN' ; nomavp = 'MOTS' 'SCAL' ; nomapp = 'MOTS' 'PN' ; vit = 'NOMC' nomavv nomapv vit ; vitm1 = 'NOMC' nomavv nomapv vitm1 ; pre = 'NOMC' nomavp nomapp pre ; forlap = 'KOPS' matlap '*' vit ; forkon = 'KOPS' matkon '*' vit ; forgrp = 'KOPS' matgrp '*' pre ; fordfdt = 'KOPS' matdfdt '*' ('-' vit vitm1) ; * forlap = ('NOMC' nomapv nomavv forlap) ; forkon = ('NOMC' nomapv nomavv forkon) ; forgrp = ('NOMC' nomapv nomavv forgrp) ; fordfdt = ('NOMC' nomapv nomavv fordfdt) ; * fort = forlap '+' forkon '+' forgrp '+' fordfdt ; * fmax = 'MAXIMUM' fort 'ABS' ; echf = '/' 0.5 fmax ; * vlap = 'VECTEUR' forlap echf UX UY VERT ; vkon = 'VECTEUR' forkon echf UX UY ROUG ; vgrp = 'VECTEUR' forgrp echf UX UY JAUN ; vdfdt = 'VECTEUR' fordfdt echf UX UY TURQ ; vf = 'VECTEUR' fort echf UX UY BLAN ; vtot = vlap 'ET' vkon 'ET' vgrp 'ET' vdfdt 'ET' vf ; * 'SI' graph ; tit = 'CHAINE' 'vert=lap;roug=kon;jaun=grp;turq=dfdt;blan=reac;fmax=' fmax ; 'TRACER' vtot _mt ('CONTOUR' _mt) 'TITRE' tit ; 'FINSI' ; * pmt = 'CHANGER' mt 'POI1' ; pcmt = 'CHANGER' cmt 'POI1' ; pint = 'DIFF' pmt pcmt ; fortint = 'REDU' fort pint ; fmaxint = 'MAXIMUM' fortint 'ABS' ; 'MESSAGE' ('CHAINE' 'fmax=' fmax) ; 'MESSAGE' ('CHAINE' 'fmaxint=' fmaxint) ; errbf = '*' ('/' fmaxint fmax) 100. ; 'MESSAGE' ('CHAINE' 'erreur bilan forces =' errbf ' %') ; forcer = 'REDU' fort cercle ; rfor = 'RESULT' forcer ; 'MESSAGE' 'Resultante des forces sur le cylindre' ; rfx = 'MAXIMUM' ('EXCO' 'UX' rfor) ; rfy = 'MAXIMUM' ('EXCO' 'UY' rfor) ; 'MESSAGE' ('CHAINE' ' en x =' rfx) ; 'MESSAGE' ('CHAINE' ' en y =' rfy) ; 'MESSAGE' ('CHAINE' 'tcpu=' tcpu) ; * * Tests * Erreur sur le bilan des forces au 08/10/2024 * 2.e-3 % pour la méthode directe * 3.5 % pour la méthode projection algébrique * 5.3e-2 % pour la méthode double projection algébrique * Valeur de la force de trainée (fx) (Linux x86) * -7.261-1 pour la méthode directe * -7.029e-1 pour la méthode projection algébrique * -7.262e-1 pour la méthode double projection algébrique * Valeur de la force de trainée (fx) (MacOS ARM) * -7.236-1 pour la méthode directe * 'SI' ('EGA' mproj 0) ; vrbf = 3.D-3 ; 'FINS' ; 'SI' ('EGA' mproj 1) ; vrbf = 4. ; 'FINS' ; 'SI' ('EGA' mproj 2) ; vrbf = 1.D-1 ; 'FINS' ; * vrftd = -7.25e-1 ; * 'SI' ('EGA' mproj 0) ; tolfx = 0.03e-1 ; 'FINS' ; 'SI' ('EGA' mproj 1) ; tolfx = 0.4e-1 ; 'FINS' ; 'SI' ('EGA' mproj 2) ; tolfx = 0.03e-1 ; 'FINS' ; * tst = '<' errbf vrbf ; 'SI' ('NON' tst) ; cherr = 'CHAINE' '!!! Erreur bilan force > ' vrbf ' %' ; 'MESSAGE' cherr ; 'SINO' ; cherr = 'CHAINE' 'OK Erreur bilan force < ' vrbf ' %' ; 'MESSAGE' cherr ; 'FINSI' ; lok = lok 'ET' tst ; * tst = '<' ('ABS' ('-' rfx vrftd)) tolfx ; 'SI' ('NON' tst) ; cherr = 'CHAINE' '!!! Erreur force trainee' ' ' rfx ' <>' ' ' vrftd ' a' ' ' tolfx ' pres' ; 'MESSAGE' cherr ; 'SINO' ; cherr = 'CHAINE' 'OK Erreur force trainee' ' ' rfx ' ~=' ' ' vrftd ' a' ' ' tolfx ' pres' ; 'MESSAGE' cherr ; 'FINSI' ; lok = lok 'ET' tst ; 'FIN' imet ; 'SAUTER' 2 'LIGNE' ; 'SI' lok ; 'MESSAGE' 'Tout sest bien passe' ; 'SINON' ; 'MESSAGE' 'Il y a eu des erreurs' ; 'FINSI' ; 'SAUTER' 2 'LIGNE' ; * 'SI' interact ; 'OPTION' 'DONN' 5 'ECHO' 1 ; 'FINSI' ; 'SI' ('NON' lok) ; 'ERREUR' 5 ; 'FINSI' ; * 'FIN' du jeu de données. fin;