Thanks for the suggestion, I tried to follow your supply code and i have implemented the above equations, But I get an error stating that :
*** Error: Unable to assemble system.
*** Reason: expected a bilinear form for a.
*** Where: This error was encountered inside SystemAssembler.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
I have attached my code below, It would be appreciable to point out the mistake:
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import timeit
#units: m, kg, s, V, K
################################## Mesh part ##########################
“”"
$PhysicalNames
2D 2 “top”
2D 3 “bottom”
3D 1 “cube”
“”"
load the mesh data from Gmsh
mesh_file = ‘cube.h5’
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), mesh_file, ‘r’)
hdf.read(mesh, ‘/mesh’, False)
cells = MeshFunction(‘size_t’, mesh, dim=3)
hdf.read(cells, ‘/subdomains’)
facets = MeshFunction(‘size_t’, mesh, dim=2)
hdf.read(facets, ‘/facets’)
JZ: use the Measure function
dx = Measure(‘dx’, domain=mesh, subdomain_data=cells)
ds = Measure(‘ds’, domain=mesh, subdomain_data=facets)
Function spaces
N = FacetNormal(mesh)
Ue = VectorElement(“CG”, mesh.ufl_cell(), 2) # displacement vector element
Ve = FiniteElement(“CG”, mesh.ufl_cell(), 1) # voltage finite element
#u,phi
W = FunctionSpace(mesh, MixedElement([Ue, Ve]))
Parameters of Piezomaterial (PIC-255)
--------------------
c_11 = 1.23e11
c_12 = 7.67e10
c_13 = 7.025e10
c_33 = 9.711e10
c_44 = 2.226e10
eps_11 = 8.234
eps_33 = 7.588
e_15 = 11.91
e_13 = -7.15
e_33 = 13.75
rho_Piezo = Constant(7800) #in kg/m3
Order of matrices similar to ANSYS notation: x, y, z, xy, yz,xz
IEEE notation: x, y, z, yz, xz, xy
Elasticity tensor c^E
C = as_tensor([[c_11, c_12, c_13, 0., 0., 0.],
[c_12, c_11, c_13, 0., 0., 0.],
[c_13, c_13, c_33, 0., 0., 0.],
[0., 0., 0., c_44, 0., 0],
[0., 0., 0., 0., c_44, 0.],
[0., 0., 0., 0., 0., (c_11-c_12)/2]])
piezoelectric coupling tensor e
et_piezo = as_tensor([[0. , 0. , 0., 0. , 0., e_15],
[0., 0., 0., 0., e_15, 0.],
[e_13, e_13, e_33, 0., 0., 0.]])
transpose form of piezo tensor
e_piezo = as_tensor([[0., 0., e_13],
[0., 0., e_13],
[0., 0., e_33],
[0., 0, 0.],
[0., e_15, 0.],
[e_15, 0., 0.]])
Permittivitatstensor epsilon^S
eps_s = as_tensor([[eps_11, 0., 0.],
[0., eps_11, 0.],
[0., 0., eps_33]])
rewrite the tensor into a vector using Voigt notation
def strain3voigt(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[0,1],ten[1,2],ten[0,2]])
rewrite a vector into a tensor using Voigt notation
def voigt3stress(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[3], vec[5]],
[vec[3], vec[1], vec[4]],
[vec[5], vec[4], vec[2]]])
first define the function for the strain tensor in fenics notation
def B(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#define the electric field vector
def E_field_vector(Ei):
return as_vector([Ei[0], Ei[1], Ei[2]])
#define the electric field in relation to potential
def E_field(v):
return -grad(v)
#definition of the different components of my Stress tensor equation (sigma=C^E:Bu-e_piezo.E)
def sigma(u, v):
return voigt3stress(dot(C, strain3voigt(B(u))) + e_piezo * E_field_vector(E_field(v)))
the electric displacements vector
def disp_D(Di):
return as_vector([Di[0], Di[1], Di[2]])
definition of the different components of my electric displacements equation (D=e^s:Bu+eps^S.E)
elect_disp_D for coupled linear elasicity and electrostatics
def D_field(u, v):
D = et_piezo * strain3voigt(B(u)) - eps_s * E_field(v)
return disp_D(D)
#Defining the “Trial” functions
w = Function(W)
(U, V) = split(w)
#Defining the test functions
(utest, vtest) = TestFunctions(W)
dirichlet boundary condition for fixed support
bcs = [DirichletBC(W.sub(0), Constant((0.,0.,0.)), facets, 2)]
Forms & matrices
U_form = inner(sigma(U, V), B(utest)) * dx
V_form = inner(D_field(U,V), grad(vtest)) * dx
K_form = U_form + V_form
M_form = rho_Piezo*inner(U, utest)*dx
L = utest[0]*dx
K, M = PETScMatrix(), PETScMatrix()
b = PETScVector()
assemble_system(K_form,L,bcs,A_tensor=K)
assemble_system(M_form,L,bcs,A_tensor=M)
[bc.zero(M) for bc in bcs]
--------------------Eigen-solver--------------------
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters[‘problem_type’] = ‘gen_hermitian’
eigensolver.parameters[‘spectral_transform’] = ‘shift-and-invert’
eigensolver.parameters[‘spectral_shift’] = 0.
N_eig = 6 # number of eigenvalues
print(“Computing {} first eigenvalues…”.format(N_eig))
eigensolver.solve(N_eig)
