Dear everyone:
To constrain the electric potential equal on the top surface, I defined the submesh based on the top surface, and set Real functionspace for this part submesh for the reason I thought the Real functionspace would only have one constant in its region. And then I used langrange multiplier to tie the
dof in submesh of the top surface to the main solid mesh.
But the potential result of the top surface is still not uniform and is the same as the result without any constraint. Is there anything wrong with what I am thinking about R0 functionspace?
I attached my result of potential contour and my code:
from __future__ import print_function
from fenics import *
from dolfin import *
from ufl import nabla_div
import numpy
# Scaled variables
# L = 1
# W = 0.2 # L is length; W is width
rho = 5550
g = 9.8
# define elasticity tensor C
c11 = 226e9
c12 = 125e9
c13 = 124e9
c33 = 216e9
c44 = 44.2e9 # * 2
c66 = (c11 - c12) / 2 # * 2
# 6x6
C = as_matrix([[c11, c12, c13, 0.0, 0.0, 0.0],
[c12, c11, c13, 0.0, 0.0, 0.0],
[c13, c13, c33, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, c44, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, c44, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, c66]])
# define piezoelectricity tensor e
e15 = 5.8
e31 = -2.2
e33 = 9.3
# 6x3
e1 = as_matrix([[0.0, 0.0, e31],
[0.0, 0.0, e31],
[0.0, 0.0, e33],
[0.0, e15, 0.0],
[e15, 0.0, 0.0],
[0.0, 0.0, 0.0]])
# define permittivity tensor p
p11 = 5.64e-9
p33 = 6.35e-9
# 3x3
p = as_matrix([[p11, 0.0, 0.0],
[0.0, p11, 0.0],
[0.0, 0.0, p33]])
# rewrite the tensor into a vector using Voigt notation
def as_voigt_vector(ten):
# FEniCS does not know anything about Voigt notation,
# so one need to access the components directly as eps[0, 0] etc.
# return as_vector([ten[0, 0], ten[1, 1], ten[2, 2], ten[1, 2], ten[0, 2], ten[0, 1]])
return as_vector([ten[0, 0], ten[1, 1], ten[2, 2], 2 * ten[1, 2], 2 * ten[0, 2], 2 * ten[0, 1]])
def as_voigt_vector2(ten):
# FEniCS does not know anything about Voigt notation,
# so one need to access the components directly as eps[0, 0] etc.
return as_vector([ten[0], ten[1], ten[2]])
# rewrite a vector into a tensor using Voigt notation
def from_voigt_vector(vec):
# FEniCS does not know anything about Voigt notation,
# so one need to rewrite the vector by hand into a tensor.
return as_tensor([[vec[0], vec[5], vec[4]],
[vec[5], vec[1], vec[3]],
[vec[4], vec[3], vec[2]]])
def from_voigt_vector2(vec):
# FEniCS does not know anything about Voigt notation,
# so one need to rewrite the vector by hand into a tensor.
return as_tensor([vec[0], vec[1], vec[2]])
# first define the function for the strain tensor in fenics notation
def epsilon(u):
return 0.5 * (nabla_grad(u) + nabla_grad(u).T)
# return sym(nabla_grad(u))
# then define the function for the static electric field in fenics notation
def e_field(phi):
return -nabla_grad(phi)
# now define the stress tensor in fencis notation
def sigma(u, phi):
# calculate sigma with C-matrix in Voigt notation
sig = C * (as_voigt_vector(epsilon(u))) - e1 * as_voigt_vector2(e_field(phi))
# return sigma as symmetric tensor
return from_voigt_vector(sig)
# then define the function for the electric displacement in fenics notation
def e_dis(u, phi):
ed = e1.T * (as_voigt_vector(epsilon(u))) + p * as_voigt_vector2(e_field(phi))
return from_voigt_vector2(ed)
# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(1, 0.2, 0.22), 50, 10, 11)
class UpperBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 0.22)
#Marking facets that are on RightBoundary and creating the corresponding submesh
boundary_marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_marker.set_all(0)
UpperBoundary().mark(boundary_marker, 1)
submesh = MeshView.create(boundary_marker, 1)
ds = Measure("ds", subdomain_data=boundary_marker)
V1 = VectorFunctionSpace(mesh, 'CG', 1) #For the displacement in main part
V2 = VectorFunctionSpace(submesh, 'CG', 1) #For the displacement in electrode part
Q1 = FunctionSpace(mesh, 'CG', 1) #For the potential in main part
Q2 = FunctionSpace(submesh, 'R', 0) #For the potential in electrode part
L_M1 = FunctionSpace(submesh, 'CG', 1) #For the Lagrange multiplier
L_M2 = VectorFunctionSpace(submesh, 'CG', 1) #For the Lagrange multiplier
#V displacement, Q1 electric potential, Q2 electric potential on top surface
Z = MixedFunctionSpace(V1, V2, Q1, Q2, L_M1, L_M2)
sol = Function(Z)
(u1, u2, phi1, phi2, LM1, LM2) = TrialFunctions(Z)
(v1, v2, q1 , q2 , lm1, lm2) = TestFunctions(Z)
# Define boundary condition
tol = 1E-14
def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < tol
def potential(x, on_boundary):
return on_boundary and x[2] < tol
bc1 = DirichletBC(V1, Constant((0, 0, 0)), clamped_boundary)
bc2 = DirichletBC(Q1, Constant(0), potential)
bc = [bc1, bc2]
# Define variational problem
dV = Measure("dx", domain=Z.sub_space(0).mesh())
dL = Measure("dx", domain=Z.sub_space(3).mesh())
d = u1.geometric_dimension() # space dimension
f = Constant((0, 0, -rho * g))
T = Constant((0, 0, 0))
a = inner(sigma(u1, phi1), epsilon(v1)) * dx - dot(e_dis(u1, phi1), e_field(q1)) * dx
#a += inner(sigma(u2, phi2), epsilon(v2)) * ds(1) - dot(e_dis(u2, phi2), e_field(q2)) * ds(1)
#a += LM1*(phi1-phi2) * ds(1) + lm1*(q1-q2) * ds(1)
#a += dot(LM2,(u1-u2))* ds(1) + dot(lm2, (v1-v2))*ds(1)
a += inner(sigma(u2, phi2), epsilon(v2)) * dL - dot(e_dis(u2, phi2), e_field(q2)) * dL
a += LM1*(phi1-phi2) * dL + lm1*(q1-q2) * dL
a += dot(LM2,(u1-u2))* dL + dot(lm2, (v1-v2))* dL
L = dot(f, v1) * dV + dot(T, v1) * ds
#solve(a == L, sol, bc)
solve(a == L, sol, bc, solver_parameters={"linear_solver": "direct"})
# solve(a == L, sol, bc, solver_parameters={"linear_solver":"minres"})
u = sol.sub(0)
phi1 = sol.sub(2)
phi2 = sol.sub(3)
L_M_p = sol.sub(4)
# Plot solution
import matplotlib.pyplot as plt
# plot(u, title='Displacement', mode='displacement')
# Plot stress
s = sigma(u, phi1) - (1. / 3) * tr(sigma(u, phi1)) * Identity(d) # deviatoric stress
von_Mises = sqrt(3. / 2 * inner(s, s))
V = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V)
# plot(von_Mises, title='Stress intensity')
# Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
# plot(u_magnitude, title='Displacement magnitude')
# print('min/max u:',
# u_magnitude.vector().array().min(),
# u_magnitude.vector().array().max())
# Compute Electric field
M = VectorFunctionSpace(mesh, 'P', 1)
ef = project(e_field(phi1), M)
# Save solution to file in VTK format
File('R0_test1/R0_test1_Electrod_displacement.pvd') << u
# File('right100v/polar_von_mises.pvd') << von_Mises
# File('right100v/polar_magnitude.pvd') << u_magnitude
File('R0_test1/R0_test1_Electrod_electric_potential.pvd') << phi1
File('R0_test1/R0_test1_Electrod_electric_field.pvd') << ef
File('R0_test1/R0_test1_Electrod_L_M_p.pvd') << L_M_p
