Why Does Increasing Final Run Time Cause Divergence in Finite Differencing?

[member 428835]

I'm hoping that you are only solving the equation over the domain from 0 to 1, and not from -1 to + 1.

Absolutely!

And I'm hoping that, at y = 0, you are solving the differential equation for y(0) using a reflective boundary condition, rather than just setting y(0) equal to y(Δz)

I'm a little unsure of what you are saying here. Can you specify please?

 

[Chestermiller]

Absolutely!
I'm a little unsure of what you are saying here. Can you specify please?

$$\left(\frac{\partial y}{\partial t}\right)_Z=Z\left(\frac{\partial y}{\partial Z}\right)_t\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{\partial ^2(y^{3/2})}{\partial Z^2}$$
At z = 0, $\frac{\partial y}{\partial Z}=0$, then $y^{3/2}(-ΔZ)=y^{3/2}(+ΔZ)$. This means that:
$$y^{3/2}(-ΔZ)-2y^{3/2}(0)+y^{3/2}(+Δz)=2[y^{3/2}(+Δz)-y^{3/2}(0)]$$
So, at Z = 0,
$$\frac{y(0,t+Δt)-y(0,t)}{Δt}=\frac{4}{3L^2}\frac{[y^{3/2}(t,+Δz)-y^{3/2}(t,0)]}{(ΔZ)^2}$$

This is the finite difference equation that should be used at Z = 0.

Chet

 

[member 428835]

I'm with you! Thanks for clearing that up! Ok, so since our thoughts have been all over the place (mainly due to my questions) I've decided to type up a latex document. I'll post the pdf when I finish it for both your review and others who may want to see the answer but don't want to look through all the conversation. Thanks a ton for your help. I'll probably be writing on this thread before the pdf is finished, as I want to ensure the FD is correct.

 

[Chestermiller]

Please note also, that, as things stand now, you are still using Forward Euler to integrate with respect to time. If stability is an issue, you can control it by taking the time step small enough. But please note that, in terms of stability, Forward Euler is the worst scheme that one can use. So, if the time step is too large, the scheme will go unstable. Stability problems can typically be removed by going to so-called implicit schemes, such as Backward Euler, in which the right hand side of the equation is evaluated at t + Δt, rather than t. But this requires the solution of a set of coupled algebraic equations at each time step. Hopefully, for your problem, Forward Euler will be adequate. But, if not, I can help.

Chet

 

[member 428835]

Chet, I have a quick question about taking a finite difference over non-linear terms. Notice the following $$ \frac{\partial (h^2)}{\partial t} = 2h\frac{\partial h}{\partial t}$$ Applying a finite difference to each yields $$\frac{\partial (h^2)}{\partial t} = \frac{h^2(t+\Delta t) - h^2(t)}{\Delta t}$$ but also $$ 2h\frac{\partial h}{\partial t} = 2 h(t) \frac{h(t+\Delta t)-h(t)}{\Delta t} $$ Clearly the finite differences are not equal. Can you explain what is happening?

 

[member 428835]

Attached is the finite difference technique along with how we arrived. Please review it, especially the end, and let me know what you think. The conclusion doesn't look right. EDIT: the limits should be from the left, but you get the idea.

 

[Chestermiller]

Chet, I have a quick question about taking a finite difference over non-linear terms. Notice the following $$ \frac{\partial (h^2)}{\partial t} = 2h\frac{\partial h}{\partial t}$$ Applying a finite difference to each yields $$\frac{\partial (h^2)}{\partial t} = \frac{h^2(t+\Delta t) - h^2(t)}{\Delta t}$$ but also $$ 2h\frac{\partial h}{\partial t} = 2 h(t) \frac{h(t+\Delta t)-h(t)}{\Delta t} $$ Clearly the finite differences are not equal. Can you explain what is happening?

These are two slightly different first order finite difference approximations to the exact same quantity, and, as such, neither one of them is favored over the other. In each case, the discretization error decreases linearly with Δt when Δt becomes small. To see that, expand h(t+Δt) in a Taylor series about t, and substitute the resulting series into each of the two finite difference approximations.

Chet

 

[Chestermiller]

Attached is the finite difference technique along with how we arrived. Please review it, especially the end, and let me know what you think. The conclusion doesn't look right. EDIT: the limits should be from the left, but you get the idea.

Corrections:

Eqn. 16 should be the derivative of h³, not h.

In Eqn. 25, there should be a 2 in the numerator, not a 4. I would also prefer to see a 4 in the denominator, but I can't strongly argue against the 4.

Eqn. 26 is incorrect. The second derivative of h³ with respect to Z has been evaluated incorrectly.

I didn't check Eqns. 28 and 29.

Josh,

I've had a huge amount of practical experience with numerical analysis, particularly on problems like this one. And that practical experience involved many situations in which a lot of money was riding on the outcome of the correct numerical solution of the equations. Based on this experience, I am very strongly urging you to solve this problem numerically in terms of y rather than in terms of h, for the following reasons:

1. There are many fewer mathematical operations involved in solving it this way
2. There is much less of a chance of making a coding error
3. The finite difference equations that result in this approach automatically conserve volume, unlike the scheme involving h, in which volume is conserved only in the limit of vanishing $\Delta Z$

If I were coding this problem, starting from the beginning of each time step (so that y is known at all grid points at time t), I would carry out the following sequence of steps:

1. Calculate y^3/2 at all the grid points
2. Calculate √y at the grid point adjacent to the right boundary. This is h at that location.
3. Calculate L²(t+Δt) and (1/L)dL/dt
4. Evaluate the right hand side of the main differential equation at each grid point
5. Execute a time step at each grid point

Chet

 

[member 428835]

These are two slightly different first order finite difference approximations to the exact same quantity, and, as such, neither one of them is favored over the other. In each case, the discretization error decreases linearly with Δt when Δt becomes small. To see that, expand h(t+Δt) in a Taylor series about t, and substitute the resulting series into each of the two finite difference approximations.
Chet

Ok, this is cool! I just did the expansion and it works. Awesome!

 

[member 428835]

Josh,

I've had a huge amount of practical experience with numerical analysis, particularly on problems like this one. And that practical experience involved many situations in which a lot of money was riding on the outcome of the correct numerical solution of the equations. Based on this experience, I am very strongly urging you to solve this problem numerically in terms of y rather than in terms of h
Chet

I'll definitely do it this way now that I see (from your previous post) that the expansion is the same! This was not immediately obvious to me, hence why I never introduced the $y$ transform. Thanks so much fro your help, Chet! Stay tuned, as I'll fix the corrections, introduce the transform, and also post a flow chart for the code (although you've seem to have already done this).

You're awesome!

 

[member 428835]

Ok, I think I corrected all the errors you pointed out. How does this attached version look?

As an aside, do we really need (12) through (16)? They seem to not really be used.

 

[Chestermiller]

You're awesome!

Not really. Just lots of experience.

 

[Chestermiller]

Ok, I think I corrected all the errors you pointed out. How does this attached version look?

As an aside, do we really need (12) through (16)? They seem to not really be used.

Yes. I would lose 12 - 16.

There's a typo in Eqn. 29: the y should be under the radical in the numerator.

Please show the derivation of Eqns. 30 and 34-36. I'm not able to follow in detail. I just want to make sure that it was done correctly. For example, I think there should be a 1/2 in Eqn. 30 instead of the factor of 2; and, in Eqn. 34, I feel like there should be a 4 in the denominator. Also, in Eqn. 36, there should be a ζ in the denominator rather than a ξ.

Chet

 

[member 428835]

How does this look?

(23) is evaluated at $t$ and yet requires $y$ to be evaluated at $t$ as well. How would you deal with this?

I also got rid of zeta and xi.

 

[Chestermiller]

How does this look?

I also got rid of zeta and xi.

I found a couple of errors in the new version:
Eqns. 17,19, 22,23,33,34,and 35 should have a 4 instead of a 2.

In Eqn. 32, you are missing an additional factor of 2 in the denominator of the first term in brackets. This is because ∂y/∂Z is evaluated numerically using a central difference over 2 grid intervals.

(23) is evaluated at $t$ and yet requires $y$ to be evaluated at $t$ as well. How would you deal with this?

In Eqns. 33-34, you need to replace the j's with (j+1) and you need to replace the (j-1)'s with j. And in Eqn. 35, you need to replace the (j-1) with j (so you have j on both sides). As far as the L_j² in the denominator of Eqn. 32 is concerned, you have already calculated that from the previous time step (Or, to the same degree of accuracy, you could also use L_j+1² or the average of the two).

 

[member 428835]

Thanks for pointing these errors out. I agree with all of your corrections, and believe this is now error free. Thanks for all the help!

How would I go about an error analysis though? It's non-linear so I don't think a von Neumann approach is valid since separation of variables doesn't work.

I'll start coding this too (in MatLab) and post a few pictures of what I'm getting in the .pdf for you're consideration.

 

[Chestermiller]

Thanks for pointing these errors out. I agree with all of your corrections, and believe this is now error free. Thanks for all the help!

How would I go about an error analysis though? It's non-linear so I don't think a von Neumann approach is valid since separation of variables doesn't work.

I'll start coding this too (in MatLab) and post a few pictures of what I'm getting in the .pdf for you're consideration.

Are you talking about error analysis, or are you talking about stability?

Chet

 

[member 428835]

Are you talking about error analysis, or are you talking about stability?

Chet

Sorry, I'm talking about stability.

 

[Chestermiller]

Sorry, I'm talking about stability.

I don't know of any methods of analyzing the stability of time-dependent schemes in non-linear problems, but, experience has shown that stability is guaranteed if you go to fully implicit (e.g., Backward Euler). It is also possible to increase stability by going partially implicit (i.e., more implicit) by handling certain quantities implicitly, while treating others explicitly.

Solving the problem implicitly means that you have to solve a set of non-linear equations at each time step. The easiest method of solving a problem like this implicitly is to use the so-called "method of lines" in conjunction with a stiff ODE software package based on the C. W. Gear method. Such a package transparently solves the set on non-linear equations at each time step. The method of lines changes the non-linear PDE into a non-linear set of coupled ODEs, by finite differencing the spatial derivatives, but not the time derivative. So there is an ODE for each grid point. The stiff package automatically integrates these stiff ODEs.

A partially implicit method can be developed for your problem that requires a tri-diagonal set of linear equations to be solved at each time step. This can be done using a tri-diagonal matrix inversion routine. However, there will be some quantities in the difference equations that will have to be handled explicitly in this approach. Still, in my judgement, it would improve stability very substantially (and allow much larger time steps).

The first thing to do is to see whether the existing scheme gives you stability when small (but acceptable) time steps are used.

Chet

 

[member 428835]

Thanks a ton! So I just finished writing the code. It's at the end of the attached pdf. I put plots in the pdf for a time value of 14 (line 8 in the code) and it ran well. However, if I set $tf = 15$ I get some errors, specifically imaginary parts. My guess is from line 37/38, when the square root is introduced. I should say that if I change line 9 to $N = (tf/10)*M^2;$ I can run for $tf$ values of 100 and more.

What is your take on this technique? Is this acceptable or should I go for an implicit scheme?

Also, for your curiosity, the percent change in volume was less than 1.5% for the code and pictures I posted. Changing line 9 as suggested and running for $tf=100$ produces a volume change of about 0.003%.

 

[Chestermiller]

Thanks a ton! So I just finished writing the code. It's at the end of the attached pdf. I put plots in the pdf for a time value of 14 (line 8 in the code) and it ran well. However, if I set $tf = 15$ I get some errors, specifically imaginary parts. My guess is from line 37/38, when the square root is introduced. I should say that if I change line 9 to $N = (tf/10)*M^2;$ I can run for $tf$ values of 100 and more.

What is your take on this technique? Is this acceptable or should I go for an implicit scheme?

Also, for your curiosity, the percent change in volume was less than 1.5% for the code and pictures I posted. Changing line 9 as suggested and running for $tf=100$ produces a volume change of about 0.003%.

I would say it looks pretty OK. I haven't checked the coding, but the solution looks like you expected, correct?

Regarding the need to go implicit, my question is "is computation time a problem?" If not, then I see no need to.

How come you are using such a large M. I would use a value of 20, and then check the accuracy of the solution by redoing the calculation at M = 40. But, I think 20 will be adequate. That should help tremendously with your stability.

What does L(t) look like?

Chet

 

[member 428835]

I haven't checked the coding, but the solution looks like you expected, correct?

Yes, the solution looks good.

Regarding the need to go implicit, my question is "is computation time a problem?" If not, then I see no need to.

Computation time is not an issue.

How come you are using such a large M. I would use a value of 20, and then check the accuracy of the solution by redoing the calculation at M = 40. But, I think 20 will be adequate. That should help tremendously with your stability.

I wanted a large M to reduce error. Having M=20 yields a final $L$ of about 8.4. Changing M=40 yields a final $L$ of about 8.75.

What does L(t) look like?
Chet

I'll attach a pdf.

When changing M=20 and N=M^2 the maximum final run time before crashing is only 20. If I go to M=21 I get complex y values. However, if I let $N=(tf/10)*M$ I can go further in time. What are your thoughts here?

 

[Chestermiller]

I'm a little confused over N. My understanding is that, when you increase N, it makes Δt small. This should increase stability. Yet when you switched from N=M² to $N=(t_f/10)M^2$, you basically reduced N by a factor of about 0.7 (tf went to 100 from 14). This should have made the calculation less stable instead of more stable. Also, when you went to smaller values of M, it should have increased stability at the same value of tf. To make this all easier to follow, you should calculate the value of $\frac{t_fM^2}{(N-1)L^2}$ for each case, and use this to try to map out the region where the solution is stable.

Chet

 

[member 428835]

Let me clarify. For $N=M^2$ and $tf=19$ the code runs well when $M=20$. However, letting $M=40$, all else constant, crashes the code. However, for $N=10*M^2$ and $M=40$ at $tf$ can reach values of 150 without crashing the code. Volume changes at this level are 9.4%.

To make this all easier to follow, you should calculate the value of $\frac{t_fM^2}{(N-1)L^2}$ for each case, and use this to try to map out the region where the solution is stable.
Chet

Why do you choose to use this value?

 

[Chestermiller]

Why do you choose to use this value?

This is Δt/(Δz)² which, aside from the diffusivity factor, is the standard stability parameter for a parabolic equation.

Chet

 

[member 428835]

Thanks. So I checked and the code is roughly stable for $\Delta t / \Delta x^2 < .6$. Additionally, $\Delta t / \Delta x^2$ values of about .65 and up crash. What do you think about this?

As an aside, could you show me where are we actually conserving volume in the pdf? I ask because I was also going to try the same differential equation and the same tip boundary condition but instead of conserving volume apply the following integral constraint $$
-2 \int_0^{tf} \left.h^2 h_z\right|_{z=0} \,d t =\int_0^{L(tf)} \left. h^2 \right|_{t=tf} \, d z$$ I can discretize this for the new $L$. What do you think? Won't this boundary condition affect the previous $L$ derived?

 

[Chestermiller]

The way we have it set up right now, it doesn't automatically conserve volume. But there is a way of setting it up so that it does. Let's start with the differential equation:
$$\frac{\partial h^2}{\partial t}=Z\left(\frac{\partial h^2}{\partial Z}\right)\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{\partial ^2h^3}{\partial Z^2}$$
If we multiply both sides by L, and re-express the first term on the right hand side using integration by parts, we obtain:
$$L\frac{\partial h^2}{\partial t}=\left[\frac{\partial (Zh^2)}{\partial Z}-h^2\right]\frac{dL}{dt}+\frac{2}{3L}\frac{\partial ^2h^3}{\partial Z^2}$$
If we next mover the second term in brackets to the left hand side of the equation, we obtain:
$$\frac{\partial (Lh^2)}{\partial t}=\frac{\partial (Zh^2)}{\partial Z}\frac{dL}{dt}+\frac{2}{3L}\frac{\partial ^2h^3}{\partial Z^2}$$
The finite difference form of this (that will automatically conserve volume) is:

$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{dL}{dt}+\frac{2}{3L}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}$
This yields:

$h^2(t+Δt,Z)=\frac{L(t)}{L(t+Δt)}h^2(t,Z)+[\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{4ΔZ}\frac{1}{L^2}\frac{dL^2}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}]Δt$

If I did this right, it should give automatic conservation of volume.

Chet

 

[member 428835]

If I did this right, it should give automatic conservation of volume.

I follow your math, but how does it now conserve volume? It seems like you have only written it differently, but have not introduced anything new. Am I missing something?

 

[Chestermiller]

I follow your math, but how does it now conserve volume? It seems like you have only written it differently, but have not introduced anything new. Am I missing something?

In my next-to-last equation, if you add the equations for all the grid points together, the right hand sides will cancel out. This will give you that the trapazoidal rule integration of h^2LdZ for the volume will have a time derivative of zero.

Chet

 

[member 428835]

In my next-to-last equation, if you add the equations for all the grid points together, the right hand sides will cancel out.
Chet

This is not immediately obvious to me. I'm happy to crunch the numbers, but I can't see the cancelation.

This will give you that the trapazoidal rule integration of h^2LdZ for the volume will have a time derivative of zero.
Chet

But h^2LdZ is not volume. Volume was h^2dZ.

 

[Chestermiller]

But h^2LdZ is not volume. Volume was h^2dZ.

No way. h^2dz=h^2LdZ is volume.

Chet

 

[Chestermiller]

This is not immediately obvious to me. I'm happy to crunch the numbers, but I can't see the cancelation.

It works kinda like this: (A-0) + (B-A) + (C-B) + (D-C) + (0 - D) = 0

 

[member 428835]

No way. h^2dz=h^2LdZ is volume.

Chet

Aaaaaaand now I feel stupid. Thanks for clarifying; this makes perfect sense now!

 

[member 428835]

$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{dL}{dt}+\frac{2}{3L}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}$
This yields:
$h^2(t+Δt,Z)=\frac{L(t)}{L(t+Δt)}h^2(t,Z)+[\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{4ΔZ}\frac{1}{L^2}\frac{dL^2}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}]Δt$

If I did this right, it should give automatic conservation of volume.

I think there is a mistake. I'll begin where you did
$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{dL}{dt}+\frac{2}{3L}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2} \implies$
$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{L Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2} \implies$
$h^2(t+Δt,Z)=\frac{L(t)}{L(t+Δt)}h^2(t,Z)+[\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{4ΔZ}\frac{1}{L^2}\frac{dL^2}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}]\frac{ΔtL}{L(t+Δt)}$

I checked back through the works but can't figure out if the $L$'s are evaluated at $t$ or $t+\Delta t$. I left it as simply $L$ since I am not sure. What do you think?

I now have the reflective boundary condition as $$y_{new}(1) = \frac{L_{old}}{L_{new}}y(1)+\left[ \frac{y(2)}{2 L^2}\frac{dL^2}{dt}-\frac{4 y^{3/2}(1)}{3 L^2 \Delta z^2}\right] \frac{L \Delta t}{L_{new}}$$ where 1 implies the first element of y. Do you agree with this? I should add that everywhere listed as simply $L$ I ran in the code as $L(t+\Delta t)$ yet volume isn't conserved. Any ideas?

 

[Chestermiller]

I think there is a mistake. I'll begin where you did
$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{dL}{dt}+\frac{2}{3L}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2} \implies$
$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{L Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2} \implies$
$h^2(t+Δt,Z)=\frac{L(t)}{L(t+Δt)}h^2(t,Z)+[\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{4ΔZ}\frac{1}{L^2}\frac{dL^2}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}]\frac{ΔtL}{L(t+Δt)}$

I checked back through the works but can't figure out if the $L$'s are evaluated at $t$ or $t+\Delta t$. I left it as simply $L$ since I am not sure. What do you think?

I now have the reflective boundary condition as $$y_{new}(1) = \frac{L_{old}}{L_{new}}y(1)+\left[ \frac{y(2)}{2 L^2}\frac{dL^2}{dt}-\frac{4 y^{3/2}(1)}{3 L^2 \Delta z^2}\right] \frac{L \Delta t}{L_{new}}$$ where 1 implies the first element of y. Do you agree with this? I should add that everywhere listed as simply $L$ I ran in the code as $L(t+\Delta t)$ yet volume isn't conserved. Any ideas?

Yes.

It doesn't matter what we do with the L's, since each of the spatial terms on the right hand side have to integrate to zero individually. The reason that they didn't in your calculation is that I used a finite difference approximation in the first term on the right hand side doesn't integrate to zero (my mistake). The approximation I should have used was:
$$\frac{\partial (Zh^2)}{\partial Z}=\frac{(Z+\frac{ΔZ}{2})\frac{(h^2(Z+ΔZ)+h^2(Z))}{2}-(Z-\frac{ΔZ}{2})\frac{(h^2(Z-ΔZ)+h^2(Z))}{2}}{ΔZ}=\frac{(Z+\frac{ΔZ}{2})(h^2(Z+ΔZ)+h^2(Z))-(Z-\frac{ΔZ}{2})(h^2(Z-ΔZ)+h^2(Z))}{2ΔZ}$$

If you try it with this, I guarantee that volume will be conserved.

Chet

 

[member 428835]

The approximation I should have used was:
$$\frac{\partial (Zh^2)}{\partial Z}=\frac{(Z+\frac{ΔZ}{2})\frac{(h^2(Z+ΔZ)+h^2(Z))}{2}-(Z-\frac{ΔZ}{2})\frac{(h^2(Z-ΔZ)+h^2(Z))}{2}}{ΔZ}=\frac{(Z+\frac{ΔZ}{2})(h^2(Z+ΔZ)+h^2(Z))-(Z-\frac{ΔZ}{2})(h^2(Z-ΔZ)+h^2(Z))}{2ΔZ}$$

If you try it with this, I guarantee that volume will be conserved.

Chet

I've redone the code with this expression and no luck. In fact, now the slope is not zero at Z=0. I've attached the code so you can see for yourself. Am I doing something terribly wrong here?

 

[Chestermiller]

I've redone the code with this expression and no luck. In fact, now the slope is not zero at Z=0. I've attached the code so you can see for yourself. Am I doing something terribly wrong here?

I'm going to look over the code and get back with you.

Chet

 

[Chestermiller]

Line 43 looks incorrect to me. I get:

ynew( 1 ) = Lold/Lnew*y ( 1 )+((y ( 2 )+y ( 1 ) ) /(4Lnew^2) * DL2+4/3*(y3( 2 )-y3 ( 1 )) /(Lnew^2 dz ^2) ) dt

My recommended corrections are in bold and italic.
Regarding line 44 (and maybe other lines), I definitely would not overwrite any of the y's until all the ynew's are calculated.

Regarding line 23, this is an incorrect application of the trapazoidal rule. The equation should read:

iv = iv+((h^2( i )+h^2( i+1) ) /2) dz

That's about all I checked so far.

Chet

 

[member 428835]

Line 43 looks incorrect to me. I get:

ynew( 1 ) = Lold/Lnew*y ( 1 )+((y ( 2 )+y ( 1 ) ) /(4Lnew^2) * DL2+4/3*(y3( 2 )-y3 ( 1 )) /(Lnew^2 dz ^2) ) dt

My recommended corrections are in bold and italic.
Regarding line 44 (and maybe other lines), I definitely would not overwrite any of the y's until all the ynew's are calculated.

Regarding line 23, this is an incorrect application of the trapazoidal rule. The equation should read:

iv = iv+((h^2( i )+h^2( i+1) ) /2) dz

That's about all I checked so far.

Chet

Awesome! Not sure how I missed this but it looks great now, and conserves volume! You've been a huge help, Chet! I'll post a final pdf when I finish it so you can see how it all looks!

 

[member 428835]

Hey Chet, I have a quick question. When transforming from $Z$ to $z$ I did the following: $Z=L(tf)z$, but now I'm not sure why I used the last value of $L$. Should $L$ be evaluated at the final time? Could you explain your reasoning?

 

[Chestermiller]

Hey Chet, I have a quick question. When transforming from $Z$ to $z$ I did the following: $Z=L(tf)z$, but now I'm not sure why I used the last value of $L$. Should $L$ be evaluated at the final time? Could you explain your reasoning?

You did it so you could always have a grid point exactly at the moving boundary, making it convenient to apply any boundary condition there. Also, you always have the same number of grid points, and they are uniformly distributed between the boundaries. You've taken the complexity out of applying the moving boundary condition and transferred it to the differential equation where it can be handled much more easily.

Chet

 

[member 428835]

Thanks! For some reason I was wondering which L I should use to transform back from Z to z but obviously it's the most recent L. Thanks Chet! You're definitely brilliant!

 

[member 428835]

Hey Chet, sorry to bother you again about this problem, but how much more effort would need to be put into the finite difference scheme to go implicit? I've never done this before and am curious about it.

 

[Chestermiller]

Hey Chet, sorry to bother you again about this problem, but how much more effort would need to be put into the finite difference scheme to go implicit? I've never done this before and am curious about it.

It depends. What you would do would be to only discretize the differential equations in the spatial direction, but not in the time direction. So then you would have a coupled set of non-linear first order ODEs in time for the y values at the M grid points (plus an additional first order ODE for L²). This approach is called the Method of Lines. To solve these equations implicitly, you would typically need a stiff ODE equation solver subroutine. The problem would be done in FORTRAN, and you would be using a commercially available stiff package like the ones available on IMSL or online (e.g., DASSL). Or, if Matlab has a stiff equation solver, you could use that.

Chet

 

[member 428835]

Thanks a ton Chet! I really appreciate your help.

 

[Strum]

Is #86 the newest code? I would like to look it over when it is fully done :)

 

[member 428835]

after chet's correction, yep!