Electrochemical phase field can not converge

zhang2022 · March 11, 2023, 6:04am

I try to solve electrochemical phase field in the paper ‘Dendrite-free Lithium Based on Lessons Learned from Lithium and Magnesium Electrodeposition Morphology Simulations’. But the it does not converge. No matter which parameter I change, the r(abs) does not change significantly.
Here is my code
#####################from fenics import *
import numpy as np
import matplotlib.pyplot as plt

tol = 1.0e-16

num_steps = int(1000)
R = 8.314 #J/mol/k
T = 298 #k
Far = Constant(96485) #As/mol
RT = Constant(R*T)
alpha = Constant(0.5)
Phie_i = -0.45
i0 = Constant(0.376) # 0.376 A/m2

delta = Constant(1.0e-4) # 1.0e-4 m
tao = 4000.0 # 4000s
C0 = 1000 # 1000 mol/m3
E0 = 1.5e8 # 1.5e7 J/m3

Cs_i = 7.692e4 # 7.692e4 mol/m3 #
D_s_i = 7.5e-14 # 7.5e-14 m2/s
D_l_i = 7.5e-10 # 7.5e-10 m2/s
Sig_s_i = 1.0e7 # 1.0e7 S/m
Sig_l_i = 1.0 # 1.0 S/m
################################Following 4 parameters influence the ‘r’ slightly######################
kappa_i = 3.68e-6 # 3.68e-06 J/m
W_i = 1.18e+6 #1.18e+06 J/m3
Leta_i = 5.1e-3 # 5.10e-03 1/s
Lsig_i = 2.50e-6 # 2.50e-06 m3/(J*s)

dt_i = 1.0e-2tao
lx_i= 2.0delta
ly_i = delta
nx = 200
ny = 100
#########################################################################################################
#Standardization
dt = Constant(dt_i/tao) # 1.0e-5
lx = Constant(lx_i/delta) #2.0
ly = Constant(ly_i/delta) #1.0
W = Constant(W_i/E0) # 0.079
kappa = Constant(kappa_i/(E0delta**2)) #2.45e-5
Cs = Constant(Cs_i/C0) # 76.92
Leta = Constant(Leta_itao) # 0.0051
Lsig = Constant(Lsig_iE0tao) # 37.5
Dc = Constant(delta2/tao)
D_s = Constant(D_s_i/Dc) # 7.5e-6
D_l = Constant(D_l_i/Dc) # 0.075
Sigmac = Constant(C0*Far2delta**2/(taoRT))
Sig_s = Constant(Sig_s_i/Sigmac) # 2.68e5
Sig_l = Constant(Sig_l_i/Sigmac) # 0.0268
Phie = Constant(Phie_i/(RT/Far)) #-17.407
#################################################################################################################

mesh = RectangleMesh(Point(0, 0), Point(lx, ly), nx,ny)
P1 = FiniteElement(‘P’,triangle, 1)
V = FunctionSpace(mesh, MixedElement(P1, P1, P1))
v1,v2,v3 = TestFunctions(V)

u = Function(V)
u_n = Function(V)
Xi, C_Li, Phi = split(u)
Xi_n, C_Li_n, Phi_n = split(u_n)

u_init = Expression((‘0.5*(1.0-1.0tanh((x[0]-0.05)/1.0e-3))‘,‘x[0] < (0.05) ? 0.0 : 1.0’,’-0.225(1.0-tanh((x[0]-0.05)/1.0e-3))/(RT/Far)’), degree=3,Far = Far, RT = RT)
u.interpolate(u_init)
u_n.interpolate(u_init)

def boundary0(x,on_boundary):
return on_boundary and near(x[0], 0)
def boundaryl(x,on_boundary):
return on_boundary and near(x[0], lx)

bc_Xi0 = DirichletBC(V.sub(0), Constant(1.0), boundary0)
bc_Xi1 = DirichletBC(V.sub(0), Constant(0.0), boundaryl)

bc_C_Li0 = DirichletBC(V.sub(1), Constant(0.0), boundary0)
bc_C_Li1 = DirichletBC(V.sub(1), Constant(1.0), boundary0)

bc_Phi0 = DirichletBC(V.sub(2), Constant(Phie), boundary0)
bc_Phi1 = DirichletBC(V.sub(2), Constant(0.0), boundaryl)

bcs = [bc_Xi0, bc_Xi1, bc_C_Li0, bc_C_Li1, bc_Phi0, bc_Phi1]

def g(xi):
return xi2*(Constant(1.0)-xi)2
def h(xi):
return xi3*(Constant(6.0)*xi2-Constant(15.0)xi+Constant(10.0))
def dg(xi):
return Constant(2.0)xi(Constant(1.0)-xi)(Constant(1.0)-Constant(2.0)xi)
def dh(xi):
return Constant(30.0)xi**2(xi-Constant(1.0))**2
def D(xi):
return D_sh(xi) + D_l*(Constant(1.0)-h(xi))
def sigma(xi):
return 1.0e-6*(Sig_sh(xi) + Sig_l(Constant(1.0)-h(xi)))

F1 = (Xi-Xi_n)/dtv1dx + kappaLsigdot(grad(Xi),grad(v1))dx + LsigWdg(Xi)v1dx + Letadh(Xi)(exp((Constant(1.0)-alpha)Phi)-C_Liexp(-alphaPhi))v1dx
F2 = (C_Li-C_Li_n)/dtv2dx + D(Xi)dot(grad(C_Li),grad(v2))dx + D(Xi)C_Lidot(grad(Phi),grad(v2))dx + Cs(Xi-Xi_n)/dtv2dx
F3 = sigma(Xi)dot(grad(Phi),grad(v3))dx + Cs(Xi-Xi_n)/dtv3*dx

F = F1 + F1 + F3

J = derivative(F,u)

t = 0
for n in range(num_steps):
# Update time
t+=dt
# Solve problem
solve(F == 0, u, bcs, J=J, solver_parameters={“newton_solver”:{“absolute_tolerance”:1.0e-6, “relative_tolerance”:1.0e-6, “maximum_iterations”:100}})
# Update previous solution
u_n.assign(u)
if (n+1)%int(100) == 0:
plot(Xi, title = ‘Xi’,mode = ‘color’)
plt.show()
plot(C_Li, title = ‘C_Li’,mode = ‘color’)
plt.show()
plot(Phi, title = ‘Phi’,mode = ‘color’)
plt.show()
###########################################################

dokken · March 11, 2023, 7:36am

Please format your code with 3x` encapsulation, and make sure that the code can be copy pasted and reproduce the error.

```python
#add code here 

```

zhang2022 · March 11, 2023, 9:56am

Thank you.

Below is the code

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

tol = 1.0e-16

num_steps = int(1000)
R = 8.314 #J/mol/k
T = 298 #k
Far = Constant(96485) #As/mol
RT = Constant(R*T)
alpha = Constant(0.5)
Phie_i = -0.45
i0 = Constant(0.376) # 0.376 A/m2

delta = Constant(1.0e-4) # 1.0e-4 m
tao  = 4000.0 # 4000s
C0 = 1000 # 1000 mol/m3
E0 = 1.5e8 # 1.5e7 J/m3

Cs_i = 7.692e4 # 7.692e4  mol/m3 #
D_s_i = 7.5e-14 # 7.5e-14 m2/s
D_l_i = 7.5e-10 # 7.5e-10 m2/s
Sig_s_i = 1.0e7 # 1.0e7 S/m
Sig_l_i = 1.0 # 1.0 S/m
################################Following 4 parameters influence the 'r' slightly######################
kappa_i = 3.68e-6 # 3.68e-06 J/m
W_i = 1.18e+6 #1.18e+06 J/m3
Leta_i = 5.1e-3 # 5.10e-03 1/s
Lsig_i = 2.50e-6 # 2.50e-06 m3/(J*s)

dt_i = 1.0e-2*tao
lx_i= 2.0*delta
ly_i = delta
nx = 200
ny = 100
#########################################################################################################
#Standardization
dt = Constant(dt_i/tao) # 1.0e-5
lx = Constant(lx_i/delta) #2.0
ly = Constant(ly_i/delta) #1.0
W = Constant(W_i/E0) # 0.079
kappa = Constant(kappa_i/(E0*delta**2)) #2.45e-5
Cs = Constant(Cs_i/C0) # 76.92
Leta = Constant(Leta_i*tao) # 0.0051
Lsig = Constant(Lsig_i*E0*tao) # 37.5
Dc = Constant(delta**2/tao)
D_s = Constant(D_s_i/Dc) # 7.5e-6
D_l = Constant(D_l_i/Dc) # 0.075
Sigmac = Constant(C0*Far**2*delta**2/(tao*RT))
Sig_s = Constant(Sig_s_i/Sigmac) # 2.68e5
Sig_l = Constant(Sig_l_i/Sigmac) # 0.0268
Phie = Constant(Phie_i/(RT/Far)) #-17.407
#################################################################################################################

mesh = RectangleMesh(Point(0, 0), Point(lx, ly), nx,ny)
P1 = FiniteElement('P',triangle, 1)
V = FunctionSpace(mesh, MixedElement(P1, P1, P1))
v1,v2,v3 = TestFunctions(V)

u = Function(V)
u_n = Function(V)
Xi, C_Li, Phi = split(u)
Xi_n, C_Li_n, Phi_n = split(u_n)

u_init = Expression(('0.5*(1.0-1.0*tanh((x[0]-0.05)/1.0e-3))','x[0] < (0.05) ? 0.0 : 1.0','-0.225*(1.0-tanh((x[0]-0.05)/1.0e-3))/(RT/Far)'), degree=3,Far = Far, RT = RT)
u.interpolate(u_init)
u_n.interpolate(u_init)

def boundary0(x,on_boundary):
    return on_boundary and near(x[0], 0)
def boundaryl(x,on_boundary):
    return on_boundary and near(x[0], lx)

bc_Xi0 = DirichletBC(V.sub(0), Constant(1.0), boundary0)
bc_Xi1 = DirichletBC(V.sub(0), Constant(0.0), boundaryl)

bc_C_Li0 = DirichletBC(V.sub(1), Constant(0.0), boundary0)
bc_C_Li1 = DirichletBC(V.sub(1), Constant(1.0), boundary0)

bc_Phi0 = DirichletBC(V.sub(2), Constant(Phie), boundary0)
bc_Phi1 = DirichletBC(V.sub(2), Constant(0.0), boundaryl)

bcs = [bc_Xi0, bc_Xi1, bc_C_Li0, bc_C_Li1, bc_Phi0, bc_Phi1]

def g(xi):
    return xi**2*(Constant(1.0)-xi)**2
def h(xi):
    return xi**3*(Constant(6.0)*xi**2-Constant(15.0)*xi+Constant(10.0))
def dg(xi):
    return Constant(2.0)*xi*(Constant(1.0)-xi)*(Constant(1.0)-Constant(2.0)*xi)
def dh(xi):
    return Constant(30.0)*xi**2*(xi-Constant(1.0))**2
def D(xi):
    return D_s*h(xi) + D_l*(Constant(1.0)-h(xi))
def sigma(xi):
    return 1.0e-6*(Sig_s*h(xi) + Sig_l*(Constant(1.0)-h(xi)))

F1 = (Xi-Xi_n)/dt*v1*dx + kappa*Lsig*dot(grad(Xi),grad(v1))*dx + Lsig*W*dg(Xi)*v1*dx + Leta*dh(Xi)*(exp((Constant(1.0)-alpha)*Phi)-C_Li*exp(-alpha*Phi))*v1*dx
F2 = (C_Li-C_Li_n)/dt*v2*dx + D(Xi)*dot(grad(C_Li),grad(v2))*dx + D(Xi)*C_Li*dot(grad(Phi),grad(v2))*dx + Cs*(Xi-Xi_n)/dt*v2*dx
F3 = sigma(Xi)*dot(grad(Phi),grad(v3))*dx + Cs*(Xi-Xi_n)/dt*v3*dx

F = F1 + F1 + F3

J = derivative(F,u)

t = 0
for n in range(num_steps):
    # Update time
    t+=dt
    # Solve problem
    solve(F == 0, u, bcs, J=J, solver_parameters={"newton_solver":{"absolute_tolerance":1.0e-6, "relative_tolerance":1.0e-6, "maximum_iterations":100}})
    # Update previous solution
    u_n.assign(u)
    if (n+1)%100 == 0:
        print('#####################################################################################################')
        print('step number :', n)
        print('#####################################################################################################')
    if (n+1)%int(100) == 0:
        plot(Xi, title = 'Xi',mode = 'color')
        plt.show()
        plot(C_Li, title = 'C_Li',mode = 'color')
        plt.show()
        plot(Phi, title = 'Phi',mode = 'color')
        plt.show()

DcOtaku · December 16, 2023, 11:15am

Hi, I’m working on a similar problem and encountering the same issue, did you manage to resolve it after?

Topic		Replies	Views
[Phase Field] Coupled Vectorized Order Parameters variational formulation	1	319	October 5, 2022
Error does not converge. As mesh size is increased, error decreases but also increases sporadically	0	436	February 16, 2019
Newtonsolver does not converge - Hyperelasticity	5	933	May 2, 2020
Advice on: Newton solver does not converge variational formulation	3	2408	August 25, 2020
Error: Newton solver did not converge Errors	1	593	June 12, 2022

Electrochemical phase field can not converge

Related Topics