## Coupled non-linear equations

Hi,

I'd really appreciate any help on this as I've spent many many hours trying to get my head around it. I am simulating the vaporization of a droplet in hot air and I have two equations to use:

$m_lc_l\frac{dT_l}{dt} = \widetilde{h}A_0(T_\infty - T_l) + \dot{m}_l(h_v-h_l)$
$\dot{m_l} = -\frac{\pi d D_{AB}p}{R_vT_m}ln \left ( \frac{p-p_{v_{\infty}}} {p-p_{v_{0}}} \right ) \left(2 + 0.6Re_d^{1/2}Sc^{1/3} \right )$

They are coupled nonlinear - I want to get the droplet diameter, mass, temperature etc. as a function of time so I can graph them. I've tried using MATLAB (ode45 function) but I think I'm missing a fundamental piece of knowledge regarding the maths.

If anyone could summarize the steps needed to obtain the variables as a function of time I would really appreciate it! I'm not asking anyone to do the work for me!

 PhysOrg.com science news on PhysOrg.com >> Hong Kong launches first electric taxis>> Morocco to harness the wind in energy hunt>> Galaxy's Ring of Fire
 Mentor Tl and ml are unknown? What happens to p? Is that time-dependent (and if yes, how?). If there is no analytic solution, a numerical simulation with small timesteps should work.
 Hey, the second equation is formed via simplification of the equation below (I think): $\frac{dm_l}{dt} = \left ( \frac{\rho _l \pi d^3}{2} \right )\frac{d(d)}{dt} + \left ( \frac{\pi d^3}{6} \right )\frac{d\rho_{l}}{dT_l}\frac{dT_l}{dt}$ I believe this equation and... $m_lc_l\frac{dT_l}{dt} = \widetilde{h}A_0(T_\infty - T_l) + \frac{dm_l}{dt}(h_v-h_l)$ ...are the coupled ones. My mistake! This shows what is and isn't variable with time. I would have initial conditions to work with, and I want to implement a loop that will then store values into an array (to plot). Many thanks.

## Coupled non-linear equations

The classic papers that you should have read:

Abramzon and Sirignano - Droplet vaporization model for spray combustion simulations, Int. J. Heat Mass Transfer 32(9) 1989
Miller, Harstad, Bellan - Evaluation of equilibrium and non-equilibrium evaporation
models for many-droplet gas-liquid Żow simulations, International Journal of Multiphase Flow 24 (1998) 1025-1055
And in case your particles are moving in a turbulent flow:

Beacuse the equations are transient, it doesn't matter (much) that they are nonlinear. Just use a standard forward Euler or Runge-Kutta scheme. Actually, for small timesteps you can assume all coefficients in the temperature equation constant and derive the exact solution. This is the exact Euler integration scheme and is much more robust than standard Euler. You can also use Runge-Kutta of course.

mass is coupled to diameter and density, so you do something like this:
Code:
// calculate change in radius dr, new radius r_p and new mass m
dr = m_dot/(4.0*M_PI*rho_l(T_s)*r_p*r_p);
r_p = r_p + dr*dt;
m = 4.0/3.0*M_PI*r_p*r_p*r_p * rho_l(T_s);

I don't see where your mass vaporization law comes from. It is usually derived by integrating the species convection-diffusion equation from the droplet surface to infinity, using the boundary conditions and then integrating again. you should get something like:
$\frac{dm_p}{dt}=-\pi d_p\rho D Sh\ln(1+B)$, which is the d-squared law because the square of the droplet diameter decreases linearly in time. B is he Spalding mass transer coefficient $B=\frac{Y_{F_s} - Y_{F_{\infty}}}{1-Y_{F_s}}$
D is the binary diffusion coefficient of the fuel vapor and air.
The fuel mass fraction at the surface of the droplet is obtained from the Clausius-Claperon equation. The Sherwood number for a sphere is known from experiments and is a function of Reynolds number and Schmidt number.

I also would like to point your attention to a recent paper I've co-authored:
Roekaerts, Jenny, Beishuizen, Modeling of turbulent dilute spray combustion, Progress in Energy and Combustion Science, Volume 38, Issue 6, December 2012, Pages 846-887
I'm assuming you're a student, so you should have access to all the papers I've mentioned. Have fun!

 Recognitions: Gold Member In repsonse #3, Dillman put his finger on your difficulty. Also, for the problem you are solving, the temperature dependence of the density is almost certainly negligible. Thus, you get the rate of change of diameter with respect to time in terms of m-dot, which is known explicitly in terms of the other parameters from your second equation. This provides closure on the problem. You have two coupled time dependent ordinary differential equations in two unknowns, T and d. These can be solved numerically. Chet
 Thank you for your response! I really appreciate it. I am getting really into this work (it's for an important project) but as it is the holidays I cannot ask for assistance from my supervisor. I have attached an image of the process I think is correct. I would really appreciate it if you could have a look at it for me. I have a few questions regarding the process. My equations are included in the attached image. Am I correct in saying m_dot is a function of T_s (m_dot(T_s))? Evaluating values at T_s,0 (initial surface temperature), and then using revised values on each loop? I have trouble evaluating B_m in this case too, but I understand the principle. Revising the value of r_s makes sense to me, I think that part is correct? Calculating the new value of T_s is where my major problem lies. Can I take m_dot to be constant here, i.e use the value obtained in the first step in the attached image, and then integrate the whole thing and use whatever value of t I am on (t = t + dt)? I would then use this new value of T_s (to revise the values of Re, Sc, Y_Fs etc.) along with new r_s to get a new m_dot and repeat the cycle. Many thanks! Attached Thumbnails

 Quote by Chestermiller In repsonse #3, Dillman put his finger on your difficulty. Also, for the problem you are solving, the temperature dependence of the density is almost certainly negligible. Thus, you get the rate of change of diameter with respect to time in terms of m-dot, which is known explicitly in terms of the other parameters from your second equation. This provides closure on the problem. You have two coupled time dependent ordinary differential equations in two unknowns, T and d. These can be solved numerically. Chet
So if I take the equations and combine density and other parameters into a single constant A, I have...

$\frac{dr_S}{dt} = -\frac{\dot{m}}{Ar_s^{2}}$

...which is the first equation. So I then use my m_dot equation to get an actual value for m_dot? I'm worried that m_dot is dynamic and changes with temperature.

I need to solve the above equation and temperature equation together using Runge Kutta or Euler, correct? Leaving dm/dt in there, would the resulting solution eliminate the need for me to revise the values in m_dot over and over?

I'm worried I have a fundamental misunderstanding of this type of problem! I am mainly confused about what happens with m_dot.

Thank you!

 Mentor If you cannot eliminate m_dot, you can keep track of m, temperature and radius at the same time. I would expect that you can express m fully as function of the radius, this gives m_dot as function of the radius (and its time-derivative) as you wrote in post 7. This can be used everywhere to substitute m_dot.

Recognitions:
Gold Member
 Quote by Dilman Thank you for your response! I really appreciate it. I am getting really into this work (it's for an important project) but as it is the holidays I cannot ask for assistance from my supervisor. I have attached an image of the process I think is correct. I would really appreciate it if you could have a look at it for me. I have a few questions regarding the process. My equations are included in the attached image. Am I correct in saying m_dot is a function of T_s (m_dot(T_s))? Evaluating values at T_s,0 (initial surface temperature), and then using revised values on each loop? I have trouble evaluating B_m in this case too, but I understand the principle. Revising the value of r_s makes sense to me, I think that part is correct? Calculating the new value of T_s is where my major problem lies. Can I take m_dot to be constant here, i.e use the value obtained in the first step in the attached image, and then integrate the whole thing and use whatever value of t I am on (t = t + dt)? I would then use this new value of T_s (to revise the values of Re, Sc, Y_Fs etc.) along with new r_s to get a new m_dot and repeat the cycle. Many thanks!
This should work, although it is not exactly how I would have done it. You can simplify things a little by writing

dm/dt = -m_dot, and integrating to get m. Then after taking a time step, calculate the new value of rs from the new value of m. Also, as I said earlier, the effect of temperature on density is not going to be important in this problem.

Work the equations into a form so that you have two coupled simultaneous first order ordinary differential equations:

dT/dt = A (m,T)
and
dm/dt = B (m,T)

The rest of the equations are algebraic, and can be used to establish the values of A(m,T) and B(m,T). The numerical method of solution you proposed is Forward Euler, and gives accurate results if dt is sufficiently small. Runge Kutta is another possibility, and can use larger values of dt. There are also automatic equation solver packages available (for free) that solve sets of coupled simultaneous first order ordinary differential equations using techniques that automatically control error tolerance. I like these quite a bit, and almost always use them. The IMSL library has a package that should suit you needs.

 Thanks for the responses. I'm not really sure which formula I need to integrate to get m(t)? I replaced dm/dt in the temperature equation so I have everything in terms of T and m: $\frac{dT_l}{dt} = \frac{{}\widetilde{h}A_0(T_\infty - T_l) -\left (2\pi \rho _gD_g\left ( \frac{3m_l}{4\pi \rho _l} \right )^{1/3}Sh^*(T_l) ln(1+B_m(T_l)) \right )(h_v-h_l)}{m_lC_l}$ I also have this which is in terms of T and m: $\frac{dm_l}{dt} = -2\pi \rho_gD_g \left (\frac{ 3m_l}{4\pi\rho_l} \right )^{1/3}Sh^*(T_l)ln(1+B_m(T_l))$ I do understand what I need, I am just a bit unsure on how to get there. Sorry if I'm being thick!
 Recognitions: Gold Member What you have is perfectly OK. m is the same as ml. Before you calculate the derivatives, you might consider first calculating m_dot and also r. Then your equation for dml/dt is dml/dt = -m_dot Then you calculate dTl/dt The next step is to take a time step in both dependent variables, say using forward Euler. Then loop back to the beginning.
 The basic form of the mass vaporization law: $\dot m_p = 4 \pi r_p \frac{\lambda_g}{Cp_g}Nu_p \ln(1+B_T)$ use $m=\rho_p \frac{4}{3}\pi r_p^3$ to get: $\frac{dm_p}{dt}=4 \pi r_p^2 \rho_p \frac{dr_p}{dt}=\frac{\pi}{2}r_p\frac{d_p^2}{dt}$ so you arrive at: $\frac{d d_p^2}{dt} = 2 \frac{\lambda_g}{\rho_p Cp_g}Nu_p\ln(1+B_T)$ and the right hand side is only implicitly dependent on $T_p$. The temperature equation: Direct application of the balance between cooling due to evaporation (with the latent heat of evaporation L) and heating due to convection leads to: $m_p Cp_p\frac{dT_p}{dt}=hA(T_g-T_p) - \dot m_p L$ with $h=\frac{Nu \lambda_g}{2r_p}$ we get: $m_p Cp_p\frac{dT_p}{dt}=Nu \frac{\lambda_g}{2r_p}\pi r_p^2(T_g-T_p) - \dot m_p L$ or: $\frac{dT_p}{dt}=Nu \frac{6 \lambda_g}{ Cp_g (2r_p)^2}(T_g-T_p) - \frac{\dot m_p}{m_p} \frac{L}{Cp_p}$ Now introduce the stokes number: $\tau_p=\frac{\rho_p (2r_p)^2}{18 \mu_g}$ And 'eliminate' the diameter for easy integration: $\frac{dT_p}{dt} = \frac{Nu}{3Pr} \frac{Cp_g}{Cp_p} \frac{T_g-T_p}{\tau_p} - \frac{\dot m_p}{m_p} \frac{L}{Cp_p}$ You can also rewrite the mass vaporization using the stokes number. Now you can do the following, although other alternatives are possible: calculate a solution on timestep t+1 for the d-squared equation to get the diameter calculate the mass calculate the change in mass using the difference in mass on timestep t and timestep t+1 calculate the temperature at timestep t+1 Hope this helps. edit: you can rewrite $\dot m_p/m_p$ in the temperature equation using an expression similar to what you have on the right hand side of the mass vaporization law and you will have a dependence only on $d_p^2$.
 Excellent help thank you all so much!
 I've gone through the steps of derivation and I feel like I understand the whole thing much more now! Can anybody tell me if Nu, B_T, L and Pr should be evaluated at a reference temperature (e.g. Hubbard 1/3 rule: T_r = T_l + (1/3)(T_g - T_l)) or if they have to be dynamic functions of temperature? (involves long winded calculations). In fact, should all the parameters (apart from d, r, m_dot, T_p) in the final two equations be evaluated and treated constant beforehand? Many many thanks!
 Recognitions: Gold Member Use the current values of T_l and T_g (at time t) in the Hubbard rule.