Piecewise terms in non-linear thermo-elastic equation

Hi Everyone

I need to incorporate to my code a non-linear piecewise function for Kappa which depend on the temperature.

My code is already running with

def kappa(T):
return alpha0*(0.8383+0.00142*T)**(-0.5853)*(2*E0/2/(1+nu0)+3*E0*nu0/((1+nu0)*(1-2*nu0)))

However if I try to use

def kappa(T):
if T<600:
alpha=alpha0*(0.8383+0.00142*T)**(-0.5853)*(2*E0/2/(1+nu0)+3*E0*nu0/((1+nu0)*(1-2*nu0)))
else:
alpha=alpha0
return alpha

I got this error:

Traceback (most recent call last):
File â€śnon_linear_thermoelasticity_fenics.pyâ€ť, line 120, in
mech_form = inner(sigma(du, dTheta), eps(u_))*dx - inner(f, u_)dx + qdot(n,u_)ds(104)
File â€śnon_linear_thermoelasticity_fenics.pyâ€ť, line 115, in sigma
return (lmbda
tr(eps(v)) - kappa(Theta)Theta)Identity(2) + 2mueps(v)
File â€śnon_linear_thermoelasticity_fenics.pyâ€ť, line 67, in kappa
if ge(T,600):
File â€ś/home/numericalmining/anaconda3/envs/fenics2018/lib/python3.7/site-packages/ufl/conditional.pyâ€ť, line 46, in bool
error(â€śUFL conditions cannot be evaluated as bool in a Python context.â€ť)
File â€ś/home/numericalmining/anaconda3/envs/fenics2018/lib/python3.7/site-packages/ufl/log.pyâ€ť, line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: UFL conditions cannot be evaluated as bool in a Python context.

Here is my full code and the mesh to run it (Dropbox - Files.zip - Simplify your life):

from dolfin import *
from mshr import *
from fenics import *
import numpy as np
if has_linear_algebra_backend("Epetra"):
parameters["linear_algebra_backend"] = "Epetra"

has_linear_algebra_backend('PETSc')
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

mesh = Mesh("dyke_mesh_super_fine.xml")
facets = MeshFunction("size_t", mesh, "dyke_mesh_super_fine_facet_region.xml")
ds = ds(subdomain_data=facets)

#caracteristica de la malla
#W=6000 #Model width m
#H=3000 #Model height m
#Prof=1500 #Mean depth magma chamber m
#w1=10 #Horizontal length of the magma chamber m
#h1=500 #Vertical length of the magma chamber m
#el (0,0) del modelo estĂˇ en el vertice superior izquierdo

Temperatura=data[0,0]
Poisson=data[0,2]
Expansiontermica=data[0,3]
Presionmagma=data[0,6]
#nmodel=len(Temperatura)

#for i in np.arange(nmodel):
tmax=24*5   #horas
num_step=10
dtime=tmax/num_step
nu0=Poisson # razon de poisson, adimensional
alpha0=Expansiontermica #coeficiente de expansion termica, 1/Â°C
p0=Presionmagma #Magma pressure, MPa
filename="sol_thermo_Tintro="+str(Tintru)+"_rho="+str(rho0)+"_E0="+str(E0)+"_nu0="+str(nu0)+"_alpha0="+str(alpha0)+"_cv0="+str(cv0)+"_k0="+str(k0)+"_p0="+str(p0)+".xdmf"
Tm=Constant(Tintru) #temperatura magma
DTdike=Expression('Tm-T0',degree=1,Tm=Tm,T0=T0)
E=Constant(E0)
nu = Constant(nu0)
lmbda = Constant(E*nu/((1+nu)*(1-2*nu)))
mu = Constant(E/2/(1+nu))
rho = Constant(rho0)     # density
alpha = Constant(alpha0)  # 8e-6 thermal expansion coefficient
#kappa = Constant(alpha*(2*mu + 3*lmbda))
cV = Constant(cv0)*rho  # specific heat per unit volume at constant strain
k = Constant(k0)  # thermal conductivity
rho_g = 0   # 0 si no hay fuerzas de cuerpo, 0.027 si las hay
f = Constant((0, -rho_g))
q=Constant(p0)

def kappa(T):
return alpha0*(0.8383+0.00142*T)**(-0.5853)*(2*E0/2/(1+nu0)+3*E0*nu0/((1+nu0)*(1-2*nu0)))

#4. Define space of solutions and boundary conditions
Vue = VectorElement('CG', mesh.ufl_cell(), 2) # displacement finite element
Vte = FiniteElement('CG', mesh.ufl_cell(), 1) # temperature finite element
V = FunctionSpace(mesh, MixedElement([Vue, Vte]))

#modelo semi
#bc1 = DirichletBC(V.sub(0).sub(1), Constant(0.), facets,102) #inferior
#bc2 = DirichletBC(V.sub(0).sub(0), Constant(0.), facets,101) #laterales
#bc3 = DirichletBC(V.sub(0).sub(0), Constant(0.), facets,106) #laterales
#bc4 = DirichletBC(V.sub(1), DTdike, facets,104) #interior
#bc5 = DirichletBC(V.sub(1), Constant(0.), facets,101) #inferior
#bc6 = DirichletBC(V.sub(1), Constant(0.), facets,103) #superficie
#bcs = [bc1, bc2, bc3, bc4, bc5, bc6]

#modelo entero
bc1 = DirichletBC(V.sub(0), Constant((0.,0.)), facets,102) #inferior
bc2 = DirichletBC(V.sub(0).sub(0), Constant(0.), facets,101) #laterales
bc4 = DirichletBC(V.sub(1), DTdike, facets,104) #interior
bc5 = DirichletBC(V.sub(1), Constant(0.), facets,102) #inferior
bc6 = DirichletBC(V.sub(1), Constant(0.), facets,103) #superficie
bc7 = DirichletBC(V.sub(1), Constant(0.), facets,101) #laterales
bcs = [bc1, bc2, bc4, bc5, bc6, bc7]

n = FacetNormal(mesh)

#5.Espacios de funciones y formulaciones debiles
U_ = TestFunction(V)
(u_, Theta_) = split(U_)
dU = Function(V)
#dU = TrialFunction(V)
(du, dTheta) = split(dU)
Uold = Function(V)
(uold, Thetaold) = split(Uold)

def eps(v):

def sigma(v, Theta):
return (lmbda*tr(eps(v)) - kappa(Theta)*Theta)*Identity(2) + 2*mu*eps(v)

dt=Constant(dtime)
T0=project(T0,V.sub(1).collapse())

mech_form = inner(sigma(du, dTheta), eps(u_))*dx - inner(f, u_)*dx + q*dot(n,u_)*ds(104)

therm_form = (cV*(dTheta-Thetaold)/dt*Theta_ +
kappa(dTheta+T0)*T0*tr(eps(du-uold))/dt*Theta_ +

form = mech_form + therm_form

U = Function(V)
file_results = XDMFFile(filename)
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
t=0

for i in np.arange(num_step):
t=t+dtime
solve(form == 0, dU, bcs, solver_parameters={"newton_solver":{"relative_tolerance": 1e-6}})
Uold.assign(dU)
u, T = dU.split()
u.rename('Displacement','label')
file_results.write(u, t)
T.rename('Temperature variation','label')
file_results.write(T,t)
Vsig = TensorFunctionSpace(mesh, "DG",degree=1)
sig = Function(Vsig, name="Cauchy_Stress")
sig.assign(project(-sigma(u,T), Vsig))
file_results.write(sig, t)
Epssig = TensorFunctionSpace(mesh, "DG",degree=1)
epsilon = Function(Epssig, name="Strain")
epsilon.assign(project(eps(u), Epssig))
file_results.write(epsilon,t)
sx = sig.sub(0)
sy = sig.sub(3)
txy = sig.sub(1)

def sigma1(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
return((sigma_x+sigma_y)/2 + sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2))

def sigma3(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
return((sigma_x+sigma_y)/2 - sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2))

S1 = sigma1(sx,sy,txy)
S3 = sigma3(sx,sy,txy)
Difstress=S1-S3
Promstress=(S1+S3)/2
Vp = FunctionSpace(mesh, "P",1)
sig1=Function(Vp,name="Sigma 1")
sig1.assign(project(S1,Vp))
sig3=Function(Vp,name="Sigma 3")
sig3.assign(project(S3,Vp))
sigdif=Function(Vp,name="Differential stress")
sigdif.assign(project(Difstress,Vp))
sigprom=Function(Vp,name="Mean Stress")
sigprom.assign(project(Promstress,Vp))
file_results.write(sig1,t)
file_results.write(sig3,t)
file_results.write(sigdif,t)
file_results.write(sigprom,t)
print('Progress: '+str(100*t/tmax)+'%')

Some idea how I should solve it?

1 Like