* fichier :  plas14.dgibi
************************************************************************
* Section : Mecanique Plastique
************************************************************************
*           Test Plas14.dgibi: Jeux de données        *
*           ---------------------------------         *
*                                                     *
graph='N';
saut page;
opti echo 1;
*------------------------------------------------------
*
*   EXEMPLE OF A CONCRETE SQUARE SECTION WITH 4 STEEL 
*                      20mm BARS 
*     CONCRETE DIMENSIONS: (0.25*0.25)m2 - ALL SECTION
*                          (0.20*0.20)m2 - CORE SECTION
*     TRANSVERSAL REINFORCEMENT: 10mm Bars / 7.0cm
*                     NOVEMBER 1993
*
*------------------------------------------------------
*    Quadrangular Elements
*------------------------------------------
*   opti dime 2 elem qua4 echo 1 ;
   opti dime 2 echo 1 ;
*------------------------------------------
*    Triangular Elements
*------------------------------------------
*  opti dime 2 elem tri3;
*------------------------------------------
*    DEFINITION OF THE STEEL GEOMETRY
*------------------------------------------
secste = 3.14159*((10.e-3)**2);
*
cars1 = ( 8.50e-2    8.50e-2);
cars2 = ( 8.50e-2   -8.50e-2);
cars3 = (-8.50e-2   -8.50e-2);
cars4 = (-8.50e-2    8.50e-2);
*
meshf = cars1 et cars2 et cars3 et cars4;
meshf = coul meshf bleu;
*------------------------------------------------------
*    Quadrangular Elements
*------------------------------------------
   opti dime 2 elem qua4 echo 1 ;
*------------------------------------------
*    DEFINITION OF THE CONCRETE GEOMETRY
*------------------------------------------
*        - UNCONFINED
*------------------------------------------
*        verd - vertical number of divisions
*        hord - horizontal number of divisions
*------------------------------------------
verd = 4;
hord = 4;
*
divu = 1;
*------------------------------------------
pc1 = -12.5d-2  12.5d-2;
pc2 = -12.5d-2  10.0d-2;
pc3 =   2.5d-2   0.0d-2;
carc1 = pc1 d (divu) pc2 tran (divu) pc3;
*
carc2 = carc1 plus ( 22.5d-2   0.0);
carc3 = carc2 plus ( 0.0  -22.5d-2);
carc4 = carc3 plus (-22.5d-2   0.0);
*
pc4 = -10.0d-2  12.5d-2;
pc5 = -10.0d-2  10.0d-2;
pc6 =  20.0d-2   0.0d-2;
carc5 = pc4 d (divu) pc5 tran (verd) pc6;
*
carc6 = carc5 plus ( 0.0  -22.5d-2);
*
pc7 = -12.5d-2  10.0d-2;
pc8 = -10.0d-2  10.0d-2;
pc9 =   0.0d-2 -20.0d-2;
carc7 = pc7 d (divu) pc8 tran (hord) pc9;
*
carc8 = carc7 plus ( 22.5d-2   0.0);
*
meshu = carc1 et carc2 et carc3 et carc4 et  carc5 et carc6 et carc7 et carc8;
meshu = coul meshu jaune;
*opti donn 5;
*------------------------------------------
*        - CONFINED
*------------------------------------------
pc10 = -10.0d-2  10.0d-2;
pc11 =  10.0d-2  10.0d-2;
pc12 =   0.0d-2 -20.0d-2;
meshc = pc10 d (verd) pc11 tran (hord) pc12;
meshc = coul meshc rouge;
*
elim (meshf et meshu et meshc) 0.001;
*
titre  'Section:blue=steel,yellow=unconfined concrete ,red=confined concrete';
*
si (ega graph 'Y');
trac (meshf et meshu et meshc);
finsi;
*
opti dime 3 elem seg2;
*
pp0 =   0.0  .0 .0;
pp1 =   1.0  .0 .0;
*
llb = pp0 d 1 pp1;
*
elim (llb et pp0 et pp1) 0.001;
*------------------------------------------
*    CARACTERIZATION OF THE STEEL AND CONCRETE MODELS
*------------------------------------------
*------------------------------------------
*    Quadrangular Elements
*------------------------------------------
   modf = MODE meshf mecanique elastique PLASTIQUE  ACIER_UNI pojs;
   modbs = MODE meshf mecanique elastique PLASTIQUE  ACIER_ANCRAGE pojs;
   modu = MODE meshu mecanique elastique PLASTIQUE  BETON_UNI quas;
   modc = MODE meshc mecanique elastique PLASTIQUE  BETON_UNI quas;
   modlab = MODE (meshc et meshu) mecanique elastique PLASTIQUE  UNILATERAL quas;
   modshea = MODE (meshc et meshu)  mecanique elastique  plastique cisail_nl quas;
   modshea2 = MODE (meshc et meshu)  mecanique elastique  plastique cisail_nl quas;
*------------------------------------------
*    Triangular Elements
*------------------------------------------
*  modf = MODE meshf mecanique elastique fibre_nl 
*         ferraille tris;
*  modu = MODE meshu mecanique elastique fibre_nl 
*         beton tris;
*  modc = MODE meshc mecanique elastique fibre_nl 
*         beton tris;
*------------------------------------------
*   Steel
*------------------------------------------
matf = MATE modf 'YOUN' 2.03e5 'NU'   0.30   'STSY' 440.0 'EPSU' .090   'STSU' 760.0  'EPSH' 0.030 'FALD' 4.375  'A6FA' 620.0 'CFAC'  0.5  'AFAC' 0.008 'ROFA' 20.0   'BFAC' 0.010 'A1FA'  18.5 'A2FA' 0.15 'RHO ' 7.8D-3 ;
carf = 'CARA' modf 'ALPY' 1.     'ALPZ' 2.  'SECT' secste;
*------------------------------------------
*   Bond slip for lap splices
*------------------------------------------
G12 =  0.1*(30000./(2.*(1+0.25)*0.05));
xs1t = 0.0006;
xs2t = 0.0020;
xs3t = 0.0060;
xt1t = 5.;
xt3t = 0.15*xt1t;
xalfa = 0.4;
*
mess G12;
mess (xt1t/xs1t);
*
matbs = MATE modbs 'YOUN' 2.03e5 'NU'   0.30   'STSY' 440.0 'EPSU' .090   'STSU' 760.0  'EPSH' 0.030 'FALD' 4.375  'A6FA' 620.0 'CFAC'  0.5  'AFAC' 0.008 'ROFA' 20.0   'BFAC' 0.010 'A1FA'  18.5 'A2FA' 0.15 'RHO ' 7.8D-3  'G12 ' G12  'S1T' xs1t 'S2T' xs2t  'S3T' xs3t 'T1T' xt1t 'T3T' xt3t  'ALFA' xalfa  'SECB' (pi*((0.020/2.)**2)) 'LANC' (10.*0.020) 'SECT' secste;
carbs = 'CARA' modbs 'ALPY' 0.     'ALPZ' 0.;

*------------------------------------------
*   Unconfined concrete
*------------------------------------------
matu = MATE modu 'YOUN' 0.30e5 'NU'   .20     'STFC' 30.0   'EZER' .002  'STFT' 3.0     'ALF1' .22687 'OME1' .32912  'ZETA' 100.0   'ST85' .0      'TRAF' 3.0 'FACL'  1. 'FAMX' 10. 'STPT' .0 'FAM1' 1.  'FAM2' 10. 'RHO ' 2.3D-3;
caru = 'CARA' modu 'ALPY' .66    'ALPZ' .00;
*------------------------------------------
*   Confined concrete
*------------------------------------------
*   Initial concrete Young modulus =
*     2 * STIFC / ( BETA * EZERO )
*------------------------------------------
matc = MATE modc 'YOUN' 0.2254e5 'NU' .25  'STFC' 30.0   'EZER' .002  'STFT' 3.0    'ALF1' .22687 'OME1' .32912 'ZETA' 0.0      'ST85' 6.0    'TRAF' 3.0 'FACL'  1. 'FAMX' 10. 'STPT' .0 'FAM1' 1. 'FAM2' 10. 'RHO ' 2.3D-3;
carc = 'CARA' modc 'ALPY' .66      'ALPZ' .00;
*
* Laborderie
*
matlab = MATE modlab 'YOUN' 0.297E+5 'NU'  .20  'YS1 ' 2.5E-4 'YS2 ' 1.5E-3 'A1  ' 5000.  'A2  '  10. 'B1  ' 1.5  'B2  ' 1.5 'BET1' 1.0  'BET2' -40. 'SIGF' 3.5  'RHO ' 7.8D-3;
carlab = 'CARA' modlab 'ALPY' .66      'ALPZ' .00;
*----------------------------------------------
* Linear shear behaviour
*-----------------------------------------------
EbetJoi = 22540.;
xnub = 0.25;
*
tultshea = (0.9*500.*pi*((10./2000.)**2))/(0.20*0.07);
*
GY = EbetJoi/(2.*(1.+XNUB));
*
UXX1=PROG  0. (0.98*tultshea/GY) (tultshea/GY) (2.*(tultshea/GY)) (10.*(tultshea/GY)) ;
SHEAR1=PROG  0. (0.99*tultshea) tultshea tultshea tultshea;
*
EP = (extr SHEAR1 2)/(extr UXX1 2);
*
DMAXP = 1. - ((extr SHEAR1 3)/(EP*(extr UXX1 3)));
DMAXN = DMAXP;
DELAP = (extr UXX1 2);
DELAN = DELAP;
E2F = (1.-DMAXP)*EP;
XNU = 0.;
*
XMONOP = PROG;
XMONON = PROG;
YMONOP = PROG;
YMONON = PROG;
*
*
J0 = 2;
  REPETER LAB2 ((DIME SHEAR1) - 2);
    J0 = J0 + 1;
    YY = (EXTR SHEAR1 J0);
    XX = ((EXTR UXX1 J0) - ((extr SHEAR1 J0)/E2F));
    XMONOP = INSE XMONOP (J0 - 2) (MAXI (PROG XX 0.));
    YMONOP = INSE YMONOP (J0 - 2) YY;
  FIN LAB2;
*
XMONON =  (XMONOP);
YMONON =  (YMONOP);
*
monop = evol manu XMONOP YMONOP;
monon = evol manu XMONON YMONON;
*
matshea = mate modshea 'YOUN' ebetjoi 'NU  ' XNUB 'DELP' delap 'DMAP' dmaxp  'DELN' delan 'DMAN' dmaxn  'BETA' 0.2 'ALFA' 0. 'TETA' 1. 'MONP' monop 'MONN' monon  'RHO ' 0. 'ALPY' 0. 'ALPZ' 1.;
*
modq = modf et modu et modc et modshea;
macq = matf et matu et matc et carf et caru et carc et matshea;
*
modq2 = modbs et modu et modc et modshea;
macq2 = matbs et matu et matc et carbs  et caru et carc  et matshea;
modq3 = modf et modlab et modshea;
macq3 = matf et matlab et carf et carlab et matshea;
*
*
*------------------------------------------
*    USE OF "MOMCUR" PROCEDURE FOR THE ANALYSIS OF THE
*               PLASTIC BEHAVIOUR SECCION
*------------------------------------------
*    CARACTERIZATION OF THE ACTION 
*   (CURVATURES ALONG OY AXIS AND CONSTANT AXIAL FORCE)
*------------------------------------------
eppl = 440.0/2.03e5;
*
cy = prog 0 pas .0005 .005 pas .005 .138;
*
ncur = dime cy;
cz = prog ncur *  .00;
fa = (prog ncur * -.25);
*------------------------------------------
*    RESOLUTION
*------------------------------------------
my mz ea moc1 = mocu cy cz fa modq macq (1.d-6*eppl) verif;
*opti donn 5;
si (ega graph 'Y');
nste = dime (moc1 . contraintes);
repete bouc nste;
toto = redu (modc et modu) (moc1 . contraintes . (&bouc - 1));
titre 'pas' &bouc '---' 'SMXX' '---'  'deformation normale =' (moc1 . deformations .(&bouc - 1));
trac (exco toto SMXX) (modc et modu) (matc et matu) (meshu et meshc) ;
toto = redu (modc et modu) (moc1 . variables_internes . (&bouc - 1));
titre 'pas' &bouc 'EPSO';
trac (exco toto EPSO) (modc et modu) (matc et matu) (meshu et meshc) ;

 fin bouc;
finsi;
my2 mz2 ea2 moc2 = mocu cy cz fa modq2 macq2 (1.d-6*eppl) verif;
si (ega graph 'Y');
nste = dime (moc2 . contraintes);
repete bouc nste;
toto = redu (modc et modu) (moc2 . contraintes . (&bouc - 1));
titre 'pas' &bouc 'SMXX';
trac (exco toto SMXX) (modc et modu) (matc et matu) (meshu et meshc) ;
 fin bouc;
finsi;
my3 mz3 ea3 = mocu cy cz fa modq3 macq3 (1.d-6*eppl);

*------------------------------------------
*    OUTPUT DIAGRAMS
*------------------------------------------
c1= evol rouge manu 'Curvature' cy 'Moment'  (my*1.d3);
c2= evol vert manu 'Curvature' cy 'Moment'  (my2*1.d3);
c3= evol bleu manu 'Curvature' cy 'Moment'  (my3*1.d3);
*------------------------------------------
*    TRILINEAR CURVE FOR A TAKEDA MODEL 
*          FOR THE SAME SECTION
*------------------------------------------
*abstak=prog 0. 2.03791E-03 1.85207E-02 1.38834E-01;
abstak=prog 0. 2.03791E-03 1.85207E-02 1.38E-01;
*ordtak=prog 0. 2.05353E+01 7.10923E+01 7.06633E+01;
ordtak=prog 0.  21.133       71.1       72.817;
albnl=evol vert manu 'Curvature' abstak 'Moment' ordtak;
*------------------------------------------
*    PLOT
*------------------------------------------
si (ega graph 'Y');
  tt = table;
  tt.1 = 'MARQ CARR';
  tt.2 = '';
*
  titre 'courbe mocu (blanc: beton uni et bleu: unilateral) , takeda (vert)  et avec lap splices (rouge)';
  dess (albnl et c1 et c2 et c3) tt;
*
finsi;
*------------------------------------------
*    ERREUR
*------------------------------------------
ordtak=ipol cy abstak ordtak;
errlis=ordtak - (my*1.d3);
errea=((ltl errlis errlis)**0.5) / (dime ordtak);
denom=((ltl ordtak ordtak)**0.5) / (dime ordtak);
errel=errea/denom;
mess 'erreur relative=' errel '(+-=3.5%)';
*
si (errel > 4.d-2); erre 5;
sinon;              erre 0;
finsi;
*------------------------------------------
*             TEST 3D
*    Modèle avec une rotule non lineaire 
*    et un element de poutre linéaire
*------------------------------------------
pp0 =   0.0  .0 .0;
pp1 =   0.10 0. 0.;
pp2 =   1.0  .0 .0;
*
llhinge = pp0 d 1 pp1;
llelast = pp1 d 2 pp2;
*
modhing = 'MODE'  llhinge mecanique elastique SECTION PLASTIQUE  SECTION TIMO;
mathing = MATE modhing MODS modq MATS macq 'VECT' (0. 1. 0.);
*
modhing2 = 'MODE'  llhinge mecanique elastique SECTION PLASTIQUE  SECTION TIMO;
mathing2 = MATE modhing2 MODS modq2 MATS macq2 'VECT' (0. 1. 0.);
*
*
*modelas = 'MODE'  llelast mecanique elastique SECTION PLASTIQUE
*                           SECTION TIMO;
*matelas = MATE modelas MODS modq MATS macq
*           'VECT' (0. 1. 0.);
XXINRZ = (0.25**4)/12.;
*
modelas = 'MODE'  llelast mecanique elastique POUT;
matelas = MATE modelas YOUN ((1./3.)*0.297E+5) NU .20 INRZ XXINRZ INRY XXINRZ TORS XXINRZ SECT (0.25*0.25) 'VECT' (0. 1. 0.) 'RHO ' 2.3D-3;
*
MODTOT = MODELAS et MODHING;
MATTOT = MATELAS et MATHING;
*
MODTOT2 = MODELAS et MODHING2;
MATTOT2 = MATELAS et MATHING2;
*
*---------------------------------------------
*  Check of the total mass 
*--------------------------------------------
MASTOT = MASS MODTOT MATTOT;
valmas = maxi (resu (MASTOT * (manu chpo  (llhinge et llelast) UX 9.81)));
valmasth =  1.41019E-03;
errel = (valmas - valmasth)/valmasth;
*
si (errel > 4.d-2); erre 5;
sinon;              erre 0;
finsi;
*
bl0 = BLOQ DEPL ROTA pp0;
bl2 = BLOQ UZ pp2;
dep2 = DEPI bl2 1.;
FV = FORC  ((-0.25) 0. 0.) pp2;
*
time = prog 0. 0.1 1.0;
tidep = prog 0. 0. 0.01;
tiforv = prog 0. 1. 1.;
timecalc = prog 0. 0.1 pas 0.05 1.;
*
evde = evol manu time tidep;
evfv =evol manu time tiforv;
chade = charg dimp dep2 evde;
chafv = charg fv evfv;
*
*------------------------------------------------------------
*   Linear shear
*------------------------------------------------------------
   TAB = TABLE ;
   TAB.'BLOCAGES_MECANIQUES' = BL0 et BL2;
   TAB.'MODELE' = MODTOT;
   TAB.'CHARGEMENT' = CHADE et CHAFV;
   TAB.'TEMPS_CALCULES' = timecalc;
   TAB.'CARACTERISTIQUES' = MATTOT;
   TAB.'MOVA'             = RIEN;
*
TMASAU=table;
tab . 'MES_SAUVEGARDES'=TMASAU;
TMASAU .'DEFTO'=VRAI;
TMASAU .'DEFIN'=VRAI;
   PASAPAS TAB ;
*
dtab1=index(tab.deplacements) ;
*
ndime=dime dtab1 ;
*
prdep = prog 0.;
prfor = prog 0.;
*
i=1 ;
*
REPETER BOU1 (ndime - 1);
*
i=i+1 ;
*
d=dtab1.i ;
*
dep0 = tab.deplacements.d ;
sig0 = tab.contraintes.d ;
var0 = tab.variables_internes.d ;
def0 = tab.deformations_inelastiques.d ;
*
prdep = prdep et (prog (extr dep0 pp2 UZ));
*
reabase = reac bl0 dep0;
*
prfor = prfor et (prog ((-1.)*(extr reabase pp0 FZ)));
*
FIN BOU1;
*
depfor = evol manu prdep prfor;
*
si (ega graph 'Y');
  titre 'courbe effort tranchant - déplacement';
  dess depfor;
finsi;
*
*
errl = abs ((maxi (extr depfor ordo) abs) - 7.72465E-02);
*
*
si (errl > 4.d-2); erre 5;
sinon;              erre 0;
finsi;
*
fin;
*