Variable Coefficients

3. Variable Coefficients#

\(\DeclareMathOperator{\opdiv}{div}\) A setup with varying heat-conductivity is modeled by the equation

\[ -\opdiv \lambda(x) \nabla u(x) = f(x), \]

where \(\lambda : \Omega \rightarrow {\mathbb R}^+\) is the heat-conductivity, relating the heat-flux

\[ q := -\lambda \nabla u \]

to the temperature gradient \(\nabla u\).

In the case of discontinuous \(\lambda\), the equation is understood in distributional sense. This includes the interface conditions: The temperature on the left and right side are equal, and the heat flux leaving the left side enters the right side:

\[\begin{eqnarray*} u_l & = & u_r \\ \lambda_l \frac{\partial u_r}{\partial n} & = & \lambda_r \frac{\partial u_r}{\partial n} \end{eqnarray*}\]

The variational form is: find \(u \in H^1(\Omega)\)

\[ \int_\Omega \lambda(x) \nabla u \nabla v \, dx = \int_\Omega f v \, dx \]

There is no issue with discontinuous coefficients. Both interface conditions are included:

  • continuity of temperature \(u\) by the continuity of the trial space

  • continuity of the heat flux in weak sense, similar as Neumann boundary conditions

from ngsolve import *
from ngsolve.webgui import Draw

Create a 2D geometric model:

from netgen.occ import *
box = Rectangle(1,1).Face()
circ1 = Circle((0.3,0.7), 0.1).Face()
circ2 = Circle((0.7,0.7), 0.1).Face()
bar = MoveTo(0.2,0.2).Rectangle(0.6,0.1).Face()
air = box-circ1-circ2-bar

circ1.faces.name = "left"
circ2.faces.name = "right"
air.faces.name = "air"
bar.faces.name = "bar"
air.edges.Min(Y).name ='b'
air.edges.Max(X).name ='r'

shape = Glue([air,circ1, circ2,bar])
Draw (shape);

geo = OCCGeometry(shape, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.05)).Curve(3)
Draw (mesh);

print (mesh.GetMaterials())

('air', 'left', 'right', 'bar')

print (mesh.GetBoundaries())

('b', 'default', 'r', 'default', 'default', 'default', 'default', 'default', 'default', 'default')

fes = H1(mesh, order=3, dirichlet="b|r")
u,v = fes.TnT()

make a coefficient function taking one constant, or even another coefficient function per material:

lam = mesh.MaterialCF({ "bar" : 5, "left" : 2.5, "right" : 2.5 }, default=1 )
Draw (lam, mesh, "lambda");

use conductivity in bilinear-form:

a = BilinearForm(lam*grad(u)*grad(v)*dx).Assemble()
f = LinearForm(50*v*dx("left")).Assemble()

reg = mesh.Materials("left|bar")
print (type(reg))
print ("volume or boundary=", reg.VB())
print ("bitmask =", reg.Mask())

<class 'ngsolve.comp.Region'>
volume or boundary= VorB.VOL
bitmask = 0: 0101

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec

Draw (gfu, mesh);

The gradient is small in regions with high conductivity:

Draw (grad(gfu), mesh, "gradient");

The heat-flux is large in highly conductive materials:

Draw (-lam*grad(gfu), mesh, vectors={"grid_size":30});

Exercise:

  • experiment with different (positive!) values for the conductivity coefficients