Recently, I'm working on the numerical simulation of Nano-second pulsed plasma. The equation to solve is the drift diffusion equation, which has a similar form to 🍔's equation. If it's a steady problem, then the equation becomes elliptical one. However, there are several challenges for me including: 1. Finite element method with weak form 2. New solver: COMSOL 3. The settings of boundary condition
So, to solve the problem, I decide to start with the simple one. I first try to set up a simple PDE solver in COMSOL for Poisson's Equation, as the following: $$ -\frac{d }{dx} \left[ \eta \frac{d T}{d x} + h \cdot T\right] = - (A \cdot x + B)$$
The equation is made up based on the analytical solution: $$ T = \frac{A}{2h} x^2 + \frac{1}{h} \left(B-\frac{\eta}{h}A\right)x + T(x=0)$$
The calculation domain is $x\in [0,1]$. Then, for Neumann boundary condition, we have:
$$ \begin{aligned} q_L & = \left. -\eta\frac{d T}{d x} - hT \right|_{x=0} \\ q_R & = \left. -\eta\frac{d T}{d x} - hT \right|_{x=1m} \\ \end{aligned} $$
Analytical Solution
Poisson's equation doesn't have a unique solution of the dependent variable, although its derivate is unique.
We can first solve for $q$:
$$ \begin{aligned} -q(x) & = \eta \frac{d T}{d x} + h \cdot T \\ & = \int_0^{x} (A\cdot x + B)\mathrm{d} x \\ & = \frac{A}{2}\cdot x^2 + B\cdot x - q_L \end{aligned} $$
🔔 Interestingly, the right boundary doesn't play a role here. Physically, it means we can only provide specific heat flux on one side of the rod. The heat flux on the other side is controlled by the heating rate of the whole rod for steady state. Otherwise, the temperature of the whole rod will keep changing. However, we can connect one side to a heat source and the other side to a heat bath, i.e. constant temperature.
Then we solve the following ODE:
$$ \eta \frac{dT}{dx} + h T = \frac{A}{2}\cdot x^2 + B\cdot x - q_L $$
Assuming there exists a function $y(x)$ satisfying $y'/y = h/\eta $, then we can get:
$$ \frac{d(T\cdot y)}{dx} = \frac{y}{\eta}\left(\frac{A}{2}\cdot x^2 + B\cdot x - q_L\right) $$
Obviously, $y(x) = \exp(h/\eta\cdot x)$ satisfies the condition. Denote:
$$ \begin{aligned} \alpha &= \frac{A}{2h} \\ \beta &= \frac{1}{h} \left(B-\frac{\eta}{h}A\right)x \\ \end{aligned} $$
We can rearrange the ODE:
$$
\begin{aligned}
\frac{d(T\cdot y)}{dx} & = \frac{y}{\eta}\left(\alpha h x^2 + (2\alpha \eta + \beta h)x - q_L\right)
& = \frac{h}{\eta}y \cdot (\alpha x^2 + \beta x - \frac{q_L}{h}) + y\cdot(2\alpha x + \beta)
\end{aligned}
$$
We can then get the solution:
$$ \begin{aligned} T(x) & = \alpha x^2 + \beta x - \frac{q_L}{h} + C\cdot \exp(-\frac{h}{\eta}\cdot x) \end{aligned} $$
Const $C$ is determined by the value of $T$ at the boundary. > The temperature doesn't have a unique solution with Neumann conditions on both boundary. Moreover, the problem can become overdefined with inappropriate Neumann conditions. To get a unique solution, the temperature on one of the boundary needs to be provided.