Hello,
I’m currently try to solve the system of coupled PDEs. More precisely, it is the model of the thermoelastic euler-bernoulli beam defined on the 2-dimensional plate.
The heat distribution is described by the heat equation:
\rho c_p \frac{\partial T(x,z,t)}{\partial t} -\lambda \Delta T(x,z,t) = -kz \frac{\partial^{2} \partial u}{\partial x^{2} \partial t} \hspace{2cm} \ \, (x,y) \in (0,L) \times (-\frac{h}{2},\frac{h}{2})
BCs: \hspace{0.5cm}\frac{\partial T(0,z,t)}{\partial x} = \frac{\partial T(L,z,t)}{\partial x} = \frac{\partial T(x,h/2,t)}{\partial x} = 0, \ \, \, T(x,-h/2,t) = T_R+T_At
ICs: \hspace{0.75cm}T(x,z,0) = 0
The equation for the bending of the plate is giving as follows:
\rho A \frac{\partial^2 u (x,t)}{\partial t^2} + E I \Delta^2 u(x,t) + \gamma \frac{\partial^{2} \hat{T}}{\partial x^{2}} = 0 \hspace{0.5cm} x \in (0,L)
BCs: \hspace{0.5cm} u(0,t) = u(L,t) = 0, \ EI \Delta u(0,t) = EI\Delta u(L,t) = 0
ICs: \hspace{0.5cm} u(z,0) = 0, \hspace{0.5cm} \dot{u}(x,0) = 0 \hspace{5.2cm}
With integral kernel defined as follows \hat{T}(x,z,t) = \int_{-\frac{h}{2}}^{\frac{h}{2}}z T(x,z,t) \, dz .
Here are \rho, A, \gamma, \kappa, c_p, E, I, \lambda arbitrary constants.
Since I am quite new in the matter of discretization of PDEs with FEM. I follow the examples of the solving the heat equation https://fenicsproject.org/pub/tutorial/html/._ftut1006.html and the biharmonic equation https://fenicsproject.org/olddocs/dolfin/1.5.0/cpp/demo/documented/biharmonic/cpp/documentation.html.
Following this examples and using the central difference for time discretization I attempted to solve this problem by employing the (I think) equivalent weak formulation for the heat equation:
a(u,w) = \rho c_p \int_{\Omega}T^{n+1}w \, d \Omega + dt \lambda \int_{\Omega} (\frac{\partial T^{n+1}}{\partial x} + \frac{\partial T^{n+1}}{\partial z})\frac{\partial w}{\partial x} d\Omega - \kappa z \int_{\Omega} \frac{\partial u^{n+1}}{\partial x} \frac{\partial w}{\partial x} d \Omega
L(v) = \rho c_p \int_{\Omega} T^n w \, d \Omega - \kappa z \int_{ \Omega}\frac{\partial u^n}{\partial x} \frac{\partial w}{\partial x} \, d \Omega
and also for the beam equation:
a(u,w) = \rho A \int_{0}^{L} u^{n+1}v \, dx + dt^2 E I \int_{0}^{L} \frac{\partial^2 u^{n+1}}{\partial x^2} \frac{\partial^2 v}{\partial x^2} \,dx - dt^2 \gamma \int_{0}^{L} \frac{\partial \hat{T}^{n+1}}{\partial x} \frac{\partial v}{\partial x} \, dx
L(v) = 2 \rho A \int_{0}^{L} u^n v\, dx - \rho A \int_{0}^{L} u^{n-1} v \,dx
Where (u^{n+1},T^{n+1}) are the searched results and the (w,v) are the test functions. Since fenics does not support a C^1 finite elements the expression
the second order derivative expression \int_{0}^{L} \frac{\partial^2 u^{n+1}}{\partial x^2} \frac{\partial^2 v}{\partial x^2} \,dx is implements using the C^0 Interior Penalty Method.
Due to the fact that the heat equation is defined on 2D and the beam equation on 1D domain I using the mixed-dimensional branch by cecile for coupling fenics-project / DOLFIN / cecile/mixed-dimensional — Bitbucket.
I attempted to implement the coupled problem, unfortunately I encountered a few problems:
-
At first the implementation of the euler-bernoulli equation with the C⁰ Interior Point Method causes an error when assembling. I think the problem lies in the implementation of the stiffnes matrix with C⁰IP Method.
-
My second problem is the formulation of the integral therm \hat{T}(x,z,t) = \int_{-\frac{h}{2}}^{\frac{h}{2}}z T(x,z,t) \, dz in fenics. In the code example i have left it blank. Does anyone have an idea how to implement it?
I’m using the docker container of the mixed-dimensional branch:latest on Linux Ubuntu.
I hope I have not made too many mistakes and expressed it understandably. For any help I am very grateful, thank you.
from fenics import *
###---------------- define mesh --------------###
N = 30
order = 2
x_start = 0
x_end = 1
z_start = -0.1
z_end = 0.1
p1 = Point(x_start,z_start)
p2 = Point(x_end,z_end)
mesh = RectangleMesh(p1, p2, N, N)
marker = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
for e in edges(mesh):
marker[e] = 0.0 - DOLFIN_EPS < e.midpoint().y() < 0.0 + DOLFIN_EPS
submesh = MeshView.create(marker, 1)
# define function space for the temperature and displacement
Vt = FunctionSpace(mesh, "CG", 1) ## 2D
Vu = FunctionSpace(submesh, "CG", 1) ## 1D
# Define the mixed function space
V = MixedFunctionSpace(Vt,Vu)
# Redefine integration measure for the beam
dxB = Measure('dx', domain=V.sub_space(1).mesh())
# Redefine interior facets integration measure for the beam
dSB = Measure('dS', domain=V.sub_space(1).mesh())
# define the Dirichlet boundary condition for the heat equation
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],z_start) and on_boundary
source = Expression("T_R + T_A*t", T_R=0, T_A=8.0, omega=pi, t=0.0, degree =2)
# impose the dirichlet boundary condition
bc1_t = DirichletBC(V.sub_space(0), source , Bottom())
bcs_t = [bc1_t]
## impose dirichlet BC's for the beam subdomain
# define the boundaries of the beam
class BeamLeft(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],x_start)
class BeamRight(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],x_end)
bc1_b = DirichletBC(V.sub_space(1), Constant(0.0), BeamLeft())
bc2_b = DirichletBC(V.sub_space(1), Constant(0.0), BeamRight())
bcs_b = [bc1_b, bc2_b]
bcs = [bc1_t, bc1_b, bc2_b]
# time parameters
Period = 20
Nsteps = 100
dt = Constant(Period/Nsteps)
# define the test and trialfunction
(w,v) = TestFunctions(V)
(T,u) = TrialFunctions(V)
# intialize the initial temperature
T_n = Function(V.sub_space(0))
# initialize the intial deflection u_n
u_n = Function(V.sub_space(1))
u_prev = Function(V.sub_space(1))
# weak formulation of the heat equation
a_heat = inner(T,w) * dx + dt * inner(nabla_grad(T), nabla_grad(w)) * dx - inner(grad(u),grad(v))*dxB
L_heat = dt * inner(nabla_grad(T_n),nabla_grad(w)) * dx + inner(div(grad(u_n)), v)*dxB
Th = Function(V, name = "Temperature")
T_hat = Function(V.sub_space(1))
#formulation of the weak form with C⁰IPM
# stiffness form
def k(u, v):
# Define normal component, mesh size and right-hand side
h = CellDiameter(submesh)
h_avg = (h('+') + h('-')) / 2.0
n = FacetNormal(submesh)
# Penalty parameter
sigma = Constant(8.0)
ex = inner(div(grad(u)), div(grad(v))) * dxB \
- inner(avg(div(grad(u))), jump(grad(v), n)) * dSB \
- inner(jump(grad(u), n), avg(div(grad(v)))) * dSB \
+ sigma('+') / h_avg * inner(jump(grad(u), n), jump(grad(v), n)) * dSB
return ex
# weak formulation for beam equation
a_beam = inner(u,v)*dxB + dt*dt*k(u,v) #+ \int_{-\frac{h}{2}}^{\frac{h}{2}}z \frac{\partial^{2} T}{\partial x^{2}} \text{dz}
L_beam = 2.0 * inner(u_n,v) * dxB - inner(u_prev,v) * dxB
# merge the heat and beam forms
a_form = a_heat + a_beam
L_form = L_heat + L_beam
# define the solution variable
U = Function(V)
# solve the heat equation for the temperature
solve(a_form==L_form, U, bcs)
And the Error:
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to define mixed linear variational problem a(u, v) == L(v) for all v.
*** Reason: Expecting the left-hand side (block 0, domain 0) to be a bilinear form (not rank 0).
*** Where: This error was encountered inside MixedLinearVariationalProblem.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: ed85093177b62be6ae55401b4ec808fa32c0a3f8
*** -------------------------------------------------------------------------