Polarization induced charges

After the strain distribution has been calculated, the strain tensor in the crystal system can be obtained by the subroutine get_strain_tensor_cxyz for each point on material grid.
On this basis, one can provide the piezoelectric polarization vector for each of these points.
The pyroelectric polarization is implemented in a similar was. However, it is independent of strain.
This is done in MODULE input_block_four:

•  subroutine get_pyroelectricity(i,PolV)
!-----------------------------------------------------------------------------
! Provides for point number i on material grid
! the polarization vector in coordinates of simulation region.
! Units: C/m^2
! 1D: PolV(1)   = component of polarization vector in simulation direction
! 2D: PolV(1:2) = component of polarization vector (Px,Py)
! 3D: PolV(1:3) = component of polarization vector (Px,Py,Pz)
!-----------------------------------------------------------------------------

    vector = sp_polar_xyz(mnu,x)        --> get the vector in xyz system
1D: compV  = mask_domain_to_xyz(1)
2D: compV  = mask_domain_to_xyz(1:2)    --> mask_domain_to_xyz says which is the
                                            simulation_x, simulation_y axis:
                                            mask_domain_to_xyz(1) --> x_axis in simulation
                                                              (2) --> y_axis in simulation
3D: compV  = mask_domain_to_xyz(1:3)    --> mask_domain_to_xyz says which is the
                                            simulation_x, simulation_y, simulation_z axis
                                            mask_domain_to_xyz(1:3) = (1 2 3)
                                            -------- CHECK -------------
1D: PolV =     vector(compV)            --> polarization vector in simulation system
  2D: PolV = (/ (vector(compV(k)) , k=1,2) /)
3D: PolV = (/ (vector(compV(k)) , k=1,3) /)

•  subroutine get_piezoelectricity(i,PolV)
!-----------------------------------------------------------------------------
! Provides for point number i on Güenther's grid
! the polarization vector in coordinates of simulation region.
! Units: C/m^2
! 1D: PolV(1)   = component of polarization vector in simulation direction
! 2D: PolV(1:2) = component of polarization vector (Px,Py)
! 3D: PolV(1:3) = component of polarization vector (Px,Py,Pz)

Polarization charges are calculated according to  - div vec P = charge.
They are stored in MODULE polarizations.

!----------------------------------------------------------------------------------
 MODULE polarizations
!----------------------------------------------------------------------------------
!
! Pyroelectricity:
!
! background_pyroV    (1..lengvecV)                --> For graded regions, e.g. ternaries with
!                                                      linear alloy profile.
!                                                      Only different from zero for single points.
!                                                      Also defined for multiple points.
!                                                      Units: C/m^3
! 1D: interface_pyro1DV(1..num_mult_points)        --> Interface charge due to
! 2D: interface_pyroV  (1..num_mult_points,1:4)        pyroelectricity in SI units [C/m^2].
! 3D: interface_pyroV  (1..num_mult_points,1:12)
!
! 1D: boundary_pyro1DV  (1:2)                      --> The same at boundary.
!     boundary_pyroV does not exist in 2D/3D.
!
!
! Piezoelectricity:
!
! background_piezoV   (1..lengvecV)                --> For graded regions, e.g. ternaries with
!                                                      linear alloy profile.
!                                                      Only different from zero for single points.
!                                                      Also defined for multiple points.
!                                                      Units: C/m^3
! 1D: interface_piezo1DV(1..num_mult_points)       --> Interface charge due to
! 2D: interface_piezoV  (1..num_mult_points,1:4)       piezoelectricity in SI units [C/m^2].
! 3D: interface_piezoV  (1..num_mult_points,1:12)
!
! 1D: boundary_piezo1DV (1:2)                      --> The same at boundary.
!     boundary_piezoV does not exist in 2D/3D.
!
!----------------------------------------------------------------------------------

 REAL(8),DIMENSION(2)           :: boundary_pyro1DV   ! 1D only
 REAL(8),DIMENSION(2)           :: boundary_piezo1DV  ! 1D only

 REAL(8),DIMENSION(:)  ,POINTER :: interface_pyro1DV  ! 1D only
 REAL(8),DIMENSION(:)  ,POINTER :: interface_piezo1DV ! 1D only

 REAL(8),DIMENSION(:,:),POINTER :: interface_pyroV    ! 2D/3D only
 REAL(8),DIMENSION(:,:),POINTER :: interface_piezoV   ! 2D/3D only

 REAL(8),DIMENSION(:)  ,POINTER :: background_pyroV   ! 1D/2D/3D
 REAL(8),DIMENSION(:)  ,POINTER :: background_piezoV  ! 1D/2D/3D

Note:

In a general 3D simulation the program gets the polarization vectors of all eight octants:
There are 12 (3D) internal interfaces between these 8 octants.
In case there is an interface, a corresponding interface charge is set: interface_pyroV,interface_piezoV
In case there is no interface a corresponding volume charge is set: background_pyroV,background_piezoV
The sum of the volume charges of ALL octants is written to  background_pyroV(adr), where adr is the position of the FIRST octant on Stefan's grid: adr=adresse_3DM(i,j,k)

Note:

On domain boundaries the box is assumed to be mirrored to the outer region, so the polarization vectors across domain boundaries are assumed to be equal, so no surface charges are set.
--> boundary_pyroV, boundary_piezoV do not exist in 2D/3D
(In 1D they exist but do not make sense due to our convention --> set to zero.)

 

Internal details:

 !------------- mapping ---------------------------------------------
 !
 ! mapping contains for each interface sigma_1..12
 ! which octants are involved mapping(i,1) --> octant1
 !                            (i,2) --> octant2
 !
 ! and which component of polarization vector is relevant
 !                     mapping(i,3) =1 --> x-component
 !                                  =2 --> y-component
 !                                  =3 --> z-component
 !-------------------------------------------------------

 ! sigma_1
  mapping( 1,1)=1            --> interface 1 is the interface between octant 1 and 5
  mapping( 1,2)=5                the normal vector of this interface is the z-axis, so
  mapping( 1,3)=3                the relevant component of the polarization vector is 3

 ! sigma_2
  mapping( 2,1)=2 ! octant2
  mapping( 2,2)=6 ! octant6
  mapping( 2,3)=3 ! z-component

 ! sigma_3
  mapping( 3,1)=3 ! octant3
  mapping( 3,2)=7 ! octant7
  mapping( 3,3)=3 ! z-component

 ! sigma_4
  mapping( 4,1)=4 ! octant4
  mapping( 4,2)=8 ! octant8
  mapping( 4,3)=3 ! z-component

 ! sigma_5
  mapping( 5,1)=1
  mapping( 5,2)=3
  mapping( 5,3)=2 ! y-component

 ! sigma_6
  mapping( 6,1)=2
  mapping( 6,2)=4
  mapping( 6,3)=2 ! y-component

 ! sigma_7
  mapping( 7,1)=5
  mapping( 7,2)=7
  mapping( 7,3)=2 ! y-component

 ! sigma_8
  mapping( 8,1)=6
  mapping( 8,2)=8
  mapping( 8,3)=2 ! y-component

 ! sigma_9
  mapping( 9,1)=1
  mapping( 9,2)=2
  mapping( 9,3)=1 ! x-component

 ! sigma_10
  mapping(10,1)=3
  mapping(10,2)=4
  mapping(10,3)=1 ! x-component

 ! sigma_11
  mapping(11,1)=5
  mapping(11,2)=6
  mapping(11,3)=1 ! x-component

 ! sigma_12
  mapping(12,1)=7
  mapping(12,2)=8
  mapping(12,3)=1 ! x-component

 !---------------------------------------------------------------------------

In 2D the mapping is as follows:

 !------------- mapping ---------------------------------------------
 !
 ! mapping contains for each interface sigma_1..4 (i=sigma_1..4)
 ! which octants are involved mapping(i,1) --> octant1
 !                            (i,2) --> octant2
 !
 ! and which component of polarization vector is relevant
 !                     mapping(i,3) = 1 --> x-component
 !                                  = 2 --> y-component
 !
 !                       |
 !                       |
 !                    sigma_3
 !                       |
 !       octant3         |       octant4
 !                       |
 !                       |
 !   ----sigma_2-----------------sigma_4----
 !                       |
 !                       |
 !       octant1         |       octant2
 !                       |
 !                    sigma_1
 !                       |
 !                       |
 !
 !-------------------------------------------------------

 ! sigma_1
  mapping(1,1)=1 ! octant1
  mapping(1,2)=2 ! octant2
  mapping(1,3)=1 ! x-component

 ! sigma_2
  mapping(2,1)=1 ! octant1
  mapping(2,2)=3 ! octant3
  mapping(2,3)=2 ! y-component

 ! sigma_3
  mapping(3,1)=3 ! octant3
  mapping(3,2)=4 ! octant4
  mapping(3,3)=1 ! x-component

 ! sigma_4
  mapping(4,1)=2 ! octant2
  mapping(4,2)=4 ! octant4
  mapping(4,3)=2 ! y-component

!---------------------------------------------------------------------------

The problem is complicated at contacts which are defined by Poisson clusters.
The polarization charges at interfaces from a pyroelectric material (say GaN) to a metal contact are very large.

At contacts, however, one has specific contact models (i.e. ohmic --> field = 0), so one does not want to have interface charges there.
 --> The subroutine interface_to_poisson_clstr removes these interface charges, in detail:

It loops through all multiple points and checks for each octant if it is contained in a Poisson cluster.
Then it loops through all 12 internal interfaces and if at least one of the adjoining octants belongs to a Poisson cluster this interface charge (pyro and piezo) is set to zero.

Note the $numeric-control specifiers (see there for details):
- piezo-charge-at-boundaries
-  pyro-charge-at-boundaries
- piezo-constants-zero
-  pyro-constants-zero