Hello, I am working on a multi-field problem with an approach of the theory of porous media / poromechanics (solid and fluid phase). The code I implemented is for large strains. For small loads (it is loaded with a stress vector) it works, but if I increase the load, the Newton solver does not converge anylonger.
The load is controlled by a ramp function. From the same point in time (= load), the solver no longer converges regardless of changes in the time step or the mesh. It seems very strange to me, because the first iteration gets more precisely from the beginning, but the second iteration becomes more and more imprecise.
When I look at the results that have been created up to then, there do not seem to be any major numerical instabilities. That is why I wonder what the problem is.
Has anyone already had a similar problem? Could it help to use a different solver? Does someone have another advice? Thank you in advance!
Below you can find an example code:
start = time.time()
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
##############################################################################
# Definition of Parameters #
##############################################################################
#fluid phase F and solid phase S
phi_F = 0.7
phi_S = 1-phi_F
#Lame-Konstanten
mu = 0.05
lambda_ = 0.05
#Permeability
K = 1.54e-3
# initial pressure
p_0 = 0.0
## Definition of time parameters
T1 = 10.0
num_steps = 100
dt = 0.1 # time step size
p_0 = Constant(p_0)
mu = Constant(mu)
K = Constant(K)
rho_FR = Constant(rho_FR)
lambda_ = Constant(lambda_)
##############################################################################
# mesh #
##############################################################################
mesh = Mesh("../mesh.xml")
##############################################################################
# Definition mixed elements #
##############################################################################
#Mixed Elements
P2 = VectorElement('Lagrange', Netz.ufl_cell(), 2)
P1 = FiniteElement('Lagrange', Netz.ufl_cell(), 1)
element = P2 * P1 #mixed element
# Define FunctionSpace
V = FunctionSpace(Netz, element) #mixed FunctionSpace für mixed element
##############################################################################
# Boundarys and conditions #
##############################################################################
# Tolerancd
tol = 1E-14
# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, "../mesh_facet_region.xml")
# Initialize mesh function for subdomains
subdomains = MeshFunction("size_t", mesh, "../mesh_physical_region.xml")
# boundaries
ds_test = Measure('ds', domain=Netz, subdomain_data=boundaries)
# normal vector
n_vec = FacetNormal(Netz)
## Stress vector - Parameters from Matlab
A_1F = -1.493395995800785e-4
A_2F = -0.004440986342513
A_3F = 0.069085083600410
A_4F = 0
A_1T = -0.006793085035899
A_2T = 0.071327392876938
A_3T = 0
t_vec3 = Expression(('0.0' , '-(A_1F * (pow(x[0],3)) + A_2F * (pow(x[0],2))+ A_3F * x[0] + A_4F) * (T_ges - m * t )/ T1 * phi_S'),degree=2, phi_S=0.45, m=-1, T_ges=0, A_1F = A_1F, A_2F = A_2F, A_3F = A_3F, A_4F = A_4F, T1=T1, t=0)
p_RF = Expression( '(p_0 + (A_1F * (pow(x[0],3)) + A_2F * (pow(x[0],2))+ A_3F * x[0] + A_4F) * (T_ges - m * t )/ T1 * phi_F )', degree=2, p_0 = p_0, phi_F=0.55, m=-1, T_ges=0, A_1F = A_1F, A_2F = A_2F, A_3F = A_3F, A_4F = A_4F, T1=T1, t=0)
p_RT = Expression( '(p_0 + (A_1T * (pow(x[0],2)) + A_2T * x[0] + A_3T) * (T_ges - m * t )/ T1 * phi_F )', degree=2, p_0 = p_0, phi_F=0.55, m=-1, T_ges=0, A_1T = A_1T, A_2T = A_2T, A_3T = A_3T, T1=T1, t=0)
# Verschiebungsrandbedingung
bc_u11 = DirichletBC(V.sub(0).sub(1), Constant(0.0), boundaries, 1)
bc_u2 = DirichletBC(V.sub(0), Constant((0.0, 0.0)), boundaries, 2)
# Druckrandbedingungen
bc_p1 = DirichletBC(V.sub(1), p_RT, boundaries, 1)
bc_p3 = DirichletBC(V.sub(1), p_RF, boundaries, 3)
bc_p4 = DirichletBC(V.sub(1), p_0, boundaries, 4)
# Class representing the intial conditions
class InitialConditions(UserExpression):
def eval(self, values, x):
values[0] = 0 # ux
values[1] = 0 # uy
values[2] = p_0 # p
def value_shape(self):
return (3,)
##############################################################################
Tensor_grad = 2
#Test- und Trial-Functions
dfunc = TrialFunction(V)
func = Function(V)
(u, p) = split(func)
(v_1, v_2) = TestFunctions(V)
x_n = Function(V)
u_n, p_n = split(x_n) # u, p at the time n
# Interpolate Initial Conditions
x_ini = InitialConditions()
x_n.interpolate(x_ini)
# Dimension (1D/2D/3D)
d = len(u_n) # Dimension
I = Identity(d)
# Kinematics
def F(u):
return I + grad((u))
def J(u):
return det(F(u))
C = F(u).T*F(u) # Right Cauchy-Green tensor
B = F(u)*F(u).T # Linker Cauchy-Green-tensor
# Invariants of deformation tensors
Ic = tr(C)
# Total-Lagrange - 2. PIOLA-Kirchhoff-Stresstensor
def T(u):
return mu*(I - inv(C)) + lambda_/2 * (J(u)**2 - 1) * inv(C)
def Piola2(u,p):
return T(u) - J(u) * inv(F(u)) * (p*I) * inv(F(u).T)
## Variational problem
Imp = inner(Piola2(u,p), dot(F(u), nabla_grad(v_1))) * dx \
- dot(dot(F(u), t_vec3), v_1) * ds_test(3) \
Vol = (J(u) - J(u_n))/dt * v_2 * dx \
+ K * J(u) * dot(dot(dot(inv(F(u)),inv(F(u).T)), nabla_grad(p) ) , nabla_grad(v_2)) * dx \
eqn = Imp + Vol
Jac = derivative(eqn, func, dfunc)
# VTK-files
vtkfile_u = File("Ergebnisse_s30f70_V5_treppen/u_{}.pvd".format(Titel))
vtkfile_p = File("Ergebnisse_s30f70_V5_treppen/p_{}.pvd".format(Titel))
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Update functions
if t <= T1:
t_vec3.t = t
p_RT.t = t
p_RF.t = t
#boundary conditions
bc_ges = [bc_u11, bc_u2, bc_p1, bc_p3, bc_p4]
# Solve variational problem for time step
solve(eqn == 0, func, bc_ges, J=Jac, solver_parameters={"newton_solver":
{"relative_tolerance": 1e-5,"absolute_tolerance": 1e-10}})
(u,p) = func.split(True)
# für Paraview
u.rename('u','u') # Das ist notwendig, damit Paraview nicht jeden
p.rename('p','p') # Zeitschritt als eine neue Variable ansieht.
# Save solution to file (VTK)
vtkfile_u << (u, t)
vtkfile_p << (p, t)
# Update previous solution
assign(x_n, [u,p])
So that you can imagine it better, the mesh with the different boundaries can be seen here:
