Merge branch '13-higher-order-thermal-expansion' into 'development'

Resolve "higher-order thermal expansion"

Closes #13

See merge request damask/DAMASK!19

src/kinematics_thermal_expansion.f90

@@ -174,8 +174,12 @@ pure function kinematics_thermal_expansion_initialStrain(ipc, ip, el)
  offset = thermalMapping(homog)%p(ip,el)
  
  kinematics_thermal_expansion_initialStrain = &
-   (temperature(homog)%p(offset) - lattice_referenceTemperature(phase)) * &
+   (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * &
-   lattice_thermalExpansion33(1:3,1:3,phase)
+   lattice_thermalExpansion33(1:3,1:3,1,phase) + &                                                  ! constant  coefficient
+   (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * &
+   lattice_thermalExpansion33(1:3,1:3,2,phase) + &                                                  ! linear    coefficient
+   (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * &
+   lattice_thermalExpansion33(1:3,1:3,3,phase)                                                      ! quadratic coefficient
   
 end function kinematics_thermal_expansion_initialStrain
 
@@ -215,9 +219,16 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc,
  TDot = temperatureRate(homog)%p(offset)
  TRef = lattice_referenceTemperature(phase)
  
- Li = TDot* &
+ Li = TDot * ( &
-      lattice_thermalExpansion33(1:3,1:3,phase)/ &
+               lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 &                           ! constant  coefficient
-      (1.0_pReal + lattice_thermalExpansion33(1:3,1:3,phase)*(T - TRef))
+             + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 &                           ! linear    coefficient
+             + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 &                           ! quadratic coefficient
+             ) / &
+      (1.0_pReal \
+            + lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. &
+            + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. &
+            + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. &
+      )
  dLi_dTstar3333 = 0.0_pReal 
   
 end subroutine kinematics_thermal_expansion_LiAndItsTangent

src/lattice.f90

@@ -1082,9 +1082,10 @@ module lattice
    lattice_nu, &
    lattice_trans_mu, &
    lattice_trans_nu
+ real(pReal),                              dimension(:,:,:,:), allocatable, public, protected :: &  ! with higher-order parameters (e.g. temperature-dependent)
+   lattice_thermalExpansion33
  real(pReal),                              dimension(:,:,:),   allocatable, public, protected :: &
    lattice_thermalConductivity33, &
-   lattice_thermalExpansion33, &
    lattice_damageDiffusion33, &
    lattice_vacancyfluxDiffusion33, &
    lattice_vacancyfluxMobility33, &
@@ -1380,8 +1381,8 @@ subroutine lattice_init
  allocate(lattice_C3333(3,3,3,3,Nphases),  source=0.0_pReal)
  allocate(lattice_trans_C66(6,6,Nphases),  source=0.0_pReal)
  allocate(lattice_trans_C3333(3,3,3,3,Nphases),  source=0.0_pReal)
+ allocate(lattice_thermalExpansion33     (3,3,3,Nphases), source=0.0_pReal)                       ! constant, linear, quadratic coefficients
  allocate(lattice_thermalConductivity33  (3,3,Nphases), source=0.0_pReal)
- allocate(lattice_thermalExpansion33     (3,3,Nphases), source=0.0_pReal)
  allocate(lattice_damageDiffusion33      (3,3,Nphases), source=0.0_pReal)
  allocate(lattice_vacancyfluxDiffusion33 (3,3,Nphases), source=0.0_pReal)
  allocate(lattice_vacancyfluxMobility33  (3,3,Nphases), source=0.0_pReal)
@@ -1545,11 +1546,17 @@ subroutine lattice_init
      case ('thermal_conductivity33')
        lattice_thermalConductivity33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
      case ('thermal_expansion11')
-       lattice_thermalExpansion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt)
+         do i = 2_pInt, min(4,chunkPos(1))                                                        ! read up to three parameters (constant, linear, quadratic with T)
+           lattice_thermalExpansion33(1,1,i-1_pInt,section) = IO_floatValue(line,chunkPos,i)
+         enddo
      case ('thermal_expansion22')
-       lattice_thermalExpansion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt)
+         do i = 2_pInt, min(4,chunkPos(1))                                                        ! read up to three parameters (constant, linear, quadratic with T)
+           lattice_thermalExpansion33(2,2,i-1_pInt,section) = IO_floatValue(line,chunkPos,i)
+         enddo
      case ('thermal_expansion33')
-       lattice_thermalExpansion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
+         do i = 2_pInt, min(4,chunkPos(1))                                                        ! read up to three parameters (constant, linear, quadratic with T)
+           lattice_thermalExpansion33(3,3,i-1_pInt,section) = IO_floatValue(line,chunkPos,i)
+         enddo
      case ('specific_heat')
        lattice_specificHeat(section) = IO_floatValue(line,chunkPos,2_pInt)
      case ('vacancyformationenergy')
@@ -1756,10 +1763,12 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc)
      end select
  end select
 
+ forall (i = 1_pInt:3_pInt) &
+   lattice_thermalExpansion33 (1:3,1:3,i,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
+                                                                         lattice_thermalExpansion33   (1:3,1:3,i,myPhase))
+
  lattice_thermalConductivity33  (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
                                                                          lattice_thermalConductivity33  (1:3,1:3,myPhase))
- lattice_thermalExpansion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
-                                                                    lattice_thermalExpansion33(1:3,1:3,myPhase))
  lattice_DamageDiffusion33      (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
                                                                          lattice_DamageDiffusion33      (1:3,1:3,myPhase))
  lattice_vacancyfluxDiffusion33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&