I would like to solve Landau-Khalatnikov (LK) equation and Poisson equation using finite element approach to simulate ferroelectric hysteresis loop, electric potential distribution and electric polarization distribution. Here, I will provide detail as much as possible. LK equation, which is also called time-dependent Ginzburg-Landau (TDGL) equation, can be simply expressed as follow:
#1: $-\frac{1}{\Gamma}\frac{\partial P}{\partial t}= \alpha P + \beta P^{3} + \gamma P^{5}-E$
Where Γ is kinetic coefficient, α, β, γ are Landau coefficients, P is polarization and E is electric field (of course, t is time). The picture is much like a ferroelectric material sandwiched by electrode, one end has fixed potential and other end is grounded to generate electric field. First, to simply test feasibility of MMA to solve such an equation, I wrote very simple code to solve LK equation in 1D:
Clear[\[Alpha], \[Beta], \[Gamma], Efieldamp, \[Omega]];
\[Alpha] = -1.35*10^9;
\[Beta] = -2.64*10^10;
\[Gamma] = 2.5*10^11;
Efieldamp = 2*10^9;
\[Omega] = 10^6*Pi;
\[CapitalGamma] = 1;
Efield[t_?NumericQ] := Efieldamp*Sin[\[Omega]*t];
SimpleLKSolver =
First@NDSolve[{-1/\[CapitalGamma]*D[P[t], t] == \[Alpha]*(P[t]*0.01) + \[Beta]*(P[t]*0.01)^3 +
\[Gamma]*(P[t]*0.01)^5 - Efield[t], P[0] == 0}, P, {t, 0, 2.5*1/\[Omega]*Pi}, Method -> Automatic, MaxSteps -> Infinity]
As shown in the code, I put E as Efieldamp * Sin[ω * t] where Efieldamp and ω are defined as above. The multiplication of 0.01 to P in rhs is for unit conversion (from C/m2 to uC/cm2). The code works nicely, and plotting P vs. V (V was calculated assuming that the distance between voltage source and ground is 10 nm (10-8 m)) gives nice hysteresis curve:
LK =
ParametricPlot[{Efield[t]*10^-8, P[t] /. SimpleLKSolver}, {t, 0, 2.5*1/\[Omega]*Pi}, PlotRange -> All, AspectRatio -> 2/3,
PlotStyle -> {Bold, Black}, Frame -> True,
FrameStyle -> {{Black, Bold}, {Black, Bold}}]
Now, I would like to solve the equation in 2D using finite element method. Actually, I would also like to simultaneously solve Poisson’s equation to see electric potential. In 2D, the equations are:
#2: $-\frac{1}{\Gamma}\frac{\partial P(x,z)}{\partial t}= \alpha P(x,z) + \beta P(x,z)^{3} + \gamma P(x,z)^{5}-g_{11}\frac{\partial^2 P(x,z)}{\partial z^2}-g_{44}\frac{\partial^2 P(x,z)}{\partial x^2}+\frac{\mathrm{d} \phi(x,z)}{\mathrm{d} x}$
#3: $-\varepsilon _{0}(\varepsilon _{x}\frac{\partial^2 \phi }{\partial x^2}+\varepsilon _{z}\frac{\partial^2 \phi }{\partial z^2})=-\frac{\mathrm{d} P}{\mathrm{d} z}$
In 2D space, I defined x as horizontal direction and z as vertical direction. The spatial second derivative terms are now added to LK (1st) equation in order to calculate spatial gradient of polarization; g11 (g44) are gradient coefficient, ϕ is electric potential, εx and εz are anisotropic dielectric permittivity in x and z direction. dϕ(x,z)/dz is electric field term in 2D, considering electric field only in Z-direction. In Poisson equation (2nd equation), -dP(x,z)/dz gives volume charge density considering polarization in ferroelectric materials. So, the workflow is basically this: first solve LK equation to obtain polarization distribution in 2D and obtain PV hysteresis curve, then obtain electric potential distribution by solving Poisson equation with obtained P value. My try is the following:
Needs["NDSolve`FEM`"]
region = Rectangle[{0, 0}, {100*10^-9, 10*10^-9}];
mesh = ToElementMesh[region, MaxCellMeasure -> 2*10^-18]
Clear[\[Alpha], \[Beta], \[Gamma], g44, g11, \[Epsilon]x, \[Epsilon]z, \[Epsilon]0, \[Omega], \[CapitalGamma]];
\[Alpha] = -1.35*10^9;
\[Beta] = -2.64*10^10;
\[Gamma] = 2.5*10^11;
g44 = 2*10^-10;
g11 = 0.5*10^-10;
\[Epsilon]x = 25;
\[Epsilon]z = 18;
\[Epsilon]0 = 8.854*10^-12;
\[Omega] = 10^6*Pi;
\[CapitalGamma] = 1;
Vapp[t_] := 20*Sin[\[Omega]*t]
PDESolution =
First@NDSolve[{Inactive[Div][{{\[Epsilon]0*\[Epsilon]x, 0}, {0, \[Epsilon]0*\[Epsilon]z}}.Inactive[Grad][\[Phi][t, x, z], {x, z}], {x, z}] ==
D[P[t, x, z], z], -1/\[CapitalGamma]*
D[P[t, x, z], t] == \[Alpha]*P[t, x, z] + \[Beta]*
P[t, x, z]^3 + \[Gamma]*P[t, x, z]^5 -
g11*D[P[t, x, z], {z, 2}] - g44*D[P[t, x, z], {x, 2}] +
D[\[Phi][t, x, z], z],
DirichletCondition[\[Phi][t, x, z] == 0,
z == 0 && 0 < x < 10*10^-8],
P[0, x, z] == RandomVariate[NormalDistribution[]], \[Phi][0, x, z] == 0,
DirichletCondition[\[Phi][t, x, z] == Vapp[t],
z == 10*10^-9 && 0 < x < 10*10^-8],
PeriodicBoundaryCondition[P[t, x, z], x == 0,
TranslationTransform[{10*10^-8, 0}]],
PeriodicBoundaryCondition[\[Phi][t, x, z], x == 0,
TranslationTransform[{10*10^-8, 0}]]}, {\[Phi],
P}, {x, z} \[Element] mesh, {t, 0, 2.5*1/\[Omega]*Pi},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement"}}}]
For boundary condition, I put ϕ = Vapp (defined above) for top surface (z = 10-8) and ϕ = 0 for bottom surface (z = 0). I also applied periodic boundary conditions (in x direction) for both ϕ and P. For initial condition, I put ϕ(t =0) = 0, whereas I put P= RamdomVariate[NormalDistribution[]] to assign random initial condition spatially (but sum is zero!) This code gives me the following error:
NDSolve::ivcon: The given initial conditions were not consistent with the differential-algebraic equations. NDSolve will attempt to correct the values.
LinearSolve::sing: Matrix <<1>> is singular.
LinearSolve::sing: Matrix SparseArray[<<1>>] is singular.
General::stop: Further output of LinearSolve::sing will be suppressed during this calculation.
With the given error message, I know that something is wrong with my initial condition of P, which is P(0,x,z)= RandomVariate[NormalDistribution[]]; therefore, I put 0 in rhs instead. Then, I obtained repeated error messages of the following lines:
LinearSolve::sing: Matrix SparseArray[Automatic,<<1>>,<<3>>,{1,{{<<1>>},{<<1>>}},{4.*10^9,-4.44089*10^-7,-4.*10^9,1.77636*10^-6,4.*10^9,-4.*10^9,<<51862>>,2.72987*10^-10,-5.8139*10^-10,-7.82473*10^-10,-1.44678*10^-9,-7.82473*10^-10,4.*10^9}}] is singular.
NDSolve::femdpop: The FEMLoadElements operator failed.
So, I realized that some singularity issue exists in the solution. To remove complication, I also tried removing gradient terms (second-order spatial derivatives), but similar error message came out (singularity error). Having confirmed that 1D solution can easily be obtained using NDSolve command, I assume that 2D solution also exists and that something is wrong with my code for solving 2D problem.
More explanation on my question
Basically, I'm trying to reproduce the results from the following paper published in arXiv: https://arxiv.org/ftp/arxiv/papers/2105/2105.04647.pdf
The major difference is that the authors here utilized "finite difference method", whereas I'm trying to utilize finite element method. Actually, the equations should be solved in "self-consistent" manner, meaning that, for each time step, polarization calculated from equation #2 (Now I labelled the equations!) are put into equation #3 to calculate electric potential, which is again put into equation #2; this iteration stops until the solutions converge. Though equation #3 is not time-dependent, it is solved statically for each P calculated from each time step. I believe how it is solved is analogous to the case of coupled Poisson-schrodinger equation, where charge distribution is calculated from the wave functions calculated using Schrodinger's equation. and this calculated charge distribution is put into Poisson equation to calculate the electric potential, and so on.
To build the code, I have read some tutorial articles in mathematica website (such as components and data structures) but I have no idea. Any hints or helps will be deeply appreciated!
D[P[t, x, z], z], while in the original model there is external electric field only. This term makes model unstable so that electric potential increases up to $Vapp/(\epsilon_0 \epsilon_z)$. $\endgroup$