1 : C LFROA SOURCE FANDEUR 11/07/19 21:16:22 7042 2 : 3 : SUBROUTINE LFROA(IPOGEO,IPMATR,IPMINT,IVAMAT, 4 : 1 IVACAR,MELE,MFR,LRE,NDDL) 5 : C 6 : C*********************************************************************** 7 : C * 8 : C Routine appelée par FRVISQ. * 9 : C * 10 : C Calcule l'amortissement de frontière pour les liquides de face LSE2 * 11 : C dans le cas 2D, ou de face LTR3 ou LQU4 dans le cas 3D. * 12 : C * 13 : C Entrées : * 14 : C -------- * 15 : C * 16 : C IPOGEO : pointeur sur le maillage de l'enveloppe des massifs, * 17 : C type MELEME * 18 : C IPMATR : pointeur sur le segment IMATRI, chapeau des rigidités * 19 : C élémentaires * 20 : C IPMINT : pointeur sur le segment d'intégration, type MINTE * 21 : C IVAMAT : pointeur sur un segment MPTVAL de données matériau * 22 : C IVACAR : pointeur sur un segment MPTVAL de caractéristiques * 23 : C (épaisseur dans le cas contraintes planes) * 24 : C MELE : numéro de l'élément fini associé à la face du massif * 25 : C MFR : numéro de la formulation * 26 : C LRE : taille de la matrice d'amortissement à construire * 27 : C NDDL : nombre de degrés de liberté * 28 : C * 29 : C Remplit le segment XMATRI pour chaque élément de la sous-zone. * 30 : C*********************************************************************** 31 : C 32 : IMPLICIT INTEGER(I-N) 33 : IMPLICIT REAL*8(A-H,O-Z) 34 : C 35 : -INC CCOPTIO 36 : -INC CCREEL 37 : 38 : -INC SMRIGID 39 : -INC SMELEME 40 : -INC SMCOORD 41 : -INC SMCHAML 42 : -INC SMINTE 43 : C 44 : SEGMENT MPTVAL 45 : INTEGER IPOS(NS),NSOF(NS) 46 : INTEGER IVAL(NCOSOU) 47 : CHARACTER*16 TYVAL(NCOSOU) 48 : ENDSEGMENT 49 : C 50 : SEGMENT,MWORK 51 : REAL*8 XE(3,NBNN) 52 : REAL*8 REL(LRE,LRE) 53 : REAL*8 SHPWRK(6,NBNN),BGENE(NDDL,LRE) 54 : REAL*8 VALMAT(NV1) 55 : REAL*8 RCLMAT(NDDL,NDDL) 56 : ENDSEGMENT 57 : C 58 : MELEME=IPOGEO 59 : c* SEGACT MELEME 60 : NBNN=NUM(/1) 61 : NBELEM=NUM(/2) 62 : C 63 : MINTE=IPMINT 64 : c* SEGACT,MINTE 65 : NBPGAU=POIGAU(/1) 66 : C 67 : DIM3=1.D0 68 : C 69 : xMATRI=IPMATR 70 : c* SEGACT,xMATRI*MOD 71 : c* NLIGRD=LRE 72 : c* NLIGRP=LRE 73 : 74 : NV1=5 75 : SEGINI,MWORK 76 : C 77 : C boucle sur les éléments 78 : C 79 : DO 1 IB=1,NBELEM 80 : C 81 : C on cherche les coordonnées de l'élément IB 82 : C 82 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 82 : CALL ZERO(REL,LRE,LRE) 82 : CALL ZERO(RCLMAT,NDDL,NDDL) 83 : C 84 : C boucle sur les points de Gauss 85 : C 86 : DO 10 IGAU=1,NBPGAU 87 : C 88 : C récupération de l'épaisseur 89 : C 90 : IF (IFOUR.EQ.-2) THEN 91 : IF (IVACAR.NE.0) THEN 92 : MPTVAL=IVACAR 93 : MELVAL=IVAL(1) 94 : IF (MELVAL.NE.0) THEN 95 : IGMN=MIN(IGAU,VELCHE(/1)) 96 : IBMN=MIN(IB,VELCHE(/2)) 97 : DIM3=VELCHE(IGMN,IBMN) 98 : ELSE 99 : DIM3=1.D0 100 : ENDIF 101 : ENDIF 102 : ENDIF 103 : C 104 : C récupération des données matériau 105 : C 106 : MPTVAL=IVAMAT 107 : DO 11 J=1,5 108 : MELVAL=IVAL(J) 109 : IGMN=MIN(IGAU,VELCHE(/1)) 110 : IBMN=MIN(IB,VELCHE(/2)) 111 : VALMAT(J)=VELCHE(IGMN,IBMN) 112 : 11 CONTINUE 113 : C 114 : RHO=VALMAT(1) 115 : CSON=VALMAT(2) 116 : ROREF=VALMAT(3) 117 : CREF=VALMAT(4) 118 : RLCAR=VALMAT(5) 119 : C 120 : C coefficient d'amortissement normalisé 121 : C 122 : RCL=(1.D0/CSON/rho)*ROREF*RLCAR*roref*cref**2/rlcar 123 : C 124 : C matrice RCLMAT 125 : C 126 : RCLMAT(1,1)=RCL 127 : RCLMAT(2,2)=RCL 128 : C 129 : C calcul de la matrice N des fonctions de forme 130 : C 131 : XDPGE=0.D0 132 : YDPGE=0.D0 132 : CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 133 : 1 DIM3,XE,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) 134 : C 135 : C mise à zéro des composantes sur p 136 : C 137 : DO 21 J=1,LRE 138 : BGENE(2,J)=0.D0 139 : 21 CONTINUE 140 : C 141 : C calcul du jacobien 142 : C 143 : IF (IDIM.EQ.3) THEN 144 : DXDQSI=0.D0 145 : DYDQSI=0.D0 146 : DZDQSI=0.D0 147 : DXDETA=0.D0 148 : DYDETA=0.D0 149 : DZDETA=0.D0 150 : DO 40 I=1,NBNN 151 : DXDQSI=DXDQSI+SHPTOT(2,I,IGAU)*XE(1,I) 152 : DYDQSI=DYDQSI+SHPTOT(2,I,IGAU)*XE(2,I) 153 : DXDETA=DXDETA+SHPTOT(3,I,IGAU)*XE(1,I) 154 : DYDETA=DYDETA+SHPTOT(3,I,IGAU)*XE(2,I) 155 : DZDQSI=DZDQSI+SHPTOT(2,I,IGAU)*XE(3,I) 156 : DZDETA=DZDETA+SHPTOT(3,I,IGAU)*XE(3,I) 157 : 40 CONTINUE 158 : DJAC=SQRT((DYDQSI*DZDETA-DYDETA*DZDQSI)**2+ 159 : 1 (DXDETA*DZDQSI-DXDQSI*DZDETA)**2+ 160 : 2 (DXDQSI*DYDETA-DXDETA*DYDQSI)**2) 161 : ELSE 162 : c* ELSE IF(IDIM.EQ.2) THEN 163 : DXDQSI=0.D0 164 : DYDQSI=0.D0 165 : DO I=1,NBNN 166 : DXDQSI=DXDQSI+SHPTOT(2,I,IGAU)*XE(1,I) 167 : DYDQSI=DYDQSI+SHPTOT(2,I,IGAU)*XE(2,I) 168 : ENDDO 169 : DJAC=SQRT(DXDQSI**2+DYDQSI**2) 170 : ENDIF 171 : C 172 : C calcul de l'élément de volume 173 : C 174 : IF (IFOUR.LT.0.OR.IFOUR.EQ.2) THEN 175 : R=1.D0 176 : IF (IFOUR.EQ.-2) R=DIM3 177 : ELSE 178 : R=0.D0 179 : DO I=1,NBNN 180 : R=R+SHPTOT(1,I,IGAU)*XE(1,I) 181 : ENDDO 182 : IF (IFOUR.EQ.0.OR.(IFOUR.EQ.1 183 : 1 .AND.NIFOUR.EQ.0)) THEN 184 : R=2*XPI*R 185 : ELSE 186 : c* ELSEIF (IFOUR.EQ.1.AND.NIFOUR.NE.0) THEN 187 : R=XPI*R 188 : ENDIF 189 : ENDIF 190 : C 191 : C construction de la matrice d'amortissement 192 : C 193 : DJAC=DJAC*POIGAU(IGAU)*R 193 : CALL BDBST(BGENE,DJAC,RCLMAT,LRE,NDDL,REL) 194 : C 195 : 10 CONTINUE 196 : C 197 : C remplissage de XMATRI 198 : C 198 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 199 : * 200 : ** la matrice calculée ci dessus serait correcte si l'on avait 201 : ** une formulation uniquement en pi. Comme on retient la formulation 202 : ** en p et pi et pour ne rien ajouter sur la ligne en pi (si non on 203 : ** ne satisfait plus la relation entre p et pi) on aboutit à une matrice 204 : ** disymetrique avec des termes uniquement sur le lignes en pi et collones 205 : ** en p. D'où la matrice suivante : 206 : 207 : do 30 i = 1, lre 208 : do 31 j = 2, lre 209 : ix = mod(j,2) 210 : if (ix.eq.0) then 211 : re(i,j,ib) = -re( i,(j-1),ib) 212 : re( i,(j-1),ib) = 0.d0 213 : endif 214 : 31 continue 215 : 30 continue 216 : 217 : 1 CONTINUE 218 : 219 : SEGSUP,MWORK 220 : c* SEGDES MELEME,MINTE,xMATRI 221 : 222 : RETURN 223 : END 224 : 225 :
© Cast3M 2003 - All rights reserved.
Disclaimer