C CFLUI3D SOURCE OF166741 25/11/04 21:15:18 12349 C CFLU3D SOURCE PV090527 23/02/07 11:58:34 11592 c A.Sellier nov 2023 SUBROUTINE cflui3d(WRK52,WRK53,WRK54,MWRKXE,WR14, # nbnb,idimb,teta13d,teta23d,nvarib,nstrsb,ifourb,dtb,trefb) C WRK52,53,54 segments déclarés dans le common DECHE C XMAT(NCOMAT) = COMPOSANTES DE MATERIAU C IVAL(NCOMAT) = INDICE DES COMPOSANTES DE MATERIAU C NCOMAT = NOMBRE DE COMPOSANTES DE MATERIAU C XCAR(ICARA) = CARACTERISTIQUES C MFR = NUMERO DE LA FORMULATION DE L'ELEMENT FINI C = 1 MASSIF C = 3 COQUE MINCE ( COQ2 , COQ3 ET DKT ) C = 5 COQUE EPAISSE ( COQ6 , COQ8 ) C = 7 POUTRE C = 9 COQUE MINCE AVEC CISAILLEMENT TRANSVERSE ( COQ4 ET DST ) C = 11 LIQUIDE C = 13 TUYAU C = 15 LINESPRING C = 17 TUYAU FISSURE C = 19 RACCORD MASSIF C = 21 RACCORD COQUE C = 23 SURFACE LIBRE C = 25 MEMBRANE C = 27 UNIAXIALE C = 29 THERMIQUE C = 31 INCOMPRESSIBLES C = 33 POREUX C = 35 JOINT C = 37 HOMOGENEISE C = 39 TUYO C = 41 TUYAU ACOUSTIQUE PURE C = 43 RACCORD TUYAU FLUIDE c* DDAUX = MATRICE DE HOOKE ELASTIQUE c* NSTRS = NBRE DE COMPOSANTES DES DEFORMATIONS c* CMATE = NOM DU MATERIAU c* VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU c* VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES c* N2EL = NBRE D ELEMENTS DANS SEGMENT DE HOOKE c* N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE c* MFR = NUMERO DE LA FORMULATION c* IFOU = type de formulation c* IB = NUMERO DE L ELEMENT COURANT c* IGAU = NUMERO DU POINT COURANT c* EPAIST= EPAISSEUR c* NBPGAU= NBRE DE POINTS DE GAUSS c* MELE = NUMERO DE L ELEMENT FINI c* NPINT = NBRE DE POINTS D INTEGRATION c* NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES c* NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES c* SECT = SECTION c* LHOOK = TAILLE DE LA MATRICE DE HOOKE c* TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI = TABLEAUX UTILISES c* POUR LE CALCUL DE LA MATRICE DE HOOKE c* matrice de Hoocke DDHOOK(LHOOK,LHOOK) C----------------------------------------------------------------------- C VARIABLES PASSEES PAR LES COMMONS COPTIO , ECOU ET NECOU C C IFOUR INDICE DU TYPE DE PROBLEME C -3 DEFORMATIONS PLANES GENERALISEES C -2 CONTRAINTES PLANES C -1 DEFORMATIONS PLANES C 0 AXISYMETRIQUE C 1 SERIE DE FOURIER C 2 TRIDIMENSIONNEL C ITYP TYPE DE FORMULATION MECANIQUE C --------------- ATTENTION --------------- C IL EST ACTIF APRES L APPEL DE VISAVI C ----------------------------------------- C ITYP=1 CAS DES ELEMENTS MASSIFS C ITYP=2 CAS DES COQUES C ITYP=3 CAS DES MEMBRANES C ITYP=4 CAS DES CABLES ET DES BARRES C ITYP=5 CAS QUELCONQUE C ITYP=6 CAS DES CONTRAINTES PLANES C ITYP=7 CAS DES COQUES A NU=0. OU CONTRAINTES PLANES C ITYP=8 CAS DES MEMBRANES A NU=0. OU CONTRAINTES PLANES C ITYP=9 CAS DES COQUES EPAISSES C ITYP=10 CAS DES JOINTS C ITYP=11 CAS DES POUTRES C ITYP=12 CAS DES TUYAUX C ITYP=13 CAS DES COQUES AVEC CISAILLEMENT TRANSVERSE C C ISTEP flag utilise pour separer les etapes dans un calcul non local C ISTEP=0 -----> calcul local C ISTEP=1 -----> calcul non local etape 1 on calcule les seuils C ISTEP=2 -----> calcul non local etape 2 on continue le calcul C a partir des seuils moyennes C ISTEP=3 ----->sous iteration de Helmholtz C ISTEP est défini par DEFINI.ESO C C EPIN0 déformations inélastiques initiales C EPINF déformations inélastiques finales C epst0 déformation totale initiale C EPSTF déformation totale finale C C----------------------------------------------------------------------- C SORTIES C SIGF(NSTRS) = CONTRAINTES FINALES C VARF(NVARI) = VARIABLES INTERNES FINALES C DEFP = DEFORMATIONS PLASTIQUES C KERRE = 0 TOUT OK C----------------------------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) c ajout sellier 2022 -INC PPARAM c fin ajout sellier 2022 -INC CCOPTIO -INC DECHE c segment de coordonnees des noeuds de l element integer nbnb,insb,idimb SEGMENT MWRKXE REAL*8 XE(3,nbnb) ENDSEGMENT * sellier/millard segment pour les Helmholtz 04 04 2020 integer NBVIA SEGMENT WR14 INTEGER INLVIA(NBVIA) ENDSEGMENT c----------------------------------------------------------------------- c nombre de noeuds max et table locale des coordonnees des c noeuds de l element fini -INC HNBRNEU c --rappel du contenu de HNBRNEU --------------------------------- C c ******* tableau local des coordonnees des noeuds de lEF ******** C integer NBNMAX3D C c nbr de noeuds maxi par element C parameter (NBNMAX3D=20) C c tableau des coordonnees des noeuds C real*8 XE3D(3,NBNMAX3D) c ---------------------------------------------------------------- c ****************************************************************** c ******** taille du pseudo vecteur des contraintes *************** c tjrs 6 en raison de son utilisation dans as3d integer nstrs3d parameter (nstrs3d=6) real*8 sig03d(nstrs3d),sigf3d(nstrs3d) real*8 epst03d(nstrs3d),epstf3d(nstrs3d) real*8 depst3d(nstrs3d) real*8 dtb,tref3d,alpha c variable logique : elasticite initialement isotrope logical iso1 c nombre de parametres de tailles suivant la formulation integer ifou11,ntail3d c dimension de l espace de travail integer idimb3d c temperatures debut et fin de pas , moyenne, pas de temps, volume rgi real*8 teta13d,teta23d,temp3d,dt3d c variables de transfert de donnees ( a declarer suivant idvar4 et idvisc) integer ierr1,mfr11 c variable pour passer le numero de l etape non locale integer istep3d,nvarib,nstrsb,ifourb c----------------------------------------------------------------------- c ****************************************************************** c * dimension des tableaux parametres materiau * c ****************************************************************** c *****nombre de renforts et de parametres par renfort************** c NB_RENF c NB_PARA_PAR_RENF c NB_PARA_RENF c NB_VARI_PAR_RENF c NB_VARI_RENF c-INC HNBRREN c adapte pour mc3d parameter(NB_RENF=0) parameter(NB_PARA_PAR_RENF=0) parameter(NB_PARA_RENF=0) parameter(NB_VARI_PAR_RENF=0) parameter(NB_VARI_RENF=0) c ****************************************************************** c ********Nombre de parametres pour variables non locales*********** c NB_HELM: Nombre de variable d etat principales c NB_PARA_PAR_HELM: Nombre de parametres par variable d etat c NB_VARI_PAR_HELM: Nombre de variable internes par variables d etat c NB_PARA_HELM: nombre de parametre total resultant c NB_VARI_HELM: nombre de variables internes totales resultantes c INLVIA3D(NB_HELM):copie du pointeur sur les numeros des vari traitee par helmholtz c NBVIA3D: recuperation du nombre de Helmholtz actives -INC HNBRHEL C c ------- Rappel du contenur de HNBRHEL ---------------------------- C c Declaration commune pour les varibales de Helmholtz C c Nombre de variable détat principales C integer NB_HELM C c Nombre de parametres par variable d etat C integer NB_PARA_PAR_HELM C c Nombre de variable internes par avriables d etat C integer NB_VARI_PAR_HELM C c nombre de variables non locales ISO et UNI C parameter (NB_HELM=8) C c nombre de parametre par variable iso et par variable uni C parameter (NB_PARA_PAR_HELM=17) C c nombre de vari par vari ISO (cf idvar4) (99 maxi) C parameter (NB_VARI_PAR_HELM=10) C c calcul du nombre maximum de variables de helmholtz C integer NB_PARA_HELM,NB_VARI_HELM C c nombre total de parametres non locaux (cf idmatr.eso) C parameter (NB_PARA_HELM=NB_HELM*NB_PARA_PAR_HELM) C c nombre de vari pour le non local C parameter (NB_VARI_HELM=NB_HELM*NB_VARI_PAR_HELM) C c tableau des numeros des variables internes non locales actives par Helmholtz C integer INLVIA3D(NB_HELM) C c nombre de Helmholtz actives C integer NBVIA3D c ------------------------------------------------------------------ c ****************************************************************** c parametres commun aux modele mecaniques elastiques isotropes integer NBELAS3D,NB_DECAL,NMAT3D c en isotrope le nombre de parametre de base est 4 = YOUN NU ALP RHO parameter (NBELAS3D=4) c decalage entre les parametres obligatoires (le modele) c et les parametres facultatifs (les parametres pour Helmholtz) parameter (NB_DECAL=3) c nombre de parametre ajoute par idmatr en fin de tableau des facultative parameter (NB_FIN=0) c Nombre de parametres partie non lineaire de la loi de comportement c NB_PARA_FLUENDO3D c NB_VARI_FLUENDO3D c NB_PARA_SUPP_FLUENDO3D c NB_VARI_SUPP_FLUENDO3D c-INC HNBRF3D c adapte pour mc3d apartir de idvis pour les parametres et idvar4 pour les vari c pas de fibre pour l instant dans mc3d parameter(NB_PARA_SUPP_FLUENDO3D=0) parameter(NB_VARI_SUPP_FLUENDO3D=0) parameter(NB_PARA_FIBRE_U=0) c ---cf affichage dans idvisc au lancement du modele mc3d ---------- parameter(NB_PARA_FLUENDO3D=52) c ---cf affichage dans idvar4 au lancement du modele mc3d ---------- parameter(NB_VARI_FLUENDO3D=137) c ***calcul des zones memoires pour les types de modeles *********** integer NMAT0,NMAT1,NMAT2,NMAT3,NMAT4 c nombre de parametres propre au modele sans les renforts ni le non local parameter (NMAT0=NBELAS3D) parameter (NMAT1=NMAT0+NB_PARA_FLUENDO3D) parameter (NMAT2=NMAT1+NB_PARA_SUPP_FLUENDO3D) parameter (NMAT3=NMAT2+NB_PARA_RENF+NB_DECAL) parameter (NMAT3D=NMAT3+NB_PARA_HELM+NB_FIN) c ****************************************************************** c dimension maxi du tableau des parametres y compris non local real*8 XMAT3D(NMAT3D) c ****************************************************************** c c ****************************************************************** c * dimension des tableaux variables internes * c ****************************************************************** c **** nombre de variables internes loi de comportement non lineaire c nbr de variables pour la loi de comportement non lineaire integer NVAR0,NVAR1,NVAR2,NVAR3,NVARI3D parameter (NVAR0=0) parameter (NVAR1=NVAR0+NB_VARI_FLUENDO3D) parameter (NVAR2=NVAR1+NB_VARI_SUPP_FLUENDO3D) parameter (NVAR3=NVAR2+NB_VARI_RENF) parameter (NVARI3D=NVAR3+NB_VARI_HELM) c ****************************************************************** c **** nombre total de variables internes ************************** c tableau local des variables internes real*8 VAR03D(NVARI3D),VARF3D(NVARI3D) c ****************************************************************** c----------------------------------------------------------------------- c*********************************************************************** c remarque : parametres elastiques en formulation poreux massif (mfr=33) c cas isotrope c en debut de xmat : c 'YOUN' : module d'Young c 'NU ' : coefficient de poisson c 'RHO ' : masse volumique (a la fin en non poreux) c 'ALPH' : coefficient de dilatation thermique c à la fin de xmat : c 'MOB ' : module de Biot c 'COB ' : coefficient de Biot c 'PERM' : perméabilité intrinsèque c 'VISC' : viscosité dynamique du fluide c 'ALPM' : coefficient de couplage pression - température c attention L ORDRE N EST PAS LE MEME SI LE MATERIAU EST ORTHOTROPE c cf deche.inc pour recuperer les noms exacts des variables c transfert des variables dans les tableaux du modele c*********************************************************************** c*****debut du transfert des donnees vers les modeles (sellier 2022)**** c ** on verifie si on a du non local via Helmholtz **************** c le decalage est RHO ALPH VISQ avant les parmetres non locaux si c il y a des non locales, alors qu'il est de 1 si on est en local if(WR14.EQ.0) then NBVIA = 0 else NBVIA=INLVIA(/1) c print*,'NBVIA = ',NBVIA c do i=1,NBVIA c print*, 'I' ,i, 'INLVIA ' ,INLVIA(i) c end do endif c suivant la formulation if(mfr.ne.33) then c ****** formulation non poreux ********************************** c print*,'ds cflui3d istep=',istep if(istep.eq.0) then c ***** formulation locale*********************************** if ((NMAT3D-NB_PARA_HELM).ne.nmatt) then print*,'pb dimension de xmat cflui3d local' print*,'NBELAS3D',NBELAS3D print*,'NMAT0',NMAT0 print*,'NB_PARA_FLUENDO3D',NB_PARA_FLUENDO3D print*,'NMAT1',NMAT1 print*,'NB_PARA_SUPP_FLUENDO3D',NB_PARA_SUPP_FLUENDO3D print*,'NMAT2',NMAT2 print*,'NB_DECAL',NB_DECAL print*,'NB_PARA_RENF',NB_PARA_RENF print*,'NMAT3',NMAT3 print*,'NB_PARA_HELM',NB_PARA_HELM print*,'NB_FIN',NB_FIN print*,'NMAT3D-NB_PARA_HELM',NMAT3D-NB_PARA_HELM print*,'NMATT',NMATT do i=1,nmatt print*,'xmat(',i,')=',xmat(i) enddo read* kerre=1 return else do i=1,(NMAT3D-NB_PARA_HELM) XMAT3D(i)=xmat(i) enddo c on met a zero les Helmholtz do i=(NMAT3D-NB_PARA_HELM+1), nmatt XMAT3D(i)=0. enddo end if else c **** formulation non locale ******************************* c ces parametres facultatifs sont definies dans IDMATR if(NMAT3D.ne.nmatt) then print*,'pb dimension de xmat cflui3d local' print*,'NBELAS3D',NBELAS3D print*,'NMAT0',NMAT0 print*,'NB_PARA_FLUENDO3D',NB_PARA_FLUENDO3D print*,'NMAT1',NMAT1 print*,'NB_PARA_SUPP_FLUENDO3D',NB_PARA_SUPP_FLUENDO3D print*,'NMAT2',NMAT2 print*,'NB_DECAL',NB_DECAL print*,'NB_PARA_RENF',NB_PARA_RENF print*,'NMAT3',NMAT3 print*,'NB_PARA_HELM',NB_PARA_HELM print*,'NB_FIN',NB_FIN print*,'NMAT3D',NMAT3D print*,'NMATT',NMATT do i=1,nmatt print*,'xmat(',i,')=',xmat(i) enddo read* kerre=1 return else do i=1,nmatt XMAT3D(i)=xmat(i) enddo end if end if else c *****formulation poreux **************************************** print*,'on est en poreux dans cflui3d' print*,'Pb affectation des caracteristiques dans cas3d' do i=1,nmatt print*,'xmat(',i,')=',xmat(i) end do print*,'voir idmatr et deche.inc pour parametres en poreux' kerre=1 return end if c********** recuperation des variables internes debut de pas *********** if(NVARI3D.ne.nvarib) then print*,'Pb dimension de var0 dans cas3d' print*,'Verifier dimension table des variables internes' print*,'et CMPATBLTCR avec idvar4 et as3d' read* kerre=1 return end if c transfert vers le tableau de fluendo3d do i=1,NVARI3D VAR03D(i)=var0(i) end do c*********** recuperation des contraintes totales debut de pas ********* if(mfr.ne.33) then c on est pas en poreux if(nstrsb.lt.nstrs3d) then do i=1,nstrsb sig03d(i)=sig0(i) sigf3d(i)=sigf(i) c ATTENTION les depst 3-6 doivent etre des gamas depst3d(i)=depst(i) epst03d(i)=epst0(i) epstf3d(i)=epstf(i) end do do i=nstrsb+1,nstrs3d sig03d(i)=0.d0 sigf3d(i)=0.d0 c ATTENTION les depst 3-6 doivent etre des gamas depst3d(i)=0.d0 epst03d(i)=0.d0 epstf3d(i)=0.d0 end do else do i=1,nstrs3d sig03d(i)=sig0(i) sigf3d(i)=sigf(i) c ATTENTION les depst 3-6 doivent etre des gamas depst3d(i)=depst(i) epst03d(i)=epst0(i) epstf3d(i)=epstf(i) end do end if else print*,'on est en poreux' print*,'Pb affectation des contraintes dans cflui3d' print*,'a terminer' do i=1,nstrsb print*,'sig(',i,')=',sig0(i) end do kerre=1 return end if c************** autres parametres a renseigner *********************** c initialisation indicateur d erreur ierr1=0 c indicateur isostropie elastique et de resistance iso1=.true. c numero de la formulation (33 pour poreux) mfr11=mfr c type de formulation ifour11=ifourb c pas de temps dt3d=dtb c numero pour le traitement non local eventuel istep3d=istep c************* recuperation des coordonnees de noeuds ****************** c boucle sur les noeuds if(nbnb.le.NBNMAX3D) then NBNB3D=nbnb if(idimb.le.3) then idimb3d=idimb else print*,'pb table de coord des neouds ds cflui3d' kerre=1 return end if else print*,'Element avec + de noeuds que capacite de cflui3d' kerre=1 return end if c chargement des coordonnees des noeuds do insb=1,NBNB3D do j=1,idimb3d c print*,'xe(',j,insb,')=',xe(j,insb) XE3D(j,insb)=xe(j,insb) end do if (idimb.lt.3) then do j=idimb+1,3 XE3D(j,insb)=0.d0 end do end if end do c on double les neouds si le probleme est 2D if ((ifourb.eq.0).or.(ifourb.eq.-1)) then c on double les noeuds avec dimension 3 recuperee ds xmat(1er param perso pour mc3d) dimension3=XMAT3D(NBELAS3D+1) NBNB3D=2*nbnb do insb=nbnb+1,NBNB3D do j=1,2 c print*,'xe(',j,insb,')=',xe(j,insb) XE3D(j,insb)=XE3D(j,(insb-nbnb)) end do XE3D(3,insb)=dimension3 end do end if c ****** recuperation de la temperature de reference ************** tref3d=trefb c print*,'verif temperature de reference fluendo3d',trefb c * intialisation des variables de sortie a leur valeur initiale * c en cas de 1ere etape non locale if(istep.eq.1) then do i=1,nvari3d varf3d(i)=var03d(i) end do else do i=1,nvari3d varf3d(i)=0.d0 end do endif c****************** fin du transferts des donnees ********************** c*************recuperation des numeros des variables Helmholtz ********* c transfert dans un tableau de fluendo3d dimensionne au max des c facultative helmholtz de IDMATR (actuellement 10) if(NBVIA.GT.NB_HELM) then print*,'Pb de dimensionnement avec les variables de Helmholtz' print*,'dans cflui3d:',NBVIA,' > ',NB_HELM print*,'Verifier NB_HELM et nbr de non locales' print*,'declarees nombre_helmholtz.h (NBRHEL.INC)' else c transfert vers les variables locales if (NBVIA.gt.0) then NBVIA3D=NBVIA do i=1,NB_HELM if(i.le.NBVIA3D) then INLVIA3D(i)=INLVIA(i) else INLVIA3D(i)=0 end if end do else NBVIA3D=NBVIA do i=1,NB_HELM INLVIA3D(i)=0 end do end if end if c************* fin de recuperation des numeros des variables Helmholtz * c************* appel au modele fluendo3d******************************* call flui3d(XMAT3D,NMAT3D,sig03d,sigf3d,depst3d,nstrs3d, #VAR03D,VARF3D,NVARI3D,teta13d,teta23d,dt3d,ierr1,iso1,mfr11, #ifour11,istep3d,epst03d,epstf3d,tref3d,NBELAS3D,NB_PARA_FLUENDO3D, #NB_PARA_SUPP_FLUENDO3D,NB_VARI_FLUENDO3D,NB_VARI_SUPP_FLUENDO3D, #NB_DECAL,NBNMAX3D,NBNB3D,IDIMB3D,XE3D,NBVIA3D,INLVIA3D, #NMAT0,NMAT1,NMAT2,NMAT3,NVAR0,NVAR1,NVAR2,NVAR3,NB_PARA_FIBRE_U) c************ re-affectation des tableaux de variables internes ******** do i=1,nvarib varf(i)=VARF3D(i) * print*,'cas3d f(',i,')=',varf(i),' i',var0(i) end do c read* c ecriture des contraintes totales fin de pas do i=1,nstrsb sigf(i)=sigf3d(i) end do c*********************************************************************** c traitement de l erreur d ecoulement if(ierr1.eq.0)then kerre=0 else kerre=1 end if return end c***********************************************************************