I would like to model materials with anisotropy (stiffness and piezoelectric tensors). I'm trying to start out with a stressed cantilever given by the following code for an isotropic material but I cannot figure out how to adapt it for an anisotropic material. If I try to multiply the displacement vector (1 by 3) by a stiffness tensor reduced to a 6 by 6 matrix the dimensions won't work out and I can't find any way to get the shear displacements.
This runs but is isotropic
start of code..............................................................................................................
from dolfin import *
d0 = 5.0
d1 = 1.0
d2 = 0.05
n0 = int(d06)
n1 = int(d16)
n2 = int(d2*6)
Mesh for cantilever beam
mesh = BoxMesh(Point(0.0, 0.0, 0.0),Point(d0, d1, d2),50,20,2)#n0, n1, n2
#
V = VectorFunctionSpace ( mesh , "Lagrange", 2)
E = 35010^9
nu = 0.25
lmbda = E * nu / ((1.0 + nu ) * (1.0 - 2.0 nu ))
mu = E / (2 * (1 + nu ))
u = TrialFunction ( V )
v = TestFunction ( V )
class Fixed ( SubDomain ):
def inside ( self , x , on_boundary ):
return ( near ( x [0] , 0.0) and on_boundary )
fixed = Fixed ()
class Load ( SubDomain ):
def inside ( self , x , on_boundary ):
return (( x [0] >= 4.88) and near ( x [2] ,0.05)
and between ( x [1] ,(0.43 ,0.57)) and on_boundary )
load = Load ()
bc_fixed = DirichletBC (V , Constant ((0 ,0 ,0)) , fixed ) #fixes end of beam
bcs = [ bc_fixed ] #fixes end of beam, only one boundary condition
boundaries = FacetFunction ("size_t", mesh )
boundaries . set_all (0)
fixed . mark ( boundaries ,1)
load . mark ( boundaries ,2)
ds = ds ( domain = mesh , subdomain_data = boundaries )
tr_load = -96000
contact_area = assemble (1* ds (2))
g = Constant ((0 ,0 , tr_load / contact_area ))
def epsilon ( v ):
return 0.5 * ( grad ( v ) + transpose ( grad ( v )))
def sigma ( v ):
return 2 * mu * epsilon ( v ) + lmbda * tr ( epsilon ( v )) * Identity (3)
a = inner ( grad ( v ) , sigma ( u ))* dx
L = inner (v , g )* ds (2)
u = Function ( V )
problem = LinearVariationalProblem (a ,L ,u , bcs )
solver = LinearVariationalSolver ( problem )
solver.parameters ["linear_solver"] = "petsc"
solver.parameters ["preconditioner"] = "jacobi"
solver.solve ()
plot ( u )
plot ( mesh )
file = File ("beam3.pvd")
file << u
end of code..............................................................................................................
I'm new to Fenics but also have a feeling this is a very hard problem, if anyone can help find the shear displacements and add them into the regular displacement vector or find some other way to allow for an anisotropic stiffness tensor it would be greatly appreciated.
