This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

differential of the unknown variable

0 votes

Hi,

I am working on a 3D scalar problem (temperature) where I have a flux applied to a surface. The problem is non-linear because the flux is dependent on the differential of the scalar unknown. Specifically, the flux term is equal to something like (w is the trial function):

$$a \int w dT$$

When I was working on the 2D (axisymmetric) version of my problem I was able to properly treat the integral by saying

$$dT = \frac{\partial T}{\partial z}dz + \frac{\partial T}{\partial x} dx$$

and I term in the x direction is zero because the flux only existed on $z$ axis and $ds = dz$ along the side of the applied flux. However, if I apply this approach to my current situation, I'm stuck because the differential of surface area $ds$ is dependent on all three directions. What is the best approach?

My code for the two 2D version was:

w = TestFunction(V)
L = cp_f * w * grad(T_n)[1] * ds(1)
T = Function(V)
solve(a == L, T, bc)

asked May 2, 2017 by dbfox FEniCS Novice (240 points)
edited May 2, 2017 by dbfox

Hi!

FEniCS Q&A is parsing LaTeX expressions through MathJax. For the first point I would recommend you to rewrite the math in the question into LaTeX, because I cannot understand it (and I wonder if anyone else can).

Second point: Include minimal working example of FEniCS code, or at least, some of your attempts to implement the problem into the FEniCS.

Thanks for the reply. I was confused on how to type mathematical questions, I didn't know you could simply use "$" and "$$". However, there seems to be an erroneous "|" character added at the end of each math mode expression.

I don't have a working example for my 3D implementation but now I've added what the flux term was when I was able to make the problem work in 2D.

Why can't you use a simple chain rule? Say, you have $dT = \partial_s T ds $ and then you are only missing the derivative with respect to $s$, which is $\partial_s T = \nabla T\cdot \partial_s x$ where $x$ are your coordinates. That depends on how you parameterized your surface, but it is something you can find. Let me know, best regards!

Thanks for the suggestion; I see your train of thinking as I can write
$$\int dT = \int \nabla T \cdot \frac{\partial x}{\partial S} dS$$.

In FEniCS the term becomes:

L = cp_f*w*dot(grad(T_n), dxds)*ds(1)

Now it seems like the trick is to come up with an expression for $\partial x/\partial s$. How do I go about getting this term? My geometry is a cube with a cylindrical hole in the center where the flux is being applied.

Then you can take 8 different domains, one for each face and two for the hole, and in those you know how the surface is parameterized:

$$ S_1(x) = x_1 $$
$$ S_2(x) = -x_1 $$
$$ S_3(x) = x_2 $$
etc
$$ S_7 = sqrt(1 - x_{-3}\cdot x_{-3}) $$
$$ S_8 = -sqrt(1 - x_{-3}\cdot x_{-3}) $$

From there you get $\partial x/\partial s$, and the inverse is an application of the inverse function theorem, so it would be literally 1 over the latter expression changing the coordinates accordingly. You would then have something like

for i in range(8):
    L += cp_f*w*dot(grad(T_n), normal(i))*ds(i)

I'm confused about why the hole is given two domains. What is "normal" in the above code? Since the flux is only applied in the hole, do you still need to loop over all 8 domains?

Whoops, you are right, you would only need to iterate over the last two. s fields are the hole parameterized, so that you get the normal with the gradient. It needs two domains because to get a global parametrization of the domain you would need spherical coordinates, which you do not have in your implementation (apparently).

I was able to get what I wanted following your suggestion of applying the chain rule and parameterizing the surface but I used a different parameterization than what you suggested. To recap, I was trying to model $\int dT$, where $T$ is the scalar field I am trying to solve. The problem is 3D and the geometry is a cube with a cylindrical hole in the center; picture a cube where one has drilled a hole going from the top face to the bottom face of the cube. The term arises from applying flux on the resulting cylindrical surface.

The first step was to apply the chain rule:
$$\int dT = \nabla T \cdot \frac{\partial x}{\partial S} dS$$

Note, that $x$ is a vector for the three direction. For my problem, I'm only concerned with the term $\partial T/ \partial z$, so my integral simplifies to
$$ \int dT = \frac{\partial T}{\partial z} \frac{\partial z}{\partial S} dS $$

The surface that I am applying the flux is cylindrical so $dS = 2\pi r dz$ and
$$\frac{\partial z}{\partial S} = \frac{1}{2 \pi r}$$

In fenics, my code looks like:

dzdS = Expression("1/(2pi)pow((x[0]x[0] + x[1]x[1]), -0.5)", degree=2, pi=pi)
L = cp_f * w * grad(T_n)[2] * dzdS * ds(1)

I was able to get the same answer as before when I used a 2D axisymmetric approach so hopefully this method would work when I model a scenario where symmetry is broken and the 3D modeling is required. Thanks a lot for your help and patience!

I was able to get what I wanted following your suggestion of applying the chain rule and parameterizing the surface but I used a different parameterization than what you suggested. To recap, I was trying to model $\int dT$, where $T$ is the scalar field I am trying to solve. The problem is 3D and the geometry is a cube with a cylindrical hole in the center; picture a cube where one has drilled a hole going from the top face to the bottom face of the cube. The term arises from applying flux on the resulting cylindrical surface.

The first step was to apply the chain rule:
$$\int dT = \int \nabla T \cdot \frac{\partial x}{\partial S} dS$$

Note, that $x$ is a vector for the three direction. For my problem, I'm only concerned with the term $\partial T/ \partial z$, so my integral simplifies to
$$ \int dT = \int \frac{\partial T}{\partial z} \frac{\partial z}{\partial S} dS $$

The surface that I am applying the flux is cylindrical so $dS = 2\pi r dz$ and
$$\frac{\partial z}{\partial S} = \frac{1}{2 \pi r}$$

In fenics, my code looks like:

dzdS = Expression("1/(2pi)pow((x[0]x[0] + x[1]x[1]), -0.5)", degree=2, pi=pi)
L = cp_f * w * grad(T_n)[2] * dzdS * ds(1)

I was able to get the same answer as before when I used a 2D axisymmetric approach so hopefully this method would work when I model a scenario where symmetry is broken and the 3D modeling is required. Thanks a lot for your help and patience!

...