Polar - axisymetric poisson formulation with PointSource at origin

variational formulation
1

Hi there,

I would like to understand if my math notation is correct for the problem I am tring to solve:

I have a steady heat source q (W/2 π) at the origin and a circular no-flow boundary at distance R from the origin. The material in the domain has a thermal conductivity K (W/mK) and a volumetric heat capacity C (J/m3/K). I want to model the temperature varying over time within the domain.When I took integration by parts of poissons equation, I set the boundary integral to zero, because of the no flow condition at R.

Therefore, I presume i’m looking to find the solution u®, in the variational formulation:

Does this make sense? I’m not sure I have set boundary conditions correctly. I want neumman boundaries but I don’t think I’m actually dictating anything here. When I try to solve for initial conditions I make the following formulation:

ql_0fde99454c934476a427191643b97192_l3

Any insight you could provide would be very helpful. So far in my code I have, just for a single time step:

from dolfin import *

T = 300 # final time
num_steps = 60 # number of time steps
dt = T / num_steps
K = 2 #w/mk
C = 2.5 # J/m3/k
Q =1.27/6.28

flow per metre is

mesh = IntervalMesh(10,0,0.2)

V = FunctionSpace(mesh, ‘P’, 1)
u = TrialFunction(V)
v = TestFunction(V)
r = Expression(‘x[0]’, degree=1)

u_n = Expression(‘0’, degree=1)
a = Cinner(u,v)rdx- dtK*(inner(grad(u),grad(v)))rdx
L = dtConstant(0)vrdx + inner(u_n,v)rdx
A, b = assemble_system(a, L)
delta = PointSource(V, Point(0.), Q)
delta.apply(b)
u = Function(V)
solve(A, u.vector(), b)

the values of u are much higher than I was expecting - any help would be appreciated.

Many thanks
C

3 
from dolfin import *

dt = 5 # seconds
D = 0.0000008 # diffusivity wm2/J
Q =1/6.28 #w/m/2pi

mesh = IntervalMesh(100,0,1)

V = FunctionSpace(mesh, "Lagrange", 3)
u = TrialFunction(V)
v = TestFunction(V)
r = Expression("x[0]", degree=1)
u_0 = Constant(0)
u_old = interpolate(u_0, V)

a = dt*D*(inner(grad(u),grad(v))*r*dx   + inner(u,v)*r*dx 
L = (dt*inner(Constant(0),v))*r*dx + inner(u_old,v)*r*dx

A = assemble(a)
u_new = Function(V)

for n in range(10):
    b = assemble(L)
    delta = PointSource(V, Point(1), Q)
    delta.apply(b)      
    solve(A, u_new.vector(), b)
    u_old.assign(u_new)
    print(u_new.vector()[0])```
4

Hi All,
Me again. I’m trying to replicate a needle probe used to determine material conductivities. Conductivity is linearly proportional to slope of the change in temperature plotted against the natural log of time.

Something is not right here. I am getting extremely high temperatures at the core even though I’m passing one watt of power in. Is it because it’s evaluating the delta function at the origin?

Any thoughts are welcome. I’m trying to replicate an idealised needle in a core of material.

from dolfin import *
import matplotlib.pyplot as plt
dt = 5 # seconds
k =2.0 # trial conductivity 
C = 2500000 # volumetric heat capacity
D = k/C # diffusivity m2/s
Q =1/6.28 #w/m/2pi

mesh = IntervalMesh(1000,0.000,1)

V = FunctionSpace(mesh, "Lagrange", 2)
u = TrialFunction(V)
v = TestFunction(V)
r = Expression("x[0]", degree=1)
u_0 = Constant(0)
u_old = interpolate(u_0, V)

a = dt*D*inner(grad(u),grad(v))*r*dx   + inner(u,v)*r*dx 
L = dt*inner(Constant(0),v)*r*dx + inner(u_old,v)*r*dx

A = assemble(a)
u_new = Function(V)
X = []

for n in range(60):
    b = assemble(L)
    delta = PointSource(V, Point(1), Q)
    delta.apply(b)      
    solve(A, u_new.vector(), b)
    print(u_new.vector()[0])  
    u_old.assign(u_new)```
5

I would recommend having a look at this 1 D axisymmetric heat conduction with neumann boundaries and the link within thus post

6

Thank you Dokken,

I believe my formulation is following that shown in the previous post. So I need something new. I am also using the weak formualtion for polar representation in the post linked to the previous post.

I notice that as I increase the mesh intervals, The value of U at r=0 increases. Is this to be expected? I’ve just noticed that the Jaeger and Carslaw solution for the line source uses exponential integrals. Maybe, I can’t actually evaulate the true value of U at r=0 then, and should focus on the changes in U, relation to time.

Would I be better representing the heat source as a constant line load over the first element? Eventually I want to distinguish the conductivities between the material being tested and the metal heating element.

Any more thoughts would be appreciated.

Thanks
C

7

To use discontinuous material parameters, have a look at either How to define different materials / importet 3D geometry (gmsh) or How to define boundary conditions on an irregular geometry?,

8

Thanks Dokken,

Can you explain why when I add the pointsource at Point(1), the spike in temperature is seen at r = 0? I tried placing the pointsource- at point(0) but then the heat seemed to flowing from the the opposite end.

9

To me, it seems strange that all your terms are scaled with r, i think there is a mistake in your derivation. Please do integration by parts in polar coordinates carefully.

10

I think my formulation is ok. I’m trying to represent the flow over a unit sector over a circle. r is multiplying d_theta which is one in this case. I’v divided the flow in at r =0, by two pi to get the amount of heat flow into the unit sector of the circle. I’m simplifying a 2d problem into a 1D problem, as heat flow will not vary with theta.