List of tutorials

commentnote:Before you proceed

All tutorials are written assuming that you are reasonably familiar with MOOSE. If you find most of the tutorials difficult to follow, please refer to the official MOOSE website for learning resources.

Tutorial 3: Ferroelectric thin film

This tutorial (and others) covers the basic usage of the ferroelectric (FE) phase field method implemented in FERRET. This specific example focuses on the thin film problem for (PTO) or (BTO) where periodicity is enforced along and corresponds to the film/substrate interface plane normal. We consider the calculation to be at room temperature.

As in Tutorial 2, we use the fully-coupled problem with electrostrictive free energy density. The total free energy density for this simulation can be written as,

for the bulk, elastic, gradient, electrostrictive, and electrostatic free energies respectively. We choose a computational geometry as follows,

[Mesh]
  [gen]
    ############################################
    ##
    ##  Type and dimension of the mesh 
    ##
    ############################################

    type = GeneratedMeshGenerator
    dim = 3

    nx = 32
    ny = 32
    nz = 30

    xmin = -16.0
    xmax = 16.0
    ymin = -16.0
    ymax = 16.0
    zmin = -10.0
    zmax = 20.0

    #############################################
    ##
    ##  FE type/order (hexahedral, tetrahedral
    ##
    #############################################

    elem_type = HEX8
  []
  [./cnode]
    input = gen

    ############################################
    ##
    ##   additional boundary sideset (one node) 
    ##   to zero one of the elastic displacement vectors 
    ##   vectors and eliminates rigid body translations 
    ##   from the degrees of freedom
    ##
    ##   NOTE: This must conform with the about
    ##         [Mesh] block settings
    ##
    ############################################

    type = ExtraNodesetGenerator
    coord = '-16.0 -16.0 -10.0'
    new_boundary = 100
  [../]

  [subdomains]
    type = SubdomainBoundingBoxGenerator
    input = cnode
    bottom_left = '-16.0 -16.0 -10.0'
    block_id = 1
    top_right = '16.0 16.0 0'
    location = INSIDE
  []
  [film_interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = subdomains
    primary_block = 0
    paired_block = 1
    new_boundary = 52
  []
  [film_surface]
    type = SideSetsFromNormalsGenerator
    input = film_interface
    normals = '0  0  1'
    fixed_normal = true
    new_boundary = '107'
  []
  [substrate_bottom]
    type = SideSetsFromNormalsGenerator
    input = film_surface
    normals = '0  0  -1'
    fixed_normal = true
    new_boundary = '108'
  []
[]

where and are chosen accordingly such that the mesh spacing nm. In general, the geometry defined in the 'Mesh' block never carries units. The length scale is introduced through Materials, Kernels, or other MOOSE objects. For this problem, the length scale is introduced through the units in the Materials objects that connect to the Kernels. Note that this discretization is around the upper bound for these types of calculations since the thickness of the PTO wall is about . We use SubdomainBoundingBoxGenerator and SideSetsBetweenSubdomainsGenerator objects from MOOSE to split the rectilinear computational volume into two blocks where the film sits on top of the substrate. For this example input file, we let the FE film have a thickness of 20 nm. This is a canonical problem in FE phase field problems (see Li et al. (2001)) as it allows one to predict the domain topology as a function of temperature and film thickness.

commentnote:Important!

We should note that the epitaxial strain coupling is not yet implemented in FERRET but this capability will be introduced in the future. As such, we just consider a free standing stress-free film.

The Variables block sets the ICs for where the amplitude is set to be random fluctuations near 0

[Variables]

  #################################
  ##
  ##  Variable definitions
  ##    P, u, phi, e^global_ij  
  ##  and their initial conditions
  ##
  #################################

  [./global_strain]
    order = SIXTH
    family = SCALAR
  [../]
  [./polar_x]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = RandomIC
      min = -1e-2
      max = 1e-2
    [../]
    block = '0'
  [../]
  [./polar_y]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = RandomIC
      min = -1e-2
      max = 1e-2
    [../]
    block = '0'
  [../]
  [./polar_z]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = RandomIC
      min = -1e-2
      max = 1e-2
    [../]
    block = '0'
  [../]

  [./potential_E_int]
    order = FIRST
    family = LAGRANGE
    block = '0 1'
  [../]

  [./u_x]
    order = FIRST
    family = LAGRANGE
    block = '0 1'
  [../]
  [./u_y]
    order = FIRST
    family = LAGRANGE
    block = '0 1'
  [../]
  [./u_z]
    order = FIRST
    family = LAGRANGE
    block = '0 1'
  [../]
[]

Other variables , and are also solved for. This tutorial problem evolves the time-dependent Landau-Ginzburg-Devonshire equation (TDLGD),

to find the ground state with (arbitrary time scale). We also solve (at every time step) the conditions for electrostatic (Poisson equation) and mechanical equilibrium (stress divergence),

The variational derivatives of the total free energy density yield residual and jacobian contributions that are computed within the Kernels block,

[Kernels]

  ###############################################
  ##
  ## Physical Kernel operators
  ## to enforce TDLGD evolution 
  ##
  ###############################################

  #Elastic problem
  [./TensorMechanics]
    use_displaced_mesh = false
    eigenstrain_names = eigenstrain
  [../]

  [./bed_x]
    type = BulkEnergyDerivativeEighth
    variable = polar_x
    component = 0
    block = '0'
  [../]
  [./bed_y]
    type = BulkEnergyDerivativeEighth
    variable = polar_y
    component = 1
    block = '0'
  [../]
  [./bed_z]
    type = BulkEnergyDerivativeEighth
    variable = polar_z
    component = 2
    block = '0'
  [../]

  [./walled_x]
    type = WallEnergyDerivative
    variable = polar_x
    component = 0
    block = '0'
  [../]
  [./walled_y]
    type = WallEnergyDerivative
    variable = polar_y
    component = 1
    block = '0'
  [../]
  [./walled_z]
    type = WallEnergyDerivative
    variable = polar_z
    component = 2
    block = '0'
  [../]

  [./walled2_x]
    type = Wall2EnergyDerivative
    variable = polar_x
    component = 0
    block = '0'
  [../]
  [./walled2_y]
    type = Wall2EnergyDerivative
    variable = polar_y
    component = 1
    block = '0'
  [../]
  [./walled2_z]
    type = Wall2EnergyDerivative
    variable = polar_z
    component = 2
    block = '0'
  [../]

  [./electrostr_ux]
    type = ElectrostrictiveCouplingDispDerivative
    variable = u_x
    component = 0
    block = '0'
  [../]
  [./electrostr_uy]
    type = ElectrostrictiveCouplingDispDerivative
    variable = u_y
    component = 1
    block = '0'
  [../]
  [./electrostr_uz]
    type = ElectrostrictiveCouplingDispDerivative
    variable = u_z
    component = 2
    block = '0'
  [../]

  [./electrostr_polar_coupled_x]
    type = ElectrostrictiveCouplingPolarDerivative
    variable = polar_x
    component = 0
    block = '0'
  [../]
  [./electrostr_polar_coupled_y]
    type = ElectrostrictiveCouplingPolarDerivative
    variable = polar_y
    component = 1
    block = '0'
  [../]
  [./electrostr_polar_coupled_z]
    type = ElectrostrictiveCouplingPolarDerivative
    variable = polar_z
    component = 2
    block = '0'
  [../]

  [./polar_x_electric_E]
    type = PolarElectricEStrong
    variable = potential_E_int
    block = '0'
  [../]
  [./FE_E_int]
    type = Electrostatics
    variable = potential_E_int
    block = '0 1'
  [../]

  [./polar_electric_px]
    type = PolarElectricPStrong
    variable = polar_x
    component = 0
    block = '0'
  [../]
  [./polar_electric_py]
    type = PolarElectricPStrong
    variable = polar_y
    component = 1
    block = '0'
  [../]
  [./polar_electric_pz]
    type = PolarElectricPStrong
    variable = polar_z
    component = 2
    block = '0'
  [../]

  [./polar_x_time]
    type = TimeDerivativeScaled
    variable = polar_x
    time_scale = 1.0
    block = '0'
  [../]
  [./polar_y_time]
    type = TimeDerivativeScaled
    variable = polar_y
    time_scale = 1.0
    block = '0'
  [../]
  [./polar_z_time]
    type = TimeDerivativeScaled
    variable = polar_z
    time_scale = 1.0
    block = '0'
  [../]

  [./u_x_time]
    type = TimeDerivativeScaled
    variable = u_x
    time_scale = 1.0
  [../]
  [./u_y_time]
    type = TimeDerivativeScaled
    variable = u_y
    time_scale = 1.0
  [../]
  [./u_z_time]
    type = TimeDerivativeScaled
    variable = u_z
    time_scale = 1.0
  [../]

[]

We refer the reader to the following hyperlinks,

TDLGD:

Poisson equation:

Mechanical equilibrium:

  • TensorMechanics (MOOSE Action) for

  • ElectrostrictiveCouplingDispDerivative

for the different objects. Next, we define the boundary conditions on , and ,

[BCs]
  [./Periodic]
    [./xy]
      auto_direction = 'x y'
      variable = 'u_x u_y u_z polar_x polar_y polar_z potential_E_int'
    [../]
  [../]

  [./boundary_interface_grounding]
    type = DirichletBC
    boundary = '52'
    variable = potential_E_int
    value = 0.0
  [../]

  # fix center point location
  [./centerfix_x]
    type = DirichletBC
    boundary = '108'
    variable = u_x
    value = 0
  [../]
  [./centerfix_y]
    type = DirichletBC
    boundary = '108'
    variable = u_y
    value = 0
  [../]
  [./centerfix_z]
    type = DirichletBC
    boundary = '108'
    variable = u_z
    value = 0
  [../]
[]

where and are the lateral directions of the FE film. We utilize the GlobalStrain system implemented in MOOSE to ensure periodicity of the strain tensor components along the periodic boundaries. This introduces a ScalarKernel,

with the computational volume. We find a set of global displacement vectors disp_x, disp_y, disp_z (see AuxKernels) such that the above condition is satisfied (see Biswas et al. (2020) for more description of the method). Note that the substrate has a boundary condition such that the elastic displacement vector goes to zero deep in the substrate which is handled by the DirichletBC. Some possible outputs of this problem with PTO Materials coefficients at room temperature,

[Materials]

  #################################################
  ##
  ## Landau coefficients from Li et al (2001)
  ##
  ##################################################

  [./Landau_P_FE]
    type = GenericConstantMaterial
    prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
    prop_values = '-0.1722883 -0.073 0.75 0.26 0.61 -3.67 0.0 0.0 0.0 0.0'
    block = '0'
  [../]

  [./Landau_P_substr]
    type = GenericConstantMaterial
    prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
    prop_values = '10.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0'
    block = '1'
  [../]

  [./Landau_G_FE]
    type = GenericConstantMaterial
    prop_names = 'G110 G11_G110 G12_G110 G44_G110 G44P_G110'
    prop_values = '0.173 0.6 0.0 0.3 0.3'
    block = '0'
  [../]

  [./mat_C_FE]
    type = GenericConstantMaterial
    prop_names = 'C11 C12 C44'
    prop_values = '175.0 79.4 111.1'
    block = '0'
  [../]
  [./mat_C_sub]
    type = GenericConstantMaterial
    prop_names = 'C11 C12 C44'
    prop_values = '220.0 34.4 161.1'
    block = '1'
  [../]

  ##################################################
  ##=
  ## NOTE: Sign convention in Ferret for the 
  ##        electrostrictive coeff. is multiplied by
  ##        an overall factor of (-1)
  ##
  ##################################################

  [./mat_Q]
    type = GenericConstantMaterial
    prop_names = 'Q11 Q12 Q44'
    prop_values = '-0.089 0.026 -0.03375'
    block = '0 1'
  [../]

  [./mat_q]
    type = GenericConstantMaterial
    prop_names = 'q11 q12 q44'
    prop_values = '-11.4 0.01438 -7.5'
  [../]

  [./eigen_strain]
    type = ComputeEigenstrain
    # eigen_base = 'exx exy exz eyx eyy eyz ezx ezy ezz'
    eigen_base = '1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0'
    eigenstrain_name = eigenstrain
    prefactor = 0.0
  [../]

  [./elasticity_tensor_1]
    type = ComputeElasticityTensor
    fill_method = symmetric9

    ###############################################
    ##
    ## symmetric9 fill_method is (default)
    ##     C11 C12 C13 C22 C23 C33 C44 C55 C66 
    ##
    ###############################################

    C_ijkl = '175.0 79.4 79.4 175.0 79.4 175.0 111.1 111.1 111.1'
  [../]
  [./strain_1]
    type = ComputeSmallStrain
    global_strain = global_strain
    eigenstrain_names = eigenstrain
  [../]

  [./stress_1]
    type = ComputeLinearElasticStress
  [../]

  [./global_strain]
    type = ComputeGlobalStrain
    scalar_global_strain = global_strain
    global_strain_uo = global_strain_uo
  [../]

  [./permitivitty_1]

    ###############################################
    ##
    ##  so-called background dielectric constant
    ##  (it encapsulates the motion of core electrons
    ##  at high frequency) = e_b*e_0 (here we use 
    ##  e_b = 10), see PRB. 74, 104014, (2006)
    ##
    ###############################################

    type = GenericConstantMaterial
    prop_names = 'permittivity'
    prop_values = '0.08854187'
  [../]
[]

are visualized below using ParaView Filters for colormaps of and white arrow Glyphs for the directors.

Figure 1: across the PTO film sample in 3D. The Glyph filter provides the arrows corresponding to the polar director.

If we use the BTO Materials coefficients at room temperature, we have the below output.

Figure 2: across the BTO film sample in 3D. The Glyph filter provides the arrows corresponding to the polar director.

In each of the calculations, the boundary condition on is open-circuit (free) at the surface of the thin film. This leads to a strong electrostatic action to form in-plane domains. Both input files are provided in the tutorial folder in the FERRET root directory.

References

  1. Sudipta Biswas, Daniel Schwen, and Jason D. Hales. Development of a finite element based strain periodicity implementation method. Finite Elements Anal. Design, 2020.[BibTeX]
  2. Y. L. Li, S. Y. Hu, Z. K. Liu, and L.-Q. Chen. Phase-field model of domain structures in ferroelectric thin films. Applied Physical Letters, 2001.[BibTeX]