Poisson's Equation

From KratosWiki

(Difference between revisions)

Jump to: navigation, search

Latest revision as of 09:02, 4 November 2009

Mathematical forms

The Poisson's Equation takes the form:

in one-dimensional Cartesian coordinates:

$\left[ \frac{d}{dx}\cdot \left( K_{x} \cdot \frac{d}{dx}\right) \right]\varphi(x) + Q(x) = 0$

in two-dimensional Cartesian coordinates:

$\left[ \frac{\partial}{\partial x}\cdot \left( K_{x} \cdot \frac{\partial}{\partial x}\right) + \frac{\partial}{\partial y}\cdot \left(K_{y} \cdot \frac{\partial}{\partial y} \right) \right]\varphi(x,y) = Q(x,y)$

in three-dimensional Cartesian coordinates:

$\left[ \frac{\partial}{\partial x}\cdot \left( K_{x} \cdot \frac{\partial}{\partial x}\right) + \frac{\partial}{\partial y}\cdot \left(K_{y} \cdot \frac{\partial}{\partial y} \right) + \frac{\partial}{\partial z}\cdot \left(K_{z} \cdot \frac{\partial}{\partial z} \right) \right]\varphi(x,y,z) = Q(x,y,z)$

for a domain

Ω

, equivalent to write:

$\nabla \cdot (K \nabla \varphi )= Q$

where K is the material property, $\varphi$ is the unknown to obtain and Q is a source term,

and an important secondary variable is $q = - K \nabla \varphi$ .

Boundary Conditions

Two different boundary conditions will be considered:

Dirichlet boundary condition

The first type boundary condition specifies the values the solution needs to take $(\varphi_{0})$ on a boundary $({\Gamma_{\varphi}})$ of the domain:

$\varphi\mid_{\Gamma_{\varphi}} = \varphi_{0}$

Neumann boundary condition

The second type boundary condition specifies the values that the derivative of a solution (q₀) is to take on the boundary (

Γ

_q) of the domain:

$K \nabla \varphi\mid_{\Gamma_{q}} = q_{0}$

Physical examples of the Poisson's equation

Heat Transfer[1]:
- $\varphi$ is the Temperature
- K is the Thermal Conductivity
- Q the Heat Source
- and q the Heat Flow
Electrostatics[2]:
- $\varphi$ is the Scalar Potential (Voltage)
- K is the Dielectric Constant
- Q the Charge Density
- q the Displacement Flux density
- and $\vec E = - \nabla \varphi$ is the Electrostatic Field
Electrostatics[3]:
- $\varphi$ is the Scalar Potential (Voltage)
- K is the Electric Conductivity ( $σ$ )
- and q is the Current Flow ( $\vec J = - \sigma \cdot \nabla \varphi$ )
Magnetostatics[4]:
- $\varphi$ is the Magnetic Potential
- K is the Permeability
- Q the Current Density
- q the Magnetic Field Intensity ( $\vec H$ )
- and $\vec B = - \nabla \varphi$ is the Magnetostatic Field or Magnetic Induction

Analytical Solution for One-Dimensional Domains

General solution using the Heat Transfer example

Consider the heat transfer without convection effects along the following bar:

Remember that the conduction phenomena refers to "the transfer of thermal energy from a region of higher temperature to a region of lower temperature through direct molecular communication within a medium or between mediums in direct physical contact without a flow of the material medium".

As boundary conditions, the temperature is fixed at the beginning of the bar, and the heat flow is given at the end of the bar:

$\varphi = \bar \varphi \mid_{x=x_0}$

$q = \bar q \mid_{x=x_L}$

Taking a differential length and by establishing the balance of heat flows, it can be written:

outflow $(q + dq) - (q) \frac{ }{ }$ inflow $= 0 \frac{ }{ }$

$dq = 0 \frac{ }{ }$

Using the Fourier law[5]:

$q = - k \cdot \frac{d \varphi}{dx}$

Therefore:

$dq = \frac{dq}{dx} dx = - \frac{d}{dx} \left(k \cdot \frac{d \varphi}{dx} \right) dx$

$\frac{d}{dx} \left(k \cdot \frac{d \varphi}{dx} \right) = 0$

If an external heat source is considered, the balance of heat flows becomes:

$(q + dq) - (q) - Q \cdot dx = 0 \frac{ }{ }$

$\frac{dq}{dx} - Q = 0$

$\frac{d}{dx} \left(k \cdot \frac{d \varphi}{dx} \right) + Q = 0$

Rewritting the boundary conditions:

$\varphi - \bar \varphi =0$ en $x=x_0 \frac{ }{ }$

$q - \bar q = - \left( k \cdot \frac{d \varphi}{dx} + \bar q \right) = 0$ en $x=x_L\frac{ }{ }$

Generic Solution:

Check a set of some specific examples of this analytical solution of the Poisson's equation for one-dimensional domains (including some figures and Matlab code you can modify).

Finite Element Method

The application of the Finite Element Method [6] (FEM) to solve the Poisson's equation consists in obtaining an equivalent integral formulation of the original partial differential equations (PDE). By dividing the whole domain in elements, the integral expression can be expressed as a sum of elementary integrals, easier to simplify as functions of lower order.

There are two basic procedures: variational formulation and residual formulation, that will be used here.

Weighted Residual Method (WRM)

The residual formulation is based on the Weighted Residual Method (WRM). The differential equation is converted in an integral equation with certain weighting functions applied to each equation.

Therefore, the Poisson's equation given by the governing PDE and its boundary conditions:

$A(\varphi) = \frac{d}{dx} \left( k \frac{d \varphi}{dx} \right) + Q = 0 ~~ in ~ \Omega$

$B(\varphi) = \begin{cases} \varphi - \overline{\varphi} = 0 & in ~ \Gamma_{\varphi} \\ k \frac{d \varphi}{dx} + \overline{q} = 0 & in ~ \Gamma_{q} \end{cases}$

can be written using the WRM as follows:

${ \int_{\Omega} W(x) A(\varphi) d\Omega + \oint_{\Gamma} \overline{W}(x) B(\varphi) d_{\Gamma}=0 }$

with $W(x) \frac{}{}$ and $\overline{W}(x)$ the weighting functions.

The above expressions use infinite functions, impossible to solve for general cases. Therefore, the solution to be found should be an approach able to be expressed in terms of finite functions. For example, an approach solution can be expressed as a lineal combination of some kind of well known functions $\frac{}{} N_i (x),$ called the basis functions:

$\varphi (x) \cong \hat \varphi (x) = \sum_{i=0}^n N_i (x) a_i$

Therefore, the integral equation to be solved become:

${ \int_{\Omega} W(x) A(\hat \varphi) d\Omega + \oint_{\Gamma} \overline{W}(x) B(\hat \varphi) d_{\Gamma}=0 }$

The above expression is the Weighted residual form of the equation because by using an approach solution $\frac{}{} \hat \varphi$ , the $\frac{}{} A(\hat \varphi)$ and $\frac{}{} B(\hat \varphi)$ expressions won't be zero:

$\frac{}{} r_{\Omega} = A(\hat \varphi) \ne 0 \quad in \quad \Omega$

$\frac{}{} r_{\Gamma} = B(\hat \varphi) \ne 0 \quad in \quad \Gamma$

and what it will be forced is:

${ \int_{\Omega} W(x) r_{\Omega} d\Omega + \oint_{\Gamma} \overline{W}(x) r_{\Gamma} d_{\Gamma}=0 }$

If the number of weighting functions ( $W(x) \frac{}{}$ and $\overline{W}(x)$ ) are the same to the basis functions ( $\frac{}{} N_i (x)$ ), the following system of equations:

${ \int_{\Omega} W_i A(\sum_j^n N_j a_j) d\Omega + \oint_{\Gamma} \overline{W_i} B(\sum_j^n N_j a_j) d_{\Gamma}=0 \qquad i=1..n }$

that can be written as:

$\frac{}{} \mathbf{K} \mathbf{a} = \mathbf{f}$

with K is a coefficients matrix that depends on the geometrical and physical properties of the problem, a is the vector with the n unknowns ( $\frac{}{} a_j$ ) to be obtained and f is a vector that depends on the source values and boundary conditions.

Using the same basis functions $\frac{}{} N_i (x)$ for the whole domain, and if they are selected to naturally accomplish the boundary conditions, the integral form can be written:

${ \int_{\Omega} W_i \left[ \frac{d}{dx} k \frac{d}{dx} \left( \sum_{j=0}^n N_j (x) a_j \right) + Q \right] d\Omega =0 \qquad i=1..n }$

and therefore:

$K_{ij}=\int_{\Omega} W_i(x) \left[ \frac{d}{dx} k \frac{d}{dx} \left( N_j (x) \right) \right] d\Omega$

$f_i = - \int_{\Omega} W_i(x) Q d\Omega \,$

Depending on the kind of weighting functions used, different weighted residual methods (WRM) can be considered:

Collocation Method:	$\frac{}{} W_i(x)=\delta(x - x_i)$ , $\frac{}{} i=1, 2,... , n$
Subdomain Method:	$\begin{cases} W_i(x)=1 & \forall x \in \Omega_i \\ W_i(x)=0 & \forall x \notin \Omega_i \end{cases}$
Galerkin Method:	$\frac{}{} W_i(x) \equiv N_i(x)$
Least Squares Method:	$W_i(x)=A(\varphi)$ and $\overline{W_i}(x)=B(\varphi)$
etc...

Check the resolution of an specific example of the Poisson's equation with the above diferents weights.

Weak form of the Weighted Residual Method

Coming back to the integral form of the Poisson's equation:

${ \int_{\Omega} W \left [ \frac{d}{dx} \left( k \frac{d \varphi}{dx} \right) + Q \right ]d\Omega + \oint_{\Gamma} \overline{W} \left [ k \frac{d \varphi}{dx} + \overline{q} \right ] d_{\Gamma}=0 }$

it should be noted that not always $\varphi(x)$ can be obtained, depending on the selected trial functions. Remember that some functions cannot be integrated, if they are not continuous.

Integrability conditions

The following pictures show that to integrate the $f^m \,$ derivative of a function, the previous derivatives ( $f(x), f'(x), f''(x)... f^{m-1}(x) \,$ ) should be continua.

That is $\int_{\Omega} \frac{d^m f(x)}{d x^m} dx$ can be computed if $\frac{d^{m-1} f(x)}{d x^{m-1}}$ is continuous.

A function f(x) has continuity of $C^m \,$ class, if f(x) and $\frac{d^{i} f(x)}{d x^{i}}$ with i=0 to m are also continuous.

C¹ means that f(x) and f'(x) are both continuous.

C⁰ means that just f(x) is continuous.

C^-1 means that f(x) is not continuous.

Applying this criterion to the integral form of the Poisson's equation shows:

$\varphi(x)$ needs to be $C^1 \,$ (and, therefore, the selected shape forms too)
k needs to be $C^0 \,$ , a clear problem for domains with different materials
W can be discontinuous

A new form of this integral can be obtained using the Integration by parts [7]. In 1D:

$\int_0^l W \left [ \frac{d}{dx} \left( k \frac{d \varphi}{dx} \right)\right ]dx = \left [W k \frac{d \varphi}{dx} \right ]_0^l - \int_0^l \frac{dW}{dx} k \frac{d \varphi}{dx} dx$

and the complete integral form of the Poisson's equation becomes:

$- \int_0^l \frac{dW}{dx} k \frac{d \varphi}{dx} dx + \int_0^l W Q dx + \left [W k \frac{d \varphi}{dx} \right ]_0^l + \left [\overline W \left ( k \frac{d \varphi}{dx} + \overline q \right) \right ]_l = 0$

that shows how the continuity requirements are now less strong than those with the original integral equation.

The following table summarises the new continuity requirements for each integral form:

Continuity requirements

	Original	Weak
$\varphi$	$C^1 \,$	$C^0 \,$
$W \,$ 1	$C^{-1} \,$	$C^0 \,$
$k \,$	$C^0 \,$	$C^{-1} \,$

Natural / Neumann boundary condition

If $\overline W = - W \,$ the integral form can be simplified as follows:

$\int_0^l \frac{dW}{dx} k \frac{d \varphi}{dx} dx = \int_0^l W Q dx - \left [W k \frac{d \varphi}{dx} \right ]_0 - \left [ W \overline q \right ]_l = \int_0^l W Q dx - \left [W q \right ]_0 - \left [ W \overline q \right ]_l$

Note that:

$\varphi \,$ is not present any more in the boundary equation where the heat flux is defined (x = l)
if the heat flux is zero in l ( $\overline q |_l= 0\,$ ), there is not any more reference to the boundary condition (Neumann condition) in the weak form, reason because it is called natural boundary condition
$q_0 \,$ is the heat flux in $x=0 \,$ due to the constraint that fixes the value of $\varphi \,$ in this point. ( $q_0 \,$ is the reaction to this constraint

Discretisation of the Weak Form

Following now the discretisation process seen in the precedent section, that is, by applying the basis functions:

$\varphi (x) \cong \hat \varphi (x) = \sum_{i=0}^n N_i (x) a_i$

${ \int_0^l \frac{d W_i}{dx} k \sum_{j=1}^n \frac{d N_j}{dx} a_j dx = \int_0^l W_i Q dx + \left[ W_i q \right ]_{x=0} - \left[ W_i \overline q \right ]_{x=l} \qquad i=1, 2, ... n }$

For each i:

$\int_0^l \frac{d W_1}{dx} k \left ( \frac{d N_1}{dx} a_1 + \frac{d N_2}{dx} a_2 + ... + \frac{d N_n}{dx} a_n \right ) dx = \int_0^l W_1 Q dx + \left[ W_1 q \right ]_{x=0} - \left[ W_1 \overline q \right ]_{x=l}$

$\int_0^l \frac{d W_2}{dx} k \left ( \frac{d N_1}{dx} a_1 + \frac{d N_2}{dx} a_2 + ... + \frac{d N_n}{dx} a_n \right ) dx = \int_0^l W_2 Q dx + \left[ W_2 q \right ]_{x=0} - \left[ W_2 \overline q \right ]_{x=l}$

...

$\int_0^l \frac{d W_n}{dx} k \left ( \frac{d N_1}{dx} a_1 + \frac{d N_2}{dx} a_2 + ... + \frac{d N_n}{dx} a_n \right ) dx = \int_0^l W_n Q dx + \left[ W_n q \right ]_{x=0} - \left[ W_n \overline q \right ]_{x=l}$

that can be written in matricial form:

$\begin{bmatrix} K_{11} & K_{12} & \cdots & K_{1n} \\ K_{21} & K_{22} & \cdots & K_{2n} \\ \vdots & ~ & \ddots & \vdots \\ K_{n1} & K_{n2} & \cdots & K_{nn} \end{bmatrix} \begin{Bmatrix} a1 \\ a2 \\ \vdots \\ an \end{Bmatrix} = \begin{Bmatrix} f1 \\ f2 \\ \vdots \\ fn \end{Bmatrix}$

$K a = f \,$

with:

$K_{ij}=\int_0^l \frac{d W_i}{dx} k \frac{d N_j}{dx} dx$

$f_{i}=\int_0^l W_i Q dx + (W_i q)_0 - (W_i \overline q)_l$

Depending on the selection of the weighting functions W, the corresponding Collocation, Subdomain, Galerkin method and so on is obtained. In the case that $W_i = N_i \,$ (Galerkin Method), the terms are as follows:

$K_{ij}=\int_0^l \frac{d N_i}{dx} k \frac{d N_j}{dx} dx$

$f_{i}=\int_0^l N_i Q dx + (N_i q)_0 - (N_i \overline q)_l$

Note that in this case K is symmetric. This is what is called the weak form of Galerkin.

Domain Discretisation: "local" Shape Functions

In the previous sections, a global approach has been used, that is, a set of functions to be applied over the whole domain. Nevertheless, the entire domain can be geometrically subdivided in elements and the definition of the set of functions be applied just for each of these portions of geometry.

These subdomains are what it is called the finite elements.

In a one-dimensional domain, for example and according to the above approaches, the unknown $\varphi(x) \,$ can be expressed as follows:

$\varphi(x) \cong \hat \varphi (x) = a_0 + a_1 x + a_2 x^2 + ... + a_n x^n =\sum_{i=0}^n a_i x^i$

by using n points where the $\varphi \,$ value is known, points which are called nodes. The constants $a_0, a_1, ..., a_n \,$ can be obtained just by using these $\varphi \,$ values in each node $x_i \,$ .

Therefore, the approach function $\hat \varphi (x) \,$ can also be expressed as follows:

$\hat \varphi (x) = N_1^{(e)}(x) \varphi_1^{(e)} + N_2^{(e)}(x) \varphi_2^{(e)} + ... + N_n^{(e)}(x) \varphi_n^{e} = \sum_{i=1}^n N_i^{(e)}(x) \varphi_i^{(e)}$

where:

$N_1^{(e)}(x) , N_2^{(e)}(x) , ... , N_n^{(e)}(x)$ are the polynomial interpolation functions, defined locally (in the element's domain);
$\varphi_i^{(e)} \,$ is the value of the unknown $\varphi \,$ in the i node;

This kind of well known functions $N_i^{(e)} \,$ are usually called Shape Functions and only affects to each element (e) with the contribution of the specific node i. For this reason is also called the Shape Function of the Node i. For more details about how to obtain some Shape Functions and their properties, please follow this link.

To be more specific, in one-dimension, the equation to be solved using the finite element method is:

$\int_0^l \frac{dW_i}{dx} k \frac{d \hat \varphi}{dx} dx = \int_0^l W_i Q dx - \left [W_i q \right ]_0 - \left [ W_i \overline q \right ]_l$

$\hat \varphi (x) = \sum_{i=1}^n N_i^{(e)}(x) \varphi_i^{(e)}$ and $\frac{d \hat \varphi (x)}{dx} = \sum_{i=1}^n \frac{d N_i^{(e)}(x)}{dx} \varphi_i^{(e)}$

by using the weak form of the WRM, which becomes:

$\int_0^l \frac{dW_i}{dx} k \sum_{j=1}^n \left( \frac{d N_j^{(e)}(x)}{dx} \varphi_j^{(e)} \right ) dx = \int_0^l W_i Q dx - \left [W_i q \right ]_0 - \left [ W_i \overline q \right ]_l$

If an element of two nodes is used, then:

$\hat \varphi (x) = \sum_{i=1}^2 N_i^{(e)}(x) \varphi_i^{(e)}$

If the equation is solved just for the element (e) and the Galerkin Method is applied, then the matricial form can be written as follows:

$\begin{bmatrix} K_{11}^{(e)} & K_{12}^{(e)} \\ K_{21}^{(e)} & K_{22}^{(e)} \end{bmatrix} \begin{Bmatrix} \varphi_1^{(e)} \\ \varphi_2^{(e)} \end{Bmatrix} = \begin{Bmatrix} f_1^{(e)} \\ f_2^{(e)} \end{Bmatrix}$

with:

$K_{ij}^{(e)}=\int_{l^{(e)}} \frac{d N_i^{(e)}(x)}{dx} k \frac{d N_j^{(e)}(x)}{dx} dx$

$f_{i}^{(e)}=\int_{l^{(e)}} N_i^{(e)}(x) Q dx$

Note that the boundary conditions have not been yet considered. In just a few lines will be treated.

In order to solve the whole system of equations, the assembly of the global matrix has to be performed, just by putting all together:

$\underbrace{ \begin{bmatrix} K_{11}^{(1)} & K_{12}^{(1)} & 0 & 0 & \cdots & 0 & 0 \\ K_{21}^{(1)} & K_{22}^{(1)} + K_{11}^{(2)} & K_{12}^{(2)} & 0 & \cdots & 0 & 0 \\ 0 & K_{21}^{(2)} & K_{22}^{(2)} + K_{11}^{(3)} & K_{12}^{(3)} & \cdots & 0 & 0 \\ \vdots & ~ & \ddots & ~ & ~ & ~ & \vdots \\ 0 & 0 & 0 & 0 & \cdots & K_{12}^{(n-1)} + K_{11}^{(n)} & K_{12}^{(n)} \\ 0 & 0 & 0 & 0 & \cdots & K_{21}^{(n)} & K_{22}^{(n)} \end{bmatrix} }_{K} \underbrace{ \begin{Bmatrix} \varphi_1 \\ \varphi_2 \\ \varphi_3 \\ \vdots \\ \varphi_n \\ \varphi_{n+1} \end{Bmatrix} }_{a} = \underbrace{ \begin{Bmatrix} f_1^{(1)} & + q_0\\ f_2^{(1)} & + f_1^{(2)} \\ f_2^{(2)} & + f_1^{(3)} \\ \vdots \\ f_2^{(n-1)} & + f_1^{(n)} \\ f_2^{(n)} & - \overline q \end{Bmatrix} }_{f}$

It should be noted that:

$a= \begin{Bmatrix} \varphi_1 & = \varphi_1^{(1)} \\ \varphi_2 & = \varphi_2^{(1)} & = \varphi_1^{(2)} \\ \varphi_3 & = \varphi_2^{(2)} & = \varphi_1^{(3)} \\ \vdots \\ \varphi_n & = \varphi_2^{(n-1)} & = \varphi_1^{(n)} \\ \varphi_{n+1} & = \varphi_2^{(n)} \end{Bmatrix}$

The boundary conditions have to be applied just to the boundaries, and for this reason they only have effect on the first node and last one.

The final system to be solved is:

$\begin{bmatrix} K_{22}^{(1)} + K_{11}^{(2)} & K_{12}^{(2)} & 0 & \cdots & 0 & 0 \\ K_{21}^{(2)} & K_{22}^{(2)} + K_{11}^{(3)} & K_{12}^{(3)} & \cdots & 0 & 0 \\ \vdots & ~ & \ddots & ~ & \vdots \\ 0 & 0 & 0 & \cdots & K_{12}^{(n-1)} + K_{11}^{(n)} & K_{12}^{(n)} \\ 0 & 0 & 0 & \cdots & K_{21}^{(n)} & K_{22}^{(n)} \end{bmatrix} \begin{Bmatrix} \varphi_2 \\ \varphi_3 \\ \vdots \\ \varphi_n \\ \varphi_{n+1} \end{Bmatrix} = \begin{Bmatrix} f_2^{(1)} & + f_1^{(2)} & - K_{21}^{(1)} \varphi_1\\ f_2^{(2)} & + f_1^{(3)} \\ \vdots \\ f_2^{(n-1)} & + f_1^{(n)} \\ f_2^{(n)} & - \overline q \end{Bmatrix}$

A particular solution using one and two elements is explained in this section. A Matlab and Python codes are also available.

Revision as of 13:50, 10 December 2008 (view source) JMora (Talk \| contribs) (→Domain Discretisation: "local" Shape Functions) ← Older edit		Latest revision as of 09:02, 4 November 2009 (view source) JMora (Talk \| contribs) (→Domain Discretisation: "local" Shape Functions)
Line 660:		Line 660:
	:A particular solution using one and two elements is explained in [[Resolution of the 1D Poisson's equation using local Shape Functions \|this section]]. A [[Matlab code of the Resolution of the 1D Poisson's equation using local Shape Functions\| '''Matlab''']] and [[Python code of the Resolution of the 1D Poisson's equation using local Shape Functions\| '''Python''']] codes are also available.		:A particular solution using one and two elements is explained in [[Resolution of the 1D Poisson's equation using local Shape Functions \|this section]]. A [[Matlab code of the Resolution of the 1D Poisson's equation using local Shape Functions\| '''Matlab''']] and [[Python code of the Resolution of the 1D Poisson's equation using local Shape Functions\| '''Python''']] codes are also available.

−	[[Category:~~Theory~~]]	+	[[Category:Poisson's Equation]]

Poisson's Equation

Latest revision as of 09:02, 4 November 2009

Contents

Mathematical forms

Boundary Conditions

Dirichlet boundary condition

Neumann boundary condition

Physical examples of the Poisson's equation

Analytical Solution for One-Dimensional Domains

Finite Element Method

Weighted Residual Method (WRM)

Weak form of the Weighted Residual Method

Integrability conditions

Natural / Neumann boundary condition

Discretisation of the Weak Form

Domain Discretisation: "local" Shape Functions

Views

Personal tools

Navigation

Search

Toolbox

Categories