Register to reply 
Coupled nonlinear equations 
Share this thread: 
#1
Dec2212, 10:55 AM

P: 7

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: [itex]m_lc_l\frac{dT_l}{dt} = \widetilde{h}A_0(T_\infty  T_l) + \dot{m}_l(h_vh_l)[/itex] [itex]\dot{m_l} = \frac{\pi d D_{AB}p}{R_vT_m}ln \left ( \frac{pp_{v_{\infty}}} {pp_{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
Dec2212, 04:30 PM

Mentor
P: 12,113

T_{l} and m_{l} are unknown?
What happens to p? Is that timedependent (and if yes, how?). If there is no analytic solution, a numerical simulation with small timesteps should work. 


#3
Dec2212, 05:33 PM

P: 7

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_vh_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. 


#4
Dec2412, 09:38 AM

P: 300

Coupled nonlinear 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 nonequilibrium evaporation models for manydroplet gasliquid Żow simulations, International Journal of Multiphase Flow 24 (1998) 10251055 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 RungeKutta 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 RungeKutta of course. mass is coupled to diameter and density, so you do something like this:
I don't see where your mass vaporization law comes from. It is usually derived by integrating the species convectiondiffusion 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 dsquared 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}}}{1Y_{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 ClausiusClaperon 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 coauthored: Roekaerts, Jenny, Beishuizen, Modeling of turbulent dilute spray combustion, Progress in Energy and Combustion Science, Volume 38, Issue 6, December 2012, Pages 846887 I'm assuming you're a student, so you should have access to all the papers I've mentioned. Have fun! 


#5
Dec2412, 02:16 PM

Mentor
P: 5,406

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 mdot, 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 


#6
Dec2412, 02:51 PM

P: 7

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! 


#7
Dec2412, 03:11 PM

P: 7

[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! 


#8
Dec2412, 03:35 PM

Mentor
P: 12,113

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 timederivative) as you wrote in post 7. This can be used everywhere to substitute m_dot. 


#9
Dec2412, 05:42 PM

Mentor
P: 5,406

dm/dt = m_dot, and integrating to get m. Then after taking a time step, calculate the new value of r_{s} 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. 


#10
Dec2412, 06:42 PM

P: 7

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_vh_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! 


#11
Dec2412, 08:36 PM

Mentor
P: 5,406

What you have is perfectly OK. m is the same as m_{l}.
Before you calculate the derivatives, you might consider first calculating m_dot and also r. Then your equation for dm_{l}/dt is dm_{l}/dt = m_dot Then you calculate dT_{l}/dt The next step is to take a time step in both dependent variables, say using forward Euler. Then loop back to the beginning. 


#12
Dec2512, 03:27 PM

P: 300

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_gT_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_gT_p)  \dot m_p L[/itex] or: [itex]\frac{dT_p}{dt}=Nu \frac{6 \lambda_g}{ Cp_g (2r_p)^2}(T_gT_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: [itex] \frac{dT_p}{dt} = \frac{Nu}{3Pr} \frac{Cp_g}{Cp_p} \frac{T_gT_p}{\tau_p}  \frac{\dot m_p}{m_p} \frac{L}{Cp_p}[/itex] 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 dsquared 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]. 


#13
Dec2712, 06:51 AM

P: 7

Excellent help thank you all so much!



#14
Dec2712, 07:59 PM

P: 7

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! 


Register to reply 
Related Discussions  
Numerical Simultaneous Solution of NonLinear Coupled Equations  General Math  14  
Coupled First Order Equations  Calculus & Beyond Homework  9  
System of second order linear homogenous differential coupled equations  Differential Equations  8  
Coupled Differential Equations  Differential Equations  2  
Coupled non linear ordinary differential equations.  Differential Equations  1 