# 1D heat-transfer in spherical annulus, encountering "Can't add expressions with different shapes"

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

The problem is that UFL distinguishes between vectors of dimension one and scalars. SpatialCoordinate always returns a vector, even in 1D, where it is a vector with only one component. Thus, in your MWE, r*v.dx(0) is a vector, while v is a scalar, and it is invalid to subtract a scalar from a vector. You can fix this simply, by defining r as the 0-th (and only) component of SpatialCoordinate(mesh):

#r = SpatialCoordinate(mesh)
r = SpatialCoordinate(mesh)[0]

1 Like

Thank you, Prof. Kamensky, for the prompt response. Your suggested solution works perfectly. I will apply this knowledge when trying to solve the same problem in 2D (axisymmetric spherical coordinates).

Following your explanation, I experimented a bit and found that

r = SpatialCoordinate(mesh)
a = (r[0]*v.dx(0)-v)*u.dx(0)*dx

also works.