C CAPAC1 SOURCE FANDEUR 11/07/19 21:15:25 7042 C======================================================================= C= C A P A C 1 = C= ----------- = C= = C= Fonction : = C= ---------- = C= Calcul de la matrice de CAPACITE CALORIFIQUE pour des elements de = C= COQue axisymetrique (COQ2) a integration semi-analytique = C= = C= Parametres : (E)=Entree (S)=Sortie = C= ------------ = C= NEF (E) Numero de l'ELEMENT FINI dans NOMTP = C= IMAIL (E) Numero du segment IMODEL dans le segment MMODEL = C= IPCHA1 (E) Pointeur sur un segment MCHEL1 de caracteristiques= C= CLAT (E) Chaleur latente du changement de phase = C= IPRIGI (E/S) Matrice de CAPACITE (RIGIDITE) resultat (ACTIF) = C= = C= P. DOWLATYARI, juin 1990 (adaptation de capanu.eso) = C======================================================================= SUBROUTINE CAPAC1 (NEF,IPMAIL,IPINTE,IVAMAT,NVAMAT,IVAPHA,NVAPHA, & IPMATR,NLIGR) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC CCOPTIO -INC CCREEL -INC CCHAMP -INC SMCHAML -INC SMCOORD -INC SMELEME -INC SMINTE -INC SMRIGID SEGMENT MMAT1 REAL*8 XE(3,NBNN),FORME(NBNN) REAL*8 CAPSS(NBNN,NBNN),CAPV(NLIGR,NLIGR) ENDSEGMENT SEGMENT MPTVAL INTEGER IPOS(NS),NSOF(NS),IVAL(NCOSOU) CHARACTER*16 TYVAL(NCOSOU) ENDSEGMENT REAL*8 VACOMP(3) C= Quelques constantes numeriques PARAMETER (X2Pi=6.283185307179586476925286766559D0) PARAMETER (X1s15=0.066666666666666666666666666667D0) PARAMETER (X2s15=0.133333333333333333333333333333D0) PARAMETER (X8s15=0.533333333333333333333333333333D0) PARAMETER (X1s30=0.033333333333333333333333333333D0) C- Element purement axisymetrique : IF (IFOMOD.NE.0) THEN CALL ERREUR(19) RETURN ENDIF C 1 - INITIALISATIONS ET VERIFICATIONS C ====================================== MELEME = IPMAIL c* SEGACT,MELEME NBNN = NUM(/1) NBELEM = NUM(/2) NBNN2 = 2*NBNN c* NBNN3 = 3*NBNN C ===== MINTE = IPINTE c* SEGACT,MINTE NBPGAU = POIGAU(/1) C- Petit test utile ? NBNO = SHPTOT(/2) IF (NBNO.NE.2) THEN CALL ERREUR(5) RETURN ENDIF C ===== MPTVAL = IVAMAT c* SEGACT,MPTVAL c*C- Verification sur la constance du champ d'epaisseur : c*C- epaisseur toujours placee en derniere position du mptval c* IPMELV = IVAL(NVAMAT) c* CALL QUELCH(IPMELV,IOK) c* IF (IOK.NE.0) THEN c* CALL ERREUR(566) c* GOTO 9990 c* ENDIF C ===== c* IF (IVAPHA.NE.0) THEN c* MPTVAL = IVAPHA c* SEGACT,MPTVAL c* ENDIF C ===== XMATRI = IPMATR c* SEGACT,XMATRI*MOD c* NLIGRP = NBNN3 = NLIGR c* NLIGRD = NBNN3 = NLIGR C ===== SEGINI,MMAT1 C 2 - BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IPMAIL C ============================================================ DO iElt = 1, NBELEM C ===== C 2.1 - Recuperation des coordonnees GLOABLES des noeuds de l'element C ===== CALL ZERO(CAPV,NLIGR,NLIGR) C ===== C 2.2 - Recuperation des coordonnees GLOABLES des noeuds de l'element C ===== CALL DOXE(XCOOR,IDIM,NBNN,NUM,iElt,XE) C ===== C 2.3 - Boucle sur les points de Gauss de l'element iElt C ===== DO iGau = 1, NBPGAU C ======= C 2.3.1 - Calcul du volume associe a ce point de Gauss (jacobien) C ======= DLX = SHPTOT(2,1,iGau)*XE(1,1)+SHPTOT(2,2,iGau)*XE(1,2) DLY = SHPTOT(2,1,iGau)*XE(2,1)+SHPTOT(2,2,iGau)*XE(2,2) DJAC = SQRT(DLX*DLX+DLY*DLY) C- Prise en compte de l'axisymetrie CALL DISTRR(XE,SHPTOT(1,1,iGau),NBNN,RR) DJAC = ABS(X2Pi*RR*DJAC) C ======= C 2.3.3 - Verification que le volume (jacobien) n'est pas nul en ce C point de Gauss --> Erreur C ======= IF (DJAC.LT.XPETIT) THEN INTERR(1) = iElt CALL ERREUR(259) GOTO 9990 ENDIF C ======= C 2.3.4 - Calcul du terme Rho.Cp.Vol en ce point de Gauss (avec prise C en compte d'un eventuel changement de phase) C ======= CALL CAPA5(iGau,iElt,IVAMAT,NVAMAT,IVAPHA,NVAPHA, VACOMP) CAPA = DJAC * POIGAU(iGau) * VACOMP(1) C ======= C 2.3.5 - Calcul de la contribution du point de Gauss a la matrice C CAPACITE elementaire pour cet element fini C ======= CALL ZERO(CAPSS,NBNN,NBNN) DO i = 1,NBNN FORME(i) = SHPTOT(1,i,iGau) ENDDO CALL NTNST(FORME,CAPA,NBNN,1,CAPSS) C ======= C 2.3.6 - Ajout de termes specifiques dus a l'integration (analytique) C suivant l'epaisseur de l'element de type COQUE C ======= EP = VACOMP(3) C- Erreur si l'epaisseur est est nulle c* IF (EP.LE.XPetit) THEN c* CALL ERREUR(517) c* GOTO 9990 c* ENDIF C1 = X2s15*EP C2 = X1s15*EP C3 = -X1s30*EP C4 = X8s15*EP C5 = C2 C6 = C1 DO j=1,NBNN j1 = j + NBNN j2 = j + NBNN2 DO i=1,NBNN i1 = i + NBNN i2 = i + NBNN2 Cte = CAPSS(i,j) CAPV( i, j) = CAPV( i, j) + C1*Cte CAPV(i1, j) = CAPV(i1, j) + C2*Cte CAPV(i2, j) = CAPV(i2, j) + C3*Cte CAPV(i1,j1) = CAPV(i1,j1) + C4*Cte CAPV(i2,j1) = CAPV(i2,j1) + C5*Cte CAPV(i2,j2) = CAPV(i2,j2) + C6*Cte ENDDO ENDDO ENDDO C ===== C 2.4 - Stockage de la matrice de CAPACITE pour cet element fini C (remplissage de XMATRI) C ===== CALL REMPMT(CAPV,NLIGR,RE(1,1,iElt)) ENDDO C 3 - MENAGE : DESACTIVATION/DESTRUCTION DE SEGMENTS C ==================================================== 9990 CONTINUE SEGSUP,MMAT1 c* SEGDES,MELEME,MINTE,XMATRI c* MPTVAL = IVAMAT c* SEGDES,MPTVAL c* IF (IVAPHA.NE.0) THEN c* MPTVAL = IVAPHA c* SEGDES,MPTVAL c* ENDIF RETURN END