FrontISTR【ヒートソースの導入】

FrontISTRのHeatにユーザー定義でヒートソースを導入しようとした取り組みです。なお、要素形状は361（六面体１次）要素限定です。

@analysis/heat/heat_mat_ass_bc_DFLUX

    !
    if( .true. ) then

      allocate( Bbak(size(hecMAT%B)) )
      vol =0.d0; Bbak=0.d0

      do icel = 1, hecMESH%n_elem

        ic_type = hecMESH%elem_type(icel)
        nn = hecmw_get_max_node(ic_type)
        is = hecMESH%elem_node_index(icel-1)

        allocate( temperature_on_node(nn) )

        do j = 1, nn
          nodLOCAL(j) = hecMESH%elem_node_item(is+j)
          xx(j) = hecMESH%node( 3*nodLOCAL(j)-2 )
          yy(j) = hecMESH%node( 3*nodLOCAL(j)-1 )
          zz(j) = hecMESH%node( 3*nodLOCAL(j)   )
          temperature_on_node(j) = fstrHEAT%TEMP( nodLOCAL(j) )
        enddo
        !
        temperature_on_gauss = get_temperature_361_on_gauss( &
        & ic_type ,nn, temperature_on_node )
        !
        ng = size( temperature_on_gauss )
        !
        val = 0.0d0
        !
        do j = 1, ng
            !
            val = val + heatsourcevalue( fstrHEAT, ctime, dtime, temperature_on_gauss(j) )
            !
        enddo
        !
        deallocate( temperature_on_node, temperature_on_gauss )

        vol = vol + VOLUME_C3(ic_type,NN,XX,YY,ZZ)

        if( ic_type.eq.361 ) then
          call heat_DFLUX_361(nn,xx,yy,zz,0,val,vect)
        else
          write(*,*) '###ERROR### : Element type not supported for heat analysis'
          write(*,*) ' ic_type = ', ic_type
          call hecmw_abort(hecmw_comm_get_comm())
        endif

        do j = 1, nn
          Bbak( nodLOCAL(j) ) = Bbak( nodLOCAL(j) ) - vect(j)
        enddo

      enddo

      if( vol>0 ) then
        Bbak = Bbak/vol
        hecMAT%B = hecMAT%B + Bbak
      endif

      deallocate( Bbak )

    endif

heat_mat_ass_bc_RADIATEなども参考に作りました。ヒートソースをガウス点上で設定し、heat_DFLUX_361を使って節点温度に返します。また、ガウス点上の温度に依存するヒートソースを設定したかったので、関数 get_temperature_361_on_gaussを作りました（後述）。

関数 heatsourcevalueではヒートソースの値を計算します。なお、温度（時刻）などに依存するヒートソースを導入する場合、そのタイムステップにおける内部反復で毎回値を入れる必要があるようなので注意です（heatsourcevalue内は省略）。

なお、heatsourcevalueを作るための、fstrHEATへの変数の追加は前記事参照です。 FrontISTR【伝熱解析での変数の追加】 - FEMとFortranが好きな人の技術のメモ

@/lib/heat_LIB_DFLUX

  function get_temperature_361_on_gauss(etype,nn,temperature) result(temp_gauss)

    use mMechGauss
    use elementInfo

    implicit none

    real(kind=kreal),ALLOCATABLE :: temp_gauss(:)
    integer(kind=kint), intent(in) :: etype                  !< element type
    integer(kind=kint), intent(in) :: nn                     !< number of elemental nodes
    real(kind=kreal), intent(in)  :: temperature(nn) !< temperature
    real(kind=kreal) :: naturalCoord(3)

    real(kind=kreal) :: func(nn), temp_i
    integer(kind=kint) :: LX

    ALLOCATE( temp_gauss( NumOfQuadPoints(etype) ), source=0.0d0 )
    !
    do LX = 1, NumOfQuadPoints(etype)
      call getQuadPoint(etype, LX, naturalCoord)
      call getShapeFunc(etype, naturalCoord, func)
      temp_i = dot_product(func, temperature)
      temp_gauss(LX) = temp_i
    end do

  end function

サブルーチンheat_conductivity_C3（＠heat_LIB_CONDUCTIVITY）を参考に作りました。361要素で接点温度からガウス点上の温度を返す関数です。