采用该子程序计算

51465

      subroutine flux(f,temflu,mibody,time)
      include 'C:\MSC.Software\Marc\2013.0.0\marc2013\common\implicit'
      include 'C:\MSC.Software\Marc\2013.0.0\marc2013\common\bclabel'
      dimension f(2),mibody(*),temflu(*)
c* * * * * *
c
c     user subroutine for non-uniform flux input.
c
c     f(1)         flux value (to be defined)
c     f(2)         derivative of flux with respect to temperature
c                  (to be defined; optional, this might improve the
c                   convergence behavior)
c
c     temflu(1)    estimated temperature
c     temflu(2)    previous volumetric flux
c     temflu(3)    temperature at beginning of increment
c     temflu(4,5,6)integration point coordinates
c     mibody(1)    element number
c     mibody(2)    flux type
c     mibody(3)    integration point number
c     mibody(4)    flux index
c     mibody(5)    not used
c     mibody(6)    =1 : heat transfer
c                  =2 : joule
c                  =3 : bearing
c                  =4 : electrostatic
c                  =5 : magnetostatic
c                  =6 : acoustic
c     mibody(8)    layer number for heat transfer shells elements
c                  and volume flux
c     time         time
c
c* * * * * *
      real a,Qs,rs,pi,Qv,H,b,rv,Q,aa,x0,y0,z0
c     Qs为面热源功率, a为面热源能量集中系数，rs为面热源作用范围
c     Qv为体热源功率，H为体热源深度，b为体热源能量衰减系数，rv为体热源有效作用半径
c     Q为热源功率，aa为热源有效吸收系数,(x0,y0,z0)为当前热源中心位置
        Q=7140765;        
        aa=0.8;        Qs=Q*aa*0.806;        Qv=Q-Qs
        a=3;        rs=6        
        H=2.78;        b=1;        rv=2.5
        x0=50;        y0=100/60*time;        z0=4
        pi=3.14        
        OPEN(UNIT=10,FILE='text.txt')
      WRITE(10,*) MIBODY(2)


        if (mibody(2).eq.3) then
            d1=(temflu(4)-x0)**2+(temflu(5)-y0)**2              
                if(d1.le.2.5**2)  then
                        f(1)=Qv/(pi*rv*rv*H)       //搅拌针体热源密度

                end if
        end if
        if (mibody(2).eq.6) then
        d1=(temflu(4)-x0)**2+(temflu(5)-y0)**2
                if(d1.gt.2.5**2.and.d1.le.6**2)  then
                        f(1)=3*Qs*sqrt(d1)/(2*pi*(rs-rv))   //轴肩面热源密度
                end if
        end if
        return
      end