%Paper: cmp-soft/9303005 %From: %Date: Thu, 17 Jun 1993 13:23:34 +0300 \title{A Spring Element, Coupling Macroscopic Stress and Strain\\ to the Boundary of any\\ Representative Axisymmetric Microscopic Unit Cell} \author{Niklas J\"arvstr\aa t \\ Hydro Aluminium a.s, R\&D Centre, Karm\o y,\\ N-4265 H\aa vik, Norway} \begin{document} \maketitle \begin{abstract} {A user defined element suitable for imposing arbitrary global stress/strain states on a representative microcell modelled by ABAQUS axisymmetric elements. Three pieces of software are included: \begin{itemize} \item FORTRAN code of the ABAQUS subroutine UEL, describing the element. \item ABAQUS input file for the case of thermal cycling ofa metal matrix composite. (SiC fibres of aspect ratio 1:4 in AA6061 matrix, aligned but randomly distributed.) \item A separate FORTRAN program that calculates pseudomacrostresses, using the ABAQUS output file `*.FIL'. I.e. the stress is averaged over part of the unit cell model only. \end{itemize} A full description of the element, and a theoretical derivation of the boundary conditions used will appear in `An Ellipsoidal Unit Cell for the Calculation of Micro-Stresses in Short Fibre Composites', N J\"arvstr\aa t, {\em Comp Mat Sci}, 1993.} \end{abstract} \newpage \section{ABAQUS subroutine UEL} The element described in the reference cited in the abstract. \begin{itemize} \item A representative unit cell must be modelled, using ABAQUS axisymmetric stress/strain elements. The origin (0,0) must be a mirror point, so by symmetry only half the cell need to be modelled. \item Define one extra node -- the global node. (Its coordinates are irrelevant and are not used in calculations.) \item Define the user defined element in the ABAQUS input file (`*.INP'). The global node is the first node, and {\em all} nodes on the external boundary of the cell are the remaining node. (I don't know if ABAQUS puts a limit on the number of nodes that can be included in a single element, but 17 worked all right) \item Apply global strain as prescribed displacements of the global node, or global stress as external forces on the global node, or any combination of stress and strain. For consistency, the global stress must be multiplied by the total volume of the unit cell before being applied at the node(See the reference cited in the abstract for an explanation). \newpage \begin{verbatim} SUBROUTINE UEL(RHS,A,SVARS,ENERGY,JLINES,XVAL, 1 NDOFEL,NRHS,NSVARS,PROPS,NPROPS,COORDS,MCRD,NNODE, 2 U,DU,V,ACCLRN,JTYPE,TIME,DTIME,KSTEP,KINC,JELEM,PARAMS, 3 NDLOAD,JDLTYP,ADLMAG,PREDEF,NPRED,LFLAGS) C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION RHS(NDOFEL,NRHS),A(NDOFEL,NDOFEL),SVARS(1), 1 ENERGY(6),JLINES(1),XVAL(1),PROPS(1),COORDS(MCRD,NNODE), 2 U(NDOFEL),DU(NDOFEL),V(NDOFEL),ACCLRN(NDOFEL), 3 TIME(7),DTIME(7),PARAMS(1),JDLTYP(NDLOAD,NRHS), 4 ADLMAG(NDLOAD,NRHS,2),PREDEF(2,NPRED,NNODE),LFLAGS(4), 5 AKK(2) C C ******************************************************* C * * C * micro-MACRO coupling element for axisymmetry * C * by Niklas Jarvstrat * C * * C * (node 1 is the MACRO node, 2-NNODE micro boundary) * C * * C ******************************************************* C C C Artificial stiffness constant AK should not be too small C AK=PROPS(1) IF (AK.LE.1E6)THEN WRITE(*,*) "ERROR: AK= ",AK STOP ENDIF C C C Calculation of element stiffness matrix A C DO JK=1,NNODE-1 DO JL=1,2 JK2L=2*JK+JL A(JL,JK2L)=-COORDS(JL,JK+1)*AK A(JL,JL) = A(JL,JL) + COORDS(JL,JK+1)**2*AK A(JK2L,JK2L)=AK A(JK2L,JL)=-COORDS(JL,JK+1)*AK ENDDO ENDDO C C Calculation of energy and residual force contribution C ENERGY(6)=0D0 DO JD=1,2 RHS(JD,1) = RHS(JD,1) - A(JD,JD)*U(JD) DO JA=1,NNODE-1 JAD=2*JA+JD RHS(JD,1) = RHS(JD,1) - A(JD,JAD)*U(JAD) RHS(JAD,1) = - A(JAD,JAD)*U(JAD) - A(JAD,JD)*U(JD) ENERGY(6) = ENERGY(6) - RHS(JAD,1)*U(JAD)/2 ENDDO ENERGY(6) = ENERGY(6) - RHS(JD,1)*U(JD)/2 ENDDO DO JAD=1,NSVARS SVARS(JAD)=RHS(JAD,1) ENDDO IF (LFLAGS(1).GT.1.OR.LFLAGS(3).GT.1) THEN WRITE(*,*) "ERROR: LFLAGS= ",LFLAGS STOP ENDIF RETURN END \end{verbatim} \newpage \section{ABAQUS input file} An ABAQUS input file (`*.INP') featuring thermal cycling of a SiC reinforced AA6061 matrix composite. The fortran program above is precompiled and linked to ABAQUS before execution. The mesh of the ellipsoidal unit cell (axis ratio 1:4, the same as the fibre aspect ratio) was made with IDEAS -- sorry about the length of the file. \newpage \begin{verbatim} *HEADING THERMAL CYCLING 910305 *NODE 999, 3.9080E-01, 1.5632E+00, 0.0000E+00 *NODE, NSET=ALLNOD 1, 0.0000E+00, 0.0000E+00, 0.0000E+00 2, 8.9130E-02, 0.0000E+00, 0.0000E+00 3, 1.4572E-01, 0.0000E+00, 0.0000E+00 4, 1.7973E-01, 0.0000E+00, 0.0000E+00 5, 2.0025E-01, 0.0000E+00, 0.0000E+00 6, 2.1506E-01, 0.0000E+00, 0.0000E+00 7, 0.0000E+00, 3.1919E-01, 0.0000E+00 8, 8.9130E-02, 3.1919E-01, 0.0000E+00 9, 1.4572E-01, 3.1919E-01, 0.0000E+00 10, 1.7973E-01, 3.1919E-01, 0.0000E+00 11, 2.0025E-01, 3.1919E-01, 0.0000E+00 12, 2.1506E-01, 3.1919E-01, 0.0000E+00 13, 0.0000E+00, 5.3755E-01, 0.0000E+00 14, 8.9130E-02, 5.3755E-01, 0.0000E+00 15, 1.4572E-01, 5.3755E-01, 0.0000E+00 16, 1.7973E-01, 5.3755E-01, 0.0000E+00 17, 2.0025E-01, 5.3755E-01, 0.0000E+00 18, 2.1506E-01, 5.3755E-01, 0.0000E+00 19, 0.0000E+00, 6.7907E-01, 0.0000E+00 20, 8.9130E-02, 6.7907E-01, 0.0000E+00 21, 1.4572E-01, 6.7907E-01, 0.0000E+00 22, 1.7973E-01, 6.7907E-01, 0.0000E+00 23, 2.0025E-01, 6.7907E-01, 0.0000E+00 24, 2.1506E-01, 6.7907E-01, 0.0000E+00 25, 0.0000E+00, 7.6640E-01, 0.0000E+00 26, 8.9130E-02, 7.6640E-01, 0.0000E+00 27, 1.4572E-01, 7.6640E-01, 0.0000E+00 28, 1.7973E-01, 7.6640E-01, 0.0000E+00 29, 2.0025E-01, 7.6640E-01, 0.0000E+00 30, 2.1506E-01, 7.6640E-01, 0.0000E+00 31, 0.0000E+00, 8.2058E-01, 0.0000E+00 32, 8.9130E-02, 8.2058E-01, 0.0000E+00 33, 1.4572E-01, 8.2058E-01, 0.0000E+00 34, 1.7973E-01, 8.2058E-01, 0.0000E+00 35, 2.0025E-01, 8.2058E-01, 0.0000E+00 36, 2.1506E-01, 8.2058E-01, 0.0000E+00 37, 0.0000E+00, 8.6025E-01, 0.0000E+00 38, 8.9130E-02, 8.6025E-01, 0.0000E+00 39, 1.4572E-01, 8.6025E-01, 0.0000E+00 40, 1.7973E-01, 8.6025E-01, 0.0000E+00 41, 2.0025E-01, 8.6025E-01, 0.0000E+00 42, 2.1506E-01, 8.6025E-01, 0.0000E+00 43, 2.3506E-01, 0.0000E+00, 0.0000E+00 44, 2.6304E-01, 0.0000E+00, 0.0000E+00 45, 3.1043E-01, 0.0000E+00, 0.0000E+00 46, 3.9080E-01, 0.0000E+00, 0.0000E+00 47, 2.3411E-01, 3.1935E-01, 0.0000E+00 48, 2.6077E-01, 3.1958E-01, 0.0000E+00 49, 3.0592E-01, 3.1996E-01, 0.0000E+00 50, 3.8249E-01, 3.2062E-01, 0.0000E+00 51, 2.3233E-01, 5.3777E-01, 0.0000E+00 52, 2.5649E-01, 5.3808E-01, 0.0000E+00 53, 2.9740E-01, 5.3860E-01, 0.0000E+00 54, 3.6678E-01, 5.3949E-01, 0.0000E+00 55, 2.3062E-01, 6.7928E-01, 0.0000E+00 56, 2.5239E-01, 6.7957E-01, 0.0000E+00 57, 2.8925E-01, 6.8007E-01, 0.0000E+00 58, 3.5177E-01, 6.8091E-01, 0.0000E+00 59, 2.2932E-01, 7.6658E-01, 0.0000E+00 60, 2.4928E-01, 7.6682E-01, 0.0000E+00 61, 2.8307E-01, 7.6723E-01, 0.0000E+00 62, 3.4039E-01, 7.6793E-01, 0.0000E+00 63, 2.2842E-01, 8.2072E-01, 0.0000E+00 64, 2.4711E-01, 8.2091E-01, 0.0000E+00 65, 2.7876E-01, 8.2124E-01, 0.0000E+00 66, 3.3243E-01, 8.2179E-01, 0.0000E+00 67, 2.2770E-01, 8.6036E-01, 0.0000E+00 68, 2.4539E-01, 8.6050E-01, 0.0000E+00 69, 2.7535E-01, 8.6075E-01, 0.0000E+00 70, 3.2615E-01, 8.6117E-01, 0.0000E+00 71, 0.0000E+00, 8.9967E-01, 0.0000E+00 72, 8.7644E-02, 8.9883E-01, 0.0000E+00 73, 1.4329E-01, 8.9830E-01, 0.0000E+00 74, 1.7674E-01, 8.9798E-01, 0.0000E+00 75, 1.9691E-01, 8.9779E-01, 0.0000E+00 76, 2.1148E-01, 8.9765E-01, 0.0000E+00 77, 0.0000E+00, 9.6567E-01, 0.0000E+00 78, 8.5155E-02, 9.6343E-01, 0.0000E+00 79, 1.3922E-01, 9.6201E-01, 0.0000E+00 81, 1.9132E-01, 9.6064E-01, 0.0000E+00 85, 1.2945E-01, 1.1150E+00, 0.0000E+00 89, 2.2361E-01, 8.9884E-01, 0.0000E+00 90, 2.4058E-01, 9.0050E-01, 0.0000E+00 91, 2.6933E-01, 9.0331E-01, 0.0000E+00 92, 3.1809E-01, 9.0808E-01, 0.0000E+00 93, 2.1659E-01, 9.6324E-01, 0.0000E+00 95, 2.5849E-01, 9.7444E-01, 0.0000E+00 96, 3.0317E-01, 9.8639E-01, 0.0000E+00 99, 2.2763E-01, 1.1442E+00, 0.0000E+00 101, 0.0000E+00, 1.4140E+00, 0.0000E+00 102, 3.6124E-02, 1.4040E+00, 0.0000E+00 103, 7.3091E-02, 1.3832E+00, 0.0000E+00 106, 0.0000E+00, 1.4886E+00, 0.0000E+00 107, 4.5434E-02, 1.4756E+00, 0.0000E+00 108, 8.9871E-02, 1.4428E+00, 0.0000E+00 109, 1.3346E-01, 1.3912E+00, 0.0000E+00 111, 0.0000E+00, 1.5632E+00, 0.0000E+00 112, 5.3338E-02, 1.5485E+00, 0.0000E+00 113, 1.0533E-01, 1.5054E+00, 0.0000E+00 114, 1.5535E-01, 1.4343E+00, 0.0000E+00 116, 0.0000E+00, 1.3135E+00, 0.0000E+00 118, 4.4000E-02, 1.3100E+00, 0.0000E+00 119, 5.7000E-02, 1.2170E+00, 0.0000E+00 140, 2.0300E-01, 9.7000E-01, 0.0000E+00 141, 2.2900E-01, 9.7900E-01, 0.0000E+00 142, 1.7000E-01, 9.7600E-01, 0.0000E+00 148, 1.5800E-01, 1.0950E+00, 0.0000E+00 150, 2.1200E-01, 1.1120E+00, 0.0000E+00 155, 2.6417E-01, 1.1520E+00, 0.0000E+00 156, 1.9700E-01, 1.0020E+00, 0.0000E+00 157, 1.8500E-01, 1.1090E+00, 0.0000E+00 158, 1.7100E-01, 1.1540E+00, 0.0000E+00 163, 8.0000E-02, 1.3050E+00, 0.0000E+00 165, 6.7000E-02, 1.1160E+00, 0.0000E+00 166, 1.2400E-01, 1.3100E+00, 0.0000E+00 167, 1.0900E-01, 1.3650E+00, 0.0000E+00 170, 1.9759E-01, 1.3487E+00, 0.0000E+00 171, 2.3581E-01, 1.2466E+00, 0.0000E+00 172, 0.0000E+00, 1.1097E+00, 0.0000E+00 173, 1.9500E-01, 1.2390E+00, 0.0000E+00 174, 1.6800E-01, 1.3210E+00, 0.0000E+00 175, 1.4500E-01, 1.2300E+00, 0.0000E+00 176, 1.0500E-01, 1.2250E+00, 0.0000E+00 177, 0.0000E+00, 1.2173E+00, 0.0000E+00 *NSET,NSET=INMAT 6,12,18,24,30,36,37,38,39,40,41,42 *NCOPY,OLD SET=INMAT,CHANGE NUMBER=200,NEW SET=INFIB,SHIFT 0,0,0 0,0,0, 1,0,0 ,0 *NSET,NSET=ALLNOD INFIB *NSET,NSET=BOUND 46,50,54,58,62,66,70,92,96,155,171,170,114,113,112,111 *NSET,NSET=SYMM2 1,2,3,4,5,6,43,44,45,46 *NSET,NSET=SYMM1 1,7,13,19,25,31,37,71,77,101,106,111,116,172,177 *USER ELEMENT,NODES=17,TYPE=U1,COORDINATES=2,PROPERTIES=1,VARIABLES=34 1,2 *UEL PROPERTY,ELSET=MICMAC 1.E15 *ELEMENT,TYPE=U1 ,ELSET=MICMAC 999, 999, 46,50,54,58,62, 66,70,92,96,155 ,171,170,114,113 112,111 *ELEMENT,TYPE=CAX4 ,ELSET=FIBER 1001, 1, 2, 8, 7 1002, 2, 3, 9, 8 1003, 3, 4, 10, 9 1004, 4, 5, 11, 10 1005, 5, 206, 212, 11 1006, 7, 8, 14, 13 1007, 8, 9, 15, 14 1008, 9, 10, 16, 15 1009, 10, 11, 17, 16 1010, 11, 212, 218, 17 1011, 13, 14, 20, 19 1012, 14, 15, 21, 20 1013, 15, 16, 22, 21 1014, 16, 17, 23, 22 1015, 17, 218, 224, 23 1016, 19, 20, 26, 25 1017, 20, 21, 27, 26 1018, 21, 22, 28, 27 1019, 22, 23, 29, 28 1020, 23, 224, 230, 29 1021, 25, 26, 32, 31 1022, 26, 27, 33, 32 1023, 27, 28, 34, 33 1024, 28, 29, 35, 34 1025, 29, 230, 236, 35 1026, 31, 32, 238, 237 1027, 32, 33, 239, 238 1028, 33, 34, 240, 239 1029, 34, 35, 241, 240 1030, 35, 236, 242, 241 *ELEMENT,TYPE=CAX4 ,ELSET=MATRIX 31 6 43 47 12 32 43 44 48 47 33 44 45 49 48 34 45 46 50 49 35 12 47 51 18 36 47 48 52 51 37 48 49 53 52 38 49 50 54 53 39 18 51 55 24 40 51 52 56 55 41 52 53 57 56 42 53 54 58 57 43 24 55 59 30 44 55 56 60 59 45 56 57 61 60 46 57 58 62 61 47 30 59 63 36 48 59 60 64 63 49 60 61 65 64 50 61 62 66 65 51 36 63 67 42 52 63 64 68 67 53 64 65 69 68 54 65 66 70 69 55 37 38 72 71 56 38 39 73 72 57 39 40 74 73 58 40 41 75 74 59 41 42 76 75 60 71 72 78 77 61 72 73 79 78 62 73 74 142 79 63 74 75 81 142 64 75 76 140 81 65 77 78 165 172 66 78 79 85 165 67 79 142 148 85 70 42 67 89 76 71 67 68 90 89 72 68 69 91 90 73 69 70 92 91 74 76 89 93 140 75 89 90 141 93 76 90 91 95 141 77 91 92 96 95 80 141 95 99 150 81 95 96 155 99 82 101 102 107 106 83 102 103 108 107 84 103 167 109 108 85 167 166 174 109 86 106 107 112 111 87 107 108 113 112 88 108 109 114 113 89 109 174 170 114 90 172 165 119 177 91 165 85 176 119 92 177 119 118 116 93 116 118 102 101 94 119 176 163 118 95 118 163 103 102 103 163 166 167 103 104 176 175 166 163 108 175 173 174 166 109 173 171 170 174 110 142 81 140 156 111 140 93 141 156 112 142 156 157 148 113 156 141 150 157 114 85 148 157 158 115 157 150 99 158 116 85 158 175 176 117 158 99 173 175 118 99 155 171 173 *ELSET,ELSET=WHOLE FIBER,MATRIX *ELSET,ELSET=PLOT1 FIBER *ELSET,ELSET=PLOT1,GENERATE 31,54,1 *ELSET,ELSET=PLOT2,GENERATE 55,67,1 *ELSET,ELSET=PLOT2,GENERATE 70,77,1 *ELSET,ELSET=PLOT2,GENERATE 80,95,1 *ELSET,ELSET=PLOT2 103, 104 *ELSET,ELSET=PLOT2,GENERATE 108,118,1 *SOLID SECTION,ELSET=MATRIX,MATERIAL=AA6061 *MATERIAL,NAME=AA6061 ** METALS HANDBOOK *ELASTIC 69.E9 ,0.33 ,25. *DENSITY 1.0 *EXPANSION 21.8E-6 ,-60 21.8E-6 ,-20 23.4E-6 ,60. 25.1E-6 ,150. 27.6E-6 ,250. 27.6E-6 ,300. *PLASTIC 160.E6 ,0 ,-195. 300.E6 ,0.3,-195. 150.E6 ,0 ,-80 260.E6 ,.26 ,-80 145.E6 ,0 ,24 241.E6 ,.25 ,24 103.E6 ,0 ,204 131.E6 ,.28 ,204 34.E6 ,0 ,260 51.E6 ,.60 ,260 19.E6 ,0 ,316 32.E6 ,.85 ,316 12.E6 ,0 ,371 24.E6 ,.95 ,371 *SOLID SECTION,ELSET=FIBER,MATERIAL=SIC *MATERIAL,NAME=SIC *ELASTIC ** Eng Prop Data on Sel Ceramics 448.2E9,0.2 *DENSITY 1.0 *EXPANSION ** Richard Warren 4.5E-6 *RESTART,WRITE,FREQUENCY=350 *MPC 5,INFIB,INMAT *BOUNDARY 1,PINNED SYMM1,XSYMM SYMM2,YSYMM *INITIAL CONDITIONS,TYPE=TEMPERATURE ALLNOD,450 *STEP ,INC=80 ,CYCLE=15 ,AMPLITUDE=RAMP ,SUBMAX ** ** Water quench ** *STATIC,PTOL=0.1 1,350,1,10 *TEMPERATURE ALLNOD,100 *EL PRINT, FREQUENCY=0,ELSET=MICMAC SDV *NODE PRINT, FREQUENCY=0 RF *EL FILE,ELSET=WHOLE,POSITION=CENTROIDAL, FREQ=0 S,E *ENDSTEP \end{verbatim} \newpage \section{FORTRAN program calculating pseudomacrostresses} This is a program that averages the stress over two separate sets of elements of the model (element numbers $\le 1000$, $>1000$). Two files are needed to run this program: \begin{enumerate} \item A file `*.FIL' containing the kinetic energy at centroidal points of all elements, obtained from a separate ABAQUS calculation, using the same geometry, but plane (strain or stress) elements and imposing a constant velocity in any direction. The density must be given in this calculation, but need only be homogeneous and can be taken as unity. This file, the `mass' file, is necessary to specify the volume of each element as needed for the calculation of averages. \item The file `*.FIL' containing stresses at centroidal points of the elements calculated using axisymmetric elements. \end{itemize} The stress is averaged over all elements with element number less than 1000. The program is easily modified for the calculation of other averaged entities by using information stored under other keys. \newpage \begin{verbatim} PROGRAM PSEUDOMACCAL C IMPLICIT DOUBLE PRECISION (A-H,K-Z), INTEGER (I,J) LOGICAL STOP CHARACTER FILNAM*11,FILNUT*7,DATAFIL*15,HEADING*80 PARAMETER(ITOTEL=500) DIMENSION ARRAY(513),JRRAY(2,513),R(ITOTEL),MASS(ITOTEL), & STRESS(3),MATSTR(3) EQUIVALENCE (ARRAY(1),JRRAY(1,1)),(ARRAY(3),HEADING) COMMON ARRAY ITYP=0 I2001=0 IELNO=0 C C C READ FILENAMES AND OPEN FILES C C WRITE(*,*) ' *******************************************************' WRITE(*,*) ' * *' WRITE(*,*) ' * ABAQUS *.FIL post processing: *' WRITE(*,*) ' * Pseudomacro stress calculation *' WRITE(*,*) ' * *' WRITE(*,*) ' * by Niklas Jarvstrat *' WRITE(*,*) ' * *' WRITE(*,*) ' *******************************************************' WRITE(*,*) '' WRITE(*,*) ' GIVE MASS FILENAME (WITHOUT .fil)' 111 FORMAT (A7) READ(*,111) FILNAM DATAFIL=FILNAM DO J=1,8 IF (DATAFIL(J:J).EQ.' ') GOTO 1111 ENDDO 1111 DATAFIL(J:J+3)='.fil' OPEN (UNIT=97, FILE=DATAFIL, STATUS='OLD',FORM='UNFORMATTED') REWIND 97 WRITE(*,*) ' GIVE STRESS FILENAME (WITHOUT .fil)' READ(*,111) FILNAM DATAFIL=FILNAM DO J=1,8 IF (DATAFIL(J:J).EQ.' ') GOTO 1112 ENDDO 1112 DATAFIL(J:J+3)='.fil' OPEN (UNIT=98, FILE=DATAFIL, STATUS='OLD',FORM='UNFORMATTED') REWIND 98 WRITE(*,*) ' GIVE FINAL STEP NUMBER' READ(*,*) IFSTEP WRITE(*,*) ' GIVE FINAL INCREMENT NUMBER' READ(*,*) IFINCR WRITE(*,*) ' GIVE OUTPUT FILENAME' READ(*,111) FILNUT OPEN(UNIT=99,FILE=FILNUT ,STATUS='NEW') JB=1 FIBMAS=0.0 MATMAS=0.0 STOP=(.FALSE.) C C Read the masses of all elements C 10 CALL READFILE(JB,STOP,97) IF (JRRAY(1,2).EQ.2000) THEN I=0 ELSEIF (JRRAY(1,2).EQ.2001) THEN IF (I2001.GT.1) GOTO 100 I2001=I2001+1 ELSEIF (JRRAY(1,2).EQ.1922) THEN WRITE(99,3) HEADING ELSEIF (JRRAY(1,2).EQ.1) THEN I=I+1 IELNO=JRRAY(1,3) ELSEIF (JRRAY(1,2).EQ.8) THEN R(I)=ARRAY(3) ITYP=8 ELSEIF (JRRAY(1,2).EQ.19) THEN C Key no 19 = kinetic energy IF (ITYP .EQ. 8) THEN I=1 ENDIF ITYP=19 MASS(I)=ARRAY(3)*6.283185308*R(I) IF (IELNO .GT. 1000) THEN WRITE(99,5) IELNO,R(I),MASS(I),ARRAY(3) FIBMAS=FIBMAS+MASS(I) ELSE WRITE(99,6) IELNO,R(I),MASS(I),ARRAY(3) MATMAS=MATMAS+MASS(I) ENDIF IMAX=I ENDIF IF (.NOT.STOP) GO TO 10 100 JB=1 I=0 I2001=0 TOMAS=FIBMAS+MATMAS WRITE(99,*) ' ' WRITE (99,4) FIBMAS,MATMAS,TOMAS,IMAX WRITE(99,*) ' ' WRITE (99,2) ' (MPa or %) ','Strans','Saxial','Sshear' WRITE(99,*) ' ' STOP=(.FALSE.) C C Read the file containing stresses C 20 CALL READFILE(JB,STOP,98) IF (JRRAY(1,2).EQ.2000) THEN I=0 I2001=0 TIME=ARRAY(3) ISTEP=JRRAY(1,8) IINCR=JRRAY(1,9) IF ( (ISTEP.GT.IFSTEP) .OR. & (ISTEP.GE.IFSTEP .AND. IINCR.GT.IFINCR) ) GO TO 200 DO 2001 IS=1,3 STRESS(IS)=0.0 2001 MATSTR(IS)=0.0 ELSEIF (JRRAY(1,2).EQ.1) THEN C key no 1 = element output: read element number I=I+1 IELNO=JRRAY(1,3) ILOCID=JRRAY(1,6) ELSEIF (JRRAY(1,2).EQ.11) THEN C key no 11 = stresses: read and accumulate stresses IF (ILOCID.EQ.1) THEN C centroidal IF (IELNO .LE. 1000) THEN CALL ACC(MATMAS,TOMAS,MASS(I),1.D6,MATSTR,STRESS) ENDIF IF (I.GT.IMAX) THEN WRITE(*,*) 'IMAX=',IMAX,' > I=',I,IELNO,JRRAY(1,2),ISTEP,IINCR ENDIF ELSEIF (ILOCID.EQ.0) THEN C integration points(four noded elements, C but less accurate than centroidal) IF (IELNO .LE. 1000) THEN CALL ACC(MATMAS,TOMAS,MASS((I+3)/4)/4,1.D6,MATSTR,STRESS) ENDIF IF (I.GT.4*IMAX) THEN WRITE(*,*) 'IMAX=',IMAX,' > I=',I,IELNO,JRRAY(1,2),ISTEP,IINCR ENDIF ENDIF ELSEIF (JRRAY(1,2).EQ.2001) THEN WRITE(*,1) ISTEP,IINCR IF (I.NE.0) THEN WRITE(99,*) WRITE(99,*) TIME,MATSTR IF (I2001.GT.2) GOTO 200 ENDIF I2001=I2001+1 ENDIF IF (.NOT.STOP) GOTO 20 200 CLOSE (97) CLOSE (98) CLOSE (99) 1 FORMAT (' ','STEP:',I3,' INC:',I2,' ',5(' ',E8.2)) 7 FORMAT (' ',A,5(' ',E8.2)) 2 FORMAT (' ',A,4(' ',A)) 3 FORMAT (A,A) 4 FORMAT (' ','FIBRE MASS= ',F8.3,' MATRIX MASS= ',F8.3, & ' TOTAL MASS= ',F8.3,' I= ',I4) 5 FORMAT (' ',I5,'Radius= ',F8.5,' FIBRE MASS= ',2F8.5) 6 FORMAT (' ',I5,'Radius= ',F8.5,' MATRIX MASS= ',2F8.5) END SUBROUTINE READFILE(JB,STOP,IFIL) C C READS FROM FILE, and converts to the format described in ABAQUS user's manual C (Tested on Apollo and IBM workstations, ABAQUS versions 4.8, 4.9) C IMPLICIT DOUBLE PRECISION (A-H,K-Z), INTEGER (I,J) LOGICAL STOP DIMENSION ARRAY(513),JRRAY(2,513),BLOCK(513) EQUIVALENCE (ARRAY(1),JRRAY(1,1)) COMMON ARRAY SAVE JRCD,BLOCK,JR C C C IF (JB .EQ.1) THEN READ(IFIL,IOSTAT=JRCD,END=20) (BLOCK(JR),JR=1,512) ENDIF 20 IF (JB .GE. 513) THEN IF (JRCD.NE.0) GO TO 200 READ(IFIL,IOSTAT=JRCD,END=21) (BLOCK(JR),JR=1,512) 21 JB=1 ENDIF ARRAY(1)=BLOCK(JB) JB=JB+1 DO J=2,JRRAY(1,1) IF (JB .GE. 513) THEN IF (JRCD.NE.0) GO TO 200 READ(IFIL,IOSTAT=JRCD,END=22) (BLOCK(JR),JR=1,512) 22 JB=1 ENDIF ARRAY(J)=BLOCK(JB) JB=JB+1 ENDDO IF (JRCD.NE.0 .AND. JB.GE.JR) THEN GOTO 200 ENDIF RETURN 200 STOP=.TRUE. RETURN END SUBROUTINE ACC(MC,MT,M,SCALE,VC,VT) C C ADD STRESS FROM ONE POINT TO THE C COMPONENT STRESS VALUE, weighted by the mass fraction C associated with the point C IMPLICIT DOUBLE PRECISION (A-H,K-Z), INTEGER (I,J) DIMENSION ARRAY(513),VC(3),VT(3) COMMON ARRAY C C X0=M/SCALE C The two transverse directions are identical(in average), so C the correct value is obtained by averaging pointwise first, C thus avoiding the effects of cylindrical coordinates. X1=(ARRAY(3)+ARRAY(5))*X0/2.0 VC(1)=VC(1)+X1/MC VT(1)=VT(1)+X1/MT DO J=2,3 VC(J)=VC(J)+ARRAY(2*J)*X0/MC VT(J)=VT(J)+ARRAY(2*J)*X0/MT ENDDO RETURN END \end{verbatim} \end{document} {}.