Solving PDE using FEM

In this lab you will familiarize yourself with how to solve PDE's using modern Finite Element methods. We will revisit some of the problems solved in lectures and lab 2, but now instead of finite difference methods, solving them using FEM methods.

Square

(Tutorial, no need to hand in)

First consider Poisson's equation on a square:

$$-\nabla^2 u(x,y) = f(x,y)$$

$$x \in [-1, 1], y \in [-1, 1]$$

With u=0 boundary conditions (homogeneous Dirichlet). The solution $u$ corresponds to the electric potential on an a membrane charged with $f(x,y)$ electrons, and the potential 0 Volt on the boundary. Such charged membranes are used in among things electrostatic loudspeakers (Quad, a British brand used in recording studios) and headphones (Stax). (To drive the membrane and generate tones, two additional plates are placed, one on each side of the membrane and an alternating current, e.g. music is fed with opposite polarity to the plates.)

To enter this problem, start ``pdetool'', Select the square tool and draw a square. To more easily select the desired size, use the options meny to set grid spacing, and then ``snap to grid''. To define the boundary conditions switch to boundary mode by clicking the $\Re \partial \Omega$ button. Double click on each boundary in your figure to enter values. Homogeneous Dirichlet is the default.

After having defined the boundary, you use the ``PDE'' button to enter the desired PDE. Poisson's is the default. Use the ``Solve'' menu to solve the PDE. Different ways to plot the resulting solution can be accessed by pressing the graph like button between ``='' and zoom. A contour plot is useful for quantitative evaluation, that e.g. symmetry is as expected. A 3D plot with color shading is useful for a quick qualitative view.

Also try to change the boundary conditions on one of the sides to say 8. Set $f(x,y)=0$ to see only the effect of the boundary condition.

Hollow square

(Tutorial, no need to hand in)

Next do a hollow square with outer boundary $[-1, 1]$ and inner boundary $[-.5, .5]$, see figure. Recall that a finite difference version of this problem was solved in the lectures. You may compare that solution with the finite element one done here.

Draw two squares, and then add formula SQ1-SQ2 on ``set formula box''. Switch to mesh mode and you should have a mesh similar in resolution to in this figure. Set boundary conditions back to homogeneous. Solve and study the contour plot of the solution. Is the solution plausible? What are possible sources of errors? Try ``refine mesh'', study the resulting mesh. What happened to each triangle? Now solve the problem again on the refined mesh. Is the result npow plausible? How many times can you refine, and still have a reasonable (seconds) computation time?

In many practical situations it is important to reduce the size of the problem. We can reduce the domain by utilizing symmetry. Across symmetry boundaries the directional derivative is 0. Hence we can reduce the domain, and use Neumann boundary conditions on the edges that represent symmetry boundaries. Consider how to reduce the domain using symmetries to 1/2, 1/4 and 1/8. For the 1/8 case input this geometry, boundary conditions, and solve the problem. Verify that solution is identical to the full domain you just did. How much time is saved?

Heat Equation

Next we consider the heat problem solved using finite differences in lab 2. You can set the workspace to a suitable domain using ``axes limits'' under options.

Use the cylindrical coordinate form $(r,\alpha,z)$ of the the Laplacian operator, $\nabla^2 T = {1 \over r} {\partial \over \partial r} (r {\partial T \over \pa... ...er r^2} {\partial^2 T \over \partial \alpha} + {\partial^2 T \over \partial z}$, to solve the heat equation:

$${\partial T \over \partial t} = \kappa \nabla^2 T + {\kappa h \over \lambda}$$

Where $\lambda$ is the thermal conductivity, and $\kappa={\lambda \over c_v \rho}$, is the thermal conduction capacity $c_v$ = specific heat capacity $\rho$ is material density, and finally $h$ is the heat addition.

For the FEM solution in our heat pump problem it is easier to model the cold addition as a boundary condition for Bore hole temp $T_b$. Most heat pumps use a glycol or brine solution in the cold circuit, and run a few degrees below freezing.

Thermal conductivity of soil in a bore hole $\lambda = 2$ (W/(m C)

Bore hole temp $T_b = -8C$

Earth temperature $T_e = 4$ C

Surface temperature $T_s = 0$ C

Bore hole radius, $r_b = 0.1$ m

We consider a heat pump installation (see figure above) with one vertical ground loop bored 45m into the ground. For numerical solution, we need to discretize it on a finite domain. A reasonable choice is a domain significantly larger than 45m. We choose a domain extending 90m horizontally, and 100m vertically.

Stationary problem (cmput 4/5/6xx)

After a sufficiently long time ${\partial T \over \partial t} = 0$ and we can write the stationary heat equation as:

$$\nabla^2 T = -{h \over \lambda}$$

To enter the cylindrical coordinate version of the Laplacian, you use the variable ``x'' in the PDE specification.

Time dependent problem (cmput 4/5/6xx)

Here we assume the heat pump is on for half the year, off the other half. Let the initial temperature in our region (before the heat pump is ever turned on) $T_0 = T_e$, and solve for time t = [1,2,3, ..., 48] months, using h=10m and h=1m given:

Specific heat capacity of (dry) soil $c_v = 840$ (J/(kg C)

Density of soil $\rho = 1555$ (kg/$\rm {m}^3$)

Bore hole temp $T_b = -8C$ for Oct-Mar, $T_b = 8C$ otherwise. (Ideal would be to switch to a Neumann boundary condition Apr-Sept, but I don't think that is possible using the pdetool.)

Surface temperature $T_s = 0$ C for Oct-Mar, 15C otherwise

Note: Above quantities give ${\partial T \over \partial t}$ in Celcius/second. In the pde tool you have to give an predetermined time discretization. (accessed under ``solve'' then parameters. Note that you have to switch the PDE type to parabolic first to get the right parameters. A reasonable time step would be days or one week. To enter the time varying boundary condition it is convenient to change the time ($t$) unit to years and use round(t-floor(t)) as a way to change the boundary condition from warm to cold.

You can rescale above given quantities to get ${\partial T \over \partial t}$ in Celcius/year.