11institutetext: Brandeis University, Waltham MA 02453, USA
11email<EMAIL_ADDRESS>22institutetext: Swarthmore
College, Swarthmore PA 19081, USA
22email<EMAIL_ADDRESS>
# Time Majority Voting, a PC-based EEG Classifier for Non-expert Users

Guangyao Dou 11 0000-0001-8011-9658 Zheng Zhou 11 0000-0001-9313-2106 Xiaodong
Qu 1122 0000-0001-7610-6475

###### Abstract

Using Machine Learning and Deep Learning to predict cognitive tasks from
electroencephalography (EEG) signals is a rapidly advancing field in Brain-
Computer Interfaces (BCI). In contrast to the fields of computer vision and
natural language processing, the data amount of these trials is still rather
tiny. Developing a PC-based machine learning technique to increase the
participation of non-expert end-users could help solve this data collection
issue. We created a novel algorithm for machine learning called Time Majority
Voting (TMV). In our experiment, TMV performed better than cutting-edge
algorithms. It can operate efficiently on personal computers for
classification tasks involving the BCI. These interpretable data also assisted
end-users and researchers in comprehending EEG tests better.

###### Keywords:

Brain-Machine Interface Machine Learning Ensemble Methods Voting Time Series
Interpretable AI

##  1 Introduction

Researchers from Computer Science, Neuroscience, and Medical fields have
applied EEG-based Brain-Computer Interaction (BCI) techniques in many
different ways [24, 15, 34, 2, 26, 22, 19], such as diagnosis of abnormal
states, evaluating the effect of the treatments, seizure detection, motor
imagery tasks [17, 23, 27, 4, 6, 5], and developing BCI-based games [14].
Previous studies have demonstrated the great potential of machine learning,
deep learning, and transfer learning algorithms [25, 39, 42, 37, 29, 20, 40,
28, 8, 21, 18, 41, 38, 1, 3, 7, 16, 12] in such clinical and non-clinical data
analysis.

However, the data size of such experiments is still relatively small compared
to the areas of computer vision or natural language processing. Thus, some
deep learning or big data approaches still struggling with the limitation of
small dataset size. Also, EEG signals have noise issues, partly because of the
contact of sensors and skin for several current non-invasive consumer-grade
devices. The outlier issue is also a concern for the EEG data because of the
difficulties subjects have in concentrating on the experimental tasks during
the entire session. Current machine learning and deep learning algorithms are
more for clinical experiments and less for the possible experiment for non-
expert user to conduct at home.

Our research questions are: Can we develop a PC-based machine learning
algorithm for non-expert end-users to do EEG classification at home? Can we
achieve reasonably high accuracy while keeping the run time in an acceptable
range? Can we make the machine learning classification results explainable to
the end-users? To answer these questions, we proposed a new machine learning
classification algorithm, Time Majority Voting (TMV). We found TMV
outperformed other state-of-the-art classifiers. Also, its run time on a PC is
still acceptable compared to the deep learning algorithms. The classification
results are adequately interpretable to the end-users.

The paper is organized as follows: section two discusses several most
frequently used classification algorithms for BCI research, then present our
new algorithm. Section three presents our experiment conducted to test the new
algorithm. Section four elaborates our result followed by sedition five, which
discusses the limitation and future work. Lastly, section six concludes the
study and summarized our answers to the research questions.

##  2 Algorithms

All of the code was run on a 2018 Macbook Pro with a 2.2GHz 6-core Intel Core
i7 processor and with 16 GB of memory. The Python version is 3.8. The scikit-
learn [31] version is 0.24.1. The PyTorch [30] version is 1.10. The code
discussed in this paper is available online
(https://github.com/GuangyaoDou/Time_Majority_Voting).

###  2.1 Existing Algorithms

u We reviewed and implemented several machine learning algorithms commonly
used in the field [9, 10, 11, 13]. For examples, Linear Classifiers, Nearest
Neighbors, Decision Trees, and Ensemble Methods.

Linear Classifiers: The Shrinkage Linear Discriminant Analysis (Shrinkage LDA)
performed adequately on EEG datasets with simple tasks. The Support Vector
Machine (SVM), effective in high dimensional spaces, performed reasonably well
based on the previous research. These algorithms are simple to implement and
are computationally efficient.

Nearest Neighbor: Such a classifier implements the K-Nearest Neighbor(KNN).
KNN performs voting to determine an unseen dataset to one of the k nearest
neighbors. The KNN performed pretty well compared to most other classifiers on
the EEG dataset.

Decision Tree: The Decision Trees classifier is easy to understand, implement,
and interpret. The decision tree creates a model that predicts the outcome of
a data point based on decision rules. The computational cost is low and can
handle both numerical and categorical data. However, it might overfit as trees
are too complex.

Ensemble Methods: We used Random Forest and boosting. These are Ensemble
Machine Learning algorithms that combine the predictions of several weaker
learners and form more robust and more accurate predictions. These have been
widely used in EEG-based experiments, and research [24, 33, 35].

Deep Learning: We implemented CNN with ReLu, and RNN, especially LSTM, mainly
using toolsets from the PyTorch platform. [37, 15, 34]. These DLs performed
very well on EEG datasets. However, we excluded these algorithms in this paper
due to their high runtime.

###  2.2 Our New Algorithms

This paper proposed a new voting approach based on the top two individual
classifiers. Ensemble methods, especially boosting, bagging, and voting, have
demonstrated excellent performance in previous research. [24, 37, 15, 36, 34]
In EEG-based BCI classification research, the following voting methods have
been investigated in several experiments [24, 33]: majority voting, weighted
voting, and time continuity voting. Here we considered the advantages of both
majority voting and time continuity voting and developed our new Time Majority
Voting (TMV) algorithm.

Figure 1: Our New Algorithm: Time Majority Voting

Time Majority Voting (TMV).
Figure 1 demonstrated the concept of the new algorithm. More details are in
the Experiment section and the Result section. There are two phases in the new
algorithm. First, we investigated the state-of-the-art machine learning
algorithms [24, 37, 15, 34] and found the top two performers on average. In
our experiments during phase 1, We tested Random Forest (RF), RBF and linear
SVM, kNN, Decision Tree, and several boosting algorithms. We found Random
Forest performed the best, and RBF SVM performed the second on average. Next,
we entered phase two. For each subject, we picked the majority task predicted
by the best performing classifier, the Random Forest classifier, for each time
interval of each task from each session.

The next step is voting. We used the Random Forest and the RBF SVM to conduct
the voting process. The best algorithm is the Random Forest, and the second
algorithm is RBF SVM. The voting details are shown in the four examples at the
bottom of Figure 1. If both the Random Forest and the RBF SVM algorithms agree
with the results, as shown in the first two rows, the results reflect both
algorithms’ results. For example, as the second row shows, if both the
classifiers predicted the task 2, then the Time Majority Voting will yield
task 2 as the result. If the two algorithms do not agree on the prediction
results, the results will be labeled as the majority of tasks already
determined by the Random Forest classifier, as shown in the last two rows. No
matter what the two algorithms predict out of the five tasks, even if they
predicted it as two different tasks other than the majority task, as the last
row shows, Task 2 and Task 3, the TMV result is still set to be Task one. This
concept is based on the temporal dependence time-series features in the
previous research [25, 36].

##  3 Experiments

Several EEG experiments focus on the high-level cognitive tasks that college
students frequently conduct, as mentioned in Table 1. In this experiment, we
used the dataset from the Think-Count-Recall (TCR) paper [33]. Scalp-EEG
signals were recorded from seventeen subjects. Each one was tested in six
sessions. Each session is five minutes long, with five tasks, each task is one
minute. Tasks were selected by the subjects together with the researchers
based on frequent tasks in study environments for students in their everyday
lives. Each subject completed six sessions over several weeks. The five tasks
are Think(T), Count(C), Recall(R), Breathe(B), and Draw(D).

(E) T | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
[32] | Math | Close-eye Relax | Read | Open-eye Relax | None
[36] | Python Passive | Math Passive | Python Active | Math Active | None
[35] | Read | Write Copy | Write Answer | Type Copy | Type Answer
[33] | Think | Count | Recall | Breathe | Draw
Table 1:  Tasks (T) in Experiments (E)

###  3.1 Data Preprocess

Data cleaning: As mentioned in [35, 33], for each task during each session,
the first 30% of the data, which is the first 18 seconds of each task and
during the transition phase, will be removed. Thus, each one-minute task only
had 42 seconds left. This has been proved reasonable during the data cleaning
phase. Some electrodes may have temporarily lost contact with the subjects’
scalp during the EEG recording. The result was that multiple sequential
spectral snapshots from one or more electrodes had the same value. In this
paper, we decide to remove such anomaly when detected for a consecutive 1.4
seconds. Such a cleaning action caused a different level of loss of the data
for each subject.

Subjects: We had a total of seventeen subjects. After the data cleaning
actions, subjects who lost more than 65% of the total data will be excluded
from the subsequent analysis. In the end, there were twelve subjects left to
continue the analysis. Moreover, for the six sessions of each subject, if a
session lost more than 65% of the data, then that session will also be
excluded for further analysis.

Time-wise Cross-Validation: We adopted time-wise cross-validation. We divided
each task into seven subsets, meaning each subset had six seconds, evenly and
continuously. Then, we created a total of seven folds. Each fold contained six
seconds of data for each task for each session. We checked any folds that lost
more than 65% of the original data in each fold and discarded these folds for
future analysis. Next, we used one fold for testing and the remaining non-
discarded folds for training, and we cross-validated them.

##  4 Results

Algorithms | Average Accuracy | Average code run-time (s)
---|---|---
Random Forest Phase 1 | 0.55 | 42.0
RBF SVM | 0.53 | 30.5
Nearest Neighbors | 0.48 | 1.9
Decision Tree | 0.44 | 0.9
Linear SVM | 0.42 | 23.2
Shrinkage LDA | 0.42 | 0.2
Adaboost Classifier | 0.39 | 47.8
RUSBoost | 0.39 | 28.2
GradientBoost | 0.31 | 24.0
Table 2:  State-of-the-art Algorithms with Accuracy and Run-time

###  4.1 Existing Algorithms

We reported the average accuracy for all subjects and the runtime of each
classifier for each subject we trained and tested during phase 1 in table 2.
As we can see, Random Forest had the highest accuracy of 0.55 in our
experiment, and the SVM with RBF kernel performed adequately on the TCR
dataset with an average accuracy of 0.53. Though Nearest Neighbors did not
perform as well as the Random Forest and the RBF SVM, it was one of the
fastest algorithms on personal computers. Other ensemble methods such as
Adaboost, RusBoost, and GradientBoosting performed relatively lower than these
top three algorithms.

Figure 2: Accuracy for Different Algorithms

The individual difference may impact the accuracy of each subject. But we can
still recognize a general pattern from Fig 2. We ordered all twelve subjects
by prediction accuracy using Random Forest. Most of the algorithms
demonstrated consistent patterns for the different algorithms. Random Forest
and RBF SVM are above most of the other algorithms. We kept the threshold of
maintaining the subject to 35% of the remaining data, as we believed that when
a subject has little data left, the high accuracy from that subject
contributes little to our research.

Figure 3: subject 3’s task 1’s RF and SVM in Phase 1

###  4.2 Identify Noisy Sessions

Figure 3 shows what the Random Forest and the RBF SVM in phase 1 predicted
during the 42 seconds of subject 3’s task 1 for all six sessions. As we can
see, both the Random Forest and the RBF SVM mainly produced results with
relatively high accuracy in sessions one, four, and six. On the other hand,
session two and three had much noise, and session five predicted task 3 for
the majority of the time. To minimize the impact of noisy datasets, we
calculated the accuracy compared to the ground truth based on the output of
Random Forest for each session for each task in phase 1.

We excluded any sessions that yielded an accuracy of less than 50%. In the
case of Figure 3, we excluded sessions two, three, and five for further
machine learning analysis. We reported to subject three and started discussing
what may have happened in these sessions. With the exclusion of noisy
sessions, the new accuracy for the Random Forest in phase 2 is referred to as
"Random Forests Phase 2" later in the paper. We also performed Time Majority
Voting on this cleaner dataset.

###  4.3 Time Majority Voting

S/T | t1 RF | t1 T | t2 RF | t2 T | t3 RF | t3 T | t4 RF | t4 T | t5 RF | t5
---|---|---|---|---|---|---|---|---|---|---
s1 | 0.686 | 0.745 | 0.669 | 0.683 | -1 | -1 | 0.569 | 0.781 | 0.798 | 0.898
s2 | -1 | -1 | -1 | -1 | -1 | -1 | 0.693 | 0.838 | 0.633 | 0.707
s3 | -1 | -1 | -1 | -1 | -1 | -1 | 0.731 | 0.824 | 0.938 | 0.969
s4 | 0.543 | 0.59 | 0.743 | 0.869 | 0.74 | 0.788 | 0.521 | 0.671 | -1 | -1
s5 | -1 | -1 | 0.645 | 0.738 | 0.588 | 0.8 | 0.593 | 0.743 | 0.662 | 0.829
s6 | 0.762 | 0.871 | -1 | -1 | 0.507 | 0.714 | 0.812 | 0.917 | 0.681 | 0.845
Average | 0.663 | 0.736 | 0.686 | 0.763 | 0.612 | 0.767 | 0.653 | 0.796 | 0.742 | 0.85
Table 3:  TCR, Accuracy of Top 1 (RF) Phase 2, and TMV(T) (subject 3, by
Session(S)/Task(T))

As shown in Table 3, the Time Majority Voting (TMV) has achieved a higher
accuracy for subject three, all the six sessions. A value of -1 means that we
excluded that session for that task, as we discussed in the previous section.
The Random Forest classifier also reached a higher accuracy after cleaning the
noisy sessions. As Figure 4 and Figure 5 shows, the pattern is consistent
across all the subjects. Figure 5 also shows a clear pattern that not only the
Random Forest but also the RBF SVM also increased accuracy in phase 2 across
all subjects. Table 4 shows the TMV achieved an 80% average accuracy with an
average 74.3 seconds run time. The runtime consists of 39.1 seconds of running
the Random Forest and 29.5 seconds of running the RBF SVM. The time for the
actual voting is, on average, 5.7 seconds. The training process is the most
time-consuming part of this analysis.

Figure 4: TMV and RF Phase 2 Algorithms | Average Accuracy | Average code run-time (s)
---|---|---
Time Majority Voting | 0.80 | 39.1 + 29.5 + 5.7 = 74.3
Random Forest Phase 2 | 0.7 | 39.1
RBF SVM Phase 2 | 0.66 | 29.5
Random Forest Phase 1 | 0.55 | 39.1
RBF SVM Phase 1 | 0.53 | 29.5
Table 4:  TMV and Random Forest with Accuracy and Run-time Figure 5: TMV, RF
Phase 2, RBF SVM Phase 2, RF phase 1, and RBF SVM Phase 1 Figure 6: subject
3’s task 1’s TMV

Figure 6 shows what the Time Majority Voting(TMV) predicted during the 42
seconds of subject 3’s task 1 for all six sessions. As you can see, sessions
2, 3, and 5 have values of -1, which means that they have been excluded.
Sessions 1, 4, and 6 have less noise than sessions 2, 3, and 6, and the
accuracy is relatively high.

##  5 Discussion

###  5.1 Accuracy and Data Remain

The innovation of this method is mainly about temporal dependency. As [24, 37,
15, 36, 34] mentioned, EEG data has a significant temporal dependency. The
signal of the same task takes about 12 to 18 seconds to switch to the next
task. Using majority voting can catch this type of time continuity effect. If
both classifiers recognize the same pattern, it is more likely to assure the
results. If both of the classifiers recognize the same pattern that is
different from the majority result, it is possible that the participants were
doing other tasks during the data collection. If only one classifier detects
some unusual behaviors, we label it as the majority task of the session. Thus
we highlight the noise and keep the remaining data to reflect more on the time
continuity nature of the EEG signals.

Figure 7: Data Remained for Training and Run-time

###  5.2 Runtime and Training Data

As Figure 7 shows, the run time directly correlated to the training data size.
After the data pre-processing, we cleaned up the noise with a plateau longer
than a threshold, as mentioned in [24, 32, 34]. We first identified more noisy
sessions during the Time Majority Voting process based on the top two
classifiers. During this step, more sessions were excluded for further
analysis. The training process was the most time-consuming step during the
coding running process. Thus the runtime changed together with the size of the
training data.

###  5.3 Interpretability

Figure 3 and Figure 6 shows two examples of the feedback results we present to
the end-users. In this figure 3, the six sessions of subject three, task one
is listed as six horizontal charts. The first chart represents session one.
Both the Random Forest and the RBF SVM identified the majority task as task
one. And for the results of the different predictions, both of the classifiers
agreed at some time spots, but not all of them. Our Time Majority Voting
algorithm favorite the majority voting results.

We started with the sessions with good prediction results when demonstrating
these figures to each subject. For example, in this figure 3, sessions one,
four, and six show pretty consistent patterns. The majority of task prediction
results were the designed task one. That implies that the subject may spend
more time on task one as planned during these sessions. Session two and three
had a lot of different prediction results from both classifiers. Thus we
suspected some unexpected reason might cause this situation. We referred back
to the experiment notes and reached out to subject three. After discussing
with the end-user, we figured out that he had many issues with the sensor
signals and was adjusting the EEG headset most of the time during the sessions
2 and 3. Thus we had more information to exclude these sessions from further
data analysis. Session five is another situation. The data implied that the
subject was doing task three, but the experiment notes were missing for that
session, and the subject did not remember the details about that session.
Thus, we left a question mark for that session, excluded the session for now,
and came up with an improvement plan to keep better experiment notes. This
type of machine learning result is explainable to the end-user.

Such interpretability could contribute to a better understanding of the
results and better design for future experiments.

##  6 Conclusion

This paper investigated the state-of-the-art machine learning algorithms that
can run on mainstream personal computers for EEG-based BCI. We then proposed a
new algorithm, Time Majority Voting (TMV). The results demonstrated that TMV
outperformed other existing classifiers. The run time for TMV is still within
the acceptable range on a PC. The interpretability of TMV can contribute to a
better understanding of the machine learning analysis and an improved design
for future experiments.

## References

  * [1] An, S., Ogras, U.Y.: Mars: mmwave-based assistive rehabilitation system for smart healthcare. ACM Transactions on Embedded Computing Systems (TECS) 20(5s), 1–22 (2021)
  * [2] Appriou, A., Cichocki, A., Lotte, F.: Modern machine-learning algorithms: for classifying cognitive and affective states from electroencephalography signals. IEEE Systems, Man, and Cybernetics Magazine 6(3), 29–38 (2020)
  * [3] Basaklar, T., Tuncel, Y., An, S., Ogras, U.: Wearable devices and low-power design for smart health applications: challenges and opportunities. In: 2021 IEEE/ACM International Symposium on Low Power Electronics and Design (ISLPED). pp. 1–1. IEEE (2021)
  * [4] Bashivan, P., Bidelman, G.M., Yeasin, M.: Spectrotemporal dynamics of the EEG during working memory encoding and maintenance predicts individual behavioral capacity. European Journal of Neuroscience 40(12), 3774–3784 (2014)
  * [5] Bashivan, P., Rish, I., Heisig, S.: Mental state recognition via wearable eeg. arXiv preprint arXiv:1602.00985 (2016)
  * [6] Bashivan, P., Rish, I., Yeasin, M., Codella, N.: Learning representations from EEG with deep recurrent-convolutional neural networks. arXiv preprint arXiv:1511.06448 (2015)
  * [7] Bhat, G., Tuncel, Y., An, S., Ogras, U.Y.: Wearable IoT Devices for Health Monitoring. In: TechConnect Briefs 2019. pp. 357–360 (2019)
  * [8] Bird, J.J., Manso, L.J., Ribeiro, E.P., Ekart, A., Faria, D.R.: A study on mental state classification using eeg-based brain-machine interface. In: 2018 International Conference on Intelligent Systems (IS). pp. 795–800. IEEE (2018)
  * [9] Breiman, L.: Bagging predictors. Machine learning 24(2), 123–140 (1996)
  * [10] Breiman, L.: Random forests. Machine learning 45(1), 5–32 (2001)
  * [11] Breiman, L.: Classification and regression trees. Routledge (2017)
  * [12] Chen, L., Lu, Y., Wu, C.T., Clarke, R., Yu, G., Van Eyk, J.E., Herrington, D.M., Wang, Y.: Data-driven detection of subtype-specific differentially expressed genes. Scientific reports 11(1), 1–12 (2021)
  * [13] Chevalier, J.A., Gramfort, A., Salmon, J., Thirion, B.: Statistical control for spatio-temporal meg/eeg source imaging with desparsified multi-task lasso. arXiv preprint arXiv:2009.14310 (2020)
  * [14] Coyle, D., Principe, J., Lotte, F., Nijholt, A.: Guest editorial: Brain/neuronal-computer game interfaces and interaction. IEEE Transactions on Computational Intelligence and AI in games 5(2), 77–81 (2013)
  * [15] Craik, A., He, Y., Contreras-Vidal, J.L.: Deep learning for electroencephalogram (eeg) classification tasks: a review. Journal of neural engineering 16(3), 031001 (2019)
  * [16] Derby, J.J., Zhang, C., Seebeck, J., Peterson, J.H., Tremsin, A.S., Perrodin, D., Bizarri, G.A., Bourret, E.D., Losko, A.S., Vogel, S.C.: Computational modeling and neutron imaging to understand interface shape and solute segregation during the vertical gradient freeze growth of babrcl: Eu. Journal of Crystal Growth 536, 125572 (2020)
  * [17] Devlaminck, D., Waegeman, W., Bauwens, B., Wyns, B., Santens, P., Otte, G.: From circular ordinal regression to multilabel classification. In: Proceedings of the 2010 Workshop on Preference Learning (European Conference on Machine Learning, ECML). p. 15 (2010)
  * [18] Gu, J., Zhao, Z., Zeng, Z., Wang, Y., Qiu, Z., Veeravalli, B., Goh, B.K.P., Bonney, G.K., Madhavan, K., Ying, C.W., et al.: Multi-phase cross-modal learning for noninvasive gene mutation prediction in hepatocellular carcinoma. In: 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). pp. 5814–5817. IEEE (2020)
  * [19] Kastrati, A., Płomecka, M.M.B., Pascual, D., Wolf, L., Gillioz, V., Wattenhofer, R., Langer, N.: Eegeyenet: a simultaneous electroencephalography and eye-tracking dataset and benchmark for eye movement prediction. arXiv preprint arXiv:2111.05100 (2021)
  * [20] Kaya, M., Binli, M.K., Ozbay, E., Yanar, H., Mishchenko, Y.: A large electroencephalographic motor imagery dataset for electroencephalographic brain computer interfaces. Scientific data 5(1), 1–16 (2018)
  * [21] Li, S., Zhao, Z., Xu, K., Zeng, Z., Guan, C.: Hierarchical consistency regularized mean teacher for semi-supervised 3d left atrium segmentation. In: 2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). pp. 3395–3398. IEEE (2021)
  * [22] Lotte, F.: A tutorial on EEG signal-processing techniques for mental-state recognition in brain–computer interfaces. In: Guide to Brain-Computer Music Interfacing, pp. 133–161. Springer (2014)
  * [23] Lotte, F.: Signal processing approaches to minimize or suppress calibration time in oscillatory activity-based brain–computer interfaces. Proceedings of the IEEE 103(6), 871–890 (2015)
  * [24] Lotte, F., Bougrain, L., Cichocki, A., Clerc, M., Congedo, M., Rakotomamonjy, A., Yger, F.: A review of classification algorithms for EEG-based brain–computer interfaces: a 10 year update. Journal of neural engineering 15(3), 031005 (2018)
  * [25] Lotte, F., Congedo, M., Lécuyer, A., Lamarche, F., Arnaldi, B.: A review of classification algorithms for EEG-based brain–computer interfaces. Journal of neural engineering 4(2), R1 (2007)
  * [26] Lotte, F., Guan, C.: Regularizing common spatial patterns to improve bci designs: unified theory and new algorithms. IEEE Transactions on biomedical Engineering 58(2), 355–362 (2010)
  * [27] Lotte, F., Jeunet, C.: Towards improved bci based on human learning principles. In: The 3rd International Winter Conference on Brain-Computer Interface. pp. 1–4. IEEE (2015)
  * [28] Lotte, F., Jeunet, C., Mladenović, J., N’Kaoua, B., Pillette, L.: A bci challenge for the signal processing community: considering the user in the loop (2018)
  * [29] Miller, K.J.: A library of human electrocorticographic data and analyses. Nature human behaviour 3(11), 1225–1235 (2019)
  * [30] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A.: Automatic differentiation in pytorch. NeurIPS (2017)
  * [31] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E.: Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825–2830 (2011)
  * [32] Qu, X., Hall, M., Sun, Y., Sekuler, R., Hickey, T.J.: A personalized reading coach using wearable EEG sensors-a pilot study of brainwave learning analytics. In: CSEDU (2). pp. 501–507 (2018)
  * [33] Qu, X., Liu, P., Li, Z., Hickey, T.: Multi-class time continuity voting for eeg classification. In: International Conference on Brain Function Assessment in Learning. pp. 24–33. Springer (2020)
  * [34] Qu, X., Liukasemsarn, S., Tu, J., Higgins, A., Hickey, T.J., Hall, M.H.: Identifying clinically and functionally distinct groups among healthy controls and first episode psychosis patients by clustering on eeg patterns. Frontiers in psychiatry p. 938 (2020)
  * [35] Qu, X., Mei, Q., Liu, P., Hickey, T.: Using eeg to distinguish between writing and typing for the same cognitive task. In: International Conference on Brain Function Assessment in Learning. pp. 66–74. Springer (2020)
  * [36] Qu, X., Sun, Y., Sekuler, R., Hickey, T.: EEG markers of stem learning. In: 2018 IEEE Frontiers in Education Conference (FIE). pp. 1–9. IEEE (2018)
  * [37] Roy, Y., Banville, H., Albuquerque, I., Gramfort, A., Falk, T.H., Faubert, J.: Deep learning-based electroencephalography analysis: a systematic review. Journal of neural engineering 16(5), 051001 (2019)
  * [38] Zeng, Z., Zhao, W., Qian, P., Zhou, Y., Zhao, Z., Chen, C., Guan, C.: Robust traffic prediction from spatial-temporal data based on conditional distribution learning. IEEE Transactions on Cybernetics (2021)
  * [39] Zhang, X., Yao, L., Wang, X., Monaghan, J.J., Mcalpine, D., Zhang, Y.: A survey on deep learning-based non-invasive brain signals: recent advances and new frontiers. Journal of Neural Engineering (2020)
  * [40] Zhao, Z., Chopra, K., Zeng, Z., Li, X.: Sea-net: Squeeze-and-excitation attention net for diabetic retinopathy grading. In: 2020 IEEE International Conference on Image Processing (ICIP). pp. 2496–2500. IEEE (2020)
  * [41] Zhao, Z., Qian, P., Hou, Y., Zeng, Z.: Adaptive mean-residue loss for robust facial age estimation. In: 2022 IEEE International Conference on Multimedia and Expo (ICME). IEEE (2022)
  * [42] Zhao, Z., Zhang, K., Hao, X., Tian, J., Chua, M.C.H., Chen, L., Xu, X.: Bira-net: Bilinear attention net for diabetic retinopathy grading. In: 2019 IEEE International Conference on Image Processing (ICIP). pp. 1385–1389. IEEE (2019)


Meshfree One-Fluid Modelling of Liquid-Vapor Phase Transitions

Pratik Suchde1,2,
Heinrich Kraus1,3,
Benjamin Bock-Marbach 1,
Jörg Kuhnert1

1Fraunhofer ITWM, Fraunhofer-Platz 1, 67663 Kaiserslautern, Germany.
2University of Luxembourg, 2 avenue de l'université, L-4365 Esch-sur-alzette, Luxembourg
3 Institut für Mathematik, Universität Kassel, Heinrich-Plett-Straße 40, 34132 Kassel, Germany

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
We introduce a meshfree collocation framework to model the phase change from liquid to vapor at or above the boiling point. While typical vaporization or boiling simulations focus on the vaporization from the bulk of the fluid, here we include the possibility of vaporization from the free surface, when a moving fluid comes into contact with a superheated surface. We present a continuum, one-fluid approach in which the liquid and vapor phases are modeled with the same constitutive equations, with different material properties. The novelty here is a monolithic approach without explicit modeling of the interface between the phases, neither in a sharp nor diffuse sense. Furthermore, no interface boundary conditions or source terms are needed between the liquid and vapor phases. Instead, the phase transition is modeled only using material properties varying with temperature. Towards this end, we also present an enrichment of strong form meshfree generalized finite difference methods (GFDM) to accurately capture derivatives in the presence of jumps in density, viscosity, and other physical properties. The numerical results show a good agreement with experimental results, and highlight the ability of our proposed framework to model phase changes with large jumps.

§ INTRODUCTION

The process of phase transition or phase change occurs everywhere around us, both in nature and in industrial processes. They have been widely studied both theoretically and experimentally. While computer simulations of phase transitions have become increasingly popular in the last few decades, there remain many challenges in modeling and simulating these processes accurately.

In the present work, we restrict our discussion to so-called first-order phase transitions that involve latent heat. Several different approaches have been considered for simulating such phase change processes. On one hand, microscopic simulations have been widely considered [7], which simulate the process at a molecular or atomic level. These methods can be classified into one of two approaches: (i) molecular dynamics (MD) simulations which integrate equations of motions on interacting particles, see, for example, [9, 25, 29, 59], and (ii) Monte Carlo (MC) methods which rely on stochastic processes,
see, for example, [2, 35, 55, 58]. While these microscopic phase change models can be very physically accurate, they are only suitable for simulations over very small length scales, typically smaller than millimeters.

On the other hand, continuum approaches are used to model the phase change process on a macroscopic level. These models are more appropriate for large length scale simulations of phase transitions. The primary difficulties in these models arise due to the volume expansion or contraction, and the evolving phase boundaries that undergo large deformations and possibly even topological changes. Based on similar notions in classical multiphase flow problems [21], continuum phase change simulations are usually classified into two categories based on how the interface between the phases is treated: (i) The classical approach is a sharp interface method, which assumes a clear and distinct interface of zero thickness between the two phases. The phase transition is modeled as interface boundary conditions for heat and mass transfer, see for example, [4, 13, 57]. (ii) An alternative approach that has become very popular is the diffusive interface approach, that postulates an interface of a small non-zero thickness. Here, the interface conditions are modeled as source terms in the governing equations across the entire thickness of the diffuse interface, see for example, [14, 30, 33, 46]. Mesoscopic approaches for phase change modeling have also been proposed, and typically use the Lattice Boltzmann method (LBM) [22].

However, for many practical applications of phase transitions, none of the above methods are suitable. The motivation behind the present work is the vaporization of cutting fluid during metalworking processes.
In these processes, chips are removed from a metal workpiece, using a cutting tool, to obtain the desired shape of the workpiece. These processes generate a large amount of heat due to friction and metal deformation [51]. Thus, a cutting fluid, also referred to as coolant or lubricant, is used to rapidly remove heat without needing to pause the production process, while also lubricating and preventing unwanted thermal expansion [49, 50]. The temperatures at the cutting region can often be much higher than the boiling temperature of the cutting fluid. As a result, as the fluid comes into contact with the workpiece or the tool in the cutting zone, it instantly vaporizes. Thus, to accurately simulate such cutting processes, the liquid-vapor transition must be modeled. However, existing phase change models are unable to represent this complex process. Firstly, microscopic models are infeasible for such large scale simulations. A macroscopic sharp interface method is not applicable due to the complicated 3D interface with changing topology between the liquid and vapor phases. Moreover, unlike typical macroscopic phase change simulations where the only free boundary is the interface between the phases, here both the liquid and the gas themselves have their own external free surfaces too. In fact, rapid boiling of the liquid occurs at and near a free surface. Diffuse interface methods would fail near the region where multiple interfaces intersect: the liquid-vapor interface, external free surfaces for both liquid and vapor phases, and solid-liquid and solid-vapor interfaces where the cutting fluid hits the workpiece and tool.

Meshfree methods have become a popular alternative to conventional mesh-based simulation methods, especially for fluid flow applications. Their advantages become the most relevant for applications with complex geometries, moving interfaces, and free boundaries, each of which holds true for our vaporization of cutting fluid application. Another advantage of a meshfree method is the ease of incorporating moving Lagrangian frameworks. While Lagrangian frameworks are considered to be very accurate in representing moving boundaries and interfaces, their use with mesh-based methods is limited to applications with small deformations [21]. This is primarily because mesh deformation as a result of moving the mesh requires an expensive remeshing step. In contrast, in meshfree methods point cloud deformation can be fixed locally and is thus much cheaper. As a result, meshfree Lagrangian methods are especially suited for applications with large deformations and topology changes of interfaces. We thus choose a meshfree approach to model the application at hand.

In this paper, we introduce a novel meshfree Lagrangian framework to simulate phase change processes. Motivated by similar ideas in cavitation modeling [15], we use a one-fluid or homogeneous approach, in which both phases of the phase transition process are represented by a single fluid, with one set of conservation equations and varying material properties representing the different phases. We note that the so-called front tracking method for phase change [18] also uses a one-fluid or a “one-field" approach. However, they explicitly track the interface and use two different grids for the flow equations and the interface. In contrast, here, we do not track the interface at all, and thus only one grid (more appropriately, one point cloud) is needed. Mass and energy transfer between the phases are represented directly by the particle-based Lagrangian nature of the discretization. This has the advantage of not needing to explicitly track or capture the interface, nor needing to determine and impose any interfacial conditions, which could have been cumbersome for interfaces with very complex topology. Using such an approach in a fully meshfree setting raises several challenges, which are highlighted and addressed in the present work. We introduce modifications to collocation-based meshfree derivative computation to enable accurate capturing of the sudden changes in the material properties. Furthermore, we introduce a framework for modeling the varying material properties that helps in attaining numerical stability while maintaining physical consistency.

Given the aforementioned motivation, we focus on liquid-to-vapor phase change throughout the present work. However, the methods introduced here could be directly used to model any first-order phase transition.

This paper is organized as follows. We first introduce the constitutive equations underlying the one-fluid model in Section <ref>, followed by explaining how the phase change is modeled by varying material properties in Section <ref>. The meshfree discretization procedure and enrichment process are then explained in Section <ref>. Section <ref> describes some details of the numerical scheme, and a short comparison with other phase change models is presented in Section <ref>. Numerical applications of the method are explained in Section <ref>, followed by brief concluding remarks in Section <ref>.

§ CONSTITUTIVE EQUATIONS

The underlying constitutive equations used to model both the liquid and the vapor phase as a single one-fluid model are the standard conservation equations of mass, momentum, and energy, written in a Lagrangian formulation

\begin{align}
\frac{D \rho}{Dt} &= -\rho \nabla \cdot \vec{v}, \label{Eq:CoMass}\\
\frac{D  \vec{v}  }{Dt} &= \frac{1}{\rho} \nabla \cdot \mathbf{S} - \frac{1}{\rho} \nabla p
+  \vec{g}, \label{Eq:CoMomentum}\\
\rho \frac{D  E }{Dt} &= \nabla \cdot (\mathbf{S}\vec{v})
- \nabla \cdot ( p \vec{v} ) +
\rho \vec{g} \cdot \vec{v} +
\nabla \cdot \left( \lambda_{\text{eff}} \nabla T \right) + q, \label{Eq:CoEnergy}
\end{align}

where the energy is given by $E = c_v T + \frac{1}{2}\vec{v} \cdot \vec{v}$. Furthermore, $\rho$ is the density, $\vec{v}$ is the velocity, $p$ is the pressure, $T$ is the temperature, while $t$ is the time, $\vec{g}$ is composed of both gravity and body forces, $c_v$ is the specific heat capacity, $\lambda_{\text{eff}}$ is the effective heat conductivity (see Section <ref>), and $q$ models the heat sources. The material or Lagrangian derivative is denoted by $\frac{D}{Dt}$. The stress tensor $\mathbf{S}$ can be composed of both viscous and non-viscous terms

\begin{equation}
\mathbf{S}  = \mathbf{S}(\vec{v})
= \mathbf{S}_{\text{visc}}(\vec{v}) + \mathbf{S}_{\text{solid}}(\vec{v}).
\end{equation}

The viscous stress tensor is given by

\begin{equation}
\mathbf{S}_{\text{visc}}(\vec{v})
= \eta_{\text{eff}} \left(
\left( \nabla \vec{v} \right) + \left( \nabla \vec{v} \right)^T  -
\frac{2}{3} \nabla \cdot \vec{v} \, \mathbf{I} \right),
\end{equation}

where $\eta_{\text{eff}}$ denotes the effective viscosity  (see Section <ref>). The non-viscous parts of the stress tensor are grouped into the $\mathbf{S}_{\text{solid}}$ term governed by the appropriate material model. Since we are focusing on Newtonian liquid-to-vapor phase transition, we have $\mathbf{S}_{\text{solid}}=\textbf{0}$ in the applications considered.

The governing equations (<ref>) require a closure relation, which is given, depending on the application, either by an equation of state or by discrete curves determined by experiments. In each case, we have a relation $\rho = \rho(T, p)$, with the temperature and pressure both being functions of space $\vec{x}$ and time $t$. In many boiling simulations in literature, an assumption of the incompressibility of the vapor phase is often imposed (example, [10, 33, 57]), to focus on the modelling of the phase transition process itself. Here, we allow the possibility of both incompressible and compressible vapor phases by using different closure relations. For actual applications, the choice of the appropriate closure model is a big challenge, as is also the case for phase change models which track the interface.

The phase change from liquid to vapor is inherently a turbulent process. We thus incorporate a standard $k$-$\varepsilon$ turbulence model [20, 53] with fluctuating dilatation and source terms omitted. For more details on the turbulence model used, the parameters, and the numerical interpretation and integration of the turbulence model, we refer to our earlier work [26].

§ MATERIAL PROPERTIES

Within the one-fluid model used here (see Section <ref>), the different phases are only distinguished by different material properties. This could lead to large jumps, or at minimum large gradients, in material properties where the phase transition occurs. We now introduce the numerical handling of these changing material properties such that both numerical stability and physical consistency are ensured.

Here, the term interface refers to the interface in the physical event of phase transition from liquid to vapor. We emphasize that this interface is not tracked explicitly in the numerical setting, neither in the sharp sense nor in a diffuse sense. Similarly, “two phases" refers only to the physical liquid and vapor phases, while in the numerical setting, they are treated as a single phase computationally.

§.§ Boiling point

An important question that arises in the one-fluid modeling of phase change concerns when the latent heat of the phase change must be taken into account. For complete physical accuracy, this should occur at a constant temperature. This would entail the latent heat of vaporization being modeled only at the normal boiling temperature or boiling point of the liquid. However, numerically, this can prove to be quite challenging, as it could mean reducing the simulation time step when a liquid nears the boiling point. To overcome this, we assume the phase change occurs over a small temperature range rather than modeling it at a single temperature.

For a given pressure, the boiling point of the liquid phase is represented by $T_p^*$. Numerically, we spread this boiling point across a small interval of size $2\Delta T^*$. Thus, the vaporization process in our model occurs in the boiling temperature range $[ T_p^* - \Delta T^*, T_p^* + \Delta T^*]$. The incorporation of the latent heat of vaporization is spread over this interval (more details in the next subsection). This provides a straightforward approach to modeling the phase change process without reducing the time step as the temperature approaches the phase change temperature. The use of a boiling temperature range here is similar in essence to the use of “mushy zones" in solidification modelling [54].

§.§ Effective specific heat capacity

In the boiling temperature range, specifically, in the semi-open intervals $[ T_p^* - \Delta T^*, T_p^* )$ and $(T_p^*,  T_p^* + \Delta T^*]$, our model must account for both the heat capacity of the material, as well as the latent heat of vaporization. For this, we introduce a numerical notion of heat capacity of vaporization, referred to by $c_v^{\text{vaporization}}$. This will combine both the specific latent heat for the phase change process and the specific heat capacity to raise the temperature without phase change, into a single entity. As a result, the isochoric structure of the specific heat capacity is taken as

\begin{equation}
c_v = \begin{cases}
c_v^{\text{liquid}},       & \text{if }   T - T_p^*   < -\Delta T^*, \\
c_v^{\text{vaporization}}, & \text{if } | T - T_p^* | <  \Delta T^*, \\
c_v^{\text{vapor}},       & \text{if }   T - T_p^*   >  \Delta T^*,
\end{cases}
\end{equation}

where $c_v^{\text{liquid}}$ and $c_v^{\text{vapor}}$ are the specific heat capacities of the liquid and vapor phases, respectively. Both of these can be either constant or given by temperature-dependent curves.

The value of $c_v^{\text{vaporization}}$ is given by ensuring physical consistency in the sense that the amount of heat required to take a unit mass of fluid from $T - T_p^*$ to $T + T_p^*$ should include the effects of both the specific latent heat and the specific heat capacity. Integrating the specific heat capacity over the phase transition region, we get

\begin{align}
\int_{T_p^* - \Delta T^*}^{T_p^* + \Delta T^*} c_v\, dT &=
\int_{T_p^* - \Delta T^*}^{T_p^* + \Delta T^*} c_v^{\text{vaporization}} \, dT \\
&= c_v^{\text{liquid}} \Delta T^* + c_v^{\text{vapor}} \Delta T^* + \Delta H_p^{\text{vap}}, \label{Eq:cVapInt}
\end{align}

where $\Delta H_p^{\text{vap}}$ is the (specific) heat of vaporization. If $c_v^{\text{liquid}}$ or $c_v^{\text{vapor}}$ are not constant in this region, Eq. (<ref>) can be changed accordingly with, for example, the term $c_v^{\text{liquid}} \Delta T^*$ replaced with $\int_{T_p^* - \Delta T^*}^{T_p^*} c_v^{\text{liquid}}\, dT$. For notational brevity, henceforth, we only consider the case of constant $c_v^{\text{liquid}}$ and $c_v^{\text{vapor}}$.

Furthermore, we consider the numerical term $c_v^{\text{vaporization}}$ to also be constant, which leads to $\int_{T_p^* - \Delta T^*}^{T_p^* + \Delta T^*} c_v^{\text{vaporization}} \, dT  = c_v^{\text{vaporization}} 2 \Delta T^*$ and thus, the numerical notion of heat capacity of vaporization is given by

\begin{equation}
c_v^{\text{vaporization}} =  \frac{1}{2}( c_v^{\text{liquid}} + c_v^{\text{vapor}} ) + \frac{\Delta H_p^{\text{vap}}}{2\Delta T^*}.
\end{equation}

§.§ Conductivity, viscosity, surface tension

The heat conductivity $\lambda$ is given by a heat conductivity for the liquid phase $\lambda^{\text{liquid}}$ at $T < T_p^* - \Delta T^*$, and a heat conductivity for the vapor phase $\lambda^{\text{vapor}}$ used at $T > T_p^* + \Delta T^*$. The heat conductivity in the boiling temperature region is taken to be a linear interpolation between these two values. The effective heat conductivity is taken as a sum of the laminar part $\lambda$ and the turbulent part, given by

\begin{equation}
\lambda^{\text{turbulent}} =
\left( c_{\mu} \frac{k^2}{\varepsilon} \right) \rho c_v.
\end{equation}

Similarly, the viscosity is given by the liquid phase viscosity $\eta^{\text{liquid}}$ and a vapor phase viscosity $\eta^{\text{vapor}}$, with a linear interpolation between the two in the boiling temperature zone. The effective viscosity is taken as a sum of the laminar part $\eta$ and the turbulent part, given by

\begin{equation}
\eta^{\text{turbulent}} =
\left( c_{\mu} \frac{k^2}{\varepsilon} \right) \rho.
\end{equation}

Similarly, the surface tension is also given by a liquid phase and vapor phase surface tension with a linear interpolation in the boiling temperature zone, without any turbulent surface tension model.

§.§ Density

As mentioned above, the density for both the liquid and vapor phase is given by the closure relation, depending on both pressure and temperature $\rho = \rho(T,p)$, where $T$ and $p$ are a function of both time $t$ and space $\vec{x}$. At a fixed pressure $p$, we use the notation $\rho(T,p) = \rho_p(T)$.

The choice of the density in the phase change temperature region can have a significant impact on the numerical stability of the complete simulation, and thus plays a much larger role than that of the materials in Section <ref>. The simplest option is a linearly varying density in the phase transition range, as done in Section <ref>. Another option is a trigonometric “smearing out", as is employed by many phase field methods for diffuse interface capturing [14]. Instead, we choose a weighted harmonic interpolation.

For this, we first define the evaporation state variable

\begin{equation}
x^{\text{vapor}} = \frac{T - (T_p^* - \Delta T^* )}{ 2 \Delta T^*}.
\end{equation}

Using this, and for a density $\rho^{\text{liquid}}$ at $T = T_p^* - \Delta T^*$, and $\rho_p^{\text{vapor}}$ (could be pressure dependent) at $T = T_p^* + \Delta T^*$, the density is determined by

\begin{equation}
\rho(T,p) =
\begin{cases}
\rho^{\text{liquid}},       & \text{if }   T - T_p^*   < -\Delta T^*, \\
\left( \rho^{\text{liquid}}\rho_p^{\text{vapor}} \right) / \left( x^{\text{vapor}} \rho^{\text{liquid}} + (1-x^{\text{vapor}}) \rho_p^{\text{vapor}} \right),                     & \text{if } | T - T_p^* | <  \Delta T^*, \\
\rho_p^{\text{vapor}},      & \text{if }   T - T_p^*   >  \Delta T^*.
\end{cases}
\end{equation}

Note that $\rho^{\text{liquid}}$ and $\rho_p^{\text{vapor}}$ can be temperature dependent.

§ MESHFREE DISCRETIZATION

In this section, we describe the meshfree point cloud based domain discretization, and the meshfree collocation approach to derivative computation. We then go on to introduce an enrichment procedure for derivative computation that is necessary to handle the large jumps and large gradients in the material properties.

§.§ Meshfree point cloud based discretization

The discretized computational domain $\Omega$, consisting of both the fluid and the vapor phases modeled as a single fluid, is composed of a cloud of $N = N(t)$ points. This includes points both in the interior of the domain and on all boundaries. The initial point cloud to discretize the domain is generated using a meshfree advancing front method [40].
These points are approximation locations or collocation nodes, and not mass-carrying particles. Each point carries all physical and material properties of the fluid system. For every point $i = 1,2,\dots, N$, at location $\vec{x}_i$, numerical approximations are done on a set of neighboring points, referred to as its support or neighborhood

\begin{equation}
S_i = \{ j \in \{1,2,\dots, N\} | \, \| \vec{x}_j - \vec{x}_i  \| \leq h_i  \},
\end{equation}

where $h_i = h(\vec{x}_i, t)$ is the smoothing length or interaction radius. Note that we have $i \in S_i$. An example of proximity-based neighborhoods is illustrated in Figure <ref>.

A proximity-based neighborhood on a meshfree point cloud. For the point marked with the black cross, the black circle encloses the neighboring points over which all derivative approximations are carried out.

§.§ Adaptivity

Since points are moving with the fluid velocity in a Lagrangian sense, the point cloud could become distorted. This would manifest in either points coming too close to each other, or hole formations in the point cloud where no points are present. In contrast to mesh distortion in mesh-based Lagrangian methods, this point cloud distortion can be easily fixed locally. This is done in two parts: (i) points that are closer than $r_{\text{min}}h$ apart are merged into a single point at the center location, and (ii) holes of size $r_{\text{max}}h$ in which no points are present are identified and filled with a point at the center of the hole. In both cases, all physical properties are interpolated at a new point location. Here, $r_{\text{min}}$ and $r_{\text{max}}$ are fixed application-independent constants that determine the number of points in the support domain. For more details about the procedures of identification of holes, the interpolation, and the thinning/filling procedures in general, we refer to [8, 37, 42]. The same distortion fixing algorithms are also used for spatial adaptivity of the discretization. For example, if large material property gradients are detected, $h$ can be reduced, which would automatically trigger the hole-filling algorithm to refine the point cloud in that region.

§.§ Spatial derivatives

Here, we compute numerical derivatives with the strong form meshfree generalized finite different method (GFDM) [1, 11, 23, 24]. GFDM is a robust method and has been used for a variety of applications, especially in modeling fluid flow, see for example [17, 26, 27, 47]. The meshfree GFDM is a collocation approach, and as the name suggests, they generalize classical finite differences to arbitrarily distributed points [23]. The derivatives of a function $u$ are approximated by a linear combination of discrete function values, $u_i$ for each $i = 1, \dots, N$, of each neighboring point

\begin{equation}
\label{Eq:GFDM_Definition}
\partial^* u(\vec x_i)\approx \partial^*_i u = \sum_{j\in S_i}c_{ij}^* u_j,
\end{equation}

where ${}^*=x,y, z, \Delta$ represents the differential operator being approximated, $\partial^*$ represents the continuous $^*$-derivative, and $\partial^*_i$ represents the discrete derivative at point $i$. For each point $i$, the stencil coefficients $c_{ij}^*$ are found using a weighted least squares approach. The numerical differential operators can be found locally at each point by solving a small linear system, independently of the operators at the rest of the point cloud. Since the point cloud moves in time, the derivative stencils need to be recomputed for each point at every time step.

For a point $i$, consider Taylor expansions around $\vec{x}_i$ up to order $2$ terms for each neighbouring point $j\in S_i$

\begin{equation}
\label{Eq:TaylorExp}
e_{ij} + u(\vec{x}_j)=u(\vec{x}_i)+\nabla u \cdot (\vec{x}_j-\vec{x}_i) + \frac12 (\vec{x}_j-\vec{x}_i)^T D (\vec{x}_j-\vec{x}_i).
\end{equation}

The unknown coefficients of $\nabla u$ and $D$ are computed by a weighted least squares method by minimizing

\begin{equation}
\label{Eq:minJ}
\text{min } J_i = \sum _{j \in S_i} W_{ij}^2e_{ij}^2,
\end{equation}

where $W$ is a weighting function used to make sure that the points closer to the central point $i$ have a larger impact than the points farther away. We note that the square of the weighting function is used only for notational convenience. The weighting function is usually taken as a Gaussian distribution

\begin{equation}
\label{Eq:WeigthingKernel}
W_{ij} = \exp\!\left(-\alpha_W \frac{\|\vec{x}_j-\vec{x}_i\|^2}{h_i^2+h_j^2}\right),
\end{equation}

where $\alpha_W$ is a positive constant usually taken in the range of $(2,8)$. Note that the weighting function is only defined on the local support $S_i$ consisting of $n(i)$ points. For the sake of brevity, we present only the case of one spatial dimension. Eq. (<ref>) leads to the following system which is solved at each point $i=1,\dots, N$

\begin{equation}
\label{Eq:GFDAsystem}
\underbrace{\left(\begin{array}{c}
e_{i1} \\
\vdots \\
\end{array}\right)}_{\vec E_i} =
\underbrace{ \left(
\begin{array}{ccc}
1 & \delta x_{i1}     & \frac12\delta x_{i1}^2    \\
\vdots & \vdots         & \vdots          \\
1 & \delta x_{in}   & \frac12\delta x_{in}^2  \\
\end{array}  \right) }_{M_i}
\underbrace{\left(\begin{array}{c}
(u_0)_i \\
(u_x)_i \\
\end{array}\right)}_{\vec a_i} -
\underbrace{\left(\begin{array}{c}
\vdots \\
\end{array}\right)}_{\vec b_i},
\end{equation}

where $\delta x_{ij} = x_j - x_i$. Or, in short form $\vec{E}_i = M_i\vec{a_i} - \vec{b}_i$. The minimization Eq. (<ref>) can be rewritten as

\begin{align}
\text{min } J_i &= \vec E_i ^{\,T} W_i^2 \vec E_i\\
&= (M_i\vec a_i-\vec b_i)^T W_i^2 (M_i\vec a_i - \vec b_i) 	\label{Eq:minJ2},
\end{align}

where $W_i$ is a diagonal matrix with entries $W_{i1},\dots,W_{in}$. A formal minimization leads to

\begin{equation}
\label{Eq:FinalSolution_GFDM}
\vec a_i = [(M_i^T W_i^2 M_i)^{-1} M_i^T W_i^2 ]\vec b_i.
\end{equation}

This leads to the differential operator stencils as in Eq. (<ref>)

\begin{align}
(u_x)_i &= \sum_{j \in S_i} c_{ij}^x u_j,\label{Eq:ux}\\
(u_{xx})_i &= \sum_{j \in S_i} c_{ij}^{xx} u_j, \label{Eq:uxx}
\end{align}

where $c_{ij}^x$ and $c_{ij}^{xx}$ represent the values in the second and third row respectively of the matrix $[(M_i^T W_i^2 M_i)^{-1}M_i^T W_i^2]$ in Eq. (<ref>).

We note that similar GFDM approaches have already been used to simulate phase change processes [32, 36, 52]. They model solidification at low velocities, while the present work models a lot more turbulent phase change process in vaporization. Furthermore, [32, 36, 52] do not take the latent heat into account directly, as done in the presented model.

§.§.§ Polynomial Formulation

We now rewrite the Taylor expansion-based derivation to a polynomial formulation that can be easily extended for our purpose. An alternative, but equivalent, way to obtain the stencil coefficients in  Eq. (<ref>) – Eq. (<ref>) is to ensure that the derivatives of monomials $m\in\mathcal{M}$ up to the order of accuracy desired are exactly reproduced, leading to the optimization problem

\begin{align}
\sum_{j\in S_i}c_{ij}^{*}m_j &= \partial^* m (\vec{x}_i),\qquad \forall m\in\mathcal{M},\label{Eq:Consistency}\\
\text{min } J_i &= \sum_{j\in S_i} \left( \frac{c_{ij}^{*}}{W_{ij}} \right)^2. \label{Eq:BasicMin}
\end{align}

where ${}^* \in \{x, y, xx, \Delta, \dots \}$.

In our earlier work [39], we have proven the equivalence of these formulations. An efficient method to compute the differential operator stencils in this formulation using a QR decomposition can also be found in [39].

§.§ Enrichment

Due to the smooth shape functions and the strong form nature of the GFDM, representing discontinuities is not straightforward. Using the polynomial formulation of Section <ref>, we now introduce an enrichment mechanism that is needed to accurately capture the sharp changes in the material properties, which is especially necessary for handling jumps and discontinuities. Motivated by our earlier work for two-dimensional manifolds [43], we extend the polynomial space of the test functions to include test functions with discontinuities.

We note that the differential operator for which discontinuities play the largest role is the $\nabla \cdot ( \alpha \nabla \phi )$ for some material property $\alpha$ and some primary variable $\phi$. For example, the term $\nabla \cdot ( \frac{1}{\rho} \nabla p )$ arises in the pressure Poisson equation,  and $\nabla \cdot ( \lambda_{\text{eff}} \nabla T )$ in the energy conservation equation. Similarly, the momentum conservation equation contains the divergence of the stress tensor which includes the term $\nabla \cdot \left( \eta_{\text{eff}} \nabla  \left( \left( \nabla \vec{v} \right) + \left( \nabla \vec{v} \right)^T  \right) \right)$. Each of these contains a diffusion operator of the form $\nabla \cdot ( \alpha \nabla \phi )$. Here, we introduce an enrichment of this numerical differential operator to handle jump discontinuities in material properties.

Using the shorthand $D^{\alpha}\phi = \nabla\cdot(\alpha\nabla\phi)$, the diffusion operator is computed as

\begin{align}
\sum_{j\in S_i}c_{ij}^{D^{\alpha}}m_j &= D^{\alpha} m (\vec{x}_i),\qquad \forall m\in\mathcal{M},\label{Eq:DEG_Consistency}\\
\text{min } J_i &= \sum_{j\in S_i} \left( \frac{c_{ij}^{D^{\alpha}}}{W_{ij}} \right)^2, \label{Eq:DEG_Min}
\end{align}

in the polynomial formulation. The space $\mathcal{M}$ of consistency conditions in Eq. (<ref>) is enhanced with discontinuous functions, in a manner similar to that done in [60, 61] in the context of crack propagation. For $\vec{s} = \frac{\nabla \alpha}{\| \nabla \alpha \|}$ denoting the normalized gradient of the material property being considered, the extra conditions added are given by

\begin{align}
\sum_{j\in S_i}c_{ij}^{D^{\alpha}} \frac{1}{\alpha_j} &=  - \Delta \log \alpha,\\
\sum_{j\in S_i}c_{ij}^{D^{\alpha}} \frac{\delta s_{ij}}{\alpha_j} &= - \frac{\partial}{\partial s} \log \alpha,\\
\sum_{j\in S_i}c_{ij}^{D^{\alpha}} \frac{\left( \delta s_{ij} \right)^2}{\alpha_j} &= 2,
\end{align}

where $\delta s_{ij} = \vec{s} \cdot \delta \vec{x}_{ij}$, and $\frac{\partial}{\partial s} = \vec{s} \cdot \nabla$.

§ NUMERICAL SCHEME

The discretization of all spatial derivatives is done as explained in Section <ref>. The time integration of the constitutive equations in Section <ref> is done using a segregated approach which forms a Chorin-type projection method [3, 5]. The integration begins with the Lagrangian movement of the point cloud, done with a second-order method [41]. Then an intermediate velocity is computed by a first-order implicit discretization of the momentum equation. This is then projected to the space of the desired velocity divergence with the help of a pressure Poisson equation, followed by a pressure update. After this, a Crank-Nicolson time-integration scheme is used for the energy conservation followed by an implicit integration of the turbulence equations [26]. For details on the entire time-integration scheme, we refer to our earlier work [8, 17, 19, 26].

§.§ Temperature smoothing

The strong form nature of the derivative computation can lead to numerical instabilities while trying to capture jumps. Since the temperature field can have jump discontinuities due to the phase change process, we numerically smooth the temperature with a Gaussian smoothing kernel [42]. The smoothed temperature field at a point $i$ is given by

\begin{equation}
\widetilde{T}_i = \frac{ \sum_{j \in S_i} K_{ij} T_j}{ \sum_{j \in S_i} K_{ij}},
\end{equation}

where the smoothing kernel is given by

\begin{equation}
K_{ij} = \exp\left( -2 \frac{\| \vec{x}_j - \vec{x}_i \|^2}{h_i^2 + h_j^2}  \right).
\end{equation}

All material properties in Section <ref> dependent on temperature are taken to be dependent on $\widetilde{T}$ numerically. As a result, this smoothed temperature field automatically has a smoothing effect on, for example, the numerical density field. Empirically, we observe that this helps ensure numerical stability.

§.§ Volume expansion

An important challenge in this macroscopic phase change modeling is to accurately capture the volume expansion (or contraction) as a result of the phase change. If a numerical point is in the liquid state at a certain time step, and subsequently in the vapor state in the next time step, this could result in a large time gradient in the density, which could lead to numerical instabilities when using the strong form GFDM. To overcome this issue, we rewrite the mass conservation in logarithmic form

\begin{align}
%	\frac{D \rho}{Dt} &= -\rho \nabla \cdot \vec{v} \,,\\
\nabla \cdot \vec{v}\, \vert^{(n+1)} &= - \frac{d}{dt} \left( \log ( \rho^{(n+1)} ) \right),\\
&\approx - \frac{1}{\Delta t} \left( \log ( \rho^{(n+1)}) - \log ( \rho^{(n)}) \right),
\end{align}

which is used as the desired divergence of velocity in the projection step and the pressure Poisson equation.

§ COMPARISON WITH INTERFACE METHODS

Verification of the present one-fluid phase change model used here can not be done using standard numerical benchmark cases used in interface capturing or tracking methods. This holds for the benchmarks used by both two-fluid methods, and one-fluid methods with explicit interface capturing or tracking, which we club together under the term interface methods. There are several modeling choices made by these interface methods (see, for example, [6, 16, 31, 48]) that do not have a one-to-one correspondence in the one-fluid model considered here.

* Interface methods typically require boundary conditions at the interface. For the temperature field, this could be a Dirichlet boundary condition or a prescribed heat flux. Occasionally, velocity boundary conditions are also imposed. This is not possible in the present one-fluid model since there is no explicit boundary between the phases.

* The interface methods often use verification cases that distinguish between liquid and vapor phases at saturation temperature. Whereas, in the present model, we only have a fluid at saturation temperature that cannot be distinguished between liquid and vapor. Furthermore, many verification tests in the literature assume an initial condition of a “super-saturated" liquid in which a liquid phase is present above the boiling point. This, once again, is not directly possible in the present approach.
We note here that since the aim of this work is to model phase change in highly dynamic fluid flow, the topic of modeling metastable super-heated liquid states is not relevant.

* A significant advantage of the present one-fluid model is that it can be used to simulate scenarios where initially only a liquid phase is present without any vapor phase. To achieve this, interface methods assume a very small vapor phase at the initial state of the simulation. However, in applications where the boiling process occurs at a fluid free surface (see the jet impinging on a hot plate test case in Section <ref>), assuming the presence of a vapor phase near an evolving free surface is not feasible.

§ NUMERICAL SIMULATIONS

All the numerical methods mentioned above have been incorporated into the in-house software suite MESHFREE [12].

The verification of the basic numerical schemes and the validation of the underlying numerical framework for fluid flow can be found in our earlier work [8, 37, 45]. In all the test cases considered below, radiative heat transfer is neglected. We further note that all simulations are full three-dimensional simulations. However, pseudo-two-dimensional clips of the results are shown for ease of visualization.

§.§ Vapor film growth around a heated sphere

We start by considering the growth of a vapor film around a heated solid sphere, shown in Figure <ref>, as also considered by [6]. A heated sphere is immersed in the middle of a fluid in a tank with an open top. The fluid tank has a square base area with a side length of 4. The tank is filled with a resting liquid of initial temperature 89 up to a height of 4. A solid sphere with a diameter of 1.5 is submerged in the middle of the tank.

Simulation parameters:
The liquid and vapor densities are taken to be constant at $\rho^{\text{liquid}} = \qty{1000}{\kilogram\per\cubic\meter}$ and $\rho^{\text{vapor}} = \qty{10}{\kilogram\per\cubic\meter}$ respectively. The phase change takes place within the vaporization interval of $\Delta T^* = \qty{5}{\degreeCelsius}$ around the boiling point of $T_p^* = \qty{100}{\degreeCelsius}$, while neglecting pressure dependence of the boiling point. The specific heat capacities are $c_{v}^{\text{liquid}} = \qty{4184}{\joule\per\kilogram\per\kelvin}$ and $c_{v}^{\text{vapor}} = \qty{2030}{\joule\per\kilogram\per\kelvin}$ for the liquid and vapor, respectively. The latent heat of vaporization is $\Delta H_p^{\text{vap}} = \qty{2501}{\kilo\joule\per\kilogram\per\kelvin}$. The heat conductivities are $\lambda^{\text{liquid}} = \qty{0.6}{\watt\per\meter\per\kelvin}$ and $\lambda^{\text{vapor}} = \qty{0.05}{\watt\per\meter\per\kelvin}$ for the liquid and vapor respectively. The liquid has a surface tension of $\sigma = \qty{0.07286}{\newton\per\meter}$ and a viscosity of $ \eta^{\text{liquid}} = \qty{0.001}{\pascal\second}$, while the vapor has a viscosity of $ \eta^{\text{vapor}} = \qty{0.00013}{\pascal\second}$, with no surface tension. The domain is resolved with a coarse resolution of $h_0$ and a fine resolution of $h_{\min} = 0.1 h_0$ near the sphere surface. The resolution of the point cloud is gradually increased radially at a rate of $dh / dr = 0.15$. At the walls and the free surface at the top, we apply a zero heat flux boundary condition. The heat transfer between the solid sphere and the fluid is governed by a heat transfer coefficient of $U_s = \qty{12500}{\watt\per\square\meter\per\kelvin}$. The temperature on the sphere surface $T_S$ is growing linearly in time from $T_S(\qty{0}{\second}) = T^* + \Delta T^*$ to $T_S(\qty{6}{\second}) =  \qty{400}{\degreeCelsius}$.

A thin vaporization area forms initially with relatively small values of $x^{\text{vapor}}$. The lower density of the vapor causes it to rise under gravity. As temperature $T_S$ increases, the values of $x^{\text{vapor}}$ rise, which results in the formation of vapor bubbles around the film, which leave the film as they rise and break at the free surface of the fluid. The vapor bubble formation and breaking are shown in Figure <ref>, with Figure <ref> showing a detailed sketch of a particular bubble formation. With rising $T_S$, the fluid in the vicinity of the sphere tends to $x^{\text{vapor}}=1$, which provides an isolation layer around the sphere, illustrated in Figure <ref>. The lower heat conductivity of the vapor phase results in a rapid drop in the heat flux from the hot sphere to the fluid. These effects are also observed in experiments and are often referred to as the Leidenfrost-effect which prevents the liquid beyond the vapor film from boiling. This effect is quantified in Figure <ref> which shows the integrated heat transfer over the surface of the sphere for different resolutions. As the temperature rises, the heat transfer is rising for all resolutions. However, for each resolution, we observe the aforementioned rapid Leidenfrost-drop of the heat transfer at different temperatures. For finer resolutions, the temperatures of the drop come closer to each other, however, it seems that the numerical resolution still does not meet the thickness of the real vapor layer as it would arise in reality. The temperature at which the Leidenfrost-effect occurs is dominated by the choice of $U_s$, so experiments around the Leidenfrost state could, in the future, be used for calibrating the heat transfer coefficient.

Vaporization around a heated sphere: Formation of unsteady vapor films, bubble growth and rising under gravity. The figures show a clip of the three-dimensional domain at increasing times. The color indicates the value of $x^{\text{vapor}}$, with the red points representing higher values and the blue points representing lower ones. The arrows show the direction of the velocity. A zoomed-in version of the bubble formation is shown in Figure <ref>.

Vaporization around a heated sphere: Bubble growth around an unsteady vapor film layer. The color shows the density. A zoomed-out version of this can be found in Figure <ref>.

Vaporization around a heated sphere: Formation of thin, but steady overheated vapor layer and rising flow under gravity. This is the Leidenfrost-state, where $x^{\text{vapor}}$ tends to 1 around the sphere surface.

Vaporization around a heated sphere: Heat flux across the sphere surface to the liquid.

As a verification, we also run this same test case in the absence of gravity. As also observed in [6], without gravity, the vapor film formed is of a uniform thickness and is symmetric, see Figure <ref>. Unlike the case with gravity, the vapor film formed here is stable throughout the simulation without the formation of bubbles, even at lower temperatures of the sphere.

Vaporization around a heated sphere without gravity: Formation of a uniform thickness vapor film, which grows slowly in time. Only a clip of the three-dimensional domain is shown. The color shows the value of $x^{\text{vapor}}$.

§.§ Cooking pot

We now simulate nucleate boiling in a cooking pot of diameter 0.1 filled with a liquid of initial temperature of 85 to a height of $H = \qty{0.15}{\meter}$. An inner circle with diameter 0.05 in the middle of the bottom surface of the pot is heated by linear ramping to $T = \qty{400}{\degreeCelsius}$ in 1. The heat transfer coefficient at the hot wall is $U_s = \qty{12500}{\watt\per\square\meter\per\kelvin}$. On other surfaces, a zero heat flux boundary condition is imposed. The material properties of the liquid are as described in the previous test case. We simulate this case with a constant smoothing length over the liquid domain. Three simulations of $h = \qty{0.002}{\meter}, \qty{0.004}{\meter}$, and $\qty{0.008}{\meter}$ respectively are considered. The resulting evolution of the boiling in the cooking pot is shown in Figure <ref>.

Boiling in a cooking pot: Snapshots of instantaneous density distribution at 1.2444 (left side), 1.3269 (middle), and 1.4177 seconds (right side) in the cooking bowl. (Top: bottom surface view, bottom: side view).

To quantify the results, we study the energy flux as a result of the fixed heat transfer coefficient. For all resolutions considered, the heat transfer is rising up to approximately $\qty{9000}{\watt}$ within the first second and then remains roughly constant. We observe temporary breakdowns of the heat transfer, especially for finer resolutions, see Figure <ref>. They might be explained by temporary vapor films covering the heating zone at the bottom and thus acting as an insulating barrier, see Figure <ref>. The vapor film is metastable (Rayleigh-Taylor-like behavior) and quickly disappears. With progress in time, the breakdowns disappear even for the simulation with the finest resolution, which is due to the increasingly chaotic nature of the flow.

Convergence study for heat flux in a cooking pot.

The energy balance shown in Figure <ref> compares the internal heat of the fluid to the energy that enters the fluid through the walls. This energy conservation is important to verify in meshfree simulations since most meshfree methods often have issues with global conservative behavior [44]. The figure shows that changes in internal energy and the energy transferred through the walls are almost the same for finer resolutions. The measured internal energy contains fluctuations, which come from the measurement of internal energy as $E_{\text{internal}}=\int_V \rho h dV \approx \sum_i h_i \rho_i V_i $, where $V_i$ is an approximation of the volume taken up by a numerical point, which inherently has fluctuations.

Energy balance in the cooking pot. The figure on the right shows a zoomed in part of the figure on the left.

§.§ Fluid jet on a hot plate

As our third test case, we consider an example where boiling occurs at the free surface of the fluid, and compare simulations with experimental results. This example is motivated by the boiling of metalworking fluid introduced in Section <ref>. The solid phase is a large circular plate made of Inconel 718 of diameter 0.14 and thickness of 0.005, which is also discretized with a meshfree method, considering only energy conservation with zero velocity. A fluid jet is simulated with an impact angle of 60 on the solid plate. The inlet is considered very close to the plate (at a distance of 0.005 above) to save simulation time. The jet inlet has a diameter 0.003 and velocity 7.0.

The fluid impinging on a hot metal plate at an angle is very similar to typical conditions for simulating cooling jets in metal cutting applications [49].
These parameters are selected to match the experimental setup described by [28], except for the distance between the nozzle exit and the plate. The nozzle exit is placed closer to the plate to reduce the number of points in the simulation and thus, save computational costs. However, since all other parameters match the experimental setup, we expect no significant change in the global jet and especially the vaporization zones and cooling areas.
For a visualization of this scenario being simulated, Figure <ref> shows the situation without the consideration of the vaporization model.

Fluid jet impinging on a hot plate: Visualisation of the scenario being simulated without the consideration of the phase change model. Side view (top row), and top view (bottom row) at different times. The solid plate is shown in gray, while the fluid (purely liquid) is shown in blue.

Simulation parameters: The vaporization interval is $\Delta T_p= \qty{5}{\degreeCelsius}$ around the boiling point at $T_p=\qty{100}{\degreeCelsius}$. The  latent heat of vaporization is $\Delta H_{p,vap} = \qty{2501}{\kilo\joule\per\kilogram\per\kelvin}$. Initially, the temperature of the plate is $T_0 = \qty{400}{\degreeCelsius}$. At the top face between the plate and the fluid, the heat transfer coefficient is $U_s = \qty{50000}{\watt\per\square\meter\per\kelvin}$. The ambient air is taken at $T = \qty{20}{\degreeCelsius}$. On all other walls of the solid, and at the inflow of the fluid phase we apply a zero heat flux boundary condition. A constant smoothing length is set at $h_l = \qty{0.00075}{\meter}$ for the liquid/vapor and at $h_s = \qty{0.0009}{\meter}$ for the solid.
The material properties are described in Table <ref>. In contrast to the experiment, the vapor density is much larger. This is to avoid abrupt changes in the simulation. Hence, the vapor layer will be thinner. However, we do not expect a generally altered distribution and flow of the liquid. For this, we presume that the time scales of the expansion due to vaporization is much smaller than the time scales of convection. For the heat transfer by heat conduction, this means a different extent of the vapor layer. To counteract this effect we apply a decreased value of the heat conductivity, based on the ratio of simulated vapor density to real vapor density (0.50489).

Phase     Density         Surface         Dynamic       Specific               Heat

tension         viscosity     heat capacity          conductivity

Liquid     998      72.86     1        4185     0.59801

Vapor      100      -         0.13     2059     0.00012405

Solid      8118     -         -         479     13.0

Material properties for the fluid jet on a hot plate test case.

The resultant vaporization process is shown in Figure <ref>. As the fluid jet impinges on the hot plate, it forms thin liquid layers on the solid surface. These are slowly heated by heat transferred from the solid. Since the inflow of the jet has a constant temperature, the thinnest regions of the fluid film which are further away from the jet impinging location reach boiling point and vaporize. As the vaporization starts and the hot plate is cooled down, the impinging jet pushes the initial vaporization regions further away. Figure <ref> also shows the vapor phase rising as a result of gravity. These observations qualitatively match experimental observations.

Fluid jet impinging on a hot plate: Top view of the evolution of the boiling process at different time steps (see Figure <ref> for a side view). The color indicates the density. Top row: initial impingement and formation of the thin fluid layer. Middle row: Start of vaporization, and the jet pushing away the vaporization regions further. Bottom row: Formation of a temporally stable ring of vaporization around the main impingement location.

Figure <ref>(a) shows a comparison between experimental and numerical results for the time history of the temperature at the center point of the rear solid surface. In the experiments [28], there is an oscillation initially when the temperature drop begins. This oscillation then stabilizes and both curves follow the same trend. In Figure <ref>(b) the temperature distribution on the rear surface is displayed for two time instants in comparison between simulation and experiment. Compared to the experiment the wetting area is qualitatively the same. However, in the experiment, there is an asymmetry and the growth of the wetting area is more pronounced. This shows some differences in the local distribution, which indicate the potential for improvement for example by more detailed modeling of the heat transfer coefficient.

xmin = 0, xmax = 0.8, ymin = 369, ymax = 401, xtick distance = 0.1, ytick distance = 4, grid = both, minor tick num = 1, major grid style = lightgray, minor grid style = lightgray!25,
width = , height = , xlabel = time in , ylabel = temperature in , legend pos=south west]

[ only marks, mark=x,mark size = .8pt, blue
] file[skip first] ./Figures/T_over_time_Experiment.dat;

[smooth, black,
] file[skip first] ./Figures/T_over_time_Simulation.dat;

(a) Time history of the temperature in the center of the bottom solid surface. (b) Comparison of temperature at the bottom solid surface between experiment (top row) and simulation (bottom row) for two time instants.

§ CONCLUSION

We presented a novel one-fluid approach to simulating phase change processes, with a focus on fluid-to-vapor phase change at the boiling point. While conventional approaches to this problem explicitly model the interface between the fluid and vapor phases, either by a sharp interface or a smooth interface approach, our proposed method does not model the interface at all. The phase change process is modeled only through varying material properties, with the latent heat of phase change smeared out over a small temperature range around the boiling point. This enables us to simulate phase changes in complex flow domains, such as vaporization at free boundaries, which is extremely challenging with conventional interface-based methods. Our numerical results show the applicability of the method to such scenarios.

Simulations were done using a Lagrangian meshfree approach using a strong form generalized finite difference method to approximate spatial derivatives. To capture the large gradients in the material properties, which is often challenging using a strong form approach, we enrich the test functions space to include discontinuous functions.

While the numerical results suggest that one-fluid phase change modeling is very promising, they also raise several questions which need to be investigated. One such question is the optimal choice of the pressure interpolation in the phase change region. Another open question concerns the appropriate modeling of the very thin vapor layers and using an adaptive refinement based on the thickness of the vapor layers. Several possible improvements of the model proposed could be considered. This includes a more detailed modeling for the heat transfer coefficient by data-driven or multi-scale modeling [34, 38, 56], the incorporation of a turbulent thermal conductivity, as well as a modeling of local time scales for vaporization rates which could, for example, describe a pressure dependency.

Numerical results for the vaporization when a jet with free boundaries impinges on a hot surface, including comparisons with experiments, suggest suitability for the present method to model the vaporization of cutting fluid in metal cutting applications. Our future work is directed towards the application of the newly developed model to vaporization during metal cutting, with more detailed validations against experiments. First simulations have also shown the ability to simulate the Leidenfrost effect. A planned extension of this work is to thoroughly examine and validate the modeling of the Leidenfrost effect using the presented one-fluid phase transition model.

§ ACKNOWLEDGEMENTS

All the authors would like to acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) under the priority program SPP 2231 “FLUSIMPRO", project number 439626733.
Pratik Suchde would also like to acknowledge partial support from  the  European  Union’s  Horizon  2020 research  and  innovation program under the Marie Skłodowska-Curie Actions grant agreement No. 892761 “SURFING".

[1]
J. Bašić, N. Degiuli, Š. Malenica, and D. Ban.
Lagrangian finite-difference method for predicting green water
Ocean engineering, 209:107533, 2020.

[2]
K. Binder, K. Vollmayr, H.-P. Deutsch, J. D. Reger, M. Scheucher, and D. P.
Monte carlo methods for first order phase transitions: some recent
International Journal of Modern Physics C, 3(05):1025–1058,

[3]
D. L. Brown, R. Cortez, and M. L. Minion.
Accurate projection methods for the incompressible Navier–Stokes
Journal of Computational Physics, 168(2):464–499, 2001.

[4]
E. Can and A. Prosperetti.
A level set method for vapor bubble dynamics.
Journal of Computational Physics, 231(4):1533–1552, 2012.

[5]
A. J. Chorin.
Numerical solution of the navier-stokes equations.
Mathematics of computation, 22(104):745–762, 1968.

[6]
A. Das and P. Das.
Modeling of liquid–vapor phase change using smoothed particle
Journal of Computational Physics, 303:125–145, 2015.

[7]
J. J. de Pablo, Q. Yan, and F. A. Escobedo.
Simulation of phase transitions in fluids.
Annual Review of Physical Chemistry, 50(1):377–411, 1999.
PMID: 15012417.

[8]
C. Drumm, S. Tiwari, J. Kuhnert, and H.-J. Bart.
Finite pointset method for simulation of the liquid - liquid flow
field in an extractor.
Computers & Chemical Engineering, 32(12):2946 – 2957, 2008.

[9]
S. L. Duncan, I. S. Dalal, and R. G. Larson.
Molecular dynamics simulation of phase transitions in model lung
surfactant monolayers.
Biochimica et Biophysica Acta (BBA) - Biomembranes,
1808(10):2450–2465, 2011.

[10]
A. Esmaeeli and G. Tryggvason.
A front tracking method for computations of boiling in complex
International Journal of Multiphase Flow, 30(7-8):1037–1050,

[11]
C.-M. Fan, C.-N. Chu, B. Šarler, and T.-H. Li.
Numerical solutions of waves-current interactions by generalized
finite difference method.
Engineering Analysis with Boundary Elements, 2018.

[12]
Fraunhofer Society.

[13]
F. Gibou, L. Chen, D. Nguyen, and S. Banerjee.
A level set based sharp interface method for the multiphase
incompressible navier–stokes equations with phase change.
Journal of Computational Physics, 222(2):536–555, 2007.

[14]
H. Gomez, M. Bures, and A. Moure.
A review on computational modelling of phase-transition problems.
Philosophical Transactions of the Royal Society A: Mathematical,
Physical and Engineering Sciences, 377(2143):20180203, 2019.

[15]
E. Goncalves and R. F. Patella.
Numerical simulation of cavitating flows with homogeneous models.
Computers & Fluids, 38(9):1682–1696, 2009.

[16]
S. Hardt and F. Wondra.
Evaporation model for interfacial flows based on a continuum-field
representation of the source terms.
Journal of Computational Physics, 227(11):5871–5895, 2008.

[17]
A. Jefferies, J. Kuhnert, L. Aschenbrenner, and U. Giffhorn.
Finite pointset method for the simulation of a vehicle travelling
through a body of water.
In M. Griebel and A. M. Schweitzer, editors, Meshfree Methods
for Partial Differential Equations VII, pages 205–221, Cham, 2015. Springer
International Publishing.

[18]
D. Juric and G. Tryggvason.
Computations of boiling flows.
International Journal of Multiphase Flow, 24(3):387–410, 1998.

[19]
J. Kuhnert.
Meshfree numerical scheme for time dependent problems in fluid and
continuum mechanics.
In S. Sundar, editor, Advances in PDE Modeling and Computation,
pages 119–136, New Delhi, 2014. Anne Books.

[20]
B. Launder and D. Spalding.
The numerical computation of turbulent flows.
Computer Methods in Applied Mechanics and Engineering, 3(2):269
– 289, 1974.

[21]
M. S. Lee, V. Aute, A. Riaz, and R. Radermacher.
A review on direct two-phase, phase change flow simulation methods
and their applications.
International Refrigeration and Air Conditioning Conference,
Paper 1289, 2012.

[22]
Q. Li, K.H. Luo, Q.J. Kang, Y.L. He, Q. Chen, and Q. Liu
Lattice Boltzmann methods for multiphase flow and phase-change heat transfer
Progress in Energy and Combustion Science, 52:62–105, 2016.

[23]
T. Liszka and J. Orkisz.
The finite difference method at arbitrary irregular grids and its
application in applied mechanics.
Computers & Structures, 11(1-2):83–95, 1980.

[24]
Y. Lu, A.-k. Hu, and Y.-c. Liu.
A finite pointset method for the numerical simulation of free surface
flow around a ship.
Journal of Marine Science and Technology, 21(2):190–202, 2016.

[25]
M. Matsumoto.
Molecular dynamics of fluid phase change.
Fluid Phase Equilibria, 144(1):307–314, 1998.

[26]
I. Michel, T. Seifarth, J. Kuhnert, and P. Suchde.
A meshfree generalized finite difference method for solution mining
Computational Particle Mechanics, 8(3):561–574, 2021.

[27]
A. Möller and J. Kuhnert.
Simulation of the glass flow inside a floating process / Simulation
de l'écoulement du verre dans le procédé float.
Revue Verre, 13(5):28–30, 2007.

[28]
K. Nabbout and L. Pasternak and M. Sommerfeld and B. Bock-Marbach and J. Kuhnert and E. Barth and E. Uhlmann
Heat Transfer Study for Oil-in-water Emulsion Jets Impinging onto hot Metal Surface
Proceedings of the 8th World Congress on Mechanical, Chemical, and Material Engineering (MCM'22) , HTFF 190, 2022.

[29]
T. D. Nguyen, C. C. Nguyen, and V. H. Tran.
Molecular dynamics study of microscopic structures, phase transitions
and dynamic crystallization in ni nanoparticles.
RSC advances, 7(41):25406–25413, 2017.

[30]
L. S. Pan, J. Lou, H. Y. Li, and C. W. Kang.
A diffuse interface model for two-phase flows with phase transition.
Physics of Fluids, 31(9):092112, 2019.

[31]
A. Rajkotwala, A. Panda, E. Peters, M. Baltussen, C. van der Geld,
J. Kuerten, and J. Kuipers.
A critical comparison of smooth and sharp interface methods for phase
International Journal of Multiphase Flow, 120:103093, 2019.

[32]
E. O. Reséndiz-Flores and F. R. Saucedo-Zendejo.
Numerical simulation of coupled fluid flow and heat transfer with
phase change using the finite pointset method.
International Journal of Thermal Sciences, 133:13–21, 2018.

[33]
N. Samkhaniani and M. R. Ansari.
The evaluation of the diffuse interface method for phase change
simulations using openfoam.
Heat Transfer—Asian Research, 46(8):1173–1203, 2017.

[34]
Y. Sato and B. Niceno
A Depletable Micro-Layer Model for Nucleate Pool Boiling
Journal of Computational physics, 300:20–52, 2015.

[35]
A. Satoh.
Monte carlo simulations on phase change in aggregate structures of
ferromagnetic spherocylinder particles.
Colloids and Surfaces A: Physicochemical and Engineering
Aspects, 504:393–399, 2016.

[36]
F. R. Saucedo-Zendejo and E. O. Reséndiz-Flores.
Transient heat transfer and solidification modelling in direct-chill
casting using a generalized finite differences method.
Journal of Mining and Metallurgy, Section B: Metallurgy,
55(1):47–54, 2019.

[37]
B. Seibold.
M-Matrices in Meshless Finite Difference Methods.
PhD thesis, University of Kaiserslautern, Germany, 2006.

[38]
G. Son, V.K. Dhir, and N. Ramanujapu
Dynamics and Heat Transfer Associated With a Single Bubble During Nucleate Boiling on a Horizontal Surface
Journal of Heat Transfer, 121(3):623–631, 1999.

[39]
P. Suchde.
Conservation and Accuracy in Meshfree Generalized Finite
Difference Methods.
PhD thesis, University of Kaiserslautern, Germany, 2018.

[40]
P. Suchde, T. Jacquemin, and O. Davydov
Point Cloud Generation for Meshfree Methods: An Overview
Archives of Computational Methods in Engineering, 30(2):889–915, 2023.

[41]
P. Suchde and J. Kuhnert.
Point cloud movement for fully lagrangian meshfree methods.
Journal of Computational and Applied Mathematics, 340:89–100,

[42]
P. Suchde and J. Kuhnert.
A fully lagrangian meshfree framework for PDEs on evolving
Journal of Computational Physics, 395:38 – 59, 2019.

[43]
P. Suchde and J. Kuhnert.
A meshfree generalized finite difference method for surface PDEs.
Computers & Mathematics with Applications, 78(8):2789 – 2805,

[44]
P. Suchde, J. Kuhnert, S. Schröder, and A. Klar.
A flux conserving meshfree method for conservation laws.
International Journal for Numerical Methods in Engineering,
112(3):238–256, 2017.

[45]
P. Suchde, J. Kuhnert, and S. Tiwari.
On meshfree GFDM solvers for the incompressible Navier–Stokes
Computers & Fluids, 165:1 – 12, 2018.

[46]
Y. Sun and C. Beckermann.
Diffuse interface modeling of two-phase flows based on averaging:
mass and momentum equations.
Physica D: Nonlinear Phenomena, 198(3-4):281–308, 2004.

[47]
A. Tramecon and J. Kuhnert.
Simulation of advanced folded airbags with VPS-PAMCRASH/FPM:
Development and validation of turbulent flow numerical simulation
techniques applied to curtain bag deployments.
In SAE Technical Paper, Warrendale, PA, USA, 2013. SAE

[48]
G. Tryggvason and J. Lu.
Direct numerical simulations of flows with phase change.
Procedia IUTAM, 15:2–13, 2015.
IUTAM Symposium on Multiphase Flows with Phase Change: Challenges and

[49]
E. Uhlmann, E. Barth, A. Quellhorst, T. Seifarth, J. Kuhnert, K. Nabbout, and
M. Sommerfeld.
Fully coupled wet cylindrical turning simulation using the
Procedia CIRP, 102:43–48, 1 2021.

[50]
E. Uhlmann, E. Barth, T. Seifarth, M. Höchel, J. Kuhnert, and A. Eisenträger.
Simulation of metal cutting with cutting fluid using the
Procedia CIRP, 101:98–101, 1 2021.

[51]
E. Uhlmann, R. Gerstenberger, and J. Kuhnert.
Cutting simulation with the meshfree finite pointset method.
Procedia CIRP, 8, 2013.

[52]
L. Veltmaat, F. Mehrens, H.-J. Endres, J. Kuhnert, and P. Suchde.
Mesh-free simulations of injection molding processes.
Physics of Fluids, 34(3):033102, 2022.

[53]
H. K. Versteeg and W. Malalasekera.
An introduction to computational fluid dynamics: the finite
volume method.
Pearson education, 2007.

[54]
V.R. Voller and C. Prakash
A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems.
International journal of heat and mass transfer, 30(8):1709–1719, 1987.

[55]
C. L. Wang, J. C. Li, M. L. Zhao, M. I. Marqués, C. Aragó, and J. A. Gonzalo.
Monte carlo simulation of first order phase transitions.
Ferroelectrics, 401(1):3–8, 2010.

[56]
P.C. Wayner, Y.K. Kao, and L.V. LaCroix
The interline heat-transfer coefficient of an evaporating wetting film
International Journal of Heat and Mass Transfer, 19(5):487–492, 1976.

[57]
S. W. Welch and J. Wilson.
A volume of fluid based method for fluid flows with phase change.
Journal of Computational Physics, 160(2):662–682, 2000.

[58]
N. B. Wilding.
Computer simulation of fluid phase transitions.
American Journal of Physics, 69(11):1147–1155, 2001.

[59]
X. Xu, C. Cheng, and I. H. Chowdhury.
Molecular dynamics study of phase change mechanisms during
femtosecond laser ablation.
J. Heat Transfer, 126(5):727–734, 2004.

[60]
Y.-C. Yoon and J.-H. Song.
Extended particle difference method for weak and strong discontinuity
problems: part i. derivation of the extended particle derivative
approximation for the representation of weak and strong discontinuities.
Computational Mechanics, 53(6):1087–1103, 2014.

[61]
Y.-C. Yoon and J.-H. Song.
Extended particle difference method for weak and strong discontinuity
problems: part ii. formulations and applications for various interfacial
singularity problems.
Computational Mechanics, 53(6):1105–1128, 2014.