Hello
I have written a code in FEniCS for the simulation of the deformation(deflection) of a clamped rectangular body in a particular problem and this deformation is time-dependent(transient).
This code works perfectly right now. Here is my complete code:
from dolfin import *
import mshr
import numpy as np
#Diffusion constant
D = 10 * (10**(-11))
#The gas constant
R = 8.31
#Electrical mobility
mu = 4.11 * (10**(-14.))
#Anion concentration
C_0 = 1200
#Dielectric permittivity
epsilon = 0.025
#Temperature
T = D / (mu*R)
#Faraday number
F = 96485
#length scale
l = 200 * (10**(-6))
#Debye screening length
lambda_d = ((epsilon*R*T)/(2*pow(F,2)*C_0))**(0.5)
#Linear force coupling
A = 0.001
eps = lambda_d / l
#Material properties
E11=10.0
E22=E11
NU12=0.3
G12=E11 / (2*(1+NU12))
NU21=NU12
#Matrix properties
Q11=E11/(1-NU12*NU21)
Q22=E22/(1-NU12*NU21)
Q12=(NU12*E22)/(1-NU12*NU21)
Q66=G12
Q_numpy = np.array([[Q11, Q12, 0.0],
[Q12, Q22, 0.0],
[0.0, 0.0, Q66 ]])
Q = as_matrix(Q_numpy)
# Define boundaries
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 4.0)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
# Domain
domain = mshr.Rectangle(Point(0,0), Point(4,1))
# Mesh
mesh = mshr.generate_mesh(domain, 60)
# Marking boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 2)
top.mark(boundaries, 1)
right.mark(boundaries, 4)
bottom.mark(boundaries, 3)
#Defining the function spaces for the concentration
V_c = FunctionSpace(mesh, 'CG', 1)
#Defining the function spaces for the voltage
V_phi = FunctionSpace(mesh, 'CG', 1)
#Defining the vectorial function spaces for the displacement
V_disp = VectorFunctionSpace(mesh, 'CG', 1)
#Defining the mixed function space
W = MixedFunctionSpace((V_c, V_phi, V_disp))
#Defining the "Trial" functions
z = Function(W)
c, phi, u = split(z)
#Defining the test functions
(v_1, v_2, v_3) = TestFunctions(W)
# Time variables
dt = 0.05; t = 0; T = 3
# Previous solution
W_const = Constant(1.0)
C_previous= interpolate(W_const, V_c)
# Define Dirichlet boundary conditions at top boundary (V=4.0 volt)
Voltage = Constant(4.0)
bc_top = DirichletBC(W.sub(1), Voltage, boundaries, 1)
# Define Dirichlet boundary conditions at bottom boundary (V=0.0 volt)
bc_bottom = DirichletBC(W.sub(1), 0.0, boundaries, 3)
# Define Dirichlet boundary conditions at the left boundary (Clamped B.C)
U_const = Constant((0.0, 0.0))
bc_left = DirichletBC(W.sub(2), U_const, boundaries, 2)
#Collecting the boundary conditions
bcs = [bc_top, bc_bottom, bc_left]
def epsi(u):
return 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)]])
def sigma(u):
return Q * epsi(u)
# Define vectorial body force
F1 = Constant(A*C_0)*(c - Constant(1.))
# Variational form For nonlinear solver
F = dot(grad(phi), grad(v_2)) * dx \
+ c*v_1 * dx \
- (1./(2 * (pow(eps,2)))) * c * v_2 * dx \
+ dt * eps * dot(grad(c),grad(v_1)) * dx \
+ dt * eps * c * dot(grad(phi),grad(v_1)) * dx\
-C_previous * v_1 * dx + (1./(2 * (pow(eps,2)))) * v_2 * dx
-inner(sigma(u),epsi(v_3)) * dx + F1*v_3[0]*dx
f = File("pion.pvd")
while t <= T:
solve(F == 0, z, bcs)
(c, phi,u) = z.split(True)
C_previous.assign(c)
t += dt
plot(u,key="c", interactive=False) #plot displacement
f << u
I have taken some screen shots of the deformation of the domain at different time steps:
Now I want to assume that this rectangular domain is inside a liquid (for example a squared-shape container of water).
This is the original domain in my code:
This is the new domain(s) I am looking at:
In the above figure, the blue domain corresponds to the water surrounding the central rectangle (green domain). I need to mention that the water does not affect the deflection of the green domain.
I guess I need to use moving mesh for this purpose but I do not know how to use it.
I want to modify my original code to add the blue domain in a way that as the green domain is deflected upward (as shown in the first figure), the blue elements near the boundaries of the green domain move or get created again along with that during the time.
Could you please help me modify the domain(s) and mesh in my code to simulate what I am looking for?
Thanks in advance for your help.
