* fichier : vibr11.dgibi
*
************************************************************************
* Section : Mecanique Dynamique
************************************************************************
*
*                        VIBR11.dgibi
*
* Objectif  : Calcul des modes propres d'un tube mince isotrope 
*             axisymetrique encastre - encastre 
* Elements  : solide-coque SHB8 et coque mince DKT + coque epaisse COQ8
* Creation  : BP, 2017-02-20
* Reference :  Cas test SDLA 01/89, Société Française des Mécaniciens, 
*              Guide de validation des progiciels de calcul des 
*              structures, AFNOR ed., 1990.
*
************************************************************************
*
* dimension dans le plan (XOY)
* 
*              solide-coque
*      Y  
*  Y=L | . . . D __ B                        
*      |        |  |                  
*      |        |  |                  
*      |        |  |                  
*      |        |  |                  
*      |        |  |                    
*      |        |  |                  
*      |        |  |                  
*      | . . .  |__|   -----> X       
*      O       C    A                 
*             R=Ri R=Re
*
*               |  |
*               <-->
*                 t
*
* la coque est definie sur Rmoy=(Re+Ri)/2
*
************************************************************************

  GRAPH = FAUX ;
  
************************************************************************
* OPTIONS ET DONNEES
************************************************************************

* options 
  OPTI DIME 3 MODE TRID ELEM CUB8;
  OPTI TRAC PSC POTR HELVETICA_16 EPTR 15;

* geometrie 
  L = 0.3048;
  t = 0.254E-3;
  Re = 0.0762;
  Ri = Re - t;
  Rmoy = (Ri + Re) / 2.;

* nombre d elements selon Z et R
  ncirc = 2**8; nz = 20; nR = 1;
*   ncirc = 2**8; nz = 40; nR = 1;
*   ncirc = 2**9;  nz = 80;  nR = 1;

* materiau
  Ey1 = 2.06E11;
  nu1 = 0.3;
  rho1 = 7850.;

  
************************************************************************
* MAILLAGE
************************************************************************

  p0   = 0. 0. 0.;
  vaxe = 0. 1. 0.;
  vx = 1. 0. 0.;
  vz = 0. 0. 1. ;


***********  MAILLAGE COQUE MINCE  ***********

  p1 = Rmoy 0. 0.;
  p2 = Rmoy L  0.;
  d12 = p1 droi nZ p2;
  OPTI ELEM TRI3;
  mesh3 = ROTA d12 360. ncirc p0 vaxe;
  ELIM mesh3 (1.E-3 * t);
  OPTI ELEM CUB8;
  
* recup  
  y3 = COOR mesh3 2;
  py1_3  = POIN y3 EGINFE      1.E-6;
  py2_3  = POIN y3 EGSUPE (L - 1.E-6);
  
  
***********  MAILLAGE SOLIDE-COQUE  ************

* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
* ATTENTION : avec les SHB8, il faut definir le maillage de telle sorte
*             que l'on ait les 4 premiers noeuds du CUB8 = face 1
*             et les 4 derniers = face 2 de la coque de maniere a 
*             respecter l'orientation de l'epaisseur
* ==> seul le maillage "C" est valide
* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

* kmail = mot 'A';  
* kmail = mot 'B';  
kmail = mot 'C';  
    
*--- facon A de construire le maillage ---*
si (ega kmail 'A');
* maillage massif lineaire 2D
  pC = Ri 0. 0.;
  pA = Re 0. 0.;
  d1 = pC droi nR pA;
  d2 = d1 plus (0. L .0);
  mesh4 = d1 regl nZ d2;
* maillage 3D
  mesh8 = VOLU 'ROTA' mesh4 ncirc 360. p0 vaxe;
  ELIM mesh8 (1.E-3*t);
finsi;

*--- facon B de construire le maillage ---*
si (ega kmail 'B');
  pC = Ri 0. 0.;
  pA = Re 0. 0.;
  d1 = pC droi nR pA;
  mesh4 = rota d1 ncirc 360. p0 vaxe;
  ELIM mesh4 (1.E-3*t);
  mesh8 = VOLU mesh4 'TRAN' nZ (0. L .0);
finsi;

*--- facon C de construire le maillage ---*
si (ega kmail 'C');
  pA = Re 0. 0.;
  pB = Re L  0.;
  pC = Ri 0. 0.;
  pD = Ri L  0.;
  de = pA droi nz pB;
  di = pC droi nz pD;
  mesh4e = de ROTA ncirc 360. p0 vaxe;
  ELIM mesh4e (1.E-3*t);
  mesh4i = di ROTA ncirc 360. p0 vaxe;
  ELIM mesh4i (1.E-3*t);
  mesh8 = VOLU mesh4i nR mesh4e;
finsi;
  
* recup  
  y8 = COOR mesh8 2;
  py1  = POIN y8 EGINFE      1.E-6;
  py2  = POIN y8 EGSUPE (L - 1.E-6);
  
  
***********  MAILLAGE COQUE EPAISSE  ***********
* rem : on met 2 fois moins d'elements que pour les maillages lineaires

OPTI ELEM QUA8;
  d12q = p1 droi (nZ/2) p2;
  meshcoq8 = ROTA d12q 360. (ncirc/2) p0 vaxe;
  ELIM meshcoq8 (1.E-3 * t);
  
* recup  
  ycoq8 = COOR meshcoq8 2;
  py1_coq8  = POIN ycoq8 EGINFE      1.E-6;
  py2_coq8  = POIN ycoq8 EGSUPE (L - 1.E-6);



***********  TRACE DES MAILLAGES  ***********
  si GRAPH;
    TRAC (0. 10. -30.) mesh3 'CACH' 'TITRE' 'maillage pour DKT';
    TRAC (0. 10. -30.) mesh8 'CACH' 'TITRE' 'maillage pour SHB8';
    TRAC (0. 10. -30.) meshcoq8 'CACH' 'TITRE' 'maillage pour COQ8';
  finsi; 
  

************************************************************************
* MODELEs ET MATERIAUx
************************************************************************

* coque mince dkt
  mod3 = MODE mesh3 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'DKT';
  mat3 = MATE mod3  'YOUN' Ey1 'NU' nu1 'RHO' rho1 'EPAI' t;
  
* solide-coque shb8
  mod8 = MODE mesh8 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'SHB8';
  mat8 = MATE mod8 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  
* coq8
  modcoq8 = MODE meshcoq8 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'COQ8';
  matcoq8 = MATE modcoq8 'YOUN' Ey1 'NU' nu1 'RHO' rho1 'EPAI' t;
  

************************************************************************
* MATRICES
************************************************************************

* coque mince dkt
  cl3 = BLOQ 'DEPL' (py1_3 et py2_3);
  K3  = RIGI mod3 mat3;
  M3  = MASS mod3 mat3;
  Ktot3 = K3 et cl3;

* solide-coque shb8
  cl8 = BLOQ 'DEPL' (py1 et py2);
  K8  = RIGI mod8 mat8;
  M8  = MASS mod8 mat8;
  K8tot = K8 et cl8;

* coq8
  clcoq8 = BLOQ 'DEPL' (py1_coq8 et py2_coq8);
  Kcoq8  = RIGI modcoq8 matcoq8;
  Mcoq8  = MASS modcoq8 matcoq8;
  Kcoq8tot = Kcoq8 et clcoq8;
  
  
************************************************************************
* ANALYSE MODALE
************************************************************************
  SI GRAPH;
    mopost = mots 'TABL' 'DEFO';
  SINO;
    mopost = mots 'TABL' ;
  FINSI;

* coque mince dkt
  TMOD_3 = VIBR 'IRAM' 300. 6  Ktot3 M3 ;
  POSTVIBR TMOD_3 mopost;
  
* solide-coque shb8
  TMOD_8 = VIBR 'IRAM' 300. 6  K8tot M8 ;
  POSTVIBR TMOD_8 mopost;
  
* coq8
  TMOD_coq8 = VIBR 'IRAM' 300. 6  Kcoq8tot Mcoq8 ;
  POSTVIBR TMOD_coq8 mopost;
   

************************************************************************
* TEST DE BON FONCTIONNEMENT
************************************************************************
* REM :  On teste seulement les valeurs propres. On pourrait aussi 
*        verifier que les harmoniques de Fourier sont 6, 5 et 7.

* ref = calcul mené le 20/02/2017 par BP avec maillage tres fin 
*       et valeurs correlé en 2D Fourier et par la Reference.
 w1_ref = (   (prog 534.62  578.44  594.11)
            + (prog 533.84  577.41  593.85) ) / 2.;
            
* tolerance = 5% , 1%  et 0.5%         
wtol_3 = 0.05;
wtol_8 = 0.01;
wtol_coq8 = 0.005;

* boucles sur les modes (tous les 2 car arrivent par pair)
I= -1 ;
REPE BMOD 3; I = I + 2;
 w1_3   =  TMOD_3 . 'MODES' . I . 'FREQUENCE';
 w1_8   =  TMOD_8 . 'MODES' . I . 'FREQUENCE';
 w1_coq8   =  TMOD_coq8 . 'MODES' . I . 'FREQUENCE';
 w1 = extr w1_ref &BMOD;
 err_3 = ABS (w1_3 - w1) / w1;
 err_8 = ABS (w1_8 - w1) / w1;
 err_coq8 = ABS (w1_coq8 - w1) / w1;
 MESS 'MODE ' &BMOD  ': DKT ' w1_3 '(' (100.*err_3) '%) '
                     '- SHB8' w1_8 '(' (100.*err_8) '%) '
                     '- COQ8' w1_coq8 '(' (100.*err_coq8) '%) ';
 si ((err_3 >eg wtol_3) ou (err_8 >eg wtol_8) ou (err_coq8 >eg wtol_8));
   ERRE 5;
 finsi;
FIN BMOD  ;

ERRE 0;
  
fin ;