* fichier :  kp_test.dgibi
************************************************************************
* Section : Mecanique Flambage
************************************************************************
************************************************************
*
*   
*  flambage d'une poutre encastree-libre sous poids propre 
*  dans un champ de pression hydrostatique
*
*  le probleme est equivalent au flambage d'une poutre sous
*  poids propre, mais avec une masse volumique 
*  (RHOstructure - RHOliquide)
*
*
*                        __ surface libre
*      __________________\/______
*
*
*
*
*
*               
*               |
*               |
*               |
*               |
*               |
*               |
*     __________|________________
*
*
* 
************************************************************
OPTI ECHO 1 ;
graph = 'N';

l = 15. ;          comm 'longueur de la poutre';
nv = 30 ;          comm 'decoupage vertical';
nh = 2 ;           comm 'decoupage horizontal';
rg = 9.81*1000;    comm 'RHO*g liquide';
**********************************************************
opti dime 3 elem qua4 ;
p1 = 0 0 l ; p2 = 1 0 l ;p3 = 1 1 l ; p4 = 0 1 l ;
l1 = p1 d nh p2 ; l2 = p2 d nh p3 ;
l3 = p3 d nh p4 ; l4 = p4 d nh p1 ;
lt = l1 et l2 et l3 et l4 ;
surf1 = tran lt nv (0 0 (-1*l)) ;
surf2 = dalle l4  l3  l2  l1 plan ;
surft = surf1 et surf2 ;
p0 = poin surf2 proc (.5 .5 l)  ;
lsol = poin surf1 plan (0 0 0) (1 0 0) (0 1 0) .01 ;
surft = orie surft poin (.5 .5 0) ;
SI ( ega graph 'O');
trac surft cach ;
finsi;
coorz = coor surft 3 ; 

*
** on tourne la structure pour mieux tester les
** differents termes des matrices
depl surft tour 30. (.5 .5 0.) p0 ;
depl surft tour  60. (.5 .5 0.) (.5 1 0) ;
SI ( ega graph 'O');
trac surft cach;
finsi;
vecx = cos 30 ; vecz = sin 30 ;
vec1 = -1*(vecx 0. vecz);  comm 'sens de la pesanteur' ;
*
modl1 = mode surft mecanique coq4 ;
mat1 = mate modl1 youn 2e11 nu .0 rho 8000 epai .04 ;
rig1 = rigi modl1 mat1 ; 
blo1 = bloq lsol depl rota ;
blo2 = bloq p0 uy ;
rigt = rig1  et blo1 et blo2 ;

*
** poids propre
ch1 = manu chpo surf1 2 ux (-9.81*vecx) uz (-9.81*vecz);
modl11 = redu modl1 surf1 ; mat11 = redu mat1 modl11 ; 
mass1 = mass modl11 mat11 ;
f0 = mass1*ch1 ;
dep0 = resou rigt f0 ;
sig0 = sigm modl1 mat1 dep0 'LINE' ;
ksig0 = ksigm modl1 sig0 mat1 flam ;

*
** pression hydrostatique
ev1 = evol manu hauteur (prog 0. l) pression 
                (prog (rg*2*l) (rg*l)) ;

chp = ev1 * coorz ;
MOP = 'MODE' SURFT 'CHARGEMENT' 'PRESSION' 'COQ4' ;
MAP = 'MATE' MOP 'EPAI' 0.04 ;
*
PR  = 'CHAN' 'CHAM' MOP chp ;
MAP = 'PRES' MOP PR ;
f1  = 'BSIG' MOP MAP ;
dep1 = reso rigt f1 ;
sig1 = sigm modl1 mat1 dep1 'LINE' ;
ksig1 = ksigm modl1 sig1 mat1 flam ;

*
** relation de corps rigide pour les sections
** afin d'eviter les flambages locaux
repeter bou1 (nv - 1 ) ;
 ll3 = lsol plus (-1*&bou1*l/nv*vec1);
 si (&bou1 ega 1) ;
   ll2 = ll3 ;
   rig2 = rela cori depl rota ll3;
 sinon ;
  ll2 = ll2 et ll3;
  rig2 = rig2 et (rela cori depl rota ll3);
 finsi;
fin bou1 ;
rig2 = rig2 et (rela cori depl rota surf2);
elim ll2 surft .01 ;
rigt = rigt et rig2;

*
** on utilse les matrices de pression disymetriques
** (leur assemblage doit etre symetrique)
kp1 = kp mop chp flam asym ;
kp2 = kp mop rg vec1 flam asym ;
*kp1 = kp modl1 chp flam  ;
*kp2 = kp modl1 rg vec1 flam  ;
kpt = kp1 et kp2 ;

*
** calcul du flambage
mod1 = vibr proche (prog 10.4) rigt 
                    (kpt et (-1*ksig1) et (-1*ksig0))  ;

fr1 = mod1.modes. 1 .frequence ;
lamda = (2*pi*fr1)**2 * (SIGN fr1); 
mess 'lamda='lamda ;
dep1 = mod1.modes.1 .deformee_modale ;
def1 = defo surft dep1 4. ;
SI ( ega graph 'O');
trac def1 cach (10 25 10);
finsi;

*
** comparaison avec la solution analytique
*
s1 = 1.*1. ; i1 = .04/6 + (.04/2) ;
rhog = (8000*4*.04/s1 - 1000.)*9.81;
lamd_th = 7.837*i1*2e11/(l**3)/rhog/s1 ;
err1 = (lamda - lamd_th)/lamd_th * 100 ;
mess 'coef. multiplicateur de (RHOsolide - RHOliquide)';
mess '    C2000               THEORIE              ERREUR';
MESS 'lamda='lamda    'lamda='lamd_th      err1'%';
MESS 'Note de PV : je soupçonne la solution analytique de ne ';
MESS  'pas l''''etre ' 'tant' ' que' ' cela';
si ((abs err1) < 10.) ; erre 0; sinon; erre 5 ; finsi;

fin ;