# An Experimental Evaluation of Transformer-based Language Models in the
Biomedical Domain

Paul Grouchy
Untether AI
Shobhit Jain∗
Manulife
Michael Liu∗
Tealbook
Kuhan Wang∗
CIBC
Max Tian∗
Adeptmind
Nidhi Arora∗
Intact
Hillary Ngai∗
University of Toronto
Faiza Khan Khattak
Manulife
Elham Dolatabadi
Vector Institute
Sedef Akinli Kocak
Vector Institute
Equal ContributionCorresponding Author<EMAIL_ADDRESS>Corresponding Author<EMAIL_ADDRESS>Author:
<EMAIL_ADDRESS>
###### Abstract

With the growing amount of text in health data, there have been rapid advances
in large pre-trained models that can be applied to a wide variety of
biomedical tasks with minimal task-specific modifications. Emphasizing the
cost of these models, which renders technical replication challenging, this
paper summarizes experiments conducted in replicating BioBERT and further pre-
training and careful fine-tuning in the biomedical domain. We also investigate
the effectiveness of domain-specific and domain-agnostic pre-trained models
across downstream biomedical NLP tasks. Our finding confirms that pre-trained
models can be impactful in some downstream NLP tasks (QA and NER) in the
biomedical domain; however, this improvement may not justify the high cost of
domain-specific pre-training.

_K_ eywords BioBERT $\cdot$ Biomedical data $\cdot$ NLP downstream tasks
$\cdot$ Transformer-based models

##  1 Introduction

There have been increased exploration and extension of transformer-based deep
learning models for NLP [1, 2, 3, 4] recently. Using transformer-based models,
one can pre-train a deep learning model on large datasets and then easily
fine-tune it to adapt to downstream NLP tasks. There are three factors that
affect the performance of these models: (a) the size of the dataset, (b) the
availability of computational resources, and (c) the expressiveness of the
model architecture [2, 1]. Because of these factors, the cost and complexity
of developing pre-trained models are rising quickly and limit the capability
of reproducing results when sufficient resources are not available.

These language models are trained on text corpora of general domains. For
example, BERT [3], Bidirectional Encoder Representations from Transformers has
been trained on Wikipedia and BooksCorpus. There has also been a rapid growth
in NLP for the biomedical domain [5] and many new methods including
transformer-based methods are being used for different biomedical tasks
involving NLP.

The performance of language models that have been trained on general domains
are not yet fully investigated in more specific domains such as biomedical,
finance or legal. Therefore, it is worth investigating if large amounts of
domain-specific data may help in getting better results, or if similar results
may be acquired by using smaller-sized data.

We therefore focus on biomedical domain in order to answer the following
questions:

  1. 1.

Does domain-specific training improve performance compared to baseline models
trained on domain-agnostic corpora?

  2. 2.

Is it possible to obtain comparable results from a domain-specific BERT model
pre-trained on smaller-sized data? While it is established fact that with
transfer learning, a model can be trained once and be reused for several tasks
but there are a few cases when the model needs to be retrained. For example,
when the data is dynamic and may change over time due to the data-shift; data
can belong to a wide variety of domains, therefore the model may need to be
retrained for the new domain data or even sub-domain; data can be confidential
to an organization hence the model needs to be retrained for that
organization-specific needs. The domain and/or organization specific datasets
can be very small. Moreover, not every organization has extensive
computational resources required to train large models. In such cases it is
helpful to know if a small domain-specific data can be used to get comparable
results.

The rest of this paper is organized as follows. We review the existing studies
related to our work in Section 2. We explain our pre-training experiments and
the results in Section 3.1 and fine-tuning experiments and results in Section
3.2. We complete our paper with discussing and conclusion in Section 4 where
we discuss our results and answering our questions.

##  2 Related Work

Large-scale pre-trained language models such as BERT[3], GPT-2 [6], RoBERTa
[2] and GPT-3 [7] have shown to outperform state-of-the-art performance in
many NLP tasks such as Named Entity Recognition (NER) and Question Answering
(QA). Moreover, several studies have used transfer learning and fine-tuning of
these models on English NLP tasks (e.g., [8, 9]). These language models are
trained on text corpora of general domains but recently there has been a trend
of training language models on the domain-specific data. For example,
financial version of BERT was introduced by Araci [10] where he only studied
sentiment classification task. Ma et al [11] fine-tuned BERT on legal
documents coming from their proprietary corpus. BioBERT [4] was introduced as
an extension of BERT that is further pre-trained on the domain-specific
biomedical corpora including PubMed and PubMed Central (PMC).

In specialized domains like biomedical, recent work has shown that using
domain specific data can provide improvement over general-domain language
models [12]. In this regard, Wang et al. [13] showed that word-embeddings
trained on biomedical corpora captured the semantics of medical terms better
than those trained on general domain corpora, but may not generalize well to
downstream biomedical NLP tasks such as biomedical information retrieval. Zhao
et al [14] showed that word2vec [15] trained on a smaller and in-domain
medical data resulted in better performance than the word2vec trained on a
large and general domain dataset. Also, [16] found that performance decreases
after 4 million distinct words of training data based on experiments with
medical data from PubMed abstracts111https://www.ncbi.nlm.nih.gov/pubmed/. In
a recent study, Gu et al [12] introduced PubMedBERT. They pre-trained the BERT
from scratch with PubMed articles and a customized vocabulary (constructed
from the PubMed articles). This study indicates that a proper vocabulary helps
the performance of downstream tasks in specific domains. However, training the
model from scratch is extremely expensive in terms of data and computation.
Researchers have specifically built adaptations of BERT that attempt to
address different domain related problems but the most effective pretraining
process remains an open research problem [12]. We replicated some of the
BioBERT original study results, and better tune the training of BioBERT for
better understanding domain specific training.

##  3 Methods and Experimental Results

We set BERTBASE as our baseline model. We started with BERTBASE, then pre-
trained it on the PubMed abstracts data (BERTbase+PM) and leveraged for
evaluation of downstream tasks: Name-entity recognition (NER, relational
extraction (RE) (Table 1) and question answering (QA) (Table 2). This allows
us to conduct a fair comparison between domain-specific pre-training and fine-
tuning. PyTorch implementation of
BERT222https://github.com/huggingface/pytorch-transformers [17] was leveraged
and the replication experiments were conducted based on the work by McDermott
et. al [18].

###  3.1 Pre-training language representations in the biomedical domain

The PubMed corpora was used for pre-training consists of paper abstract from
millions of samples of biomedical text. While the original BioBERT study
considers combined pre-training on PubMed, PMC, and Pubmed+PMC together, our
model was pre-trained only on PubMed [19] to check the performance based on
smaller data and meet the shared computing resources.

The PubMed data was processed into a format amenable for pre-training. The raw
data consisted of approximately 200 million sentences in 30 GBs. The raw
sentence data was batch processed into 111 chunks of ready to consume input
data for BERT pre-training333This implies that the sampling of the second
sentence for the Next Sentence Prediction (NSP) training task is not from the
entire corpus but from within its respective chunk. As each chunk is still
close to 2 million separate sentences, this was judged to be an acceptable
compromise.. Technically, the NSP task was lower-cased sentences with a
maximum length of 512 and masked at the sub-token level analogous to the
original BERTBASE(uncased). The original BERT models were trained on data with
maximum sequence lengths of 128 for the initial 90% of the training and 512
for the remainder. The full training loss is shown in Figure 1.

Figure 1: Training loss. The effect of chunking the data into shards is seen
in the fluctuation of the training loss. The vertical line indicates
separation between batch size 28 and 14 in the training paradigm. The top
(bottom) x-axis indicate calendar date (global step).

Generally, a longer sequence length is preferable if the corpora tends to have
longer passages but the amount of data and time to train is also increased. In
this experiment, the training was maintained at the maximum length of 512
tokens throughout.
The pre-training experiment was designed to accommodate some realistic
computing resource limitations. In a time-shared computing cluster GPU
resources were allocated on limited priority basis by Slurm Workload Manager.
In order to facilitate stable training and to accommodate a hardware
environment where resources may be reallocated to higher priority users at any
time and without warning, the following training features were implemented:
(1) storing model weights and gradients at regular time intervals during
training; (2) querying and automating job submission from system task
scheduler; (3) automatically restoring training from data chunk and step as
GPU resources became available.

###  3.2 Fine-tuning language representations in the biomedical domain

Following previous work, our domain specific pre-trained model (BERTbase+PM)
as well as BERTBASE(uncased) were evaluated on common biomedical NLP tasks as
described below.

####  3.2.1 Named-entity recognition (NER)

The biomedical NER task involves the extraction of biomedical named entities
such as genes, proteins, diseases and species from unstructured biomedical
text. This is a challenging task because of the unique characteristics of
biomedical named entities such as sharing of head nouns, spelling forms per
entity, ambiguous abbreviations, descriptive naming convention, and nested
names [20]. In the original BioBERT study nine different biomedical NER
datasets were tested and their best model (BERTBASE+PM+PMC) was able to
outperform state-of-the-art on the majority. In this experiment, three NER
datasets including NCBI Disease [21], JNLPBA [22], and Species-800 [23] were
used. The model was trained on 50 tokens long sentences with a learning rate
of $\mathit{5x10^{-5}}$. These experiments were repeated several times and F1
scores corresponding to the epoch with the minimum validation loss are
reported in Table 1. The standard deviation in F1 scores across different runs
make definitive comparisons with the results reported in BioBERT study, but a
surprising outcome is that the scores for the BERTBASE starting point models
are higher on average than the ones fine-tuned with our BERTbase+PM. This runs
counter to the notion that pre-training on biomedical domains corpora improves
downstream NER.

####  3.2.2 Relation Extraction (RE)

The GAD dataset [24] was used for the RE experiment. The goal of RE on the GAD
dataset is to correctly classify whether a gene and disease are related, based
on a given sentence. Two models were fine-tuned on GAD, one starting with
BERTBASE and the other starting with BERTbase+PM. A learning rate of
$\mathit{5x10^{-5}}$ was used for all runs. The best validation F1 score over
several runs for each fold is reported. Our BERTbase+PM outperformed other
BERT models but the differences between models are not significant (see Table
1).

| BioBERT [4] | Our Experiments
---|---|---
Tasks | Datasets | Metrics | BERTbase | +PubMed | BERTbase | BERTbase+PM
NER | NCBI Disease | F1 | 85.63 | 89.71 | 84.08$\pm$2.07 | 80.33$\pm$2.40
JNLPBA | F1 | 74.94 | 77.49 | 76.43$\pm$2.29 | 75.16$\pm$3.63
Species 800 | F1 | 71.63 | 74.06 | 78.38$\pm$7.87 | 74.59$\pm$5.28
RE | GAD | F1 | 79.30 | 79.83 | 79.60 | 81.50
Table 1: NER and RE Performance Results. Our experiments including fine-tuning
of the BERTBASE(uncased) model and the BERTBASE(uncased) model further pre-
trained on PubMed abstracts (+PM). For comparison, the results of BioBERT
study using BERTBASE (i.e. Wiki+Books) and the BioBERT version of the PubMed
trained model (+PubMed) are also shown. Best scores are shown in bold.

####  3.2.3 Question Answering (QA)

The goal of QA tasks is to automatically find the answer to a question posed
in human language, usually from a context paragraph. In this study, we
explored the performance of BERT and BioBERT on two common biomedical QA tasks
including BioASQ [25] and PubmedQA [26]. For this task, three versions of BERT
model (BERTlarge, BERTbase, and our BERTbase+PM) were fine-tuned on the QA
datasets.

BioASQ: As was suggested in the BioBERT study [4], all BERT models were
initially fine-tuned on the SQuAD [27] dataset (with intermediate
evaluations), and then on the BioASQ training set before finally evaluating on
the BioASQ test sets (Table 2). It has been also shown in [28] that pre-
training BERT on SQuAD 1.1 generated better results when fine-tuned on BioASQ
in comparison to the model pre-trained on SQuAD 2.0. For tuning the model on
SQuAD 1.1 and SQuAD 2.0, we took inspiration from the training schemes
outlined in [2] and [29] to adjust the hyperparameters, namely the learning
rate, $\beta_{2}$, learning rate schedule, batch size, and number of training
epochs. We found that training with a cosine learning rate schedule with no
warm-up steps with $\beta_{2}=0.98$ consistently resulted in the best
performance, and thus carried over the heuristics generated from fine-tuning
on SQuAD to further fine-tune our models on BioASQ. Training on domain-
specific data: Overall, there is some evidence (Table 2) that pre-training on
biomedical domain corpus improves performance on the downstream BioASQ QA
task. However, the improvement is not so large as to be entirely convincing
and carefully fine-tuned BERT models (such as BERTlarge in our case) can
perform comparably to BioBERT.
Zero-shot setting: We also conducted an experiment in order to evaluate the
performance of the BERT models on BioASQ datasets in a zero-shot setting where
BERTBASE was fine-tuned on SQuAD and _not_ on BioASQ training data. We
averaged the scores over the 5 test sets within Strict Accuracy (S), Lenient
Accuracy (L), and Mean Reciprocal Rank (M) metrics and the results were 31.37,
46.8, and 37.16, respectively. As we can see the zero-shot evaluation results
on test sets are worse than the results indicated in Table 2 which emphasizes
the effects of fine-tuning on BioASQ.
Number of epochs: Additionally, we found that fine-tuning BERTlarge model on
BioASQ for more epochs than the typically recommended 1-3 epochs resulted in
much better results. Our best results across the three evaluation metrics on
BioASQ came from fine-tuning BERTlarge for 20 epochs and BERTbase for 4 epochs
with initial learning rate of $5x10^{-6}$.

| BioBERT [4] | Our Experiments
---|---|---
Datasets | Metrics | BERTbase | +PubMed | BERTlarge∗ | BERTbase | BERTbase+PM
BioASQ (4b) | S | 27.33 | 27.95 | 31.8 | 31.2 | 31.54
L | 44.72 | 44.10 | 51.8 | 44.58 | 48.36
M | 33.77 | 34.72 | 40.0 | 36.25 | 38.39
BioASQ (5b) | S | 39.33 | 46.00 | 43.0 | 41.61 | 46.05
L | 52.67 | 60.00 | 55.8 | 56.8 | 57.94
M | 44.27 | 51.64 | 48.2 | 47.86 | 50.54
BioASQ (6b) | S | 33.54 | 42.86 | 35.8 | 36.55 | 37.57
L | 51.55 | 57.77 | 54.4 | 59.3 | 58.87
M | 40.88 | 48.43 | 43.4 | 47.59 | 46.5
Average | S | 33.4 | 38.93 | 36.86 | 37.48 | 38.39
over | L | 49.65 | 52.68 | 55.12 | 53.56 | 55.06
4b,5b, and 6b | M | 39.6 | 44.93 | 43.86 | 43.9 | 45.14
PubMedQA | K-Fold Acc | — | 57.28 | 56.52 | 55.20 | 56.20
K-Fold F1 | — | 28.70 | 26.14 | 23.71 | 23.98
Table 2:  The performance results of three versions of BERT model (BERTlarge,
BERTbase and BERTbase+PM) on two distinct QA datasets. For comparison, the
results of BioBERT study using BERTBASE (i.e. Wiki+Books) and the BioBERT
version of the PubMed trained model (+PubMed) are also shown. Since, we
already explored latest versions of BERTlarge and BERTbase, we didn’t examine
the BERTbase version used in the BioBERT study. Best scores are shown in bold.
*Used BERTlarge as other versions were not providing good results.

PubMedQA: All three versions of the BERT models were fine-tuned on PQA-L i.e.,
1k expert-annotated generated QA instances of the PubMedQA dataset with the
results shown in Table 2. Since PubMedQA hasn’t been experimented in the
BioBERT study, in order to make a fair comparison between our study and the
BioBERT study [4], we ran an additional experiment where we fine-tuned PubMed
trained model (+PubMed) from the BioBERT study on the PubMedQA. 10-fold cross-
validation was performed with only 450 training instances in each fold of
validation. As seen in Table 2, the BioBERT version of the PubMed trained
model (+PubMed) has the highest accuracy 57.28 and F1 Score 28.27. This
suggests that pre-training on biomedical corpus improves the performance of
the downstream PubMedQA task. However, the performance improvement compared to
BERTlarge, BERTbase, and BERTbase+PM is minimal.

####  3.2.4 Text Summarization

Text summarization refers to automatic generation of summary of a given text.
Extractive summarization is done by extracting the most important sentences
from the document that summarize the whole document. Abstractive summarization
refers to condensing the document into shorter versions while preserving its
meaning [30].

| CNN/DailyMail | XSum | BioASQ 7b (Our experiments)
---|---|---|---
$\\#$ docs (train/val/test) | 196,961/ /12,148/10,397 | 204,045/11,332/11,334 | 22,462/4,814/4,800
| ROUGE-1 | ROUGE-2 | ROUGE-3 | ROUGE-1 | ROUGE-2 | ROUGE-3 | ROUGE-1 | ROUGE-2 | ROUGE-3
Abs | 41.72 | 19.39 | 38.76 | 38.81 | 16.33 | 31.15 | 38.20 | 27.12 | 34.15
ExtAbs | 42.13 | 19.6 | 39.18 | 38.76 | 16.50 | 31.27 | 38.89 | 29.29 | 36.01
Ext | 43.25 | 20.24 | 39.63 | $\dagger$ | $\dagger$ | $\dagger$ | 33.30 | 26.50 | 32.20
Table 3:  Text summarization F1 score. $\dagger$ Results were poor hence were
not reported by the authors [31]

BERTSum [32, 31] is a simple variant of BERT, for text summarization that
produces abstractive and extractive summaries. BERTSumAbs is based on an
encoder-decoder architecture where the encoder is pre-trained whereas decoder
is a randomly initialized transformer which is trained from scratch. There are
three versions of BERTSum. BERTSumAbs produces abstractive summary, BERTSumExt
outputs extractive summary, and BERTSumExtAbs also produces abstractive
summary but is based on two-stage approach where the encoder is fine-tuned
twice, first with an extractive objective followed by an abstractive one.

We explored the performance of these methods on BioASQ
7b444http://bioasq.org/participate/challenges_year_7 dataset consisting of
links to research papers as well as summaries are embedded in the json file
itself. BERT was trained from scratch on BioASQ 7b. The results are shared in
Table 3. For CNN/DailyMail dataset all three methods outperform on BioASQ
dataset according to ROUGE-2 score, while produce comparable results in all
other cases. For XSum dataset [33], we have comparable results according to
ROUGE-1, while BioASQ 7b performs better according to ROUGE-2 and ROUGE-3 (for
Ext summary results were not reported by the authors [31]). This maybe due to
the fact that BioASQ dataset is much smaller in size as compared to the other
two datasets. This also shows that smaller-sized dataset can be used for text-
summarization.

##  4 Discussion and Conclusion

In this study, we present our findings for experiments conducted in
conjunction with further pre-training of BERT model in the biomedical domain
as well as evaluating both the domain specific and domain agnostic pre-trained
models across downstream biomedical NLP tasks. Our experiments also included
considering data and optimization related factors. Further pre-training was
conducted on PubMed abstracts only to evaluate the performance of the models
trained on smaller-size datasets and also the effects of learning rate was
carefully evaluated for fine-tuning tasks.

(1) Does domain-specific training improves performance? This study confirms
that unsupervised pre-training in general could improve the performance on
fine-tuning tasks. However, the effectiveness of domain specific pre-training
as a way of further improving the performance of supervised downstream tasks
does not significantly outperform the effectiveness of domain agnostic pre-
trained models considering the high cost of domain-specific pre-training which
makes it challenging for most of researchers and NLP developers. In the
biomedical domain, however, this conclusion may not be wholly substantiated
owing to a lack of consistent evidence particularly in downstream NER, QA, and
Text Summarization tasks. For the SQuAD task, it has been shown that fine-
tuning results depend on the size and duration of training [34] [35]. However,
given the small size of the biomedical QA (BioASQ) datasets, it is not
possible to run the same experiments, as the current pre-trained models easily
overfit. Therefore, it should be beneficial for the biomedical community to
curate and expand the magnitude of benchmark datasets. Understandably it is
difficult and expensive, but with the ever increasing size of the deep-
learning models and significant advances in the development of pre-training
language representations, it is necessary in order to facilitate reproducible
research for health.

(2) Is it possible to obtain comparable results using BERT model pre-trained
on smaller-sized data? We present results of the experiments that were
conducted using the models pre-trained on the PubMed abstracts only. These
results were comparable to the results produced by the model trained on the
PubMed, PMC, and Pubmed+PMC together (please check [4]). Although it requires
further investigation but empirically it shows that small-sized datasets may
be used as a surrogate. Using small-sized dataset can be especially useful
when the models need to retrained instead of using pre-trained publicly
available models due to (a) data-shift, (b) wide variety of data-domains, (c)
confidential data not publicly available to train the model on, (d) small-size
of the domain specific data available, and (e) a lack of computing resources.

We position these experiments as complementary to existing literature on the
applications of transformer-based models for biomedical NLP. Our results
provide some evidence for the validity and the limitations of existing
language representations for pre-training $\&$ fine-tuning in biomedical
domain.

We would like to expand our experiments on other variants of BERT, with more
domain specific datasets and a variety of downstream tasks. This would also
allow us to use our experiments for other medical
corpora555http://www.nactem.ac.uk/resources.php for external validation and
generalization of the results.

## Acknowledgements

We want to thank Vector Institute industry sponsors, researchers and technical
staff who participated in the Vector Institute’s NLP Project
(https://vectorinstitute.ai/wp-content/uploads/2020/12/nlp-report-final.pdf).

## References

  * [1] Zhilin Yang, Zihang Dai, Yiming Yang, Jaime Carbonell, Ruslan Salakhutdinov, and Quoc V Le.  Xlnet: Generalized autoregressive pretraining for language understanding.  arXiv preprint arXiv:1906.08237, 2019.
  * [2] Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov.  Roberta: A robustly optimized bert pretraining approach.  arXiv preprint arXiv:1907.11692, 2019.
  * [3] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  arXiv preprint arXiv:1810.04805, 2018.
  * [4] Jinhyuk Lee, Wonjin Yoon, Sungdong Kim, Donghyeon Kim, Sunkyu Kim, Chan Ho So, and Jaewoo Kang.  Biobert: a pre-trained biomedical language representation model for biomedical text mining.  Bioinformatics, 36(4):1234–1240, 2020.
  * [5] Faiza Khan Khattak, Serena Jeblee, Chloé Pou-Prom, Mohamed Abdalla, Christopher Meaney, and Frank Rudzicz.  A survey of word embeddings for clinical text.  Journal of Biomedical Informatics: X, page 100057, 2019.
  * [6] Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever.  Language models are unsupervised multitask learners.  OpenAI blog, 1(8):9, 2019.
  * [7] Tom B Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al.  Language models are few-shot learners.  arXiv preprint arXiv:2005.14165, 2020.
  * [8] Matthew E Peters, Sebastian Ruder, and Noah A Smith.  To tune or not to tune? adapting pretrained representations to diverse tasks.  arXiv preprint arXiv:1903.05987, 2019.
  * [9] Wei Yang, Yuqing Xie, Aileen Lin, Xingyu Li, Luchen Tan, Kun Xiong, Ming Li, and Jimmy Lin.  End-to-end open-domain question answering with bertserini.  arXiv preprint arXiv:1902.01718, 2019.
  * [10] Dogu Araci.  Finbert: Financial sentiment analysis with pre-trained language models.  arXiv preprint arXiv:1908.10063, 2019.
  * [11] Xiaofei Ma, Peng Xu, Zhiguo Wang, Ramesh Nallapati, and Bing Xiang.  Domain adaptation with bert-based domain classification and data selection.  In Proceedings of the 2nd Workshop on Deep Learning Approaches for Low-Resource NLP (DeepLo 2019), pages 76–83, 2019.
  * [12] Yu Gu, Robert Tinn, Hao Cheng, Michael Lucas, Naoto Usuyama, Xiaodong Liu, Tristan Naumann, Jianfeng Gao, and Hoifung Poon.  Domain-specific language model pretraining for biomedical natural language processing.  arXiv preprint arXiv:2007.15779v3, 2020.
  * [13] Yanshan Wang, Sijia Liu, Naveed Afzal, Majid Rastegar-Mojarad, Liwei Wang, Feichen Shen, Paul Kingsbury, and Hongfang Liu.  A comparison of word embeddings for the biomedical natural language processing.  Journal of biomedical informatics, 87:12–20, 2018.
  * [14] Mengnan Zhao, Aaron J Masino, and Christopher C Yang.  A framework for developing and evaluating word embeddings of drug-named entity.  In Proceedings of the BioNLP 2018 workshop, pages 156–160, 2018\.
  * [15] Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean.  Efficient estimation of word representations in vector space.  arXiv preprint arXiv:1301.3781, 2013.
  * [16] Yongjun Zhu, Erjia Yan, and Fei Wang.  Semantic relatedness and similarity of biomedical terms: examining the effects of recency, size, and section of biomedical publications on the performance of word2vec.  BMC medical informatics and decision making, 17(1):95, 2017.
  * [17] Hugging Face.  A library of state-of-the-art pretrained models for natural language processing (nlp).  https://github.com/huggingface/pytorch-transformers, 2019.
  * [18] Matthew McDermott, Shirly Wang, Nikki Marinsek, Rajesh Ranganath, Marzyeh Ghassemi, and Luca Foschini.  Reproducibility in machine learning for health.  arXiv preprint arXiv:1907.01463, 2019.
  * [19] US National Library of Medicine National Institute of Health.  Pubmed, pubmed central, (accessed June, 2019).  https://www.ncbi.nlm.nih.gov/pubmed/.
  * [20] Amy Neustein, S Sagar Imambi, Mário Rodrigues, António Teixeira, and Liliana Ferreira.  Application of text mining to biomedical knowledge extraction: analyzing clinical narratives and medical literature.  Text mining of web-based medical content. De Gruyter, Berlin, 50, 2014.
  * [21] Rezarta Islamaj Doğan, Robert Leaman, and Zhiyong Lu.  Ncbi disease corpus: a resource for disease name recognition and concept normalization.  Journal of biomedical informatics, 47:1–10, 2014.
  * [22] Jin-Dong Kim, Tomoko Ohta, Yoshimasa Tsuruoka, Yuka Tateisi, and Nigel Collier.  Introduction to the bio-entity recognition task at jnlpba.  In Proceedings of the international joint workshop on natural language processing in biomedicine and its applications, pages 70–75. Citeseer, 2004.
  * [23] Evangelos Pafilis, Sune P Frankild, Lucia Fanini, Sarah Faulwetter, Christina Pavloudi, Aikaterini Vasileiadou, Christos Arvanitidis, and Lars Juhl Jensen.  The species and organisms resources for fast and accurate identification of taxonomic names in text.  PLoS One, 8(6):e65390, 2013.
  * [24] Kevin G Becker, Kathleen C Barnes, Tiffani J Bright, and S Alex Wang.  The genetic association database.  Nature genetics, 36(5):431, 2004.
  * [25] George Tsatsaronis, Georgios Balikas, Prodromos Malakasiotis, Ioannis Partalas, Matthias Zschunke, Michael R Alvers, Dirk Weissenborn, Anastasia Krithara, Sergios Petridis, Dimitris Polychronopoulos, et al.  An overview of the bioasq large-scale biomedical semantic indexing and question answering competition.  BMC bioinformatics, 16(1):138, 2015.
  * [26] Qiao Jin, Bhuwan Dhingra, Zhengping Liu, William W. Cohen, and Xinghua Lu.  Pubmedqa: A dataset for biomedical research question answering.  arXiv preprint arXiv:1909.06146v1, 2019.
  * [27] Pranav Rajpurkar, Jian Zhang, Konstantin Lopyrev, and Percy Liang.  Squad: 100,000+ questions for machine comprehension of text.  arXiv preprint arXiv:1606.05250, 2016.
  * [28] Wonjin Yoon, Jinhyuk Lee, Donghyeon Kim, Minbyul Jeong, and Jaewoo Kang.  Pre-trained language model for biomedical question answering.  arXiv preprint arXiv:1909.08229, 2019.
  * [29] Zhenzhong Lan, Mingda Chen, Sebastian Goodman, Kevin Gimpel, Piyush Sharma, and Radu Soricut.  Albert: A lite bert for self-supervised learning of language representations.  arXiv preprint arXiv:1909.11942, 2019.
  * [30] Vishal Gupta and Gurpreet Singh Lehal.  A survey of text summarization extractive techniques.  Journal of emerging technologies in web intelligence, 2(3):258–268, 2010.
  * [31] Yang Liu and Mirella Lapata.  Text summarization with pretrained encoders.  arXiv preprint arXiv:1908.08345, 2019.
  * [32] Yang Liu.  Fine-tune bert for extractive summarization.  arXiv preprint arXiv:1903.10318, 2019.
  * [33] Shashi Narayan, Shay B Cohen, and Mirella Lapata.  Don’t give me the details, just the summary! topic-aware convolutional neural networks for extreme summarization.  arXiv preprint arXiv:1808.08745, 2018.
  * [34] Alon Talmor and Jonathan Berant.  Multiqa: An empirical investigation of generalization and transfer in reading comprehension.  CoRR, abs/1905.13453, 2019.
  * [35] Georg Wiese, Dirk Weissenborn, and Mariana L. Neves.  Neural domain adaptation for biomedical question answering.  CoRR, abs/1706.03610, 2017.


# Block Spectral Stresses (BSS) estimation for shock-capturing and turbulent
modeling

Matteo Ruggeri
Department of Aeronautics & Astronautics Engineering
Purdue University
<EMAIL_ADDRESS>
&Victor C. B. Sousa
Department of Mechanical Engineering
Purdue University
&Carlo Scalo
Department of Mechanical Engineering
Purdue University
Department of Aeronautics & Astronautics Engineering by courtesy

###### Abstract

A new combined sub-filter scale turbulence and shock-capturing model is
developed for high-order finite volume numerics, extending previous work [1]
to unstructured solvers. Block Spectral Stresses (BSS) method relies on the
spectra of the velocity gradients to estimate the subfilter scale stresses,
heat-flux, and pressure-work based on the resolved field. The method is able
to capture shocks with numerical order up to 25 and in a shock-vortex
interaction simulation is able to capture the shock and not interfere with the
vortex structure. In turbulence calculations the new method is compared with
Smagorinsky, dynamic Smagorinsky, and Vreman methods adapted to a block
spectral code. In the simulations of homogeneous isotropic turbulence, the new
model is worse than the others when on coarse meshes and better on finer ones.
Instead, for supersonic and hypersonic channel flow the case is the opposite
because as expected the sub-filter terms are mostly depend on the numerical
order and not the mesh resolution.

_Keywords_ Turbulence modeling $\cdot$ Shock-capturing $\cdot$ Flux
reconstruction $\cdot$ Block spectral method $\cdot$ Wall bounded flow

##  1 Introduction

The last years have seen the proliferation of pseudospectral methods which
merge the benefits of unstructured mesh to the high-order numerics which allow
a better resolution than low order methods with a smaller number of degrees of
freedom, given their faster convergence rate [2, 3, 4, 5, 6]. Indeed,
pseudospectral methods show a convergence rate closer to spectral methods than
finite difference ones. But at the same time they shows the same limitation
[7, 8, 9] of the high-order method when it has to work with steep gradients as
in the case of shocks in compressible flow. In addition to shocks in gas
dynamic flow another phenomenon that can cause numerical instability and leads
to the blow-up of the simulation if not properly modeled is turbulence. Both
turbulence and shock formation have a similar energy cascade from large to
small scales because of nonlinear interaction [10, 11]. Differently from a
shock that has an infinite spectra, turbulence has a finite spectra therefore
standard shock-capturing methods usually tend to remove a part of the spectra
of the turbulence and therefore reduce the quality of the simulation. Instead,
Eddy viscosity models–which allow to resolve turbulence on a coarser mesh than
Direct Numerical Simulation quality (DNS) thanks to the modeling of small
scales turbulence–are not able to fully dissipate the shock. This work aims to
address these two problems with a unified model that is able to capture the
shock and the same time does not remove a broad spectra of the turbulence from
the simulation.

Block spectral methods can be implemented with three different approaches: the
discontinuous Galerkin (DG) [12, 13], the staggered grid approach (or spectral
differences) [14], and the flux reconstruction (FR) method [4]. The flux
reconstruction (FR) method is the latest to be developed and it relies on
correction functions that are used to interpolate back to the solution points
the fluxes from their updated value at the interfaces, this makes them
continuous among elements. Hyunh [4] pointed out that the flux reconstruction
method is depend to the functions that are use for correction and hence it can
be made more or less smooth, stable, or accurate than other numerical methods.
Asthana et al. [15] proved that even when the computational domain is refined
the instabilities caused by aliasing errors related to the collocation
projection of the flux at the jump conditions are not completely reduced. In
the same paper, they mathematically proved that if an artificial viscosity is
added to flux reconstruction method, it can be stabilized with all numerical
order after a certain grid refinement. They also showed how this allows to run
1D shock simulation up to numerical order 120 and even without the sensor that
detects the shock the rate of convergence is not much different from when is
active (the same order of accuracy was achieved by Sousa and Scalo [1] for 1D
Sod shock tube).

Tonnicello et al. [16] addressed the problem of preserving entropy in a shock
when applying the shock-capturing method in a high-order spectral difference
numerics code. The entropy profile across a shock is not monotonic (also shown
by Gupta and Scalo [17]) but it presents a peak inside the propagation wave
this result has been shown by Morduchow and Libby [18] and later validated by
other works [19, 20, 21]. The problem is that most of the shock-capturing
methods are not able to capture this phenomenon and make the entropy monotonic
across the shock. They compare the classical Laplacian viscosity to a physical
artificial viscosity method, showing that the latter is able to have a better
prediction of the entropy even when the Laplacian approach is able to resolve
the shock well. The main drawback of this approach is that it also removes the
spectra of the acoustic waves, which in the case of supersonic and moreover in
hypersonic flow have great importance because phenomena like second mode waves
lead the transition to turbulence in hypersonic flow. The authors pointed out
that the suppression of acoustic waves can be caused by the discontinuity-
based sensor use to the detect and shock and added that a divergence-based
sensor may have better performance. Nevertheless, the authors presented works
[22, 23] that use a divergence-based sensor and they still have the removal of
acoustic waves when compare to DNS simulations.

The current work aims to address this challenge in developing a combined SFS
turbulence shock-capturing model that is able to capture the shock, while
retaining the spectrum of turbulence. Section 2 describes the code used for
this investigation and the numerics used to solve the compressible Navier-
Stokes equations. In the second part of the section the Block Spectral
Stresses (BSS) model is introduced together with a description of the
Smagorinsky and dynamic Smagorinsky models implemented in a block spectral
code and later use for comparison with BSS. Section 3 presents how the model
performs in 1D and 2D flow with shocks. Then the model is tested to act as an
LES model on a Taylor-Green Vortex (TGV) flow and turbulent wall-bounded flow
at supersonic and hypersonic flow speed in section 4.

##  2 Models description

The code used to perform the simulations is H3AMR [24, 25, 26, 27] (HySonic,
High-Order, Hybrid Adaptive Mesh Refinement developed by HySonic Technology,
LLC) an unstructured block-spectral code [28] for compressible flow based on
flux reconstruction [4] to compute the gradients. In the current section a
description of the code 2.1 and of the LES models as implemented in the code
are provided: Smagorinsky [29] (Section 2.2.1), Dynamic Smagorinsky [30]
(Section 2.2.2), Vreman’s model [31] (Section 2.2.3), and Block Spectral
Stresses (BSS) (Section 2.2.4).

###  2.1 Code description (H3AMR)

The vector of conserved quantities reads:

$\mathbf{Q}=\begin{bmatrix}\rho\\\ \rho u_{1}\\\ \rho u_{2}\\\ \rho u_{3}\\\
\rho E\end{bmatrix}$ (1)

where $\rho$ is the density, $u_{i}$ the velocity in the $i^{th}$ direction,
and $E=e+\frac{u_{i}u_{i}}{2}$ is the specific total energy
($e=\frac{p}{\rho(\gamma-1)}$ is the specific internal energy in the case of
ideal calorically perfect gas). The flux vectors read:

$\displaystyle\mathbf{F}=$ $\displaystyle\mathbf{G}=$
$\displaystyle\mathbf{H}=$ (2) $\displaystyle\begin{bmatrix}\rho u_{1}\\\ \rho
u_{1}u_{1}+p-\mu\sigma_{11}\\\ \rho u_{1}u_{2}-\mu\sigma_{12}\\\ \rho
u_{1}u_{3}-\mu\sigma_{13}\\\ \rho(\rho e+p)u_{1}-k\frac{\partial T}{\partial
x_{1}}+\\\ \dots-\mu\sigma_{i1}u_{i}\end{bmatrix}$
$\displaystyle\begin{bmatrix}\rho u_{2}\\\ \rho u_{2}u_{1}-\mu\sigma_{21}\\\
\rho u_{2}u_{2}+p-\mu\sigma_{22}\\\ \rho u_{2}u_{3}-\mu\sigma_{23}\\\
\rho(\rho e+p)u_{2}-k\frac{\partial T}{\partial x_{2}}+\\\
\dots-\mu\sigma_{i2}u_{i}\end{bmatrix}$ $\displaystyle\begin{bmatrix}\rho
u_{3}\\\ \rho u_{3}u_{1}-\mu\sigma_{31}\\\ \rho u_{3}u_{2}-\mu\sigma_{32}\\\
\rho u_{3}u_{3}+p-\mu\sigma_{33}\\\ \rho(\rho e+p)u_{3}-k\frac{\partial
T}{\partial x_{3}}+\\\ \dots-\mu\sigma_{i3}u_{i}\end{bmatrix}$

where $\mu$ is the dynamic viscosity, $p$ the pressure, $T$ the temperature,
and $\sigma_{ij}=\frac{\partial u_{i}}{\partial x_{j}}$ the strain-rate
tensor. Given these vectors, the Navier-Stokes equations can be written as:

$\frac{\partial\mathbf{Q}}{\partial t}+\frac{\partial\mathbf{F}}{\partial
x_{1}}+\frac{\partial\mathbf{G}}{\partial
x_{2}}+\frac{\partial\mathbf{H}}{\partial x_{3}}=0$ (3)

and they are still in the physical space. Considering that the numeric is
specifically designed for a computational space the conserved quantities
vector and the flux vectors can be recast in computational space, with
coordinates $\xi_{i}$ through:

$\displaystyle Q=J\mathbf{Q}$ (4) $\displaystyle
F=J\left(\frac{\partial\xi_{1}}{\partial
x_{1}}\mathrm{F}+\frac{\partial\xi_{1}}{\partial
x_{2}}\mathrm{G}+\frac{\partial\xi_{1}}{\partial x_{3}}\mathrm{H}\right)$
$\displaystyle G=J\left(\frac{\partial\xi_{2}}{\partial
x_{1}}\mathrm{F}+\frac{\partial\xi_{2}}{\partial
x_{2}}\mathrm{G}+\frac{\partial\xi_{2}}{\partial x_{3}}\mathrm{H}\right)$
$\displaystyle H=J\left(\frac{\partial\xi_{3}}{\partial
x_{1}}\mathrm{F}+\frac{\partial\xi_{3}}{\partial
x_{2}}\mathrm{G}+\frac{\partial\xi_{3}}{\partial x_{3}}\mathrm{H}\right)$

where the italics characters are considered in computational space,
$\frac{\partial\xi_{i}}{\partial x_{j}}$ is the linear transformation matrix
from the physical $x_{j}$ to computational space $\xi_{i}$ and $J$ the
Jacobian of the transformation matrix. Therefore the Navier-Stokes equations
can be rewritten as:

$\frac{\partial\mathbf{Q}}{\partial t}+\frac{1}{J}\frac{\partial
F}{\partial\xi_{1}}+\frac{1}{J}\frac{\partial
G}{\partial\xi_{2}}+\frac{1}{J}\frac{\partial H}{\partial\xi_{3}}=0$ (5)

which are the equations solved by H3AMR.

Figure 1: Example of 1D elements which are interpolate to the interface (black
lines) and then corrected back to the solution points with their updated
values at the interface (red line).

The computational domain is a mesh divided into several elements, in each
element, the conserved quantities and the fluxes are stored in the solution
points which are the Gauss-Legendre quadrature points (N+1 points given a
polynomial reconstruction of order N). From now on every reference to
numerical order is equal to the number of solution points in each cell and
hence the polynomial order is the numerical order minus one. For the
simplicity of exposure, all the following numerics will be described in a 1D
formulation, but it is enough to use an outer product to extend it to a 3D
domain. Given a set of values of conserved quantities and fluxes in the
solution points $f(\xi_{j})$ their values are interpolated to the interfaces
of the element ($f(-1)$ for the left interface and $f(+1)$ for the right
interface) and then updated to their new value using the Rusanov [32] method:

$\mathrm{F}^{new}=\frac{1}{2}(\mathrm{F}_{R}+\mathrm{F}_{L})-\frac{S^{+}}{2}(\mathrm{Q}_{R}-\mathrm{Q}_{L})$
(6)

where $F^{new}$ is the new value of the flux at the interface, $F_{R}$ and
$F_{L}$ are the initial values of the fluxes at the right and left of the
interface, $Q_{R}$ and $Q_{L}$ are the initial values of the conserved
quantities at the right and left of the interface, and $S^{+}$ is the maximum
between the sum of the speed of sound $a$ and velocity $u$ at the left and
right of the interface. From the updated values the derivatives can be
computed using the flux reconstruction [4] method:

$\frac{\partial
f(\xi_{j})}{\partial\xi}=D_{ij}f(\xi_{j})+\left(f^{un}_{-\frac{1}{2}}-f(-1)\right)g^{\prime}_{L}(\xi_{i})+\left(f^{un}_{+\frac{1}{2}}-f(+1)\right)g^{\prime}_{R}(\xi_{i})$
(7)

where $f^{un}_{-\frac{1}{2}}$ and $f^{un}_{+\frac{1}{2}}$ are the updated
values at the left and right interface of the element, $D_{ij}$ is the matrix
representation of the first-order derivative, and $g^{\prime}_{L}(\xi_{i})$
and $g^{\prime}_{R}(\xi_{i})$ are the derivatives of the correction functions
on the solution points. The correction functions have the property of being
symmetric $g_{R}(\xi_{j})=g_{L}(-\xi_{j})$ and they can be computed from the
Radau polynomials:

$R_{R,N+1}(\xi_{j})=\frac{(-1)^{N+1}}{2}\left(L_{N+1}(\xi_{j})-L_{N}(\xi_{j})\right),$
(8)

where $R_{R,N+1}(\xi_{j})$ and $L_{N+1}(\xi_{j})$ are the Radau and Legendre
polynomials of polynomial order $N+1$ at $\xi_{j}$. From the Radau
polynomials, the correction functions can be computed as:

$g_{N+1}=\frac{N+1}{2N+1}R_{R,N+1}+\frac{N}{2N+1}R_{R,N}.$ (9)

Figure 1 shows how this process works, where the black line shows the
interpolation to the interface of the element. Then, the values at the
interface are updated to a new common value which is the red dot at the
interface. The updated quantity is interpolated back to the solution points
making the reconstruction inside the element $C^{0}$ continuous among
elements.

###  2.2 Filtered Navier-Stokes formalism

The Navier-Stokes equations in physical space can be filtered with a spatial
filter that can commute with the derivative operation, a variable with the bar
on top is considered filtered $\bar{f}$. For compressible flows it is
suggested to use a Favre-based filter:

$\check{f}=\frac{\overline{\rho f}}{\bar{\rho}},$ (10)

which leads to the Favre-filtered Navier-Stokes equations:

$\frac{\partial\overline{\rho}}{\partial
t}+\frac{\partial\overline{\rho}\check{u}_{j}}{\partial x_{j}}=0$ (11)

$\frac{\partial\overline{\rho}\check{u}_{i}}{\partial
t}+\frac{\partial}{\partial
x_{j}}\left(\overline{\rho}\check{u}_{i}\check{u}_{j}+\overline{p}\delta_{ij}-\mu\check{\sigma}_{ij}+\overline{\rho}\tau_{ij}\right)=0,$
(12)

$\frac{\partial\overline{E}}{\partial t}+\frac{\partial}{\partial
x_{j}}\left((\overline{E}+\overline{p})\check{u}_{j}-k\frac{\partial\check{T}}{\partial
x_{j}}-\mu\check{\sigma}_{ij}\check{u}_{i}+\frac{1}{2}\left(\frac{\gamma\pi_{j}}{\gamma-1}+\overline{\rho}C_{p}q_{j}\right)+\frac{1}{2}\overline{\rho}\zeta_{j}\right)=\mu\epsilon,$
(13)

$\frac{\overline{p}}{\gamma-1}=\overline{E}-\frac{1}{2}\overline{\rho}\check{u}_{i}\check{u}_{i}-\frac{1}{2}\overline{\rho}\tau_{ii}.$
(14)

The last is the state equation in the case of an ideal gas, $k$ and $C_{p}$
are the thermal conductivity and heat capacity at a constant pressure of the
gas. Among the nonlinear terms ($\tau_{ij}$, $q_{j}$, $\pi_{j}$, $\zeta_{j}$,
and $\epsilon$) the SFS kinetic energy advection $\zeta_{j}$ and the SFS heat
dissipation $\epsilon$ can be considered negligible. The Smagorinsky and
Dynamic Smagorisky models assume negligible also the SFS pressure-work
$\pi_{j}$, which is instead modeled by the BSS method because leads to better
results when the model is used for shock-capturing. The quantities can be
analytically computed as:

$\tau_{ij}=\widecheck{u_{i}u_{j}}-\check{u}_{i}\check{u}_{j}$ (15)

for the SFS stress tensor,

$q_{j}=\widecheck{Tu_{j}}-\check{T}\check{u}_{j}$ (16)

for the SFS heat flux, and

$\pi_{j}=\overline{pu_{j}}-\overline{p}\check{u}_{j}$ (17)

for the SFS pressure-work. Considering that BSS is a purely numerical method
all the gradients are computed in computational space and to have a similar
comparison also the other models use the same gradients. The Smagorinsky
models are adapted to the computational gradients and considering that the
computational domain is perfectly orthogonal in this specific case the
predicted SFS quantities would be the same in physical and computational
space. The SFS quantities are then added to the physical flux vectors:

$\mathbf{F}^{\mathrm{SFS}}_{j}=\begin{bmatrix}0\\\ \bar{\rho}\tau_{11}\\\
\bar{\rho}\tau_{12}\\\ \bar{\rho}\tau_{13}\\\
\frac{1}{2}\left(\frac{\gamma\pi_{1}}{\gamma-1}+\overline{\rho}C_{p}q_{1}\right)\end{bmatrix}\mathbf{G}^{\mathrm{SFS}}_{j}=\begin{bmatrix}0\\\
\bar{\rho}\tau_{21}\\\ \bar{\rho}\tau_{22}\\\ \bar{\rho}\tau_{23}\\\
\frac{1}{2}\left(\frac{\gamma\pi_{2}}{\gamma-1}+\overline{\rho}C_{p}q_{2}\right)\end{bmatrix}\mathbf{H}^{\mathrm{SFS}}_{j}=\begin{bmatrix}0\\\
\bar{\rho}\tau_{31}\\\ \bar{\rho}\tau_{32}\\\ \bar{\rho}\tau_{33}\\\
\frac{1}{2}\left(\frac{\gamma\pi_{3}}{\gamma-1}+\overline{\rho}C_{p}q_{3}\right)\end{bmatrix}$
(18)

and then the code continues to work as described in 2.1.

####  2.2.1 Smagorinsky

Smagorisky [29] proposed a model to close the filtered Navier-Stokes equations
using the strain-rate tensor of the flow. The model was developed for
incompressible flow and hence closing only $\tau_{ij}$, but it can be extended
to compressible flow with the closure of $q_{j}$. The model relies on the eddy
viscosity which is computed as:

$\nu_{t}=2C\Delta^{2}|\check{S}|$ (19)

where $C=0.0256$ is the Smagorinsky constant, $\Delta$ is the spatial length
scale (for our case in computational space),
$|\check{S}|=\sqrt{2\check{S}_{ij}\check{S}_{ij}}$ is the norm of the strain-
rate tensor, and $\check{S}_{ij}$ is the strain-rate tensor. The strain-rate
tensor is computed using the computational space derivatives and hence
$\check{S}_{ij}=\frac{1}{2}\left(\frac{\partial\check{u}_{i}}{\partial\xi_{j}}+\frac{\partial\check{u}_{j}}{\partial\xi_{i}}\right)$.
To be consistent with the unit measure also the spatial length scale $\Delta$
needs to be computed in computational space, which can be computed using the
Gauss-Legendre quadrature weights $\Delta=\sqrt[3]{w_{1}w_{2}w_{3}}$ where the
numbering is the direction. When the mesh is perfectly orthogonal this
formulation is the same as if the derivatives are computed in physical space
and the spatial length scale is the physical one.

From the eddy viscosity the SFS stress tensor can be computed as:

$\tau_{ij}=-\nu_{t}(\check{S}_{ij}-\frac{1}{3}\check{S}_{kk}\delta_{ij})$ (20)

where $\delta_{ij}$ is the Kronecker delta. The SFS heat flux can be computed
as:

$q_{j}=-\frac{\nu_{t}}{Pr_{t}}\frac{\partial\check{T}}{\partial\xi_{j}}$ (21)

where $Pr_{t}=0.9$ is the turbulent Prandtl number and it is assumed to be
constant in this formulation.

####  2.2.2 Dynamic Smagorinsky

Given the limitation of the Smagorinsky model of being bound to the same model
for all kinds of turbulence, Moin et al. [30] developed a dynamic model which
is able to adapt to different kinds of turbulence. The original model was
proposed for incompressible flow by Germano et al. [33] and then extended to
compressible flow by Moin et al.. The Dynamic Smagorinsky model is based on
the standard Smagorinsky model but instead of assuming $C$ and $Pr_{t}$
constant, they are computed dynamically with equations–15 and 19 in
[30]–derived from algebraic identities, which involves a lot of filtering
operations. The model was developed to work in fully spectral code and to the
best of my knowledge has never been applied to a pseudo-spectral code as
H3AMR. Hence, we choose to apply the model to each element. Therefore, each
element is treated on its own and it has a specific value of $C$ and $Pr_{t}$.
Because of the small domain covered by the elements it can happen that the
constants are negative, therefore the code imposes a lower limit of 0 for $C$
and $0.01$ for $Pr_{t}$.

####  2.2.3 Vreman’s model

On the base of the Smagorinsky model Vreman [31] developed an Eddy viscosity
model specifically for turbulent shear flow. The model estimates $\nu_{t}$ to
close the filtered Navier-Stokes equation but in this case the dissipation is
small for transitional and near-wall region, making it particularly suited to
be applied in wall bounded or boundary layer flows. The eddy viscosity is
computed as:

$\nu_{t}=C_{Vr}\sqrt{\frac{\beta_{11}\beta_{22}-\beta_{12}^{2}+\beta_{11}\beta_{33}-\beta_{13}^{2}+\beta_{22}\beta_{33}-\beta_{23}^{2}}{\alpha_{ij}\alpha_{ij}}}$
(22)

where

$\beta_{ij}=\Delta_{m}^{2}\alpha_{mi}\alpha_{mj}$ (23)

and

$\alpha_{ij}=\frac{\partial\check{u}_{j}}{\partial\xi_{i}}.$ (24)

$C_{Vr}$ is a constant that in this case is set to 0.07 and $\Delta_{m}$ is
the grid spacing along the $m^{th}$ direction, in this to be consistent with
our formulation the grid spacing is computed in computatinal space.

####  2.2.4 Block Spectral Stresses (BSS)

The Block Spectral Stresses (BSS) method is based on the work of Sousa and
Scalo [1] and specifically redesigned for pseudo-spectral code. This change of
approach is due to the fact that LSV was only validated as a shock capturing
method and does not work as LES method because it makes the simulation blow
up, even though it is based on a model Quasi Spectral Viscosity (QSV) [10],
which is designed for spectral method. The new model computes the estimated
cutoff kinetic energy block by block using the gradient of the velocity in
each solution point:

$E_{N}^{i}(\check{u}_{j})=\left|\frac{1}{\gamma_{N}}\sum_{k=0}^{N}\left(\ell_{j}\frac{\partial\check{u}_{j}(\xi_{k}^{i})}{\partial\xi_{i}}\right)^{2}L_{N}(\xi_{k}^{i})w_{k}\right|$
(25)

where $\frac{\partial\check{u}_{j}(\xi_{k}^{i})}{\partial x_{i}}$ is the
computational derivative of $\check{u}_{j}$ along the $j$ direction computed
with the correction polynomials, $L_{N}(\xi_{k}^{i})$ is the value of the
Legendre polynomial at $\xi_{k}^{i}$ ($i$ is the direction in which the value
is filtered), $w_{k}$ is the Gauss Legendre quadrature weight, $\ell^{j}$ is
the computational length scale of the solution point (which in the
computational space is equal to $w_{j}$), $\gamma_{N}=\frac{2}{2N+1}$, and $N$
is the polynomial order of the function used to reconstruct the solution
inside the element. BSS method has better shock-capturing capability than LSV
and made the SFS quantities more physical, considering that for LSV were not
even symmetrical on a 2D vortex in a freestream flow, because the cutoff
energy was computed from $\check{u}_{j}(\xi_{k}^{i})^{2}$ instead of using the
derivative. Indeed the goal of taking only the last mode of the kinetic energy
was to remove any freestream interference but considering that the velocity is
squared the freestream component is not removed. Instead, the gradient of the
velocity removes the freestream interference because it considers the rate of
change and not the actual value.

From the estimated cutoff energy, it is possible to compute the subfilter
velocity scale
$v_{i}(\check{u}_{j})=\sqrt{\frac{N}{2}E_{N}^{i}(\check{u}_{j})}$, where the
factor $N/2$ is the average grid spacing of the element in computational
space. From that it is possible to estimate the dissipation needed by the
model:

$\mathcal{D}_{ij}=v_{i}(\check{u}_{j})\ell_{j}.$ (26)

From the dissipation the SFS quantities can be computed as:

$\tau_{ij}=-\frac{1}{2}\left(\mathcal{D}_{ij}\frac{\partial\check{u}_{i}}{\partial\xi_{j}}+\mathcal{D}_{ji}\frac{\partial\check{u}_{j}}{\partial\xi_{i}}\right)$
(27)

$q_{j}=-\mathrm{diag}(\mathcal{D})_{j}\frac{\partial\check{T}}{\partial\xi_{j}}$
(28)

$\pi_{j}=-\mathrm{diag}(\mathcal{D})_{j}\frac{\partial\bar{p}}{\partial\xi_{j}}$
(29)

where $\mathrm{diag}(\mathcal{D})$ is the vector of the diagonal of the matrix
and all the indices are not used in Einstein’s notation but as the actual
index of the vector or matrix. Important consideration is that unlike the
Smagorinsky models the diagonal of the strain-rate tensor is not removed to
compute the SFS stress tensor because it has an important role when the model
is used as shock capturing. Moreover, differently than LSV were the SFS
quantities are modulated for this model they are not.

##  3 BSS shock capturing capability

Given the fact that BSS is both a shock-capturing and LES method in this
section we are going to show its capability as shock-capturing method. The
first case use for the investigation is the Sod shock tube (Section 3.1) and
then the investigation moved to a 2D shock in the shock vortex interaction
case (Section 3.2).

###  3.1 Sod shock tube

Sod shock tube [34] is a standard case widely used to test shock-capturing
methods [35, 36, 1]. The case simulates a tube with a gas with high pressure
and density on the left and lower values on the right separated by a membrane.
When the membrane breaks a shock and an expansion wave form propagating
respectively to the right and to the left. The problem is described by the
Euler equations and as initial conditions:

$[u_{1},p,\rho]=\begin{cases}[0,1,1]\ \ \ \mathrm{if}\ x\leq 0\\\ \
[0,0.1,0.125]\ \ \ \ \ \ \mathrm{if}\ \ \ x>0\end{cases}$ (30)

where $u_{1}$ is the velocity, $p$ the pressure, and $\rho$ the density. The
biggest constraint is given by the shock that is difficult to simulate with
high order polynomial and hence the need of a shock-capturing method.

Figure 2: Sod shock tube [34] velocity $u_{1}$, density $\rho$ and pressure
$p$ at $t=0.264$ s. The first row has a domain of 10 elements and the second
of 20 elements. The profiles in all plots from bottom to top used numerical
order inside each element 10, 15, and 20 (which correspond to 100, 150, and
200 DOF on the first row and 200, 300, and 400 DOF on the second row).

Figure 2 compares the velocity, pressure, and density profiles at $t=0.264$ s
of the simulated case with BSS and the analytical solution on different mesh
and numerical order $N$ resolution. On the top row the results are obtained
using 10 elements and in each plot the numerical order is 10, 15, and 20 from
the bottom to the top. Most of the oscillation around the shock are visible in
the velocity profile instead the pressure looks mostly smooth. As expected
increasing the numerical order allows to have a better resolution of the shock
and reduces the oscillation after the shock but on the same time it increases
them before the shock. However this phenomenon is not depended on the shock
but on the fact that shock is closer to the left face of the element, indeed
also on the right face it is possible to see some oscillation. On the lowest
row of Figure 2 are presented the results obtained on a mesh with 20 elements.
As before the numerical order used is 10, 15, and 20 from the bottom to the
top. The most visible outcome is the fact that the oscillation are reduced
when compare to the same numerical order on a lower elements mesh and also
allows a better resolution of the shock. Therefore, it is possible to conclude
that increasing the number of degrees of freedom (DOF) reduces the oscillation
cased by the shock, but to have a better resolution of the shock it is better
to increase the number of elements and not only the order.

Figure 3: Sod shock tube [34] velocity $u_{1}$, density $\rho$, pressure $p$,
SFS shear stress $\tau_{11}$, SFS heat flux $q_{1}$, and SFS pressure-work
$\pi_{1}$ at $t=0.264$ s, for a 10 elements mesh with numerical order 20. The
SFS are computed using equations 27-29. The last row of plots presents the
spectral content–computed using the Legendre coefficients–of the three SFS
quantities computed in the element with the shock.

Figure 3 shows the velocity, pressure, density, SFS shear stress, SFS heat
flux, and SFS pressure-work for the 10 elements and numerical order 20 case.
The exact SFS quantities are computed using equations 27-29 and the BSS one
are the a posteriori results of a simulation. As expected considering that the
shock has an infinite spectra the SFS quantities have a lot of oscillations
and the modeled quantities are able to precisely match the analytical
quantities. The last row of plots presents the spectral content computed using
the Legendre coefficients and then purely real of the SFS shear stress, SFS
heat flux, and SFS pressure-work in the element with shock. For all cases, the
a posteriori spectra match the behavior of the exact solution spectra, with
the oscillation of sign and similar magnitude. There are some differences in
the spectra of the SFS shear stress but it is possible to see that this
difference is due to the discrete nature of the investigation, indeed if the
results were continuous the exact solution spectra would oscillate further.

###  3.2 Shock vortex interaction

The shock vortex interaction is a case that investigates the interaction
between a vortex without circulation and a standing shock. For this
investigation the domain is $\left[0,2L\right]\times\left[0,L\right]$, the
freestream velocity before the shock is $V_{0}=1.5\sqrt{\gamma
p_{0}/\rho_{0}}$, the shock location is at $x_{s}/L=1/2$, and the vortex
center at $x_{cv}/L=1/4$ and $y_{cv}/L=1/2$. The vortex is initialized with
the tangential velocity that scales linearly respect to the radius in the
inner part and a combination of linear and inverse in the external part. The
inner part radius is $a/L=0.075$ and the outer part $b/L=0.175$ the velocity
of the vortex can be written as
$\mathbf{u}=u_{\theta}(r)\hat{\mathbf{e}}_{\theta}+V_{0}\hat{\mathbf{e}}_{x}$
where:

$\frac{u_{\theta}(r)}{u_{\theta}(a)}=\begin{cases}\frac{r}{a},&\text{if }r\leq
a\\\ \frac{\eta}{2}\left(\frac{r}{b}-\frac{b}{r}\right),&\text{if }a<r\leq
b\\\ 0,&\text{otherwise,}\end{cases}$ (31)

where $\eta=2(b/a)/[1-(b/a)^{2}]$ and the maximum tangential velocity is
$u_{\theta}(a)=0.9V_{0}$. The pressure is assumed for ideal gas and isentropic
compression and it is initial so that the gradient balance the centrifugal
forces:

$\frac{\partial P}{\partial r}=\rho\frac{{u}_{\theta}^{2}(r)}{r},\quad P=\rho
RT,\quad\frac{P}{P_{0}}=\left(\frac{\rho}{\rho_{0}}\right)^{\gamma}.$ (32)

Other works [37, 38, 16] run simulations on the configuration and reported
that the interaction of the vortex with shock leads to the creation of two
smaller vortices with the upper one in a forward position respect to the
lower. Figure 4 shows the evolution of the vortex passing through the shock
deforming the shock into an S-shape. The vortex is splitted in two smaller
vortices that continue to travel to the right and as described in previous
works the upper one is in a forward position. This test is important because
it tests the shock-capturing capability also on an oblique shock and the same
time it does not dissipate the vortex.

Figure 4: Shock-vortex interaction evolution using numerical order 4 and a
mesh 256x128 blocks (524,768 DOF).

Figure 5 tabulates a comparison of of different numerical order (in the raws)
and the their refinementment (in the columns). As previously seen in the Sod
shock tube investigation when the number of degrees of freedom is the same,
the simulation with the lower numerical order has the best result. This is
reported in the 131,072 DOF case where with numerical order 16 the upper
vortex is completely dissipate, instead with order 4 is fully resolved. With
524,288 DOF the resolution of all numerical order reach convergence and they
show similar results. Therefore, it is possible to conclude that the method
works better when using more elements, this was shown in the case of a shock
when looking at the Sod shock tube simulation and now it is fatherly validate
also in the case of vortices.

Figure 5: Shock-vortex interaction at $tV_{0}/L=1.65$ with different numerical
order and mesh resolution. Each row have different degrees of freedom in the
simulation, respectively 32768, 131072, and 524288. The meshes have in x and
y: a) 16x8, b) 32x16, c) 64x32, d) 128x64, and e) 256x128.

##  4 LES comparison

###  4.1 Taylor-Green Vortex (TGV)

The Taylor-Green Vortex (TGV) is a 3D flow in a cube in which the flow tends
to a state of homogeneous turbulence. Even though there exists an analytical
solution for the incompressible case, considering that in this investigation
the code used simulates compressible flow the results are compared with a DNS
of the flow. The flow starts at $t=0$ s and ends at $t=200$ s in a cubed box
$[-\pi L,\pi L]^{3}$ with the following initial conditions:

$\rho(\mathbf{x},0)/\rho_{0}=1,$ (33)

$\frac{u_{1}(\mathbf{x},0)}{V_{0}}=\mathrm{sin}\left(\frac{x_{1}}{L}\right)\mathrm{cos}\left(\frac{x_{2}}{L}\right)\mathrm{cos}\left(\frac{x_{3}}{L}\right),$
(34)

$\frac{u_{2}(\mathbf{x},0)}{V_{0}}=-\mathrm{cos}\left(\frac{x_{1}}{L}\right)\mathrm{sin}\left(\frac{x_{2}}{L}\right)\mathrm{cos}\left(\frac{x_{3}}{L}\right),$
(35)

$\frac{u_{3}(\mathbf{x},0)}{V_{0}}=0,$ (36)

$\frac{p(\mathbf{x},0)}{\rho_{0}V^{2}_{0}}=\frac{p_{0}}{\rho_{0}V^{2}_{0}}+\frac{1}{16}\left[\mathrm{cos}\left(\frac{2x_{1}}{L}\right)+\mathrm{cos}\left(\frac{2x_{2}}{L}\right)\right]\left[\mathrm{cos}\left(\frac{2x_{3}}{L}\right)+2\right]$
(37)

where $\rho_{0}=1$, $V_{0}=0.1$, $p_{0}=1/\gamma$, and $L=1$ are the non-
dimensionalized density, velocity, pressure, and length. The Mach number of
the flow is $M_{0}=\frac{V_{0}}{\sqrt{\gamma p_{0}/\rho_{0}}}=0.1$, which can
be considered incompressible but the solution completely diverges from the
analytical incompressible solution. The Reynolds number is
$Re=\frac{\rho_{0}V_{0}L}{\mu_{0}}=5000$ where the viscosity is considered to
be constant at $\mu_{0}=2\cdot 10^{-5}$. All the simulations use a third-order
Runge-Kutta time advancement method with $CFL=0.1$. Grid convergence is
achieved a $32^{3}$ elements mesh running with numerical order 9.

Figure 6: Normalized time and dissipation $\epsilon_{kin}$ with different
degrees of freedom (DOF) and numerical order. From left to right the meshes
use for the $96^{3}$ DOF case have respectively $32^{3}$, $16^{3}$, and
$8^{3}$ elements. For the $192^{3}$ DOF respectively $64^{3}$, $32^{3}$, and
$16^{3}$ elements.

Similarly to the shock-vortex interaction case, figure 6 reports how different
methods behaves with different degrees of freedom (DOF) and numerical orders.
On the y-axis there is the normalized dissipation computed as
$\epsilon_{kin}=-\frac{\partial E_{kin}}{\partial t}$ and on the x-axis the
normalized time. When using numerical order 3 only BSS is able to resolve the
turbulence instead the order method starts to oscillate and then blow up.
Differently from spectral method where the blow up is caused by the
dissipation becoming positive, in block spectral code the blow up is caused by
oscillations in the single block which in some cases can be seen in the
general statics as the dissipation. When running with $96^{3}$ DOF the BSS is
significantly more dissipative than the other method when the turbulence
starts to break down and consequently the peak dissipation is smaller than the
other methods. In the $192^{3}$ DOF case BSS is the only method able to
capture the plateau after the peak instead the other methods have a single
peak or another peak after the first one where the plateau should be. The
dynamic Smagorinsky method seems a filtered version of the case without any
model, this is caused by the fact that the model is applied by block by block
so the effect on the simulation is small but strong enough to reduce the
oscillation of not having a model.

Figure 7: Dissipation convergence of different mesh and numerical resolution
when using the BSS method. Solid lines indicate cases with $96^{3}$ DOF and
dashed lines $192^{3}$ DOF. Red, black, and blue color indicate simulation
respectively with numerical order 12, 6, and 3.

Figure 7 plots the dissipation of the different mesh and numerical order so to
highly how the method perform with different setup. All the cases with
$96^{3}$ DOF have very different dissipation indeed in the $32^{3}$ order 3
the peak dissipation happens sooner and is stronger than what it should be.
The order 12 case is the one with the best performance due to the fact that in
this case the model is able to be more active and is able to better capture
the turbulence. The cases with $192^{3}$ DOF are closer to the DNS resolution
and for this reason the mesh or the order do have a big different has it has
been seen in the shock-vortex interaction case. The biggest difference is the
peak that is stronger for the order 12 simulation thanks to a bigger spectra
available to the model.

Table 1: Time, number of iterations, and the second per iteration all normalized respect to the no model case for TGV simulation with no model, Smagorinsky, Dynamic Smagorinsky, and Legendre Spectral Viscosity (BSS) methods. The results are obtained on a $16^{3}$ elements mesh and numerical order 6. | Time | # iterations | Time per iteration
---|---|---|---
No model | 1.000 | 1.000 | 1.000
Smagorinsky | 1.009 | 0.987 | 1.024
Dynamic Smagorinsky | 0.940 | 0.801 | 1.173
BSS | 1.067 | 0.990 | 1.078

Table 1 reports the computational time, the number of iterations, and the time
per iteration by all LES models and without any model, on a $16^{3}$ elements
mesh with numerical order 6. As expected the Dynamic Smagorisky model is the
slowest due to the big number of filtered quantities that it needs to estimate
the coefficients. The standard Smagorisky is the fastest considering that it
only needs to multiply a constant with the gradients of the velocity and the
temperature. BSS is placed in the middle of the two Smagorisky models being 6%
slower than the standard one. From the number of iterations it is possible to
conclude that all models lead to a lower time constraint induced by the SFS
stresses because all of them require fewer iterations than without a model.

###  4.2 Supersonic channel flow

Figure 8: y-z section of channel flow mach 1.5 and no-slip isothermal boundary
condition at the wall. On the left column the velocity, density, and
temperature contours and on the right one the SFS shear stress, pressure-
strain, and heat flux. The first row of each quantity is on a 32x32x24 mesh
and numerical order 5. The second is on a 64x64x48 mesh with numerical order
3.

A compressible turbulent channel flow has been used to check BSS performance
in wall bounded flow, specifically in a supersonic case Mach 1.5 and an
hypersonic case Mach 6.0. The flow condition match the one used by Chen and
Scalo [39], all the flow parameters are normalized using the channel’s half-
width ($\delta$), the speed of sound at the wall ($a_{w}$), the wall
temperature ($T_{w}$), and the bulk density $\rho_{w}$. The results are
plotted using the Trettel-Larsson transformation [40], where the mean velocity
is transformed as:

$U^{+}_{TL}=\int_{0}^{\check{u}/u_{\tau}}\sqrt{\frac{\overline{\rho}}{\rho_{w}}}\left(1+\frac{1}{2}\frac{1}{\overline{\rho}}\frac{\partial\overline{\rho}}{\partial
y}y-\frac{1}{\overline{\mu}}\frac{\partial\overline{\mu}}{\partial
y}y\right)d\left(\frac{\check{u}}{u_{\tau}}\right)$ (38)

where $u_{\tau}=\sqrt{\overline{\tau_{w}}/\overline{\rho_{w}}}$, plotted
versus:

$y^{*}=\frac{\overline{\rho}(y)u^{*}_{\tau}}{\overline{\mu}(y)}y$ (39)

where $u_{\tau}=\sqrt{\overline{\tau_{w}}/\overline{\rho}(y)}$. The reference
log-law profile is computed as:

$U^{+}_{TL}=\frac{1}{\kappa}\mathrm{ln}(y^{*})+C$ (40)

where $\kappa=0.41$ and $C=5.5$.

Figure 9: TL-transformed mean velocity profile for Mach 1.5 compressible
channel flow. On the left using a 32x32x24 elements mesh numerical order 5 and
on the right a 64x48x48 mesh numerical order 3.

Figure 8 shows on the left column the velocity, density, and temperature
contours and on the right one the SFS shear stress, pressure-strain, and heat
flux for a 64x48x48 mesh order 3 and 32x32x24 mesh order 5. It is possible to
see that the finer mesh is able to better resolve the boundary layer and the
flow looks more continuous, instead in the coarse mesh the boundary layer is
thicker and sometimes the flow looks discontinuous cell by cell but this
problem is mainly due to Paraview interpolation and partially to the block
spectral nature of the solver. The SFS quantities are not as continuous as the
flow also for the finer mesh case in this instance the effect is because BSS
computes these quantities block by block and does not enforce continuity. It
is interesting to notice that the SFS quantities have a stronger magnitude on
the finer mesh case where the numerical order is lower and they are weaker
with high numerical order. This result agrees with the fact that a high
numerical order is able to better resolve the dynamic properties of the flow.

Figure 9 presents more quantitative results of the same case with the plots of
the TL-transformed mean velocity. On the left has been used a 32x32x24
elements mesh with numerical order 5, it is possible to see that the BSS model
has the best performance together with Vreman’s model indeed the Smagorinsky
model is not able to reduce the velocity as much as BSS. On the right has been
used a finer mesh with 64x48x48 elements and numerical order 3. In this case
Smagorinsky and Vreman model precisely match the DNS results instead BSS do
not fully resolve the turbulence in the log-law layer. The fact that BSS is
not able to be as good as the other models when using a finer mesh with lower
numerical order can be also explained by the SFS terms shown in Figure 8 where
the SFS quantities are stronger when using lower numerical order. Indeed BSS
is more related to the numerical order then the other two methods because to
compute the SFS energy dissipation it needs the energy spectrum. Instead
Smagorinsky and Vreman models depend on the numerical order only by the grid
spacing.

Figure 10: TL-transformed mean velocity profile for Mach 6.0 compressible
channel flow. On the left using a 96x75x48 elements mesh numerical order 3.

Figure 10 plots the results for the Mach 6.0 hypersonic case. The mesh used
for the simulation has 96x75x48 elements and uses numerical order 3. For this
case the Smagorinsky and BSS model leads to a similar solution that match the
DNS results, instead Vreman’s model predictis a lower velocity closer to the
log-law profile. In both the Mach 1.5 and 6.0 the Dynamic Smagorinsky did not
work probably because the averaging is performed element by element and not by
layer as suggested by the authors for this kind of simulation. Considering
that the averaging have to variate case by case it would be difficult both to
tune and to code in a CFD solver and hence we can conclude that Smagorisky,
Vreman and BSS are versatile enough to work with different flow condition in a
block spectral code.

##  5 Conclusion

In this paper, we presented a new shock-capturing and LES model for high-order
finite volume methods called Block Spectral Stresses (BSS) estimation. The
method relies on the last mode of the velocity gradient inside each element to
estimate the subfilter scale stresses, heat-flux, and pressure-work to close
the filtered Navier-Stokes equations. The main idea is that both turbulence
and shock formation have a cascade of energy from large to small scale and
hence a closure for the turbulence model should work also for shock-capturing.
Specifically, we noticed that in our model the subfilter scale pressure-work,
which is rarely modeled in LES closure, is really important to capture the
shock.

BSS has been implemented in a block spectral unstructured code called H3AMR
and validated on simulations with shocks and turbulence. In the Sod shock tube
case, the model allows to simulate up to numerical order 25, we were not able
to reach higher order because of the limitation but the model can potentially
allow higher orders. From this investigation, we also notice how the shock is
better predicted when the mesh is refined and not as much when the numerical
order is increased. Another shock investigation was on the shock-vortex
interaction case, where the model demonstrated to be able to capture the shock
without altering the shape of the vortex. Also in this case, the mesh
refinement leads to faster convergence than when increasing the numerical
order.

For the LES closure, BSS has been compared to the Smagorinsky, dynamic
Smagorinsky, and Vreman models, with the second changed to work in a block
spectral code. The Taylor-Green vortex was used to test the model in a
homogeneous isotropic simulation. In the case of numerical order 3 BSS is the
only model able to estimate the small scale turbulence and avoid the
simulation blow-up. For numerical order 6 and 12 and coarse mesh BSS is more
dissipative than the other methods but on finer meshes, it is able to reach
similar results to the other models. Another turbulence simulation was
performed for supersonic and hypersonic channel flow; in this case, dynamic
Smagorinsky does not work because of the averaging chosen for the block
spectral implementation, instead BSS, the standard Smagorinsky, and Vreman
were able to simulate the flow. Differently, from TGV in the supersonic case,
BSS achieved better performance at predicting the velocity profile when using
a coarser mesh with higher numerical order because of a better prediction of
the SFS terms. For the hypersonic case, Smagorinsky is slightly closer to the
DNS result than BSS but the difference is minimal, instead Vreman’s model
estimates a lower velocity closer to the log-law profile instead of the DNS.

Therefore, we can conclude that even though BSS is not the best model in every
kind of flow is able to work with every flow without big changes in the code,
as in the case of dynamic Smagorinsky which requires different averaging for
different kinds of flows. Currently, the model has been used with flux
reconstruction numerics a Legendre base but we expect that the model can work
also with different numerics and base on high-order finite volume. This last
point requires further investigation.

## References

  * [1] Victor C. B. Sousa and Carlo Scalo.  A legendre spectral viscosity (lsv) method applied to shock capturing for high-order flux reconstruction schemes.  Journal of Computational Physics, 460:111157, 2022.
  * [2] Sergei K Godunov and I Bohachevsky.  Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics.  Matematičeskij sbornik, 47(3):271–306, 1959.
  * [3] Eleuterio F Toro.  Riemann solvers and numerical methods for fluid dynamics: a practical introduction.  Springer Science & Business Media, 2013.
  * [4] Hung T Huynh.  A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods.  In 18th AIAA computational fluid dynamics conference, page 4079, 2007.
  * [5] Hung T Huynh.  A reconstruction approach to high-order schemnes including discontinuous Galerkin for diffusion.  In 47th AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition, page 403, 2009.
  * [6] HT Huynh, Zhi J Wang, and Peter E Vincent.  High-order methods for computational fluid dynamics: A brief review of compact differential formulations on unstructured grids.  Computers & fluids, 98:209–220, 2014.
  * [7] David Gottlieb, Chi-Wang Shu, Alex Solomonoff, and Hervé Vandeven.  On the Gibbs phenomenon I: Recovering exponential accuracy from the Fourier partial sum of a nonperiodic analytic function.  Journal of Computational and Applied Mathematics, 43(1-2):81–98, 1992.
  * [8] Cornelius Lanczos and John Boyd.  Discourse on Fourier series.  SIAM, 2016.
  * [9] David Gottlieb and Steven A Orszag.  Numerical analysis of spectral methods: theory and applications.  SIAM, 1977.
  * [10] Victor C. B. Sousa and Carlo Scalo.  A unified quasi-spectral viscosity (qsv) approach to shock capturing and large-eddy simulation.  Journal of Computational Physics, 459:111139, 2022.
  * [11] Prateek Gupta, Guido Lodato, and Carlo Scalo.  Spectral energy cascade in thermoacoustic shock waves.  Journal of Fluid Mechanics, 831:358–393, 2017.
  * [12] F. Bassi and S. Rebay.  A high-order accurate discontinuous finite element method for the numerical solution of the compressible navier–stokes equations.  Journal of Computational Physics, 131(2):267–279, 1997.
  * [13] B. Cockburn and C.W. Shu.  Tvb runge-kutta local projection discontinuous galerkin finite element method for conservation laws. ii. general framework.  Mathematics of Computation, 52:411–435, 1989.
  * [14] David A. Kopriva and John H. Kolias.  A conservative staggered-grid chebyshev multidomain method for compressible flows.  Journal of Computational Physics, 125(1):244–261, 1996.
  * [15] Kartikey Asthana, Manuel R. López-Morales, and Antony Jameson.  Non-linear stabilization of high-order flux reconstruction schemes via fourier-spectral filtering.  Journal of Computational Physics, 303:269–294, 2015.
  * [16] Niccolò Tonicello, Guido Lodato, and Luc Vervisch.  Entropy preserving low dissipative shock capturing with wave-characteristic based sensor for high-order methods.  Computers & Fluids, 197:104357, 2020.
  * [17] Prateek Gupta and Carlo Scalo.  Spectral energy cascade and decay in nonlinear acoustic waves.  Physical Review E, 98(3):033117, 2018.
  * [18] Morduchow M and Libby P.  On a complete solution of the one-dimensional flow equations of a viscous, heat-conducting, compressible gas.  J Aeronaut Sci, 16(11):674–84, 1949.
  * [19] Joel Smoller.  Shock waves and reaction—diffusion equations, volume 258.  Springer Science & Business Media, 2012.
  * [20] Manuel D Salas and Angelo Iollo.  Entropy jump across an inviscid shock wave.  Theoretical and computational fluid dynamics, 8(5):365–375, 1996.
  * [21] JF Colombeau and AY Le Roux.  Multiplications of distributions in elasticity and hydrodynamics.  Journal of Mathematical Physics, 29(2):315–319, 1988.
  * [22] Pablo Fernandez, Cuong Nguyen, and Jaime Peraire.  A physics-based shock capturing method for unsteady laminar and turbulent flows.  In 2018 AIAA Aerospace Sciences Meeting, page 0062, 2018.
  * [23] F Ducros, P Comte, and M Lesieur.  Large-eddy simulation of a spatially growing boundary layer over an adiabatic flat plate at low mach number.  International journal of heat and fluid flow, 16(5):341–348, 1995.
  * [24] Carson L. Running, Benjamin L. Bemis, J. Luke Hill, Matthew P. Borg, Joel J. Redmond, Karl Jantze, and Carlo Scalo.  Attenuation of hypersonic second-mode boundary-layer instability with an ultrasonically absorptive silicon-carbide foam.  Experiments in Fluids, 64(4), 2023.
  * [25] Benjamin L. Bemis, John L. Brun, C. Taber Wanstall, Jonathan L. Hill, Matthew P. Borg, Joel J. Redmond, Matteo Ruggeri, Karl Jantze, Carlo Scalo, and Carson L. Running.  Ultrasonically absorptive silicon-carbide foam for boundary-layer control.  In AIAA SCITECH 2023 Forum, 2023.
  * [26] Benjamin L. Bemis, Megan C. Sieve, Jonathan L. Hill, Matthew P. Borg, Joel J. Redmond, Matteo Ruggeri, Karl Jantze, Carlo Scalo, and Carson L. Running.  Effect of porosity on the ability of silicon-carbide foams to attenuate the second-mode boundary-layer instability.  In AIAA Aviation 2023 Forum, 2023.
  * [27] Carson L. Running, Benjamin L. Bemis, J. Luke Hill, Matthew P. Borg, Joel J. Redmond, Karl Jantze, and Carlo Scalo.  Attenuation of hypersonic second-mode boundary-layer instability with an ultrasonically absorptive silicon-carbide foam.  Experiments in Fluids, 64(4):79, 2023.
  * [28] David A. Kopriva and John H. Kolias.  A conservative staggered-grid chebyshev multidomain method for compressible flows.  Journal of Computational Physics, 125(1):244–261, 1996.
  * [29] J. Smagorinsky.  General circulation experiments with the primitive equation i the basic experiment.  Monthly Weather Review, 91:99–164, 1963.
  * [30] Parviz Moin, Kyle Squires, W Cabot, and Sangsan Lee.  A dynamic subgrid-scale model for compressible turbulence and scalar transport.  Physics of Fluids A: Fluid Dynamics, 3(11):2746–2757, 1991.
  * [31] A. W. Vremana.  An eddy-viscosity subgrid-scale model for turbulent shear flow : Algebraic theory and applications.  2004\.
  * [32] V.V. Rusanov.  Calculation of interaction of non–steady shock waves with obstacles.  USSR Computational Mathematics and Mathematical Physics, 1:267–279, 1961.
  * [33] Massimo Germano, Ugo Piomelli, Parviz Moin, and William H Cabot.  A dynamic subgrid-scale eddy viscosity model.  Physics of Fluids A: Fluid Dynamics, 3(7):1760–1765, 1991.
  * [34] A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws.  Journal of Computational Physics, 27(1):1–31, 1978.
  * [35] Takanori Haga and Soshi Kawai.  On a robust and accurate localized artificial diffusivity scheme for the high-order flux-reconstruction method.  Journal of Computational Physics, 376:534–563, 2019.
  * [36] Hiroshi Terashima, Soshi Kawai, and Mitsuo Koshi.  Consistent numerical diffusion terms for simulating compressible multicomponent flows.  Computers & Fluids, 88:484–495, 2013.
  * [37] Janet L. Ellzey, Michael R. Henneke, J. Michael Picone, and Elaine S. Oran.  The interaction of a shock with a vortex: Shock distortion and the production of acoustic waves.  Physics of Fluids, 7:172–184, 1995.
  * [38] Audrey Rault, Guillaume Chiavassa, and Rosa Donat.  Shock-vortex interactions at high mach numbers.  Journal of Scientific Computing, 19:347–371, 2003.
  * [39] Yongkai Chen and Carlo Scalo.  Trapped waves in supersonic and hypersonic turbulent channel flow over porous walls.  Journal of Fluid Mechanics, 920:A24, 2021.
  * [40] Andrew Trettel and Johan Larsson.  Mean velocity scaling for compressible wall turbulence with heat transfer.  Physics of Fluids, 28(2), 2016.