Tutorial/Wave equation with Dirichlet boundary control
Contents
The considered problem
In this problem, we consider a Wave equation with a Dirichlet control on a part of the boundary, and homogeneous Dirichlet or Neumann condition on the other part. To model this in GetDP, we will introduce a "Constraint" with "TimeFunction". We assume that the reader has already studied the parabolic example on the Heat equation.
We considered the annulus $\Omega = \mathcal D (0, 1) \setminus \mathcal D (0, 0.2)$ and we seek the functions $w_N$ and $w_D$, respective solutions of the following problems
\begin{equation} \begin{cases}\label{eq:problemNeu} \dfrac{\partial^2}{\partial t^2} w_N(x,t) -\Delta w_N(x,t) = 0 & \forall x \in \Omega, t\geq0\\ \dfrac{\partial w_N}{\partial \mathbf{n}}(x,t) = 0 & \forall x \in \Gamma_0, t \geq0\\ w_N(x,t) = u(x,t) & \forall x \in \Gamma_1, t \geq0. \end{cases} \end{equation}
\begin{equation} \begin{cases}\label{eq:problemDir} \dfrac{\partial^2}{\partial t^2} w_D(x,t) -\Delta w_D(x,t) = 0 & \forall x \in \Omega, t\geq0\\ w_D(x,t) = 0 & \forall x \in \Gamma_0, t \geq0\\ w_D(x,t) = u(x,t) & \forall x \in \Gamma_1, t \geq0. \end{cases} \end{equation}
where $\Gamma_0 = \partial \mathcal D (0, 1)$ is the boundary of $\mathcal D (0, 1)$, $\Gamma_1 = \partial \mathcal D (0, 0.2)$ is the boundary of $\mathcal D (0, 0.2)$ and $\displaystyle{\Delta = \frac{\partial^2}{\partial x_1^2} + \frac{\partial^2}{\partial x_2^2} }$ is the Laplace operator.
To simplify, we suppose that the control $u$ is defined by $$ \forall (x,y)\in \Gamma_1, t\geq0 \qquad u(x,y,t) = u(t) = 5\sin[2\pi t]. $$
Physically, these problems model the vibration of a membrane on the annulus over time. The Dirichlet boundary condition models the fact that the membrane is fixed on the boundary $\Gamma_0$ at $0$, as the Neumann boundary condition models the fact that the membrane is let free on the boundary $\Gamma_0$ to move up and down. Furthermore, the control $u$ corresponds to an up and down motion of the boundary $\Gamma_1$ of the membrane oscillating between -5 and 5.
To solve problems \eqref{eq:problemNeu} and \eqref{eq:problemDir} using the finite element method, we write their weak formulations:
TO FINISH HERE (+TIME)
\begin{equation}\label{eq:WeakFormulationNeu} \left\{\begin{array}{l} \text{Find } v_N\in \mathcal H = \{w \in H^1(\Omega) \mid w = u \text{ on } \Gamma_1 \} \text{ such that, }\\ \displaystyle{\forall \varphi\in \mathcal H, \qquad \int_{\Omega}\dfrac{d}{dt}v_N \varphi \;{\rm d}\Omega + \int_{\Omega} \nabla v_N\cdot\nabla \varphi \;{\rm d}\Omega = 0}, \end{array}\right. \end{equation}
\begin{equation}\label{eq:WeakFormulationDir} \left\{\begin{array}{l} \text{Find } v_D\in \mathcal H = \{w \in H^1(\Omega) \mid w = 0 \text{ on } \Gamma_0, \; w = u \text{ on } \Gamma_1 \} \text{ such that, }\\ \displaystyle{\forall \varphi\in \mathcal H, \qquad \int_{\Omega}\dfrac{d}{dt}v_D \varphi \;{\rm d}\Omega + \int_{\Omega} \nabla v_D\cdot\nabla \varphi \;{\rm d}\Omega = 0}, \end{array}\right. \end{equation} where the functions $\varphi$ are the test functions.
There is an additional difficulty in this case than in the Heat equation. Indeed, here we are dealing with a second-order PDE, which implies that the time scheme needs two steps to be correctly initialised. In other words, we have to initialised both position and velocity of the membrane.
We have to take that into account in the Resolution group in GetDP. If the initial velocity is equals to zero, then we just have to double the initialisation InitSolution[SysWave];, this will put the initial state in the two first steps of the time discretisation. However, if we want to solve a system with non-zero initial velocity, we have to use a trick. The idea is to initialise the velocity using a weak formulation, which has to be zero when the Time Step is greater than 2. We invite the reader to read carefully the Formulation group of the Wave.pro file to clarify what we mean here. Then, we solve separetely the second time step and the following one, as one can read in the Resolution group of the Wave.pro file.
Outline of the program
The proposed solution is composed by 2 different files:
- Annulus.geo
- GMSH file, used to build both the square $\Omega$ and its boundary $\Gamma_0$ and $\Gamma_1$.
- Wave.pro
- GetDP file, contains the weak formulations \eqref{eq:WeakFormulationNeu} and \eqref{eq:WeakFormulationDir} of the problems \eqref{eq:problemNeu} and \eqref{eq:problemDir}.
Annulus.geo: creation of the geometry with GMSH
/* ----------------------- File Annulus.geo ----------------------- */ // Mesh step h : h = 0.025; /*==============*/ /* Large Circle */ /*==============*/ // Points : Point(1) = {0,0,0,h}; Point(2) = {1,0,0,h}; Point(3) = {0,1,0,h}; Point(5) = {-1,0,0,h}; Point(6) = {0,-1,0,h}; // Arcs : Circle(1) = {2,1,3}; Circle(2) = {3,1,5}; Circle(3) = {5,1,6}; Circle(4) = {6,1,2}; // Circle : Line Loop(1) = {1,2,3,4}; /*===============*/ /* Little Circle */ /*===============*/ // Points : Point(200) = {0.2,0,0,h}; Point(300) = {0,0.2,0,h}; Point(500) = {-0.2,0,0,h}; Point(600) = {0,-0.2,0,h}; // Arcs : Circle(100) = {200,1,300}; Circle(200) = {300,1,500}; Circle(300) = {500,1,600}; Circle(400) = {600,1,200}; // Circle : Line Loop(100) = {100,200,300,400}; /*=======*/ /* Omega */ /*=======*/ // Annulus : The surface between the little circle and the large circle Plane Surface(1) = {1, 100}; /*===================*/ /* Physical Entities */ /*===================*/ // Gama0 Physical Line(1001) = {1,2,3,4}; // Gama1 Physical Line(1002) = {100,200,300,400}; // Omega Physical Surface(1005) = {1};
Direct link to file `Heat/GMSH_GETDP/Annulus.geo'
Heat.pro: weak formulations
Direct link to file `Heat/GMSH_GETDP/Wave.pro'
How to launch
The procedure is the same as described in this example, but we recall it here. Let us recall that all files (.geo and .pro) must be located in the same directory. We just describe the procedure using a terminal.
- Meshing the domain
In the right directory, type:
- gmsh Annulus.geo -algo iso -2
(again, the algorithm could be another one) After that, a file Annulus.msh should have been created.
- Solving the problem
Type the following command in the terminal:
- getdp Wave.pro -solve '#1' -pos '#1' -mesh Annulus.msh
- Displaying the result
Simply open the two .pos files with GMSH, and compare them to see the influence of the control on $\Gamma_1$ and the Neumann and Boundary conditions on $\Gamma_0$.
Result
Here is the result obtained with GMSH/GetDP. The first and second pictures show respectively the solution $w_N$ and $w_D$.
All files
A zipfile containing all files (.geo, .pro, .mesh and .pos) can be downloaded here.