Hello Everyone,
I want to study piezoelectric effect in a 2D beam. The boundary conditions are: 1) the left most bottom corner of the beam is fixed in x and y direction and the rest of the nodes on the bottom edge are restricted only in y-direction 2) load is applied to its top edge 3) potential of 0.0 is given for the bottom edge 4) potential of 0.2 applied to its top edge
The variational formulation is as follows:
where:
fB = body forces
fN = traction forces
qB = volume electric charge given by Div(eletric_displacement)
qN = surface electric charge
My implementation is as follows:
from __future__ import print_function
from dolfin import *
mesh = RectangleMesh(Point(0.0, 0.0), Point(50.0,5.0), 50,10,"crossed")
piezo_Mat=as_matrix([[0.,0.,0.,0.,15,0.],
[0.,0.,0.,15,0.,0.],
[30,30,35,0.,0.,0.]])
permit_Mat=as_matrix([[0.5,0.,0.],
[0.,0.5,0.],
[0.,0.,0.6]])
V_mech=VectorElement("CG",mesh.ufl_cell(),2)
V_elect=FiniteElement("CG",mesh.ufl_cell(),2)
Mix_ele=MixedElement([V_mech,V_elect])
Mixed_FS=FunctionSpace(mesh,Mix_ele)
(u,phi)=TrialFunctions(Mixed_FS)
(v,si) = TestFunctions(Mixed_FS)
V_ele=FunctionSpace(mesh,"DG",0)
tol=1E-14
def bottom_fixed_boundary(x,on_boundary):
return near(x[0],0,tol) and near(x[1],0,tol)
Fixed_left=Constant((0.,0.))
bc_Fixed=DirichletBC(Mixed_FS.sub(0),Fixed_left,bottom_fixed_boundary,method='pointwise')
def bottom_right_boundary(x,on_boundary):
return near(x[1],0.0,tol) and x[0]>0.0
Roller_right=Constant(0)
bc_roller=DirichletBC(Mixed_FS.sub(0).sub(1),Roller_right,bottom_right_boundary)
ground_phi=Constant(0)
bc_bottom_ground=DirichletBC(Mixed_FS.sub(1),ground_phi,bottom_fixed_boundary)
bcs_all=[bc_Fixed,bc_roller,bc_bottom_ground]
nu=0.3
E0=3000
alpha_val=100
rho0=10
rho_val=[]
for cell_s in cells(mesh):
rho_val.append(rho0)
class Top(SubDomain):
tol=1E-14
def inside(self,x,on_boundary):
return near(x[1],5.0,tol)
top=Top()
boundries=MeshFunction('size_t',mesh,1)
boundries.set_all(0)
top.mark(boundries,1)
ds=Measure('ds',domain=mesh,subdomain_data=boundries)
def calculate_E(rho):
E_updated=Function(V_ele)
E_array=E_updated.vector().get_local()
for i, rho in enumerate (rho):
E_array[i]=E0
E_updated.vector().set_local(E_array)
return E_updated
E0=calculate_E(rho_val)
def calculate_alpha(rho_values):
alpha=Function(V_ele)
alpha_array=alpha.vector().get_local()
for i, rho in enumerate (rho_values):
alpha_array[i]=alpha_val
alpha.vector().set_local(alpha_array)
return alpha
alpha0=calculate_alpha(rho_val)
def calculate_Lame_coefficients(E_val):
mu_val=(E_val)/(2*(1+nu))
lmbda_val=(E_val*nu)/((1+nu)*(1-2*(nu)))
return mu_val,lmbda_val
mu, lmbda=calculate_Lame_coefficients(E0)
def epsilon(u):
strain_u=0.5*(grad(u)+grad(u).T)
return strain_u
def strain_voigt(u):
strain_3_1=as_vector([u[i].dx(i) for i in range(2)] +[u[i].dx(j) + u[j].dx(i) for (i,j) in [(0,1)]])
strain_6_1=as_vector([strain_3_1[0], strain_3_1[1], 0, 0, 0, strain_3_1[2]])
return strain_6_1
def sigma(u,mu,lmbda):
eps=epsilon(u)
stress_u=lmbda*div(u)*Identity(d)+2*mu*eps
return stress_u
def elect_field_E(phi):
E_field_val=-grad(phi)
E_field_val_3_1=as_vector([E_field_val[0], E_field_val[1], 0])
return E_field_val_3_1
def elect_disp_D(u,phi):
elect_disp=alpha0*((piezo_Mat*strain_voigt(u)) + (permit_Mat*elect_field_E(phi)))
return elect_disp
fB=Constant((0.0,0.0))
fN=Constant(10.0)
qB=elect_disp_D(v,si)
qN=Constant(0.2)
d=u.geometric_dimension()
a = 2*mu*inner(epsilon(u),epsilon(v))*dx + lmbda*dot(tr(epsilon(u)),tr(epsilon(v)))*dx + alpha0*div((permit_Mat*elect_field_E(phi))-(piezo_Mat*strain_voigt(u)))*div(si)*dx
L=dot(fB,v)*dx + v[1]*fN*ds(1) + div(qB)*si*dx + qN*si*ds(1)
w=Function(Mixed_FS)
(u,phi)= split(w)
solve(a==L,w,bcs_all)
(u,phi)=w.split(True)
print("Displacement: ",u.vector().get_local())
print("Potential: ",phi.vector().get_local())
I am getting an error that says “Nonmatching dimensions for free indices with same ids!”
I am not sure if my implementation of the variation formulations is right and what is causing the error. Can I please get some hint or guidance?
Thank you very much.