* fichier : nonconforme.dgibi
*
* utilisation de maillage non conforme crée par raff
*
* mettre igraph='O' ; pour avoir les dessins
*
igraph = mot 'N';
opti elem qua4 mode plan defo dime 2;
dens 2.;
* mesh 10x6
pa= 0 0; pb= 10 0;pc= 20 0;
liab= pa droi pb;libc= pb droi pc;
su = (liab et libc) trans ( 0 12);
si ( ega igraph 'O'); trac su; finsi;
* definition of density
x y = coor su;
distance = ((x - 10 ) * ( x- 10) + ( y * y)) ** 0.5;
den = 0.3 + (0.18*distance);
si ( ega igraph 'O'); trac su den; finsi;
* new mesh
su2= raff su den;
hh = elem su2 SURE ;
sureal = su2 diff hh ;
* use of this mesh
* definition of model and caracteristic
mo= model su2 mecanique elastique isotrope ;
ma = matr mo YOUN 2.e5 NU 0.3 ;
* loads
psupe = su2 point droite ( 0 12) ( 10 12) 0.1;
lisupe = elem ( contour su2) appu stric psupe;
ff = pres ( redu mo sureal) massif lisupe -1.;
* displacements conditions
py0= point su2 droit pa pc 0.01;
liy0= elem ( contou su2) appu strict py0;
li2bc = liy0 elem compris pb pc;
cl1= bloqu li2bc UY;
cl2= bloq UX pb;
cltot= cl1 et cl2;
* compute elastic solution
ri = rigi mo ma;
displa = reso ( ri et cltot) ff;
stre = sigma displa mo ma;
vm = vmis stre mo ma;
si ( ega igraph 'O'); trac su2 vm mo ma; finsi;
* compute stress intensity factor
gt = table;
gt.'OBJECTIF' = MOT 'J';
lifis = liy0 elem compris pa pb;
gt.'LEVRE_SUPERIEURE' = lifis;
gt.'FRONT_FISSURE' = Pb;
gt.'CARACTERISTIQUES' = ma;
gt.'MODELE' = mo;
gt.'SOLUTION_RESO' = displa;
rea = reaction ri displa;
gt.'CHARGEMENTS_MECANIQUES'=rea;
naa = 5;
gmoy=0.;
repe no naa;
gt.'COUCHE' = &no;
G_THETA gt;
si ( &no ega 1) ; g2=prog gt.resultats;sinon;
g2 = g2 et ( prog gt.resultats); finsi;
gmoy=gmoy + gt.resultats;
fin no;
gmoy = gmoy / naa;
irate=0;
repe noo naa;
g1= extr g2 &noo;
si (abs ( ( g1 - gmoy) / gmoy) > 0.03); IRATE=1;FINSI;
fin noo;
xx = prog 1 pas 1 naa;
ev= evol manu 'nb of rows' xx 'G ' g2;
ttt=table;
ttt.1 = mot 'MARQ CROI';
tt2=table;
tt2.1= ' G ';
ttt.'TITRE'=tt2;
si ( ega igraph 'O'); dess ev lege ttt; finsi;
* use of the refined mesh for thermal computation
mmt = model su2 thermique isotrope;
mmm= mate mmt k 2.5 'C' 2.7;
cod = condu mmt mmm;
pc3= poin su2 droite pc ( pc plus ( 0 1));
pc1= poin su2 droit pa ( pa plus ( 0 1 ));
bl1 = bloqu T pc1; bl2 = bloq T pc3;ft=depi bl2 12.;
te= resou ( cod et bl1 et bl2 ) ft;
si ( ega igraph 'O'); trac te su2;finsi;
tpb = redu te pb;tpbv= maxi tpb;
si(abs ( ( tpbv - 6. ) / 6. ) > 0.001); irate=1;finsi;
*
* resultats
*
mess ' valeur moyenne de G ' gmoy;
mess ' valeurs trouvées ' ; list g2;
mess ' valeur de la temperature au point PB ' tpbv;
mess ' valeur attendue au point pb est 6.';
** trois points non conforme reliés par une relation 42 53
** et 95 au milieu
p42 = su2 point proc (16 6); t42= redu te p42; t42 = maxi t42;
p53 = su2 point proc (16 8); t53 = redu te p53 ; t53 = maxi t53;
p95 = su2 point proc (16 7); t95 = redu te p95 ; t95 = maxi t95;
tmo= t42 + t53 / 2.; eca = tmo - t95 ;
er_rela= abs (eca / t95);
si (er_rela > 1.e-8);
mess ' t42 t53 t95 ' t4é t53 t95 ; erreur 5;
finsi;
si (irate ega 1);
erreur 5;
finsi;
fin;
© Cast3M 2003 - All rights reserved.
Disclaimer