HI,
I am solving a FSI problem. I want to vary the x and y components of displacement in the solid (rectangle) domain linearly. So something like d^2u_y/dy^2=0 and d^2u_x/dy^2=0. Is there a way to do that? If so, please let me know.
Thank you!
#Define Subdomain
class SolidDomainClosure(SubDomain):
def inside(self, x, on_boundary):
return (x[1] > SOLID_BOTTOM - DOLFIN_EPS
and x[0] > SOLID_LEFT - DOLFIN_EPS
and x[0] < SOLID_RIGHT + DOLFIN_EPS)
#Defined Function Space
cell = mesh.ufl_cell()
Ve = VectorElement(“CG”, cell, 1)
Qe = FiniteElement(“CG”, cell, 1)
VQe = MixedElement((Ve,Qe))
W = FunctionSpace(mesh,VQe)
V = FunctionSpace(mesh,Ve)
Mesh motion functions:
uhat = Function(V)
uhat_old = Function(V)
du = TestFunction(V)
vhat = (uhat-uhat_old)/Dt
Fluid–structure functions:
(dv, dp) = TestFunctions(W)
w = Function(W)
v,p = split(w)
p_mean = Function(W)
w_old = Function(W)
v_old, p_old = split(w_old)
dv_dr = (v - v_old)/Dt
u = uhat_old + Dt*v
This will need to be updated to match u, for purposes of setting the boundary condition on the mesh motion subproblem.
u_func = Function(V)
bc_m_struct = DirichletBC(V, u_func, SolidDomainClosure())
u_func.assign(project(u,V))
####### Formulation of mesh motion subproblem #######
Residual for mesh, which satisfies a fictitious elastic problem:
F_m = grad_y(uhat) + I
E_m = 0.5*(F_m.T*F_m - I)
m_jac_stiff_pow = Constant(3)
Artificially stiffen the mesh where it is getting crushed:
K_m = Constant(1)/pow(det(F_m),m_jac_stiff_pow)
mu_m = Constant(1)/pow(det(F_m),m_jac_stiff_pow)
S_m = K_mtr(E_m)I + 2.0mu_m(E_m - tr(E_m)I/3.0)
res_m = (inner(F_mS_m,grad_y(du)))*dy
Dres_m = derivative(res_m, uhat)