Dear All
I am trying to solve the heat-transfer problem in a spherical annulus, by assuming that the temperature varies only along the radial direction, and has neither \theta nor \phi dependence. I am considering Dirichlet boundary conditions, with one end maintained at a dimensionless temperature of 1.0, and the other end at 0.0. I will first lay out the mathematical setup before presenting the code I wrote and the error which I encounter upon executing the code.
1. Governing equation and variational formulation
The governing equation is:
\nabla^2 u=0;\,u(r=r_{\text{i}})=1.0,\,u(r=r_{\text{o}})=0.0
with
\nabla^2 u = \dfrac{1}{r^2}\dfrac{\partial}{\partial r}\left(r^2\dfrac{\partial u}{\partial r}\right)+\dfrac{1}{r^2\sin\theta}\dfrac{\partial}{\partial \theta}\left(\sin\theta\dfrac{\partial u}{\partial \theta}\right)
The \theta-dependence is shown here only to indicate the choice of the integration domain, which is d\Omega=rdrd\theta. I will henceforth drop the \theta-dependence and write the variational form as:
\langle\nabla^2u,v\rangle_{\Omega}=0 where \langle\nabla^2u,v\rangle_{\Omega}=\int_{r=r_i}^{r=r_o}\left[ \dfrac{1}{r^2}\dfrac{d}{d r}\left(r^2\dfrac{d u}{d r}\right)\right]v(rdr)
which may be simplified, after integration by parts,
\int_{r=r_i}^{r=r_o}\left[ \dfrac{1}{r^2}\dfrac{d}{d r}\left(r^2\dfrac{d u}{d r}\right)\right]v(rdr)=\int\left(\dfrac{v}{r}\right)\dfrac{d}{d r}\left(r^2\dfrac{d u}{d r}\right)=\cancel{vr\dfrac{du}{dr}\Biggr|_{r_i}^{r_o}}-\int\left[\dfrac{d}{dr}\left(\dfrac{v}{r}\right)\right]r^2\dfrac{du}{dr}dr
to give
\int_{r=r_i}^{r=r_o}\left(r\dfrac{dv}{dr}-v\right)\left(\dfrac{du}{dr}\right)dr=0
as the variational problem to be solved.
2. FEniCS script that produces shape-related error upon execution
I am assuming that I have obtained the variational form correctly. Here is the FEniCS script that I am using to implement the math:
from dolfin import *
from mshr import *
### defining mesh
r_i=1 # radius of inner sphere
r_o=5 # radius of outer sphere
mesh = IntervalMesh(20, r_i, r_o)
r = SpatialCoordinate(mesh)
### defining Function space on this mesh using Lagrange polynoimals of degree 2.
V = FunctionSpace(mesh, "CG", 2)
### defining boundary values
u_i = Constant(1)
u_o = Constant(0)
### this functions checks whether the input x is on the boundary or not.
def inner_boundary(x, on_boundary):
tol = 1e-14
return on_boundary and near(x[0],r_i,tol)
def outer_boundary(x, on_boundary):
tol = 1e-14
return on_boundary and near(x[0],r_o,tol)
### Enforcing BCs at inner and outer ends
inner_bc = DirichletBC(V, u_i, inner_boundary)
outer_bc = DirichletBC(V, u_o, outer_boundary)
bcs = [inner_bc, outer_bc]
### Setting up the variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = (r*v.dx(0)-v)*u.dx(0)*dx
L = f*v*dx
### solving the variational problem.
u = Function(V)
solve( a == L, u, bcs)
When I execute the above code, I obtain the error:
Can’t add expressions with different shapes.
Traceback (most recent call last):
File “spher_1d_radial.py”, line 40, in
a = (r*v.dx(0)-v)*u.dx(0)*dx
Could you please help me with the error? I am new to FEniCS and would highly appreciate the help.
Thanks in advance