c prestress.for c calculates wire overtensioning during stringing of the drift chamber c Roark & Young, table 24, case 2, with a step implicit real*8 (a-h), (o-z) common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99), . holes(99), stresscon(99), n, rmid(99), ihole(99), . layer(99) dimension amat(4,4), cvec(4), vbeg(4), oswire(3,99), area(99), . twire(3), Ewire(3), dwire(3), EAwire(3), fwire(3), . loadfac(3,99), ncelltot(99), nwire(3), dhole(3), . dzf(2,99), dzs(2,99), dznow(2,99) c bi = inner radius of region i data bi / 0.236, ! inner radius of plate . 0.2530, 0.2662, 0.2781, 0.2899, 0.3016, ! superlayer 1 . 0.3124, 0.3243, 0.3361, 0.3479, 0.3596, ! superlayer 2 . 0.3646, 0.3764, 0.3882, 0.4000, 0.4118, ! superlayer 3 . 0.4168, 0.4286, 0.4404, 0.4522, 0.4639, ! superlayer 4 . 0.4747, 0.4867, 0.4984, 0.5101, 0.5219, ! superlayer 5 . 0.5269, 0.5387, 0.5505, 0.5623, 0.5741, ! superlayer 6 . 0.5791, 0.5909, 0.6027, 0.6145, 0.6262, ! superlayer 7 . 0.6370, 0.6489, 0.6606, 0.6724, 0.6842, ! superlayer 8 . 0.6892, 0.7010, 0.7128, 0.7246, 0.7364, ! superlayer 9 . 0.7414, 0.7532, 0.7650, 0.7768, 0.7900, ! superlayer 10 . 0.8090, 47*0.0 / ! outer radius of plate c ihole = 1 -- region has holes data ihole / 0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1, . 0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1, 49*0 / c layer = drift chamber layer number data layer / 0, 1, 2, 3, 4, 0, 5, 6, 7, 8, . 0, 9, 10, 11, 12, 0, 13, 14, 15, 16, . 0, 17, 18, 19, 20, 0, 21, 22, 23, 24, . 0, 25, 26, 27, 28, 0, 29, 30, 31, 32, . 0, 33, 34, 35, 36, 0, 37, 38, 39, 40, . 49*0 / c desired wire tensions in grams (sense, clear, field) data twire / 34., 86., 182. / c modulus of wire (Pa) data Ewire / 3.55e11, 2*7.0e10 / c diameter of wires (m) data dwire / 0.00002, 0.00008, 0.00012 / c diameter of holes (m) data dhole / 0.0045, 2*0.0025 / c cells per layer data ncelltot / 0, 4*96, 0, 4*112, 0, 4*128, 0, 4*144, 0, 4*176, . 0, 4*192, 0, 4*208, 0, 4*224, 0, 4*240, 0, 4*256, . 49*0 / c number of wires per cell (sense, clear, field) data loadfac / . 3*0, 1,2,3, 1,0,2, 1,0,2, 1,2,2, ! 1 extra field wire . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, ! in first and last layer . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,3, . 147*0/ pi = 4.*atan(1.) c set up wire parameters do i = 1,3 fwire(i) = twire(i)*9.8/1000. ! tension in N EAwire(i) = Ewire(i)*pi*dwire(i)**2/4. nwire(i) = 0 enddo ! i zlen = 2.764 ! wire length (m) c setup for plate divided into 51 regions n = 51 t1r = 0.024 ! thickness of rear endplate rstep = 0.47 ! radius of step t1f = 0.014 ! thickness for r < rstep t2f = 0.014 ! thickness for r > rstep Eal = 7.0e10 ! Young's modulus for aluminum anu = 0.3 ! Poissons' ratio c values chosen to match Hodgson's results holefac = 1./1.5 ! modulus reduction factor for holes stressfac = 4.1 ! stress concentration factor around holes aload = 0. do i = 1,51 rmid(i) = 0.5*(bi(i+1) + bi(i)) area(i) = pi*(bi(i+1)**2 - bi(i)**2) E(i) = Eal anui(i) = anu t(1,i) = t1f ! front endplate if (bi(i) .gt. rstep) t(1,i) = t2f t(2,i) = t1r ! rear endplate qi(i) = 0. if (ihole(i). eq. 1) then stresscon(i) = stressfac awire = 0. do j = 1,3 loadfac(j,i) = loadfac(j,i)*ncelltot(i) ! convert to wires per layer awire = awire + loadfac(j,i)*pi*dhole(j)**2/4. nwire(j) = nwire(j) + loadfac(j,i) oswire(j,i) = twire(j) ! oswire stores wire tension in grams temp = oswire(j,i)*loadfac(j,i)/1000. ! wire tension in N aload = aload + temp ! load in kg qi(i) = qi(i) + temp*9.8/area(i) ! pressure in Pa c write(1,9898) j, i, loadfac(j,i), oswire(j,i), temp 9898 format (' j, i',2i3,' loadfac, oswire, temp ',i5,2e11.3) enddo ! j c holes(i) = 1. - awire/area(i) holes(i) = holefac else holes(i) = 1. stresscon(i) = 1. endif enddo ! i write (1,1) zlen, n 1 format (' Stress analysis of a chamber ',f5.3,' m long, with ', . i3,' regions'/) write (1,2) aload, nwire 2 format (' wire load = ',f5.0,' kg for ', i5,' sense, ',i5, . ' clearing, and ',i5,' field wires'// . ' Parameters for the regions:'// . ' region, b(cm), P(kPa), E(Gpa), nu, holefac, stress tf(mm),' . ' Df, tr, Dr') do i = 1,n do j = 1,2 Di(j,i) = E(i)*holes(i)*t(j,i)**3/12./(1. - anui(i)**2) ! plate constant enddo ! j write (1,11) i, bi(i)*100., qi(i)/1000., . E(i)*1.e-9, anui(i), holes(i), stresscon(i), . t(1,i)*1000., Di(1,i), t(2,i)*1000., Di(2,i) 11 format (i4,f10.1,f8.2,f8.2,f6.2,f7.3,f8.2,f7.1,f7.0,f7.1,f7.0) 10 format ( . ' region',t30,i3/ . ' inner radius',t30,1pg11.3,' m'/ . ' thickness, b < r < c',t30,g11.3,' m'/ . ' pressure',t30,g11.3,' Pascals'/ . ' Youngs modulus',t30,1pg11.3,' Pascals'/ . ' Poissons ratio',t30,g11.3/ . ' Modulus loss to holes',t30,g11.3/ . ' Stress conc. at holes',t30,g11.3/ . ' D = plate constant',t30,g11.3/) enddo ! i write (1,12) bi(n+1) 12 format ( . ' outer radius',t30,g11.3,' m'/) c front endplate c build the transformation from inner to outer regions call build(amat,cvec,1) c solve for values at inner radius, case 2c call case2c(vbeg,amat,cvec) c print final values, and store final displacements in dzf call output(1,vbeg,0,1,dzf) c rear endplate call build(amat,cvec,2) call case2c(vbeg,amat,cvec) call output(1,vbeg,0,2,dzf) c string one region at a time, starting from the inside. write(1,105) 105 format (/' Results for wire overstressing'// . ' layer dz(desired) dz(strung) dz(final) ', . 'tension (20, 80, 120 micron)') c initialize array dznow do j = 1,n do k = 1,2 dznow(k,j) = 0. enddo ! k enddo ! j do i = 1,n if (ihole(i) .eq. 1) then ! iterate only on regions that have wires c store deflection of plate at layer i at moment of stringing dzs(1,i) = dznow(1,i) dzs(2,i) = dznow(2,i) c calculate the overtensioning needed for layer i (in grams) do k = 1,3 oswire(k,i) = . (fwire(k)*(zlen - dznow(1,i) - dznow(2,i)) + Eawire(k)* . (dzf(1,i) + dzf(2,i) - dznow(1,i) - dznow(2,i)))/ . (zlen - dzf(1,i) - dzf(2,i))/.0098 enddo ! k write(1,1234) i, (oswire(k,i),k=1,3) 1234 format (' layer ',i3,' overtension ',3e11.3) c calculate new wire load, taking into account relaxation do j = 1,n qi(j) = 0. ! wire load only for regions inside region i if (j .le. i .and. ihole(j) .eq. 1) then do k = 1,3 tense = oswire(k,j)*.0098 ! convert to N relax = tense - (tense + Eawire(k))* . (dznow(1,j) + dznow(2,j) - dzs(1,j) - dzs(2,j))/ . (zlen - dzs(1,j) - dzs(2,j)) qi(j) = qi(j) + relax*loadfac(k,j)/area(j) enddo ! k write(1,1235) j, qi(j) 1235 format (' layer ',i3,' pressure ',e11.3) endif ! j, ihole enddo ! j c front endplate call build(amat,cvec,1) call case2c(vbeg,amat,cvec) call output(1,vbeg,i,1,dznow) ! store displacement of layer i at time of stringing c rear endplate call build(amat,cvec,2) call case2c(vbeg,amat,cvec) call output(1,vbeg,i,2,dznow) endif ! ihole enddo ! i do i = 1,n if (ihole(i) .eq. 1) then write(1,110) layer(i), (dzf(j,i)*1000.,j=1,2), . (dzs(j,i)*1000.,j=1,2), (dznow(j,i)*1000.,j=1,2), . (oswire(j,i),j=1,3) 110 format (i3,6f8.4,3f7.2) c output on unit 2 for TeX table write(2,111) layer(i), (dzf(j,i)*1000.,j=1,2), . (loadfac(j,i),oswire(j,i),j=1,3) 111 format (i3,' & ',f5.3,' & ',f5.3,' & ',i3,' & ',f5.1,' & ', . i3,' & ',f5.1,' & ',i3,' & ',f5.1,' \s') endif ! ihole enddo ! i stop end subroutine build(amat,cvec,ifr) c build transformation between inner and outer radii c vec_j(1) -- y at inner edge of region j c vec_j(2) -- theta c vec_j(3) -- radial bending moment, M_r c vec_j(4) -- shear force density, Q c form of transformation of vec_j from the inner edge of region j to c vec_i at the outer edge of region j: c vec_i(k) = sum_l amat_j(k,l)*vec_j(l) + cvec_j(k), c where amat_j and cvec_j are obtained from routine smatrix c we also want to transform from the inner edge of region 1 to the c outer edge of region j: c vec_i(k) = sum_l amat_i(k,l)*vec_j(l) + cvec_i(k) c therefore c cvec_i(k) = sum_l amat_j(k,l)*cvec_i-1(l) + cvec_j(k), c amat_i(k,l) = sum_m amat_j(k,m)*amat_i-1(m,l). c I store cvec_i-1 in cvec, and amat_i-1 in amat c ifr = 1 -- front, = 2 -- rear implicit real*8 (a-h), (o-z) common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99), . holes(99), stresscon(99), n, rmid(99), ihole(99), . layer(99) dimension amat(4,4), cvec(4), amatj(4,4), cvecj(4), amati(4,4), . cveci(4) do k = 1,4 cvec(k) = 0. do l = 1,4 amat(k,l) = 0. if (k .eq. l) amat(k,l) = 1. enddo ! l enddo ! k do j = 1,n ! get coeffs at outer edge of region call smatrix(j,bi(j+1),amatj,cvecj,ifr) c write (1,93) j, bi(j+1), cvecj, ((amatj(k,l),l=1,4),k=1,4) 93 format (' region ',i3,' outer radius ',f8.4/ . ' cvec',t10,4e11.3/' amat ',(t10,4e11.3)) do k = 1,4 cveci(k) = cvecj(k) do l = 1,4 cveci(k) = cveci(k) + amatj(k,l)*cvec(l) amati(k,l) = 0. do m = 1,4 amati(k,l) = amati(k,l) + amatj(k,m)*amat(m,l) enddo ! m enddo ! l enddo ! k do k = 1,4 cvec(k) = cveci(k) do l = 1,4 amat(k,l) = amati(k,l) enddo ! l enddo ! k enddo ! j return end subroutine case2c(vbeg,amat,cvec) implicit real*8 (a-h), (o-z) common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99), . holes(99), stresscon(99), n, rmid(99), ihole(99), . layer(99) dimension amat(4,4), cvec(4), vbeg(4), vfin(4) c apply B.C.'s at inner and outer radii -- depending on case c case 2c, simply supported at both a and b c height = 0 at inner and outer radii vbeg(1) = 0. vfin(1) = 0. c bend moment = 0 at inner and outer radii vbeg(3) = 0. vfin(3) = 0. c solve for vbeg(2) and veg(4) using c vfin(1) = 0 = amat(1,2)*vbeg(2) + amat(1,4)*vbeg(4) + cvec(1) c vfin(3) = 0 = amat(3,2)*vbeg(2) + amat(3,4)*vbeg(4) + cvec(3) denom = amat(1,2)*amat(3,4) - amat(3,2)*amat(1,4) vbeg(2) = -(cvec(1)*amat(3,4) - cvec(3)*amat(1,4))/denom vbeg(4) = -(amat(1,2)*cvec(3) - amat(3,2)*cvec(1))/denom c write (1,92) amat(1,2), amat(1,4), cvec(1), c . amat(3,2), amat(3,4), cvec(3) 92 format (' M12 ',e11.3,' M14 ',e11.3,' C1 ',e11.3/ . ' M32 ',e11.3,' M34 ',e11.3,' C3 ',e11.3/) return end subroutine output(icase,vbeg,iout,ifr,dz) implicit real*8 (a-h), (o-z) common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99), . holes(99), stresscon(99), n, rmid(99), ihole(99), . layer(99) dimension vbeg(4), vbegi(4,99), vfin(4), dz(2,99) character cases(2) / 'c', 'd' / pi = 4.*atan(1.) if (iout .eq. 0 .and. ifr. eq. 1) write(1,2) 2 format (' Results for Front Endplate') if (iout .eq. 0 .and. ifr. eq. 2) write(1,3) 3 format (' Results for Rear Endplate') if (iout .eq. 0) write (1,1) cases(icase) 1 format (/' case 2',a1// . ' Values at inner edges of regions' // . ' region, y(mm), theta(rad), M_r, Load(kg)') c find values at beginning of each region do j = 1,4 vbegi(j,1) = vbeg(j) enddo ! j do i = 1,n+1 if (i. gt. 1) . call values(vbegi(1,i),vbegi(1,i-1),bi(i),i-1,ifr) c write (1,91) (vbegi(j,i-1),j=1,4), (vbegi(j,i),j=1,4) 91 format (' vbeg ',4e11.3/' vfin ',4e11.3) alb = 2.*pi*bi(i)*vbegi(4,i)/9.8 if (iout .eq. 0) write (1,11) i, vbegi(1,i)*1000., vbegi(2,i), . vbegi(3,i), alb 11 format (i4,f10.3,f9.4,f9.1,f8.0) enddo ! i c print values at middles of regions that have wires if (iout .eq. 0) write (1,29) 29 format (/' Values at middle of regions that have wires'/) if (iout. eq. 0) write (1,21) 21 format(' layer, r(cm), y(mm), theta(rad), M_r(N-m/m), M_t,' . ' Q(N), sigma_r(MPa), sigma_t') do i = 1,n if (ihole(i). eq. 1) then r = rmid(i) call values(vfin,vbegi(1,i),r,i,ifr) c store displacement; change sign!!!! dz(ifr,i) = -vfin(1) amt = vfin(2)*Di(ifr,i)*(1. - anui(i)**2)/r + . anui(i)*vfin(3) sigmar = 6.*vfin(3)/t(ifr,i)**2*1.e-6*stresscon(i) ! include stress sigmat = 6.*amt/t(ifr,i)**2*1.e-6*stresscon(i) ! concentration factor if (iout .eq. 0) write (1,30) layer(i), r*100., . vfin(1)*1000., vfin(2), . vfin(3), amt, vfin(4), sigmar, sigmat endif ! ihole enddo ! i return 9999 continue c print values for 100 radii do i = 1,101 r = bi(1) + 0.01*(i- 1)*(bi(n+1) - bi(1)) ! find which region we are in j = 0 100 j = j + 1 if (r .gt. bi(j+1)) go to 100 if (j .gt. n) j = n call values(vfin,vbegi(1,j),r,j,ifr) amt = vfin(2)*Di(ifr,j)*(1. - anui(j)**2)/r + anui(j)*vfin(3) sigmar = 6.*vfin(3)/t(ifr,j)**2*1.e-6*stresscon(j) ! include stress sigmat = 6.*amt/t(ifr,j)**2*1.e-6*stresscon(j) ! concentration factor write (1,30) j, r*100., vfin(1)*1000., vfin(2), vfin(3), . amt, vfin(4), sigmar, sigmat 30 format (1x,i3, f7.2,',',f8.3,',',f9.4,',',f9.1,',', . f9.1',',f9.0,',',f9.2,',',f9.2) c 30 format (i3,8e9.1) enddo ! i return end subroutine values(vfin,vbeg,r,i,ifr) implicit real*8 (a-h), (o-z) dimension vfin(4), vbeg(4), amat(4,4), cvec(4) call smatrix(i,r,amat,cvec,ifr) do j = 1,4 vfin(j) = cvec(j) do k = 1,4 vfin(j) = vfin(j) + amat(j,k)*vbeg(k) enddo !k enddo ! j return end subroutine smatrix(i,r,amat,cvec,ifr) implicit real*8 (a-h), (o-z) common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99), . holes(99), stresscon(99), n, rmid(99), ihole(99), . layer(99) common/nu/anu dimension amat(4,4), cvec(4) anu = anui(i) b = bi(i) D = Di(ifr,i) q = qi(i) amat(1,1) = 1. amat(1,2) = r*f1(r,b) amat(1,3) = r**2/D*f2(r,b) amat(1,4) = r**3/D*f3(r,b) amat(2,1) = 0. amat(2,2) = f4(r,b) amat(2,3) = r/D*f5(r,b) amat(2,4) = r**2/D*f6(r,b) amat(3,1) = 0. amat(3,2) = D/r*f7(r,b) amat(3,3) = f8(r,b) amat(3,4) = r*f9(r,b) amat(4,1) = 0. amat(4,2) = 0. amat(4,3) = 0. amat(4,4) = b/r cvec(1) = -q*r**4/D*f11(r,b) cvec(2) = -q*r**3/D*f14(r,b) cvec(3) = -q*r**2*f17(r,b) cvec(4) = -q*(r**2 - b**2)/2./r return end real*8 function f1(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f1 = 0.5*(1. + anu)*b/r*log(r/b) + . 0.25*(1. - anu)*(r/b - b/r) return end real*8 function f2(r,b) implicit real*8 (a-h), (o-z) f2 = 0.25*(1. - (b/r)**2*(1. + 2.*log(r/b))) return end real*8 function f3(r,b) implicit real*8 (a-h), (o-z) f3 = 0.25*b/r*((1. + (b/r)**2)*log(r/b) + (b/r)**2 - 1.) return end real*8 function f4(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f4 = 0.5*((1. + anu)*b/r + (1. - anu)*r/b) return end real*8 function f5(r,b) implicit real*8 (a-h), (o-z) f5 = 0.5*(1. - (b/r)**2) return end real*8 function f6(r,b) implicit real*8 (a-h), (o-z) f6 = 0.25*b/r*((b/r)**2 - 1. + 2.*log(r/b)) return end real*8 function f7(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f7 = 0.5*(1. - anu**2)*(r/b - b/r) return end real*8 function f8(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f8 = 0.5*(1. + anu + (1. - anu)*(b/r)**2) return end real*8 function f9(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f9 = b/r*(0.5*(1. + anu)*log(r/b) + 0.25*(1. - anu)* . (1. - (b/r)**2)) return end real*8 function f11(r,b) implicit real*8 (a-h), (o-z) f11 = (1. + 4.*(b/r)**2 - 5.*(b/r)**4 - 4.*(b/r)**2* . (2. + (b/r)**2)*log(r/b))/64. return end real*8 function f14(r,b) implicit real*8 (a-h), (o-z) f14 = (1. - (b/r)**4 - 4.*(b/r)**2*log(r/b))/16. return end real*8 function f17(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f17 = 0.25*(1. - 0.25*(1. - anu)*(1. - (b/r)**4) - . (b/r)**2*(1. + (1. + anu)*log(r/b))) return end