Difference between revisions of "Electrostatics template"

From ONELAB
Jump to: navigation, search
(→‎Weak formulation)
Line 38: Line 38:
 
[[File:omega.png|thumb|right|alt=Domain and its boundary.|Domain and its boundary.]]
 
[[File:omega.png|thumb|right|alt=Domain and its boundary.|Domain and its boundary.]]
  
−
Let us skip all the mathematical details related to the conditions for the existence and uniqueness of the solution to this problem and focus instead on reformulating of the problem in order to solve it using the finite element method with [http://geuz.org/getdp GetDP]. In order to do this, we need to transform the original partial differential equation into a so-called ''[http://en.wikipedia.org/wiki/Weak_formulation weak formulation]''. This weak formulation is obtained by ''testing'' equation  \eqref{eq:elesta} with a series of functions, which means taking the scalar product of the equation with these ''test functions''. Mathematically, one can show that \eqref{eq:elesta} is equivalent to finding the solution $v\in H^1(\Omega)$ (we will come back to what this means shortly) such that  
+
Let us skip all the mathematical details related to the conditions for the existence and uniqueness of the solution to this problem and focus instead on reformulating of the problem in order to solve it using the finite element method with [http://geuz.org/getdp GetDP]. In order to do this, we need to transform the original partial differential equation into a so-called ''[http://en.wikipedia.org/wiki/Weak_formulation weak formulation]''. This weak formulation is obtained by ''testing'' equation  \eqref{eq:elesta} with a series of functions, which means taking the scalar product of the equation with these ''test functions''. Mathematically, one can show that \eqref{eq:elesta} is equivalent to finding the solution $v$ such that  
 
\begin{equation}\label{eq:elesta_weak}
 
\begin{equation}\label{eq:elesta_weak}
−
- \int_\Omega (\Div{\epsilon\,\Grad{v}}) \,v' \,d\Omega = \int_\Omega \rho  \,v' \,d\Omega ,\quad \forall v' \in H^1_0(\Omega).
+
- \int_\Omega (\Div{\epsilon\,\Grad{v}}) \,v' \,d\Omega = \int_\Omega \rho  \,v' \,d\Omega
 
\end{equation}
 
\end{equation}
−
This is the weak formulation of \eqref{eq:elesta}, where $H^1(\Omega)$ and $H^1_0(\Omega)$ are [http://en.wikipedia.org/wiki/Sobolev_space Sobolev spaces] of square integrable functions whose derivatives are also square integrable. In addition, functions in $H^1_0(\Omega)$ have zero value on the boundary of the domain $\Omega$. The scalar product used here is the standard $L^2$ scalar product
+
holds for all test functions $v'$ in an appropriate function space.
 +
This is the weak formulation of \eqref{eq:elesta}, where the ''appropriate function spaces'' for $v$ and $v'$ are [http://en.wikipedia.org/wiki/Sobolev_space Sobolev spaces] of square integrable functions whose gradients are also square integrable, the precise definition of which depends on the applied boundary conditions. The scalar product is the standard $L^2$ scalar product
 
\begin{equation}
 
\begin{equation}
 
(u,v)_\Omega = \int_\Omega u v \, d\Omega
 
(u,v)_\Omega = \int_\Omega u v \, d\Omega
 
\end{equation}
 
\end{equation}
−
and we will usually thus write \eqref{eq:elesta_weak} in a more compact form as:
+
and we can write the weak formulation in a more compact form as: Find $v$ such that
 
\begin{equation}\label{eq:elesta_weak2}
 
\begin{equation}\label{eq:elesta_weak2}
−
- (\Div{\epsilon\,\Grad{v}} , v')_\Omega = (\rho , v')_\Omega ,\quad \forall v' \in H^1_0(\Omega).
+
- (\Div{\epsilon\,\Grad{v}} , v')_\Omega = (\rho , v')_\Omega, \forall v'.
 
\end{equation}
 
\end{equation}
−
 
+
It is advantageous to integrate the left-hand side of \eqref{eq:elesta_weak2} by parts, which symmetrizes it and reduces the continuity requirements on $\epsilon$ and $v$. Applying [http://en.wikipedia.org/wiki/Green%27s_identities Green's first identity], the weak formulation then becomes: Find $v$ such that:
−
Setting aside the technical details about Sobolev spaces, what \eqref{eq:elesta_weak2} means is: ''solving the original equation is equivalent to solving an infinity of scalar equations obtained by taking the scalar product of the original equation with appropriately chosen test functions''. We do have an infinity of scalar equations since we must test with each $v' \in H^1_0(\Omega)$, and $H^1_0(\Omega)$ is a vector space of infinite dimension.
 
−
It is advantageous to integrate the left-hand side of \eqref{eq:elesta_weak2} by parts, which symmetrizes it and reduces the order of the differential operator applied to $v$. Applying [http://en.wikipedia.org/wiki/Green%27s_identities Green's first identity], our weak formulation then becomes: Find $v\in H^1(\Omega)$ such that:
 
 
\begin{equation}\label{eq:elesta_weak3}
 
\begin{equation}\label{eq:elesta_weak3}
−
(\epsilon\,\Grad{v} , \Grad{v'} )_\Omega = (\rho , v')_\Omega ,\quad \forall v' \in H^1_0(\Omega).
+
(\epsilon\,\Grad{v} , \Grad{v'} )_\Omega + (\epsilon \vec{n}\cdot\Grad{v} , v' )_{\partial\Omega} = (\rho , v')_\Omega ,\quad \forall v'.
 
\end{equation}
 
\end{equation}
−
Note that the boundary term $(\vec{n}\cdot\Grad{v} , v' )_{\partial\Omega}$ that arises through integration by parts vanishes here since we chose $v'\in H^1_0(\Omega)$, which implies that $v'\equiv 0$ on the boundary $\partial\Omega$ of the domain $\Omega$.
 
  
 
== Implementation in GetDP ==
 
== Implementation in GetDP ==

Revision as of 18:03, 5 March 2016

\(\renewcommand{\vec}[1]{\mathbf{#1}} \newcommand{\Grad}[1]{\mathbf{\text{grad}}\,{#1}} \newcommand{\Curl}[1]{\mathbf{\text{curl}}\,{#1}} \newcommand{\Div}[1]{\text{div}\,{#1}} \newcommand{\Real}[1]{\text{Re}({#1})} \newcommand{\Imag}[1]{\text{Im}({#1})} \newcommand{\pvec}[2]{{#1}\times{#2}} \newcommand{\psca}[2]{{#1}\cdot{#2}} \newcommand{\E}[1]{\,10^{#1}} \newcommand{\Ethree}{{\mathbb{E}^3}} \newcommand{\Etwo}{{\mathbb{E}^2}} \newcommand{\Units}[1]{[\mathrm{#1}]} \)The Electrostatics.pro template focuses on one of the simplest electromagnetic problem: the determination of the distribution of the electric field $\vec{e}$ $\Units{V/m}$ due to static (non-varying, in both space and time) electric charge densities $\rho$ $\Units{C/m^3}$ and/or levels of electric potential $v$ $\text[V]$. This is called an electrostatic problem, or simply electrostatics.

  • Download the templates (they are also directly available in the ONELAB bundle)
  • Create a geometry with Gmsh, or open an existing geometry
  • File/Merge the Electrostatics.pro template
  • For each physical group in the model, you will be asked to select the material, sources or boundary conditions

Every time you run your model, the ONELAB database is saved to disk (you can later reload it using the gear menu). An flattened problem file containing the GetDP definitions is also saved in case you want to edit it manually later on.

Governing equations

The equations describing electrostatics are obtained by dropping the terms with time derivatives (since we only consider static phenomena) in Maxwell's equations. The equations involving the electric and magnetic quantities (the electric field $\vec{e}$ and the electric displacement $\vec{d}$ on one side, the magnetic field $\vec{h}$ and the magnetic induction $\vec{b}$ on the other) then become decoupled and in order to determine the electric field we can thus only consider: \begin{align} \Curl{\vec{e}} &= 0 , \label{eq:curle}\\ \Div{\vec{d}} &= \rho , \label{eq:divd} \end{align} together with the dielectric constitutive law \begin{equation} \vec{d} = \epsilon \vec{e} + \vec{p}_e , \label{eq:de} \end{equation} where the charge density $\rho$ and the source polarization $\vec{p}_e$ are given (see Maxwell's equations for additional details).

Scalar electric potential

Since $\Curl{\vec{e}} = 0$, the electric field can be derived from a scalar potential (the electric scalar potential $v$, whose MKSA unit is the Volt $\Units{V}$) such that \begin{equation} \vec{e} = - \Grad{v} . \label{eq:v} \end{equation}

Introducing \eqref{eq:v} and \eqref{eq:de} in \eqref{eq:divd}, we obtain \begin{equation}\label{eq:elesta} -\Div{\epsilon\,\Grad{v}} = \rho . \end{equation} Our goal is to solve this equation in a bounded open set $\Omega$ of the oriented Euclidean space $\Ethree$. Note that if $\epsilon$ is a constant, \eqref{eq:elesta} reduces to Poisson's equation: $\Delta\,v=-\rho/\epsilon$.

Weak formulation

Domain and its boundary.
Domain and its boundary.

Let us skip all the mathematical details related to the conditions for the existence and uniqueness of the solution to this problem and focus instead on reformulating of the problem in order to solve it using the finite element method with GetDP. In order to do this, we need to transform the original partial differential equation into a so-called weak formulation. This weak formulation is obtained by testing equation \eqref{eq:elesta} with a series of functions, which means taking the scalar product of the equation with these test functions. Mathematically, one can show that \eqref{eq:elesta} is equivalent to finding the solution $v$ such that \begin{equation}\label{eq:elesta_weak} - \int_\Omega (\Div{\epsilon\,\Grad{v}}) \,v' \,d\Omega = \int_\Omega \rho \,v' \,d\Omega \end{equation} holds for all test functions $v'$ in an appropriate function space. This is the weak formulation of \eqref{eq:elesta}, where the appropriate function spaces for $v$ and $v'$ are Sobolev spaces of square integrable functions whose gradients are also square integrable, the precise definition of which depends on the applied boundary conditions. The scalar product is the standard $L^2$ scalar product \begin{equation} (u,v)_\Omega = \int_\Omega u v \, d\Omega \end{equation} and we can write the weak formulation in a more compact form as: Find $v$ such that \begin{equation}\label{eq:elesta_weak2} - (\Div{\epsilon\,\Grad{v}} , v')_\Omega = (\rho , v')_\Omega, \forall v'. \end{equation} It is advantageous to integrate the left-hand side of \eqref{eq:elesta_weak2} by parts, which symmetrizes it and reduces the continuity requirements on $\epsilon$ and $v$. Applying Green's first identity, the weak formulation then becomes: Find $v$ such that: \begin{equation}\label{eq:elesta_weak3} (\epsilon\,\Grad{v} , \Grad{v'} )_\Omega + (\epsilon \vec{n}\cdot\Grad{v} , v' )_{\partial\Omega} = (\rho , v')_\Omega ,\quad \forall v'. \end{equation}

Implementation in GetDP

GetDP uses weak formulations such as \eqref{eq:elesta_weak3} to define the problems it can solve. Solving an electrostatic problem in GetDP thus amounts to transcribe \eqref{eq:elesta_weak3} into GetDP's own language, which is what is done in Electrostatics.pro

Retrieved from "http://onelab.info/onelab_wiki/index.php?title=Electrostatics_template&oldid=8054"