Merge branch '170-temperature-dependent-heat-capacity' into 'development'

Temperature-dependent heat capacity and thermal conductivity

See merge request damask/DAMASK!1028

Author: eisenlohr

PRIVATE

@@ -11,6 +11,6 @@ K_11: 12.86
 K_11,T: 1.618e-2
 
 C_p: 508.7
-C_p,T: 1.347e-1 # not implemented
+C_p,T: 1.347e-1
 
 T_ref: 293.15

examples/config/phase/thermal/X2CrNi18-10.yaml

@@ -11,6 +11,6 @@ K_11: 13.85
 K_11,T: 1.571e-2
 
 C_p: 497.7
-C_p,T: 1.327e-1 # not implemented
+C_p,T: 1.327e-1
 
 T_ref: 293.15

examples/config/phase/thermal/X2CrNiMo17-12-2.yaml

@@ -4,8 +4,8 @@
 submodule(phase) thermal
 
   type :: tThermalParameters
-    real(pREAL) :: C_p = 0.0_pREAL                                                                  !< heat capacity
-    real(pREAL), dimension(3,3) :: K = 0.0_pREAL                                                    !< thermal conductivity
+    type(tpolynomial) :: C_p                                                                  !< heat capacity
+    type(tpolynomial) :: K_11, K_33                                                           !< thermal conductivity
     character(len=pSTRLEN), allocatable, dimension(:) :: output
   end type tThermalParameters
 
@@ -105,22 +105,23 @@ module subroutine thermal_init(phases)
       thermal => phase%get_dict('thermal',defaultVal=emptyDict)
     end if
 
-    ! ToDo: temperature dependency of K and C_p
     if (thermal%length > 0) then
       thermal_active = .true.
 
       print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph)
       refs = config_listReferences(thermal,indent=3)
       if (len(refs) > 0) print'(/,1x,a)', refs
 
-      param(ph)%C_p = thermal%get_asReal('C_p')
-      param(ph)%K(1,1) = thermal%get_asReal('K_11')
-      if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asReal('K_33')
-      param(ph)%K = crystal_symmetrize_33(param(ph)%K,phase_lattice(ph))
+      associate(prm => param(ph))
+
+        prm%C_p = polynomial(thermal,'C_p','T')
+        prm%K_11 = polynomial(thermal,'K_11','T')
+
+        if (any(phase_lattice(ph) == ['hP','tI'])) prm%K_33 = polynomial(thermal,'K_33','T')
+
+      end associate
 
       ! sanity checks
-      if (    param(ph)%C_p <= 0.0_pREAL )  extmsg = trim(extmsg)//' C_p'
-      if (any(param(ph)%K   <  0.0_pREAL))  extmsg = trim(extmsg)//' K'
       if (    phase_rho(ph) <= 0.0_pREAL )  extmsg = trim(extmsg)//' rho'
       if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg))
 
@@ -223,9 +224,17 @@ module function phase_mu_T(co,ce) result(mu)
   integer, intent(in) :: co, ce
   real(pREAL) :: mu
 
+  real(pREAL) :: T
+
+
+  associate(ph => material_ID_phase(co,ce), &
+            en => material_entry_phase(co,ce))
+
+    T = current(ph)%T(en)
+    mu = phase_rho(ph) &
+       * param(ph)%C_p%at(T)
 
-  mu = phase_rho(material_ID_phase(co,ce)) &
-     * param(material_ID_phase(co,ce))%C_p
+  end associate
 
 end function phase_mu_T
 
@@ -238,8 +247,23 @@ module function phase_K_T(co,ce) result(K)
   integer, intent(in) :: co, ce
   real(pREAL), dimension(3,3) :: K
 
+  real(pREAL) :: T
+
+
+  associate(ph => material_ID_phase(co,ce), &
+            en => material_entry_phase(co,ce))
+
+    T = current(ph)%T(en)
+
+    K = 0.0_pREAL
+    K(1,1) = param(ph)%K_11%at(T)
+    if (any(phase_lattice(ph) == ['hP','tI'])) K(3,3) = param(ph)%K_33%at(T)
+
+    K = crystal_symmetrize_33(K,phase_lattice(ph))
+    K = crystallite_push33ToRef(co,ce,K)
+
+  end associate
 
-  K = crystallite_push33ToRef(co,ce,param(material_ID_phase(co,ce))%K)
 
 end function phase_K_T