Simulating a crack in ice-substrate interface in 3D

Hello LAMMPS mailing list,

Please find the details of an issue am encountering below:

Objective: I am trying to use the 2D crack example as a template to simulate a delaminating crack in a 3D structure of Ice Ih-Aluminium(FCC111) at exactly the interface between ice and aluminium. This is to replicate an adhesion test procedure that is used in the laboratory.

Model: I have created an Ice Ih crystal and a FCC111 surface of Aluminium and combined them using ASE Python. I minimize and equilibrate it and then use “pair_style hybrid” to assign TIP4P and EAM potentials to ice and Aluminium respectively. I have used LJ potential to define cross-interaction between Aluminium and Oxygen atoms.

Issue: Since I use TIP4P and it requires a PPPM solver, I cannot use the boundary dimensions as “s s p” since the boundaries should be periodic for use with this solver. In the 2D crack example I see that the boundaries are defined as “s s p”. I’m not sure if this a mandatory requirement for simulating the crack properly.

Query: I was wondering if there was any way around this. One solution I’m aware of is that we can use coarse-grained modelling and assign mW potential to Ice and that would negate the use of the PPPM solver. But I’m wondering if there is a way it can be done by keeping the TIP4P potential intact for ice.

Any insight is appreciated. Thank you for your time.

Kind regards,

That is a pretty tall order and more importantly, I would say that your approach is not a good idea. What you want to model is quite different from what the crack example is about.

Using a hybrid pair style the way you are doing is not such a straightforward choice since you are looking at water, which is quite polar and thus has a significant polarizing impact on the metal (depending on distance) and that has quite an impact on the details of what you can learn from your simulation. I just wrote a more detailed discussion on hybrid potential models in a different topic. It should be applicable to your case, too: Using already existing potentials with GROMOS or OPLS simultaneously - #2 by akohlmey
In your case, however, the metal-water interaction is extremely important and I have serious doubts that a simple Lennard-Jones function is suitable to represent this with sufficient accuracy.

More importantly, you also need to keep in mind that your experiment will likely have very different length and time scales than what is accessible for modeling with MD simulations. So it is questionable, if you can extract useful knowledge from mimicking your experiment but then run it with much smaller dimensions and at a much shorter time.

So rather than making up a model, you should search the published literature for how people have previously extracted useful information from simulations for comparison to similar experiments and what level of agreement they have reached. Please keep in mind that there is a difference between a simulation (that represents real physics) and an animation (that looks like what you expect it should look like). Sometimes, the difference is small, sometimes less so.

Thank you for the insight and the detailed explanation on this, Dr. Axel.

I suspect this is possibly a misguided effort as you suggest since I have not seen any published literature in this area use crack propagation to delaminate ice from the substrate and extract any useful information. Usually its a different fix like smd, spring, deform, addforce, move etc.

I have previously used Steered Molecular Dynamics (SMD) with a coarse-grained ice-substrate model to obtain the peak delamination force and then convert it to adhesion strength for ice. However, as you state, the loading rates and timescales are not comparable to the real experiments and thus apart from some trends in adhesion strength, I have been unable to extract anything worthwhile to report.

I will definitely rethink this. Appreciate your response.