search

murphyqm / pytesimal

Model the conductive cooling of small planetary bodies with temperature-dependent material properties
MIT License
1 stars 0 forks source link

Reproduce Haack et al. (1990)? #44

Open andreww opened 3 years ago

andreww commented 3 years ago

It would be nice, wouldn't it. I spent a bit of time seeing if I could reproduce the 'reference' model of Haack et al. (1990). It looks like this is non-trivial but could be done (if we allow some API breaking changes). (Not so) brief notes below.

The 'reference' model is quite complex. One issue is that it has source terms in the mantle. In principle these are easy to add (it's just another term to add to the temperature change here: https://github.com/murphyqm/pytesimal/blob/d166aba7adf0a1c9b700f6f398858477be3b159f/pytesimal/numerical_methods.py#L380 Assuming this defaults to zero and we can add to the API without loosing backwards compatibility this would be okay.

The values to feed into the source term are a bit more complex. There are two parts. A radiogenic component (which is different in the 'crust' and the 'mantle' - the only difference between the two regions) and a latent heat component (the model starts off above the solids for the mantle). Ultimately this is just a function of time, depth and temperature (which would return a heat generation rate, that we would divide by the heat capacity and multiply by the time step before adding this to the change in temperature).

A second issue is that we'll need a more complex core model. It is not totally clear to me from looking at the paper if the core is isothermal or if it is part of a single conduction problem. I think it is the latter. This wouldn't be a problem (it's just a different set of material properties and a different source term for the latent heat). This time there is no radiogenic heating, just a different latent heat term. I'm not sure how to handle this once it cools to the eutectic temperature (the model is cooling through a two phase region). Maybe this would need a separate core (with its own conductive solution and a zero flux boundary condition in the middle). Grrr.

Then the third issue (that actually made me stop) is that as far as I can see the density is no specified anywhere in the paper.

andreww commented 3 years ago

Anyway, what I've done lives at https://github.com/murphyqm/pytesimal/tree/amw-examples. No point merging this.

andreww commented 3 years ago

Maybe the thing to do for now is work out if we could do this in a way that is backwards compatible (the thing I am unsure about is essentially coupling two discretisation's together) and find an easer to reproduce example for now. Bryson et al.?

murphyqm commented 3 years ago

One thing that could be useful/easier to implement would be a solid core - applying a zero flux bc at the centre - this would mean the code could model bodies with solid rocky cores (where k isn't high enough for isothermal to be a good approximation anymore) - it should just mean copying across mantle discretisation and squishing it into an object shaped box that interacts with the mantle in the same way (even if instantiation will look quite different)...

murphyqm commented 3 years ago

A possible example to reproduce - inward crystallizing core: https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2015JE004843

andreww commented 3 years ago

I think that should be essentially the same as the implementation in the Springwater preprint, so putting that in ought to do the job.

andreww commented 3 years ago

(We should probably also make some notes somewhere about internal heat sources - just so we know that could be implemented (in a way that does not break the APIs).