* fichier : ktest_lump_dkt.dgibi
*
*
* ktest pour tester l'opérateur LUMP
* elem DKT3
*
*
*
*
graph = faux;
ep1 = 1. ;
exc1 = 0.5 ;
a1 = 8 ** 0.5 ;
b1 = 18 ** 0.5 ;
ro1 = 5.d0 ;
'OPTI' 'ECHO' 0 ;
'OPTI' 'DIME' 3 'ELEM' 'SEG2' ;
P1 = -1. 1. 1. ;
P2 = 1. -1. 1. ;
P3 = 3. 3. 1. ;
* cote du triangle
*==== TRI3
m_el1 = 'MANU' 'TRI3' P1 P2 P3 ;
mo_el1 = 'MODE' m_el1 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'DKT';
ma_el1 = 'MATE' mo_el1 'YOUN' 7.d0 'NU' 0.1 'RHO' ro1 'EPAI' ep1
'EXCE' 0. ;
mas_el1 = 'LUMP' mo_el1 ma_el1 ;
* on teste les valeurs propres de la matrice
n1 = 'DIMN' mas_el1 ;
n2 = 'DIAG' mas_el1 ;
'MESS' 'Dimension du noyau' n1 ;
'MESS' 'Valeur propres négatives' n2 ;
err1 = ( n1 'NEG' 0 ) 'OU' ( n2 'NEG' 0 ) ;
ma_el1 = 'MATE' mo_el1 'YOUN' 7.d0 'NU' 0.1 'RHO' ro1 'EPAI' ep1
'EXCE' exc1 ;
mas_el1 = 'LUMP' mo_el1 ma_el1 ;
* on teste les valeurs propres de la matrice
n1 = 'DIMN' mas_el1 ;
n2 = 'DIAG' mas_el1 ;
'MESS' 'Dimension du noyau' n1 ;
'MESS' 'Valeur propres négatives' n2 ;
err11 = ( n1 'NEG' 0 ) 'OU' ( n2 'NEG' 0 ) ;
* test sur le calcul de la taille de l'élément
* calculé avec mesu et avec la massse
* cela revient à considerer la conservation de l'energie
* cinétique en translation
s_el1 = 'MESU' m_el1 ;
un1 = 'MANU' 'CHPO' m_el1 3 'UX' 1. 'UY' 1. 'UZ' 1.;
s2_el1 = ( 'XTMX' mas_el1 un1 ) '/' ep1 / ro1 / 3. ;
dif2 = 'ABS' (( s_el1 - s2_el1 ) / s2_el1 ) ;
*mess s2_el1 ;
err2 = dif2 '>EG' 1.d-6 ;
'MESS' 'Erreur relative' ( dif2 * 100 ) '%' ;
'SI' ( err1 'OU' err2 'OU' err11 ) ;
'MESS' 'problème de conservation de la masse' ;
'ERREUR' 5 ;
'FINSI' ;
*====================================================================
* Conservation de l'energie cinétique en rotation
*
* rotation de la coque autour de (p1 p2)
x1 = 'COOR' 1 m_el1 ;
x2 = 'COOR' 2 m_el1 ;
* mouvement suivant UZ
deu2 = ( 'MANU' 'CHPO' m_el1 6 'UX' 0. 'UY' 0. 'UZ' 0.
'RX' (2 ** -0.5 ) 'RY' (2 ** -0.5 * -.1 ) 'RZ' 0. )
'+' ( 'NOMC' ( x1 + x2 / ( 2 ** 0.5) ) 'UZ' ) ;
deu2 = ( 'MANU' 'CHPO' m_el1 6 'UX' 0. 'UY' 0. 'UZ' 0.
'RX' 0. 'RY' 0. 'RZ' 0. )
'+' ( 'NOMC' ( x1 + x2 / ( 2 ** 0.5) ) 'UZ' ) ;
mv_el1 = ( 'XTMX' mas_el1 deu2 ) ;
* energie cinétique de translation suivant z
* theorie a * b^3 / 12 * ep1 *rho
* a = 8^0.5 b = 8^0.5 pour l'energie de translation suivant UZ
mvtheo = ( a1 * ( b1 ** 3. ) / 12. * ep1 * ro1 ) ;
*
* energie cinétique de transaltion selon le plan de la coque
* on peut aussi considérer l'énergie de rotation des fibres
* c'est à dire le mouvement suivant UX et UY induit par la rotation
* on ne tient compte que de l'excentrement
mvtheo2 = mvtheo + ( ro1 * ep1 * ( exc1 ** 2. ) * a1 * b1 / 2.) ;
mvtheo3 = mvtheo2 + ( ro1 * ep1 * ( ep1 ** 2. / 12 )
* a1 * b1 / 2.) ;
'MESS' 'Qte de mvt :' mv_el1 ' Théorique' mvtheo
'Avec inertie' mvtheo2 ;
* ======= on augmente le nombre d'éléments
nelem = 100 ;
l1 = d p1 nelem p2 ;
l2 = d p2 nelem p3 ;
l3 = d p3 nelem p1 ;
su1 = 'SURF' 'PLAN' ( l1 et l2 et l3 ) ;
mo_el1 = 'MODE' su1 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'DKT';
ma_el1 = 'MATE' mo_el1 'YOUN' 7.d0 'NU' 0.1 'RHO' ro1 'EPAI' ep1
'EXCE' exc1 ;
mas_el1 = 'LUMP' mo_el1 ma_el1 ;
*mas_el1 = 'MASS' mo_el1 ma_el1 ;
* rotation de la coque autour de (p1 p2)
x1 = 'COOR' 1 su1 ;
x2 = 'COOR' 2 su1 ;
* mouvement suivant UZ
deu2 = ( 'MANU' 'CHPO' su1 6 'UX' 0. 'UY' 0. 'UZ' 0.
'RX' (2 ** -0.5) 'RY' (2 ** -0.5 * -.1 ) 'RZ' 0. )
'+' ( 'NOMC' ( x1 + x2 / ( 2 ** 0.5) ) 'UZ' ) ;
mv_el1 = ( 'XTMX' mas_el1 deu2 ) ;
'MESS' 'Qte de mvt :' mv_el1 ' Théorique' mvtheo
'Avec inertie' mvtheo2 ;
dif3 = mv_el1 - mvtheo2 / mvtheo2 ;
err3 = dif3 '>EG' 1.d-1 ;
'MESS' 'Erreur relative' ( dif3 * 100 ) '%' ;
'SI' ( err3 ) ;
'MESS' 'problème de conservation de la masse en rotation' ;
'ERREUR' 5 ;
'FINSI' ;
*====================================================================
*
* test sur les fréquences propres
ro1 = 5. ;
'OPTI' 'DIME' 3 'ELEM' 'SEG2' ;
P1 = -1. 1. 1. ;
P2 = 1. -1. 1. ;
p3 = 1. 1. 1. ;
a1 = 2. ** 0.5 * 2 ;
b1 = 2. ** 0.5 ;
ep1 = 0.01 ;
you1 = 7. ;
nelem = 9 ;
l1 = d p1 nelem p2 ;
l2 = d p2 nelem p3 ;
l3 = d p3 nelem p1 ;
su1 = 'SURF' 'PLAN' ( l1 et l2 et l3 ) ;
mo_el1 = 'MODE' su1 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'DKT';
ma_el1 = 'MATE' mo_el1 'YOUN' you1 'NU' 0.3 'RHO' ro1 'EPAI' ep1 ;
rig1 = 'RIGI' mo_el1 ma_el1 ;
cl1 = 'BLOQ' ( l1 et l2 et l3 ) DEPL ROTA ;
* theorie Dr Robert Blevins
* Formulas for natural frequency and mode shape
* Van Nostrand Rheinhold Company
* p280
fith1 = 1. / 2. / pi / b1 / b1 * ( ( ep1 ** 2. * you1 / 12. / ro1
/ ( 1. - 0.09 ) ) ** 0.5 ) ;
mas_el1 = 'LUMP' mo_el1 ma_el1 ;
tab1a = 'VIBR' 'INTERVALLE' ( 10. * fith1 ) ( 100. * fith1 )
( rig1 et cl1 ) mas_el1
'BASSE' 1 'TBAS' ;
fex1 = tab1a . 'MODES' . 1 . 'FREQUENCE' ;
def1a = defo su1 ( tab1a . modes . 1 . deformee_modale ) 0.5 bleu ;
* avec la masse consistante
mas_el1 = 'MASS' mo_el1 ma_el1 ;
tab1b = 'VIBR' 'INTERVALLE' (10. * fith1) ( 100. * fith1 )
( rig1 et cl1 ) mas_el1
'BASSE' 1 'TBAS' ;
fex1b = tab1b . 'MODES' . 1 . 'FREQUENCE' ;
def1b = defo su1 ( tab1b . modes . 1 . deformee_modale ) 0.5 jaune ;
'MESS' 'Calcul de la 1 ere frequence propre' ;
'MESS' 'Lump : ' ( fex1 / fith1 ) 'Mass :' ( fex1b / fith1 )
'theo ' 46,8 ;
dif5 = 'ABS' ( fex1 / fith1 - 46.8 / 46.8 ) ;
'MESS' 'Erreur relative ' ( dif5 * 100. ) '%' ;
err5 = dif5 '>EG' 0.05 ;
si graph ;
'TRAC' ( def1a et def1b ) ;
finsi;
'SI' ( err5 ) ;
'MESS' 'problème du calcul de la frequence propre' ;
'ERREUR' 5 ;
'FINSI' ;
'MESS' 'DKT test positif' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer