Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Coupled non-linear equations

  1. Dec 22, 2012 #1

    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:

    [itex]m_lc_l\frac{dT_l}{dt} = \widetilde{h}A_0(T_\infty - T_l) + \dot{m}_l(h_v-h_l)[/itex]
    [itex]\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 )[/itex]

    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!
  2. jcsd
  3. Dec 22, 2012 #2


    User Avatar
    2017 Award

    Staff: 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.
  4. Dec 22, 2012 #3
    Hey, the second equation is formed via simplification of the equation below (I think):

    [itex]\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}[/itex]

    I believe this equation and...

    [itex]m_lc_l\frac{dT_l}{dt} = \widetilde{h}A_0(T_\infty - T_l) + \frac{dm_l}{dt}(h_v-h_l)[/itex]

    ...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.
    Last edited: Dec 22, 2012
  5. Dec 24, 2012 #4
    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 (Text):
    // 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:
    [itex]\frac{dm_p}{dt}=-\pi d_p\rho D Sh\ln(1+B)[/itex], which is the d-squared law because the square of the droplet diameter decreases linearly in time. B is he Spalding mass transer coefficient [itex]B=\frac{Y_{F_s} - Y_{F_{\infty}}}{1-Y_{F_s}}[/itex]
    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!
  6. Dec 24, 2012 #5
    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.

  7. Dec 24, 2012 #6
    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 Files:

    • evap.jpg
      File size:
      52.9 KB
  8. Dec 24, 2012 #7
    So if I take the equations and combine density and other parameters into a single constant A, I have...

    [itex]\frac{dr_S}{dt} = -\frac{\dot{m}}{Ar_s^{2}}[/itex]

    ...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!
  9. Dec 24, 2012 #8


    User Avatar
    2017 Award

    Staff: 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.
  10. Dec 24, 2012 #9
    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)
    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.
  11. Dec 24, 2012 #10
    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:

    [itex]\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}[/itex]

    I also have this which is in terms of T and m:

    [itex]\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))[/itex]

    I do understand what I need, I am just a bit unsure on how to get there.

    Sorry if I'm being thick!
  12. Dec 24, 2012 #11
    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
    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.
  13. Dec 25, 2012 #12
    The basic form of the mass vaporization law:
    [itex]\dot m_p = 4 \pi r_p \frac{\lambda_g}{Cp_g}Nu_p \ln(1+B_T)[/itex]
    use [itex]m=\rho_p \frac{4}{3}\pi r_p^3[/itex] to get:
    [itex]\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}[/itex]

    so you arrive at:
    [itex]\frac{d d_p^2}{dt} = 2 \frac{\lambda_g}{\rho_p Cp_g}Nu_p\ln(1+B_T)[/itex]
    and the right hand side is only implicitly dependent on [itex]T_p[/itex].

    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:
    [itex]m_p Cp_p\frac{dT_p}{dt}=hA(T_g-T_p) - \dot m_p L[/itex]
    with [itex]h=\frac{Nu \lambda_g}{2r_p}[/itex] we get:
    [itex]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[/itex]
    [itex]\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}[/itex]

    Now introduce the stokes number:
    [itex]\tau_p=\frac{\rho_p (2r_p)^2}{18 \mu_g}[/itex]

    And 'eliminate' the diameter for easy integration:
    \frac{\dot m_p}{m_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 [itex]\dot m_p/m_p[/itex] 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 [itex]d_p^2[/itex].
    Last edited: Dec 25, 2012
  14. Dec 27, 2012 #13
    Excellent help thank you all so much!
  15. Dec 27, 2012 #14
    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!
  16. Dec 27, 2012 #15
    Use the current values of T_l and T_g (at time t) in the Hubbard rule.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook