subroutine plotv(v,s,sp,etot,eplas,ecreep,t,mel,nip,layer, * ndj,nsh,jpltcd) c c CONTACT: ap@marc.de c c c* * * * * * c c PURPOSE: c This version of plotv will calculate the magnetostatic energy and c coenergy by integrating the b-h relation and summing over the whole c model. It also computes the magnetostatic energy density (J/m^3) c and writes it to the post file. c c energy and coenergy are defined as follows: c energy = Integral (HdB) ds (i.e. area of H vs. B curve c integrated over the whole model). c coenergy = Integral (H*B - HdB) ds (i.e. complementary area of c H vs. B curve integrated over the whole model). c c For 2D it will also write the magnetic potential (planar: A and c axisymmteric: r*A) to the post file. With this one can display c the field lines of B (the magnetic induction). c c USAGE: c The MARC input file requires post code -1 for the energy c density, total energy and coenergy and -2 for the magnetic c potential in the post option. c c NOTE: c This routine will work for 2D planar, axisymmetric and 3D models. c However, 3D will need further testing. c The routine will work for isotropic materials only. c c OUTPUT: c v variable to be written to post file c jpltcd=1: v = magnetic energy density c jpltcd=2: v = magnetic potential (2D only) c when planar: v = A c when axisymmetric: v = r*A c when 3D: v = 0 c c INPUT: c s b,h array c sp not used in magnetostatic analysis c etot not used in magnetostatic analysis c eplas not used in magnetostatic analysis c ecreep not used in magnetostatic analysis c t magnetic potential (2D only) c mel(1) user element number c mel(2) internal element number c nip integration point number c layer layer number c ndj (3) not used in magnetostatic analysis c nsh (3) not used in magnetostatic analysis c jpltcd post code (in absolute value) c c* * * * * * c implicit real*8 (a-h,o-z) dp dimension s(1),etot(1),eplas(1),ecreep(1),sp(1),mel(1) c include '../common/strnen' include '../common/arrays' include '../common/array4' include '../common/heat' include '../common/dimen' include '../common/elmcom' include '../common/matdat' include '../common/space' include '../common/lass' include '../common/concom' include '../common/maters' c c xip(1..3) integration point coordinates dimension xip(3) c c NPMAX maximum number of points in b-h relation parameter (NPMAX=200) dimension bhrel(2,NPMAX),bhsum(NPMAX) common/bh_data/bhrel,bhsum,energy,coener,ielid c bhrel array containg b-h relation for current material c bhrel(1,i)=magnetic field intensity at point i c bhrel(2,i)=magnetic induction at point i c bhsum integral of HdB at break points in b-h relation c bhsum(i)=integral under b-h relation up to point i c energy magnetostatic energy c coener magnetostatic coenergy c ielid element id as given by order of appearance in element loop c data ielid/0/ data energy,coener /0.0,0.0/ c c... return if layer number not equal to one if (layer.gt.1) return c c... compute area under b-h curve at breakpoints and store in bhsum if (nip.eq.1) then c... get b-h relation for this material from MARC internal storage do i1=1,noesl la1 = iesl + 2*(mats-1)*noslps +2*(i1-1) c... get h-value bhrel(1,i1) = vars(la1) c... get b-value bhrel(2,i1) = vars(la1+1) enddo c... integrate under curve up to the break points bhsum(1) = 0d0 do i1=2,noesl bhsum(i1) = bhsum(i1-1) + 0.5d0* * (bhrel(1,i1)+bhrel(1,i1-1))* * (bhrel(2,i1)-bhrel(2,i1-1)) enddo endif c c... compute norm of magnetic induction if (ncrd.eq.2) bn = dsqrt(s(1)*s(1)+s(2)*s(2)) if (ncrd.eq.3) bn = dsqrt(s(1)*s(1)+s(2)*s(2)+s(3)*s(3)) c... compute norm of magnetic field intensity if (ncrd.eq.2) hn = dsqrt(s(3)*s(3)+s(4)*s(4)) if (ncrd.eq.3) hn = dsqrt(s(4)*s(4)+s(5)*s(5)+s(6)*s(6)) c c... energy is computed for post code -1 if (jpltcd.eq.1) then c c... increment the element count when processing integration point 1 if (nip.eq.1) ielid = ielid+1 c c... perform integration under b-h curve for current b, h values if (bn.gt.bhrel(2,noesl)) then write(6,1000) mel(1),nip,bn,hn 1000 format('warning: b-h values out of specified range: el' * ,i5,' ip',i5,' b=',e13.5,' h=',e13.5) i1 = noesl else i1=0 10 continue i1 = i1+1 if (bn.gt.bhrel(2,i1)) goto 10 i1 = i1-1 endif c... compute energy density [J/m^3] c... i1=0 denotes a material without a b-h relation if (i1.eq.0) then v = 0.5d0*hn*bn else v = bhsum(i1) + 0.5d0*(hn+bhrel(1,i1))*(bn-bhrel(2,i1)) endif c c... data base offset (here use internal element number) lofr = (n-1)*nelstr c c... dv = the jacobian (the integration point volume factor) la1 = ijacoc + lofr + nn - 1 dv = vars(la1) c c... sum the energy and coenergy contributions energy = energy + v*dv coener = coener + (bn*hn-v)*dv c c... write energy and coenergy to output file when processing c last element and last integration point if (ielid.eq.numel.and.nip.eq.intel) then ielid = 0 write(6,5000) write(6,2000) energy write(6,3000) coener write(6,5000) 2000 format(30x,'Energy ',e13.5) 3000 format(30x,'Coenergy ',e13.5) 5000 format(30x,'**********************') energy = 0d0 coener = 0d0 endif else if (jpltcd.eq.2) then c... magnetic potential is computed for post code -2 c... for the axisymmetric case multiply with radius c... in 2D we now get to see the field lines of B if (lclass.eq.7) then c... for axisymmetric case get integration point coordinates la1 = icrxpt + (n-1)*nelstr + (nn-1)*ncrd do i1=1,ncrd xip(i1) = vars(la1+i1-1) enddo if (ncrd.eq.2) xip(3) = 0d0 v = xip(2)*t else if (ncrd.eq.2) then c... 2D planar case v = t else c... 3D case v = 0d0 endif endif c return c end