0
$\begingroup$

I have the following 1D (in space, that is) system of equations I would like to solve:

\begin{equation} \rho_{fs}\frac{\partial x_{fs}}{\partial t} = h_m\left(W_a - W_{fs}\right) - D_{eff}\left(\frac{\partial \rho_v(T_{fr})}{\partial x} \right) \end{equation} \begin{equation} \rho_{fr}c_{p,fr}\frac{\partial T}{\partial t} = \frac{\partial}{\partial x} \left( k_{fr} \frac{\partial T}{\partial x} \right) + L_{sub} \frac{\partial \rho_{fr}}{\partial x} \end{equation} \begin{equation} \frac{\partial \rho_{fr}}{\partial t} = \frac{\partial}{\partial x}\left( D_{eff} \frac{\partial \rho_v(T_{fr})}{\partial x} \right) \end{equation}

This set of equations describes the time evolution of a frost layer's height, $x_{fs}$, (local) temperature, $T_{fr}$, and (local) density, $\rho_{fr}$.

I have been struggling with trying to solve this system. I have constructed a system of ODEs following the method of lines and attempted to solve it in Matlab using ode45, but this resulted in matrices filled with NaN values.

Following this, I am now trying to solve it through a 1st order upwind scheme for the time derivatives, solving one equation at a time and trying to iterate until the differences are small enough:

\begin{equation} \frac{\left(x_{fs}^{n+1} - x_{fs}^n\right)}{\Delta t} = \frac{1}{\rho_{fs}}\left(h_m\left(W_a - W_{fs}\right) - D_{eff}\frac{\rho_v\left(T^n_{i+1}\right) - 2\rho_v\left(T^n_{i}\right) + \rho_v\left(T^n_{i-1}\right)}{(\Delta x)^2} \right) \end{equation} \begin{equation} \frac{\left(T_i^{n+1} - T_i^n\right)}{\Delta t} = \frac{1}{\rho_{fr,i}c_{p,fr}}\left( k_{fr}\frac{T^n_{i+1} - 2T^n_{i} + T^n_{i-1}}{(\Delta x)^2} + L_{sub}D_{eff}\frac{\rho_v\left(T^n_{i+1}\right) - 2\rho_v\left(T^n_{i}\right) + \rho_v\left(T^n_{i-1}\right)}{(\Delta x)^2} \right) \end{equation} \begin{equation} \frac{\left(\rho_i^{n+1} - \rho_i^n\right)}{\Delta t} = D_{eff}\frac{\rho_v\left(T^n_{i+1}\right) - 2\rho_v\left(T^n_{i}\right) + \rho_v\left(T^n_{i-1}\right)}{(\Delta x)^2} \end{equation}

In order to adhere to boundary conditions, the temperature is calculated on the $n_x$ number of nodes, and the density is calculated inbetween, in the $n_x - 1$ number of cells. As for initial and boundary conditions, the temperature at the frost surface and at the bottom of the layer are known, and an initial uniform density profile is assumed due to the very small initial frost height mentioned before.

This approach has not yielded good results yet either, although I am now starting to think it might have to do with my choice of $\Delta t$. The spatial step size $\Delta x$ is set by the layer height $x_{fs}$ and the number of nodes. The starting layer height is very small, in the order of $10^{-4}$. If I understand correctly, this means my $\Delta t$ has to be very small as well in order to satisfy the Courant condition? This would be a bit impractical though, as the simulated time is in the order of tens of minutes or even hours.

I am wondering if there is a better approach to take to solve this system of equations, or if I am making a mistake somewhere. I would appreciate any insights. Apologies if I have made any mistakes in posing my question here.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Have you tried cautiously debugging your initial fully- mcoupled implementation? Maybe there was only a small mistake... Aso, regarding time integration, a good way to check that your model behaves correctly is to use adaptive time stepping. This way, the integrator will automatically find stable time steps. If the solution fails, then it is likely that your implemented system is singular ou too stiff. By the way, in 1D, you can very efficiently use implicit integrators to lift stability issues. $\endgroup$
    – Laurent90
    Dec 16, 2022 at 7:09
  • $\begingroup$ Thank you for your reply. I have gone over my implementation multiple times at this point, but I will give it another shot. One thing I am not too sure about is defining the initial temperature field. So far I have mostly used a linearl profile between the frost surface and wall temperatures, which of course reduces the second order approximation to zero, but I was hoping the non-linear $\rho_v$ field would compensate for it. Furthermore, could you elaborate on the implicit integration? I am not sure how I would go about implementing something like that. Thank you. $\endgroup$
    – HVdB
    Dec 16, 2022 at 10:07
  • $\begingroup$ NaNs are actually quite nice, better than wrong values! That's because you can step through your program and find the point where they are appearing, and then start reasoning why that is so. $\endgroup$
    – Wolfgang Bangerth
    Dec 17, 2022 at 2:30
  • $\begingroup$ Separately, take a look here: arxiv.org/abs/2209.04198 $\endgroup$
    – Wolfgang Bangerth
    Dec 17, 2022 at 2:31
  • $\begingroup$ you have discretized the last term in the 2nd eqn. incorrectly. also, there is no D_eff coefficient in that term. $\endgroup$
    – Sthavishtha Bhopalam
    Dec 22, 2022 at 20:16

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.