Maria Al Jreidy
email: mariajreidy@hotmail.com

Joseph Constantin

Fadi Dornaika

Denis Hamad

Maria Al

Joseph Prefix

Denis Prefix

Graph Convolution Networks for Unsupervised Learning

Keywords: separated by '-') Graph convolution network -Unsupervised learning -Kernel k-means -Spectral clustering -Regularization term Graph convolution network, Unsupervised learning, Kernel k-means, Spectral clustering, Regularization term

HAL is

Introduction

Graphs are powerful forms of describing and modeling complex systems. In recent years, the problem of graph clustering has been well studied [START_REF] Shah | Deep continuous clustering[END_REF][START_REF] Aljalbout | Clustering with deep learning: taxonomy and new methods[END_REF][START_REF] Wang | Unsupervised learning for community detection in attributed networks based on graph convolutional network[END_REF]. Several methods have been proposed to achieve this goal, some of them rely only on the structure of the graph in the learning process, [START_REF] Perozzi | Deepwalk: online learning of social representations[END_REF][START_REF] Cao | Deep neural networks for learning graph representations[END_REF], while others just rely on data features.

One of the well-known methods in the clustering process is k-means clustering [START_REF] Moradi Fard | Jointly clustering with k-means and learning representations, pattern recognition letters[END_REF][START_REF] Gao | Deep Clustering with Concrete kmeans[END_REF] where the number of clusters is initially needed to know. k-means clusters graph by considering node features only. After that, the spectral clustering algorithm [START_REF] Von Luxburg | A tutorial on spectral clustering[END_REF] was introduced where the graph clustering is done based on the structure of the graph built from data features. Then, other methods combining auto-encoders with the clustering algorithm [START_REF] Yang | Deep spectral clustering using dual autoencoder network[END_REF][START_REF] Affeldt | Spectral clustering via ensemble deep autoencoder learning (SC-EDAE)[END_REF] were proposed to give better solutions with high-dimensional databases. For example, adversarially regularized variational graph autoencoder (ARVGE) [START_REF] Pan | Adversarially regularized graph autoencoder for graph embedding[END_REF], that learns node embedding by graph autoencoder and graph variational autoencoder. The spectral embedding network (SENet) [START_REF] Zhang | Spectral embedding network for attributed graph clustering[END_REF] where the graph structure and node feature information via higher-order graph convolution are used to train the model with a spectral clustering loss function.

All these methods can be considered as unsupervised methods as they do not impose prior knowledge of the fundamental truth label in their learning phase.

However, there are other graph-based methods which require a small part of their inputs to be labelled, we can mention the Graph Convolution Network (GCN) [START_REF] Kipf | Semi-supervised classification with convolutional networks[END_REF] which has a very important advantage where the structure of the graph and the features of the data are convolved together in order to obtain the final outputs, namely the soft labels. In this sense, a GCN processes nodes in a graph based on their structure as well as their feature similarities. Yet, as said before a GCN usually relies on prior knowledge of some labels that are often not available for multiple databases, as data labeling is a tedious task to perform.

To keep the advantage of the GCN, and to be able to use it with unlabeled data, unsupervised graph clustering using GCN model where proposed. In [START_REF] Jin | Community detection via joint graph convolutional network embedding in attribute network[END_REF] an unsupervised model for clustering graphs via a joint GCN has been introduced. In [START_REF] Wang | Unsupervised learning for community detection in attributed networks based on graph convolutional network[END_REF] a label sampling method was introduced and used before the GCN model. However, these methods are considered complex as they require labelling generation techniques followed by the GCN to use it in an unsupervised way.

In this paper, we propose a new unsupervised learning framework based on GCN only. The main contribution of our work is as follows: Without any prior knowledge of node labels, we propose a GCN model trained with the kernel kmeans objective function where the kernel matrix K represents the data features, integrated with a regularization term reflecting the graph structure. In other words, our proposed loss has two clustering-friendly terms. In this case, if two nodes share a larger proportion of neighbors and features, they have a higher probability of belonging to the same cluster. Therefore, both network topology and data features are considered during model training to obtain consistent clustering.

The rest of this paper is organized as following. In Sect. 2, we introduce some related work. In Sect. 3, we present some preliminary notations and formulations related to this work. Section 4 gives a detailed description of the proposed method. We present experimental results on various real-world networks to validate the performance of the method proposed in Sect. 5, we also compare our approach with other state-of-the-art algorithms. Finally, in Sect. 6, we give the conclusion.

Related Work

Early graph clustering techniques merely grouped nodes based on the graph topology. Such as, DeepWalk [START_REF] Perozzi | Deepwalk: online learning of social representations[END_REF] and node2vec [START_REF] Grover | Node2vec: scalable feature learning for networks[END_REF], which produce node sequences using reduced random walks before obtaining node embeddings using the Skip-gram model [START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF]. Different k-order proximities between nodes are captured by GraRep [START_REF] Cao | Grarep: learning graph representations with global structural information[END_REF]. Then, models based on autoencoders were created in order to capture the highly nonlinear graph structures [START_REF] Wang | Structural deep network embedding[END_REF]. In order to increase the model strength and generalizability, ANE [START_REF] Dai | Adversarial network embedding[END_REF], a clustering technique that uses adversarial learning to learn node representations, attempts to regularize the embedding learning process.

The developed node-clustering techniques that use both node properties and graph structure are referred to as the attributed methods. Some techniques employ spectral clustering [START_REF] Von Luxburg | A tutorial on spectral clustering[END_REF], random walks [START_REF] Zhou | Clustering large attributed graphs: an efficient incremental approach[END_REF], matrix factorization, and Bayesian models [START_REF] Xu | A model-based approach to attributed graph clustering[END_REF] on them to make use of both structure and feature information. Other approaches [START_REF] Zhou | Graph clustering based on structural/attribute similarities[END_REF] develop a trade-off distance metric between them.

Graph convolutional network (GCN) is a semi-supervised approach that was recently presented for graph classification problems [START_REF] Kipf | Semi-supervised classification with convolutional networks[END_REF]. GCN-based algorithms help to merge network topology and attribute information, but they depend on a large number of node labels to classify unlabeled nodes, in contrast to most semisupervised methods that concentrate on maintaining network structure. Sun et al. [START_REF] Sun | Network embedding for community detection in attributed networks[END_REF] suggested a network embedding framework for node clustering based on a graph convolutional autoencoder. Additionally, few unsupervised techniques have been recently developed. An unsupervised model for community detection using GCN embedding was proposed by Jin et al. [START_REF] Jin | Community detection via joint graph convolutional network embedding in attribute network[END_REF]. A spectral embedding network for attributed graph clustering (SENet) was suggested in [START_REF] Zhang | Spectral embedding network for attributed graph clustering[END_REF] and uses a spectral clustering loss with GCN to learn node embeddings while also enhancing graph structure. In this work,the GCN model formed with two loss functions commonly used in clustering problems, which manages both data structure and data attributes.

Preliminary

The section briefly presents the preliminary knowledge of this work, and the architecture of convolution networks of graphs.

Weighted Kernel K-Means and Spectral Clustering

Weighted kernel k-means is a k-means clustering algorithm enhanced by the use of a kernel function [START_REF] Dhillon | Kernel k-means, spectral clustering and normalized cuts[END_REF]. The kernel weighted k-means objective function is defined as follows:

D({π j } k j=1 ) = k j=1 vi∈πj w(x i ) φ(x i ) -m j 2 (1)
with w(x i ) is the weight of data x i , π j denote the clusters, k is the number of clusters, φ a non-linear function and m j =

x i ∈π j w(xi)φ(xi)

x i ∈π j w(xi)
the i-th cluster center.

According to [START_REF] Dhillon | Kernel k-means, spectral clustering and normalized cuts[END_REF], there is a direct relationship between the trace maximization of the normalized cut in spectral clustering and the kernel k-means problems. Therefore, the objective function represented in (1) can be written as follows:

D({π j } k j=1 ) = trace(W 1 2 φ T φW 1 2 ) -trace(Y T W 1 2 φ T φW 1 2 Y ) (2)
where, W is the diagonal matrix of all the weights. Y is the N * k orthonormal cluster assignment matrix, i.e., Y T Y = I.

Considering that φ T φ is simply the kernel matrix K of the data, and trace(W

1 2 φ T φW 1 2 ) is a constant term, the minimization of the objective function in (2) is equivalent to the minimization of -trace(Y T W 1 2 KW 1 2 Y ).
In situations where computing the spectral clustering algorithm is difficult, the weighted-kernel k-means algorithm is particularly useful as an alternative algorithm.

Graph Convolution Network Overview

A graph is represented by G = (V, E, X) where V is a set of N nodes, and E a set of edges such as an edge (v i , v j ) is a link connecting the nodes v i and v j . X = [x 1 ; x 2 ; ...; x N ], ∈ R N * d is a node feature matrix, where x i ∈ R d denotes a feature vector of node v i [START_REF] Zhang | Spectral embedding network for attributed graph clustering[END_REF]. The structure of a graph G can be represented by two principal matrices: the adjacency matrix A ∈ R N * N where a ij = 1 if there is an edge between nodes v i , v j , and a ij = 0 otherwise, and the similarity matrix S which is a square matrix of N * N and symmetric (s ij = s ji ), s ij represents the weight of the edge (v i , v j ).

The GCN presented in [START_REF] Kipf | Semi-supervised classification with convolutional networks[END_REF] uses a semi-supervised learning approach on data structured in graphs. It handles the challenge of labeling nodes in a network where labels are only known for a limited portion of nodes by estimating labels for unlabeled data, for that it is referred to as semi-supervised learning technique. A GCN model learns hidden layer representations that encode both local graph structure and node characteristics.

The layer-wise propagation rule of the GCN for semi-supervised learning is:

H i = σ( AH i-1 W i ) (3)
H i : output matrix of the i-th inner layer of the GCN, H i ∈ R N * hi with h i the number of features of the inner-layer i. H i-1 : the input of the i-th inner layer of the GCN,

H i-1 ∈ R N * hi-1 .
For the first layer, H i-1 is the features matrix X of the graph. A: indicates the adjacency matrix of the graph. A its normalized matrix calculated as follows:

A = D -1 2 (A + I) D -1 2
where, D = i (A + I) ij , and I the identity matrix. W i : the weight of the layer i, W i ∈ R di-1 * di . σ is an activation function of internal layers such as rectified linear unit (ReLu) [START_REF] Kipf | Semi-supervised classification with convolutional networks[END_REF].

Then, the Softmax activation function, defined in [START_REF] Kipf | Semi-supervised classification with convolutional networks[END_REF] is applied in order to obtain the predicted label matrix F ∈ R N * k where k is the number of classes.

The GCN network is usually trained with the cross-entropy loss function which is typically used in classification problems. In Kejani, Dornaika & Talebi (2020) [START_REF] Kejani | Graph Convolution Networks with manifold regularization for semi-supervised learning[END_REF], multiple regularization loss was added to the cross-entropy loss function to give better clustering results. The loss equation then becomes:

Loss = - l i=1 k c=1 y ic log(F ic ) + λ trace(F T LF ) (4)
where y ic represents the true label of the dataset, F ic is the estimated probability of the i-th sample to be in class c. l denotes the total number of labels samples, λ is a hyper-parameter, and L denotes the Laplacian matrix of the graph. The model parameters W i of GCN can be trained by minimizing Eq. ( 4) with gradient descent.

The Proposed Method

Given a graph G(V, E, X), and its adjacency matrix A, our goal is to partition the nodes of this graph into k different clusters using GCN in an unsupervised way.

Model Architecture

To achieve our goal, we design a three-layer GCN as shown in Fig. 1. We first feed the node feature matrix X ∈ R N * d into the network, then learn the first propagation layer output H 1 ∈ R N * h1 , and H 2 ∈ R N * h2 for the second layer by:

H 1 = ReLu( AXW 1 ) (5)
H 2 = ReLu( AH 1 W 2 ) (6)
where ReLu() is an activation function, A the normalized adjacency matrix of the graph, W 1 ∈ R d * h1 , and W 2 ∈ R h1 * h2 are the trainable weight matrices of the inner layers and h 1 , h 2 the number of features in the hidden layers.

In the output layer H 3 , we map H 2 to a k-dimensional space, where k indicates the number of clusters:

H 3 = H 2 W 3 (7)
W 3 ∈ R h2 * k is a learnable model parameter matrix. Then, to be able to use the output of the GCN in an unsupervised learning strategy, we perform a cholesky decomposition [START_REF] Zhang | Spectral embedding network for attributed graph clustering[END_REF] on (H 3 ) T H 3 , i.e. (H 3 ) T H 3 =Q T Q, where Q ∈ R k * k is a lower triangular matrix. Finally, we obtain the orthogonal form of the GCN output which also represents the soft clustering assignment matrix output of the model:

H = H 3 (Q -1 ) T (8)
After getting the orthogonalized output, our clustering loss function will be applied to it in order to train the model parameters.

Model Loss Function and Clustering

As mentioned before, our main goal is to train the GCN in an unsupervised manner, i.

Loss = -trace(H T KH) + λ trace(H T LH) ( 9
)
where H is the model orthogonalized output, L the Laplacian graph matrix, λ is a hyper-parameter used to balance between the two terms, usely λ ∈ [0, 1]. K is the kernel matrix, in our case, we choose K the Gaussian kernel matrix given by: K = exp(

-xi-xj 2 s 2
), x i , x j ∈ X and s 2 is equal to the average value of the distance between the data points. Using this loss function, the main advantage is that no prior knowledge is needed to train the parameters of our model. Therefore, our GCN model is an unsupervised model. Moreover, the use of kernel k-means loss function allows us also to acquire the advantages of spectral clustering as already mentioned in paragraph Sect. 3.1. The objective function of spectral clustering and kernel k-means are similar. Furthermore, during the feature propagation phase and the learning phase, the model adapts both the structure and the features of the graph.

The parameters of the model W 1 , W 2 , W 3 are trained by minimizing equation ( 9) with gradient descent. After finishing the model training process, we perform k-means on its output H to obtain the partitions of nodes.

With the proposed loss function [START_REF] Yang | Deep spectral clustering using dual autoencoder network[END_REF], two interesting and friendly clustering types of constraints are imposed on the unknown representation matrix H: (i) the kernel K-means loss enforces compactness in the space induced by the kernel matrix (information derived from features) and (ii) the spectral clustering term enforces smoothness of the representation over the graph (information derived from data structures).

Experiments

We evaluate our model on three datasets widely used for the analysis of attributed graphs: Cora, Citeseer, Pubmed [START_REF] Kipf | Variational graph auto-encoders[END_REF].

These datasets are citation networks where the nodes correspond to publications and are connected by an edge if one cites the other. Each node is associated with a feature vector, which represents the content of the document. Informations of these datasets are summarized in Table 1.

Baselines

We compare our methods with two different types of baselines:

The first type is methods that use only node features or graph structure: k-means partition nodes using the graph features matrix [START_REF] Gao | Deep Clustering with Concrete kmeans[END_REF], spectral clustering [START_REF] Von Luxburg | A tutorial on spectral clustering[END_REF] which takes the adjacency matrix of nodes hence the structure of graph.

The second type is methods adopting both node characteristics and graph structure: such as variational graph autoencoder (ARVGE) [START_REF] Pan | Adversarially regularized graph autoencoder for graph embedding[END_REF], and SENet [START_REF] Zhang | Spectral embedding network for attributed graph clustering[END_REF] method which improves graph structure by leveraging the information of shared neighbors, and learns node embedding with the help of a spectral clustering loss.

Evaluation Metrics and Experimental Setups

To assess the clustering performance of our method, we use three performance measures [START_REF] Aggarwal | Data clustering: algorithms and applications[END_REF]: clustering accuracy (Acc), normalized mutual information (NMI), and adjusted rand index (ARI). Acc is the ratio of correctly predicted data. NMI assesses cluster quality by measuring the match between true and predicted labels. ARI measures the separation ability between clusters and the recognition ability for each cluster.

For our model, we build the network with randomly initialized weights with 16 hidden units and two internal layers as in [START_REF] Zhang | Spectral embedding network for attributed graph clustering[END_REF]. We adopt the Adam optimizer [START_REF] Kingma | Adam: a method for stochastic optimization[END_REF] with a 0.001 learning rate and run experiments on Pytorch. The λ hyperparameter in the proposed loss function is fixed to 0.5.

Result Analysis

Tables 2, 3, and 4 show the clustering and standard deviation results obtained after running each method 10 times on each dataset. The best results are in bold. The interpretations are as follows: U-GCN, the algorithm proposed in this paper, outperforms clustering approaches that just use node attributes or graph structure such as k-means and spectral clustering, as it effectively incorporates and captures both kinds of information.

Furthermore, U-GNC improves the results of methods that use the attributes and structure of graph at the same time, since U-GCN adapts the structure and attributes of the graph in the propagation of features in the model and in the training process not as ARVGE as it only uses them both during the propagation phase. In addition, the use of the weighted kernel k-means objective loss function which is similar to the spectral clustering objective loss function facilitates the transmission of information about the cluster structure at each layer.

Conclusion

In this work, we proposed an unsupervised method based on GCN for the partitioning problem in graphs without prior knowledge. Inspired by the mechanism of GCN and the relationship between kernel weighted k-means and spectral clustering, we integrated two objective functions to achieve clustering of nodes by merging structural and feature information. The proposed loss function is suitable for training the GCN model when there is no information available on the labels. This approach is also a new solution for clustering graphs in an unsupervised way. The experimental results demonstrate that the proposed method achieves better performance on three citation networks, compared to state-ofthe-art algorithms. In future work, we aim to study the effect of changing the trade-off hyper-parameter λ.

Fig. 1 .

 1 Fig. 1. The model architecture of the graph convolution neural network.

2 1 2Y

 21 e. the training process does not include any prior information about the ground truth labels. To achieve it, we have merged the rectified version of the kernel k-means objective function given by-trace(Y T W 1 KW ) with the classic spectral loss that use the graph. IN the rectified version of the kernel kmeans loss we replace W by I (Identity matrix) and Y by H the model output. the model output. In this case, our proposed model is trained by minimizing the following loss function:

Table 1 .

 1 Datasets Information

	Datasets Nodes Edges Features Clusters
	Cora	2708 5429 1433	7
	Citeseer	3327 4732 3703	6
	Pubmed 19717 44338 500	3

Table 2 .

 2 Clustering results on Cora

	Methods	ACC	NMI	ARI
	k-means	0.3465 ± 0.0136 0.1673 ± 0.0154 0.0958 ± 0.0107
	Spectral	0.3419 ± 0.0203 0.1949 ± 0.0194 0.0181 ± 0.0186
	ARVGE [11]	0.6380 ± 0.0096 0.4500 ± 0.0083 0.3740 ± 0.0071
	SENet [12]	0.7192 ± 0.0066 0.5508 ± 0.0065 0.4896 ± 0.0109
	U-GCN(this paper) 0.73 ± 0.00872 0.56 ± 0.00569 0.495 ± 0.00859

Table 3 .

 3 Clustering results on Citeseer

	Methods	ACC	NMI	ARI
	k-means	0.3849 ± 0.0237 0.1702 ± 0.0206	0.1243 ± 0.0192
	Spectral	0.2591 ± 0.0109 0.1184 ± 0.0168	0.0012 ± 0.0137
	ARVGE [11]	0.5440 ± 0.0139 0.2610 ± 0.0172	0.2450 ± 0.0164
	SENet [12]	0.6752 ± 0.0075 0.417 ± 0.0082 0.4237 ± 0.0097
	U-GCN(this paper) 0.68 ± 0.00617 0.417 ± 0.00632 0.43 ± 0.00859

Table 4 .

 4 Clustering results on Pubmed

	Methods	ACC	NMI	ARI
	k-means	0.5732 ± 0.0381 0.2912 ± 0.0352	0.2505 ± 0.0346
	Spectral	0.3974 ± 0.0267 0.0346 ± ± 0.0309 0.0018 ± 0.0291
	ARVGE [11]	0.5822 ± 0.0100 0.2062 ± 0.0057	0.2045 ± 0.0065
	SENet [12]	0.6759 ± 0.0062 0.3061 ± 0.0145	0.2966 ± 0.0113

U-GCN(this paper) 0.68 ± 0.00527 0.3105 ± 0.0102 0.31 ± 0.0010 Author Proof

Author Queries

Chapter 3 Query Refs.

Details Required

Author's response AQ1 This is to inform you that corresponding author and email address have been identified as per the information available in the Copyright form.

AQ2

This is to inform you that as the Institutional email address of the corresponding author is not available in the manuscript, we are displaying the private email address in the PDF and SpringerLink. Do you agree with the inclusion of your private e-mail address in the final publication?


Falko Ziebert

Igor M Kulić

How Influenza's Spike Motor Works

While often believed to be a passive agent that merely exploits its host's metabolism, influenza virus has recently been shown to actively move across glycan-coated surfaces. This form of enzymatically driven surface motility is currently not well understood and has been loosely linked to burnt-bridge Brownian ratchet mechanisms. Starting from known properties of influenza's spike proteins, we develop a physical model that quantitatively describes the observed motility. It predicts a collectively emerging dynamics of spike proteins and surface bound ligands that combined with the virus' geometry give rise to a self-organized rolling propulsion. We show that in contrast to a Brownian ratchet, the rotary spike drive is not fluctuation driven but operates optimally as a macroscopic engine in the deterministic regime. The mechanism also applies to relatives of influenza and to man-made analogues like DNA-monowheels and should give guidelines for their optimization.

One of humanity's greatest inventions is the wheel. While reflecting about why nature overlooked wheeled propulsion, it caught us by surprise that the wheel was rolling in nature eons ago: the common influenza virus uses its whole capsid as a motorized surface rolling machine [1,2], see Fig. 1A. The reason for this fundamental discovery of Sakai et al. staying almost unnoticed by a broader audience (with few exceptions [3,4]) is possibly rooted in our lack of understanding of its underlying physical mechanism.

Being such an omnipresent molecular adversary, the influenza virus (IV) and its proteins have been extensively characterized [5][6][7][8]. The two spike proteins responsible for IV A's interaction with the host membrane are Hemaglutinin (HA) and Neuraminidase (NA), see Fig. 1B. While in certain influenza subtypes (like IV C) these two proteins are fused together [9], in general they are distinct ∼ 10 nm sized entities performing two distinct and mutually competing functions: HA binds sialic acid residues of glyco-peptides and lipids coating our cells (Fig. 1B,C). NA in turn, acts antagonistically and degrades the contacts with the glycan substrate by hydrolytically cutting the same sialic acid residue that HA binds to, see Fig. 1D. The residue, located at the very ends of branched glycans, can be either bound by one HA or by one NA molecule, but for steric reasons not to both at the same moment. Inhibition of HA abolishes virus binding to glycans (as exploited for diagnosis, e.g., in the classical hemagglutination assay [START_REF] Spackman | Hemagglutination Inhibition Assay[END_REF]) while NA inhibition abolishes its motility [2,3].

While most textbooks depict influenza as a spheroidal virus, its aspect ratio is in fact highly polymorphic and during infection of human hosts the majority of the virus mass comes in filamentous form [11,12]. The reason for their filamentous shape with lengths from 1-300 µm [11] is debated since most evolutionary arguments favor the sphere (e.g. volume to surface ratio and stability [13] or uptake dynamics [14]). Although previously speculated [15] that elongated viruses could self-propel like manmade actively rolling fibers [16], it was only the work of Sakai et al. [2] that found the direct evidence and suggests that the elongated form is in fact advantageous for robust directionally persistent motion.

Beyond IV and related viruses, interestingly DNA nanotechnology has developed synthetic rollers termed DNA-RNA monowheels [17][18][19], that use a similar de- sign principle -namely linkers and their digestion. So how does influenza, and the motif of "bridging and cutting" in general, generate the force necessary for the rolling motion? Starting from basic known aspects we show here that it is not a simple burnt-bridge fluctuation driven mechanism and that virus motion is deterministic in nature, macroscopically robust and in fact close-to inevitable.

The mechanism. Consider the interface where the viral capsid and the glycan coated substrate surface meet. In this nanoscopic region, glycan chains at a high concentration G 0 (well in excess to spike proteins, see the estimates below), are constantly binding to and unbinding from the HA proteins and in turn elongating to a length l in that process. Once bound (with dissociation constant K d ) they gain a free energy k B T ln G0 K d . On the other hand they pay a stretching energy cost E el (φ) = S 2 l 2 ≈ R 2 S 8 φ 4 which, due to the curvature of the capsid, depends on the angle φ measured from the virus symmetry axis and its radius R ( 50 nm for IV), see Fig. 2A top panel. Chains bound to NA are short lived and neglected here for simplicity. We consider the glycan chain as an ideal linear spring with spring constant S ∼ 0.01-1 k B T /nm 2 (a typical range for polymers of few nm length). The balance of the two energy terms then sets the angular size

φ c = ln G0 K d 8kT SR 2 1/4 of the contact zone φ ∈ [-φ c , φ c ].
The stretching force F el = -∂ ∂l E el resulting from a single stretched linker gives rise to a torque ∝ SR 2 (1 -cos φ) sin φ ≈ 1 2 SR 2 φ 3 . The linkers have an angular density ρ HA b (φ) given by the product of the angular density of HA spikes, ρ HA , and the angular probability density of each linker being bound, b (φ), and the total torque acting on the capsid is the integral over all bound linkers

m = -m 0 +φc -φc b (φ) φ 3 dφ = 0 . (1)
Here m 0 = 1 2 SR 2 ρ HA is the characteristic torque scale. At typical densities of linkers and typical angular speeds ω ∼ 1 s -1 [2], the linker torque dominates all other torques acting on the virus including hydrodynamic dissipation. Therefore the torque balance m = 0 holds. Denoting the concentration of HA spikes on the virus with H 0 and the initial concentration of glycans on the (cell's) surface as G 0 , we have to determine the evolution of the concentration of bound HA-glycan links B(φ, t) and the glycan concentration G(φ, t) both as functions of time t and the angle φ ∈ [-φ c , φ c ]. In addition to this binding kinetics, the NA spike enzyme progressively digests the glycan in its vicinity with a catalytic velocity V cut and a Michaelis constant K M . Combining these effects and assuming the virus to roll with angular velocity ω we have

∂ t B + ω∂ φ B = k on G (H 0 -B) -k of f B (2)
∂ t G + ω∂ φ G = -k on G (H 0 -B) + k of f B - V cut G K M + G .(3)
Here the terms ∝ ω on the l.h.s. represent the advection of concentrations in the virus fixed-frame due to its rotation. The first terms on the r.h.s. are the on/off-kinetics of glycan binding, with the kinetic constants satisfying

k of f kon = K d .
Although the off-rate is stretching force- [20] and hence angle-dependent, we neglect this effect assuming a small size of the contact interval φ 1 where the elastic energy dependence E el ∝ φ 4 is weak.

Finally, the last term in Eq. ( 3) represents the Michaelis-Menten-like degradation of the free glycans by NA, with a velocity V cut = k cut N set by the cutting rate k cut , see Fig. 1D, and the enzyme concentration N . These two equations, together with b (φ) = B/H 0 satisfying the torque balance, Eq. ( 1), completely determine the dynamics and the question is now, whether the enzymatic activity can sustain solutions with non-zero ω.

Passive frictional torque. In a fist step, it is instructive to consider the passive case, i.e. in absence of catalytic activity (V cut = 0), and to assume that the virus is forced by a weak external force/torque to roll with a steady-state angular velocity ω [21]. Related situations have been investigated when modeling cells rolling in external shear flow [22,23] and contraction/sliding motion induced by stochastic linkers [24,25]. In the steady state, Eqs. (2, 3) imply the conservation law (B + G) = 0 or G = G 0 -B (for a homogeneous G 0 ), allowing us to reduce the problem to

ω B = k on (G 0 -B) (H 0 -B) -k of f B. (4)
For the initial condition B(-φ 0 ) = 0 (rolling to the left) the exact solution is given by

B (φ) = C 0 -C 1 2 - C 1 C0+C1 C0-C1 e C 1 kon ω (φ+φ0) -1 (5)
with constants

C 0 = H 0 +G 0 +K d , C 1 = C 2 0 -4H 0 G 0 .
This profile, cf. the sketch in the bottom panel of Fig. 2A, implies an increase of the bound HA-glycan links leveling to a plateau B pl = C0-C1 2 [26]. For simplicity, we approximate the exact profile by two lines: first, in the region of its rapid increase, B is approximated by the slope at the front: for φ ∈ [-φ c , φ pl ], B (φ) = α ω (φ + φ c ) with φ pl the angle where in this approximation the plateau is reached and the linker binding velocity

α = k on H 0 G 0 . (6)
Note that the faster the rolling, the shallower the spatial gradient becomes, since the build-up of HA-glycan links needs time. And second, for φ ∈ [φ pl , φ c ] we approximate B by its plateau value, i.e. B (φ) = α ω (φ pl + φ c ) = B pl . With this slope-plateau approximation, for slow rotation (implying a steep increase of B to the plateau), we can evaluate Eq. ( 1) to get the torque-velocity relation for passive rolling

m diss (ω) = -ξ diss ω , ξ diss = m 0 H 0 B 2 pl α φ 3 c 2 . ( 7
)
Hence this is a frictional torque, acting against the motion and linear in ω. The friction constant, ξ diss , is determined by both the slope and the plateau of the distribution of B as well as the size of the contact interval, which themselves contain all system parameters. Enzyme activity induces active torque. The effect of the NA activity can be captured perturbatively and leads to an additional contribution to the torque. If the enzyme activity is slow compared to the binding kinetics, = V cut /α is a small parameter. Expanding B = B (0) + B (1) and G = G (0) + G (1) in powers of yields the leading order correction

B (1) = -αf φ+φc ω with f = H 0 -B pl K d + G pl G pl K M + G pl (8)
a dimensionless ratio of all concentrations/kinetic constants and G pl = G 0 -B pl the plateau of the glycan distribution (in the passive case). The enzymatic activity hence leads to a negative slope β = V cut f instead of the plateau, cf. Fig. 2B, and insertion into Eq. ( 1) yields the active torque in "two-slope" approximation:

m act = p act ω , p act = m 0 H 0 f 2φ 5 c 5 V cut , (9)
where p act is the power injected by NA operation. The active torque is positive (since B pl < H 0 ), it is proportional to V cut and has a 1/ω dependence, unlike the passive one which is linear in ω.

Taken together, the passive and active torques yield the torque balance m diss + m act = 0 = -ξ diss ω + p act /ω, implying a pitchfork bifurcation for the steady-state rolling velocity

ω = ± p act ξ diss ∝ φ c f √ αV cut B pl . ( 10
)
Note that the torque scale, m 0 = 1 2 SR 2 ρ HA , cancels out, but the parameters S, R are still present due to the dependence in φ c .

We can compare to the experiments [2] by inserting typical parameters for IVs: using R = 50 nm and S 0.1 k B T /nm 2 implies φ c 0.5. Typical concentrations are G 0 = 10 mM, H 0 = 1-5 (we use 2) mM [START_REF]The volume concentrations G0, H0 and N (in M = mol/l) can be estimated from typical volumes per molecule V = Aδ, with A ∼ 100 -500 nm 2 (for glycans and spike proteins respectively) the area per molecule and δ ∼ 1 nm its typical height fluctuation from the surface[END_REF]. The HA on-off kinetics has been characterized [START_REF] Sauter | [END_REF]29] yielding K d = 1-5 (2) mM, k of f = 10 -1 -1 (1) s -1 , k on = 0.01-1 (0.5) mM -1 s -1 and NA's enzymatic activity [30] to yield K M = 14.3 mM, k cat 15 s -1 , implying with a typical NA concentration of N = 1 mM a V cut = k cat N = 15 mMs -1 . Using these values in Eq. ( 10) we get ω 0.4 s -1 , which compares well to Ref. [2] where virus speeds of v = 10 -30 nm/s were reported, corresponding to angular velocities ω = v R 0.2 -0.6 s -1 . Numerical study including stochasticity. To scrutinize the robustness of the mechanism, so far described on the continuum level via concentrations, we implemented the stochastic reaction kinetics using the Gillespie algorithm [31]. For the latter, the virus cross-section was assumed to present a number of n vir discrete binding sites per (angular) contact interval [-φ c , φ c ]. Larger n vir correspond to a more elongated virus with more linkers per angle, with n vir → ∞ being the deterministic limit. The virus position is updated in each step in accordance to the vanishing torque condition. Fig. 3 compares the approximate two-slope theory (black curve) to a numerical solution of the continuum model Eqs. (2, 3) (red) and stochastic simulations (symbols). Shown is the angular frequency ω as a function of the NA concentration. The inset displays the profiles B(φ) for stationary rolling in the three cases. While in the stochastic implementation the virus inverts its rolling direction occasionally, from Fig. 3 it is evident that the mechanism is robust against finite number of binding site effects. It nevertheless works best (i.e. rolling is fastest) for the "macroscopic" i.e. continuum case, in contrast to the classical burnt-bridge mechanism, as discussed below.

Fig. 4 shows stochastic simulations highlighting the robustness of the mechanism against perturbations in the glycan distribution on the surface: the upper panel in A shows a snapshot of a virus rolling over a surface displaying small, almost glycan-depleted regions and the lower Rolling velocity and bound linker profile in the continuum model vs. the stochastic implementation. Shown is the rolling frequency ω as a function of enzymatic activity, i.e. concentration of NA, comparing the approximate theory (black), a numerical solution of the continuum model (red) and the stochastic implementation (symbols) for different numbers of discrete binding sites, nvir. The inset shows the profiles B(φ) for stationary rolling in the three cases. Here N = 0.1 mM and the stochastic profile was averaged over 50000 Gillespie moves.

panel the trajectory. Although the virus is slowed down in the depleted zones (cf. the green lines in the trajectory panel), its motion persists. Fig. 4B displays two superimposed snapshots of a virus moving on a modulated glycan distribution, the lower panel again showing the trajectory. This proves that the mechanism allows rolling on both uphill and downhill glycan gradients.

The generic double-gradient mechanism employed by the IV is the interplay of binding and digestion of linkers, the enzymatic reaction transforming the plateau of the linker distribution (cf. Fig. 2A) into a negative slope (cf. Fig. 2B). One can imagine other reaction pathways to result into a similar generic "polarization" with two different slopes for the bound linker distribution, and in fact this was recently realized [17]:

A synthetic relative: the DNA-RNA monowheel. Interestingly, without knowing about the mechanism of IV rolling, a synthetic variant of the generic mechanism was recently implemented using DNA nanotechnology [17,18]. There, the surface of silica particles was covered with DNA sequences that form heteroduplexes (of size ∼ 15 base pairs) with complementary RNA strands coated on a surface. Motion of the particles (termed monowheels in [17]) was initiated by adding RNase H, which selectively hydrolyses the hybridized RNA (the bound linkers, B in our notation), but not single-stranded RNAs (i.e. free linkers, G).

This variant can be easily cast into our theoretical framework: as the enzyme RNase H does not destroy the free linkers but instead the bound ones, the enzymatic term (∝ V cut ) is absent in Eq. ( 3), but its ana-

logue -VcutB

KM+B has to be added to Eq. ( 2). This only slightly modified model can be analyzed along the same lines leading to an active torque as in Eq. ( 9) but with

f = B pl K M + B pl . (11)
We can again insert numbers, namely R = 2.5 µm [17], S 0.1 k B T /nm 2 implying φ c 0.1, and for the binding parameters G 0 = 1-2 mM, H 0 = 2-4 mM (estimated from [17]), K d = 10 -5 mM [18,32], k of f = 10 -5 s -1 [32,[START_REF] Strunz | Proc. Natl. Acad. Sci[END_REF], implying k on = 1 mM -1 s -1 . RNase H kinetics is known [START_REF] Fang | [END_REF]: K M = 0.3 • 10 -3 mM, k cat 1 s -1 and its typical concentration used in [17]

is R H = 1.44 • 10 -4 mM, implying V cut = k cat R H = 1.4 • 10 -4 mMs -1 .
For these values we get B pl K M hence f 1 and, as the passive torque has the same dependence as before,

ω = ± 2 √ 5 √ αV cut B pl φ c . (12)
Inserting numbers yields ω 10 -2 s -1 which once again fits well to the experimentally observed velocity of 30 nm/s [17] implying ω 10 -2 s -1 .

Discussion. The double-gradient mechanism described here is very robust and gives rise to large propulsion speeds, ω ∝ √ V cut , even for weak enzymatic activity. Importantly, it is not a simple bridge-burning as recently hypothesized for both IV and the DNA-wheel [2,18]. Although bridge-burning can be operative for certain IV A strains that show phase separation of NA and HA spikes [35], these lack rolling and move much slower than reported in [1,2]. In burnt-bridge Brownian models [36,37] the random walker destroys the "bridges" it walks on, and only the prohibited back-stepping leads to the directed motion implying self-avoiding paths. Such models, like other members of the class of Brownian motors [38,39], are inherently fluctuation driven and increasing the angular density of linkers implies a slow-down of the motion. In contrast, the mechanism discussed here relies on the self-organized, "internal" polarization of the linker distribution within the contact zone and works optimally in the macroscopic regime. It is robust (cf. Fig. 4) but not self-avoiding/unidirectional, since the roller can run in reverse direction even if its trail (the glycan or RNA distribution) is substantially depleted behind. It responds to existing surface gradients, but much less than a burntbridge walker, giving the virus -controlling its contact zone -a higher motile autonomy in the evolutionary race with its host, controlling the rest of the substrate. How the "delicate balance" [4] of NA vs. HA and its adaptation orchestrates the mechanism in detail should be explored more in the future.

In conclusion, it seems that with influenza we are facing an underestimated, smart adversary that in contrast to classical virology dogmas displays a "metabolism" at its interface, providing it with a motile organelle (engulfing its whole body) that emerges from geometry and the self-organization of its spike proteins. The mechanism should also apply to many relatives of influenza bearing enzymatic spike proteins, including toro-virus and some beta-corona-viruses [40]. However, we can flip the coin and turn the adversary into an ally, by learning from its workings: The triplet of a binding molecule, a cutting enzyme and the spherical or cylindrical geometry should be abundant and was in fact already used to propel DNAcoated beads [17,18], which could now be optimized using the understanding of the mechanism.

FIG. 1 .

 1 FIG.1. Rolling Influenza, its surface structure and activity. (A) Superimposed snapshots of a rolling elongated IV C (taken from[2]), the green arrows indicating the rolling direction. (B) An IV cross-section showing its surface covered with two kinds of spike proteins: Hemaglutinin (HA, blue) and Neuraminidase (NA, red). The substrate (cell's surface) is covered with glycans exposing a sialic acid residue (green, both not to scale). (C) HA binds to and unbinds from glycan via the sialic acid with rates kon, k of f . (D) NA transiently binds, with a Michaelis constant KM = (k-1 + kcut)/k1, and hydrolytically cuts the sialic residues with a rate kcut, making the glycans inactive for HA binding.

  arXiv:2010.02182v1 [cond-mat.soft] 5 Oct 2020

FIG. 2 .

 2 FIG. 2. The double-gradient mechanism of motion, maintained by the distribution of bound HA-glycans. (A) Top: A virus rolling with constant angular frequency ω due to an externally applied force fext. Bottom: The distribution of bound HA-glycans, B(φ), within the contact interval [-φc, φc] has two regions: a sharp increase (with slope α/ω) in the rolling direction followed by a plateau B = B pl . (B) Top: self-rolling due to the enzymatic activity of NA, cutting away the sialic acid residues. Bottom: In this case, the distribution B(φ) has a negative slope (β/ω) instead of a plateau in the second region, i.e. at the rear.

  FIG. 3.Rolling velocity and bound linker profile in the continuum model vs. the stochastic implementation. Shown is the rolling frequency ω as a function of enzymatic activity, i.e. concentration of NA, comparing the approximate theory (black), a numerical solution of the continuum model (red) and the stochastic implementation (symbols) for different numbers of discrete binding sites, nvir. The inset shows the profiles B(φ) for stationary rolling in the three cases. Here N = 0.1 mM and the stochastic profile was averaged over 50000 Gillespie moves.

FIG. 4 .

 4 FIG. 4. Stochastic simulations for modulated glycan distributions. (A)The upper panel shows a snapshot of the bound linker (in blue) and the glycan (green) distributions and the lower panel the trajectory. The surface-bound glycans were locally depleted (marked also in the trajectory by the green vertical lines), but the rolling motion persists. (B) Two superimposed snapshots of a virus rolling on a surface presenting uphill and downhill glycan gradients and the corresponding trajectory below, again showing the robustness of the rolling mechanism.

Acknowledgements. We thank Jens-Uwe Sommer, Ulrich Schwarz and Felix Frey for discussions. IMK thanks the Leibniz-Institute IPF Dresden for hospitality.


Yi Sun

Sihong Liu

Yun Jia
email: yun.jia@polytech-lille.fr

Anlong Zhu

Liujiang Wang

Wen He

A cross-scale finite element analysis of concrete-faced rockfill dam 1

Keywords: concrete-faced rockfill dam, cross-scale finite element method, arbitrary-node hexahedron element, discontinuous mesh

In this paper, a new cross-scale finite element method is developed to study the mechanical response of concrete-faced rockfill dams, in which an arbitrary-node hexahedron element is used to deal with the mesh transition between the concrete slab and dam body. By using the proposed method, the mesh sizes of concrete slab and rockfill may be very different. Moreover, since the cross-scale method doesn't change the finite element calculation framework, the proposed method is easy to use. In the simulated CFRD, the mechanical behavior of the concrete slab is described by an elastic-plastic constitutive model, and a nonlinear elastic model is applied to the rockfill materials. The computational accuracy and performance of the proposed method are analyzed. The numerical results and analysis provides an effective approach to study the local behavior of concrete face slabs.

Introduction

The concrete-faced rockfill dam (CFRD) is one of the most common dam types. Particularly, due to its advantages in reliability, economic cost, and environmental adaptability, the concrete-faced rockfill dam is becoming an increasingly attractive type of dam [START_REF] Fell | Geotechnical engineering of dams[END_REF][START_REF] Cruz | Concrete face rockfill dams[END_REF][START_REF] Kan | Application of advanced bounding surface plasticity model in static and seismic analyses of Zipingpu Dam[END_REF]. As the key seepage-proof structure of a CFRD, the safety and reliable operation of the face slab remains the major research topic of the dam engineering community. In the classical finite element simulation, the concrete slab is simulated using coarse mesh as the dam body which is composed of rockfill materials. Since the concrete face slab is a thin-slab structure, the element dimensions on the slab surface will be much larger than that along the thickness direction. As a result, it is difficult to guarantee the reliability of numerical results. In practice, it is generally believed that the size of slab mesh should be controlled within a few meters. However, the size of the slab mesh should be controlled within a few centimeters to study the local failure and crack propagation of the slab [START_REF] Wei | Application of displacement multi-point constraint refinement method in simulation of concrete-faced rockfill dams[END_REF]. However, a fine mesh of the slab leads to a significant increase in the element number in the dam body, and the computational costs are inevitably increased. Given this, numerical estimation of the mechanical behavior of CFRDs requires developing a numerical method to deal with the mesh transition between the concrete slab and the rockfill body.

In practice, the mesh transition between the dam body and the face slab is realized by gradually increasing the element density in the dam body near the face slab [START_REF] Xu | Dynamic damage evaluation on the slabs of the concrete faced rockfill dam with the plastic-damage model[END_REF]. However, the previous method increases the modeling complexity and difficulty while the numerical accuracy may be also disturbed by the eventual creation of distorted elements [START_REF] Wang | Finite element method[END_REF]. Moreover, although the transition method has been successfully applied in two-dimensional conditions, it is still difficult to be used in three-dimensional conditions [START_REF] Qu | A novel interface element with asymmetric nodes and its application on concrete-faced rockfill dam[END_REF] due to the presence of complex element shapes in the transition zone. Therefore, developing a reasonable and efficient method for the mesh transition between the face slab and rockfill mass is necessary for the numerical prediction of deformation behavior of CFRDs. So far, some scholars have applied transition methods to the refined simulation of the face slab. Zhou et al. (2016a[START_REF] Zhou | Numerical evaluation of soft inter-slab joint in concrete-faced rockfill dam with dual mortar finite element method[END_REF] used the dual mortar approach to realize the load transfer between the fine mesh of the concrete slab and the coarse mesh of the dam body. Nevertheless, this method is a contact mechanics approach, including important nonlinear characteristics. This method presents the high convergence insufficiency in the application of engineering structures with complex geometry and under seismic loading conditions. As a result, the previous method can't be used to study the deformation developed at the interface between the face slab and dam body under seismic loading conditions. After that, [START_REF] Qu | A novel interface element with asymmetric nodes and its application on concrete-faced rockfill dam[END_REF] developed a plane interface element with asymmetric nodes to describe the interface characteristics. However, the developed interface element is not suitable for local mesh optimization and three-dimensional problem. Recently, [START_REF] Wei | Application of displacement multi-point constraint refinement method in simulation of concrete-faced rockfill dams[END_REF] applied the multi-point constraint method to study a concrete-faced rockfill dam with a refined mesh for the face slab. [START_REF] Chen | Global concurrent cross-scale nonlinear analysis approach of complex CFRD systems considering dynamic impervious panel-rockfill material-foundation interactions[END_REF], [START_REF] Xu | A nonlinear analysis of dynamic interactions of CFRD-compressible reservoir system based on FEM-SBFEM[END_REF], [START_REF] Gong | The simulation of high compressive stress and extrusion phenomenon for concrete face slabs in CFRDs under strong seismic loads[END_REF] and [START_REF] Qu | Three-dimensional refined analysis of seismic cracking and anti-seismic measures performance of concrete face slab in CFRDs[END_REF] applied the polyhedral scaled boundary finite element method to connect the fine slab mesh with the coarse rockfill mesh. Using the above method, the element number is significantly decreased and detailed analyses of the concrete slab can be also obtained. However, due to the complexity of such a semi-analytical numerical method, this method has only been applied in several concrete-faced rockfill dam simulations.

In view of this, a cross-scale finite element method is proposed in this study to improve the accuracy and efficiency of the numerical simulation of a concrete-faced rockfill dam. In this method, the arbitrary-node hexahedron element is used as the transition zone between two different mesh: the fine mesh of the concrete slab and the coarse mesh of the dam body. In the numerical simulation, the mechanical behaviors of rockfill materials and concretes are described by a nonlinear elastic model and 4 an elasto-plastic model, respectively. Comparing the numerical results of the proposed method and the classical finite element method, it is observed that the load transfer of discontinuous mesh is precisely predicted by the proposed method. The obtained 3D FE numerical results of the CFRD indicate that the cross-scale finite element method is capable of satisfactorily reproducing the stress gradient and the distribution of plastic hardening in the face slab. Thanks to its ease of use, the proposed method can be used in engineering structures with complex geometry and different material components.

Cross-scale finite element method

Arbitrary-node hexahedron element

In the finite element calculation, eight-node hexahedron elements are used in this study. The mesh transition between the fine element of the face slab and the relatively coarse element of the rockfill mass is performed via the arbitrary-node hexahedral element [START_REF] Wang | A multiscale extended finite element method for modeling three-dimensional crack problems[END_REF], which has also been applied in three-dimensional crack problems [START_REF] Wang | Finite element method[END_REF]Yu, 2014,2016). To the best knowledge of the authors, it is the first time that the arbitrary-node hexahedral element is applied in the numerical analysis of CFRDs. With the application of the arbitrary-node hexahedral element, the cross-scale numerical simulation of a concrete-faced rockfill dam can be easily realized in the finite element calculation framework. As the displacement mode of the arbitrary-node hexahedral element is quite different from that of the eight-node hexahedron element, the arbitrary-node hexahedral element is presented in the following.

Under three-dimensional conditions, for a given real displacement field () u ξ , its approximated displacement h () u ξ is given as follows:

    1 () M h T i i i     u ξ N ξ u a p ξ (1)
where M represents the number of sampling points in the point interpolation. For a given node i, i N is the shape function and its nodal displacement i u is a column vector with the size of 31  .   (2)

where

  T 1M ...  U u u
is the nodal displacement matrix, and

  1M ...  q p p
is the polynomial basis matrix at the node.

Substituting Eq. (2) into Eq. (1), the shape functions are then obtained:

    -1 T 1 M N ... N  qpξ (3)
According to Eq. (3), the shape functions of arbitrary-node hexahedral elements can be obtained via the point interpolation of a series of special bases. Considering a set of piecewise linear bases, the arbitrary-node hexahedral element with seven types of nodes is given in Fig. 1. The first type of node is eight nodes at corners; the second type of node is j additional nodes on the edge of The polynomial basis function of the arbitrary-node hexahedral element is expressed as

 ( ) 1 , , , , , , , ,
             p ξ 9 9 9 ( ( ) ) ( ( ) ) , , s ig n s ig n          L 8 8 8 ( ( ) ) ( ( ) ) , j j j s ig n s ig n             8 1 8 1 8 1 ( ( ) ) ( ( ) ) , , j j j s ig n s ig n                L 8 8 8 ( ( ) ) ( ( ) ) , j k j k j k s ig n s ig n                8 1 8 1 8 1 ( ( ) ) ( ( ) ) , , j k j k j k s ig n s ig n                   L 8 8 8 ( ( ) ) ( ( ) ) , j k l j k l j k l s ig n s ig n                   8 1 8 1 8 1 ( ( ) ) , , j k l j k l j k l s ig n                      L 8 8 8 ( ( ) ) , j k l p j k l p j k l p s ig n                      8 1 8 1 8 1 ( ( ) ) , , j k l p j k l p j k l p s ig n                         L 8 8 8 ( ( ) ) , j k l p q j k l p q j k l p q s ig n                         8 1 8 1 8 1 ( ( ) ) , , j k l p q j k l p q j k l p q s ig n                            L 8 8 8 ( ( ) ) j k l p q r j k l p q r j k l p q r s ig n                              (4)
The above polynomial basis is illustrated in detail as follows. When the element contains only the first two nodes, the polynomial basis function of the element becomes

 9 9 9 ( ) 1 , , , , , , , , ( ( ))
( ( )) , , s ig n s ig n                       p ξ L 8 8 8 ( ( ) ) ( ( ) ) , , ( ( ) ) ( ( )
)

i i i j j j s ig n s ig n s ig n s ig n                        L (5) where i  , i  , i 
is the local coordinate of the second type of nodes.

( ( ))( ( )) i i i s ig n s ig n         
is the basis added in the function ()

p ξ
when the second type of nodes appears in the element. The basis is written as

( 1)( 1) i       
when the second type of nodes is only located on the edge of =1

 , =1  , 1
  , and is equal to zero on other edges of 1   .

Therefore, under the influence of the basis ( 1)( 1

) i       
, the displacement distribution is piecewise linear along the edge of =1

) i       
is discontinuous, it is necessary to modify the numerical integration.

To overcome the discontinuity of the differential in the numerical integration, the 2 × 2 Gauss integration per subdomain is used [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Lim | Variable-node element families for mesh connection and adaptive mesh computation[END_REF]. According to additional nodes, the element is divided into some eight-node hexahedron sub-elements, and the shape functions still are trilinear interpolation within a sub-element despite the slope discontinuity at the element level.

Numerical implementation

To implement the arbitrary-node hexahedral element in the finite element code, an important issue is to determine the local coordinates of the additional nodes. In Fig. 2, a mesh composed of two elements (called A and B) is established using a uniform scale mesh. On the contact surface of element A and B, ijmn denotes four nodes of this contact surface. After that, element B is divided into four elements, denoted as B1-B4. Meanwhile, hanging nodes appear on the ijmn surface of element A. With the combination of element A and hanging nodes, an arbitrary-node hexahedral element is created. In order to calculate a polynomial basis function of the arbitrary-node hexahedral element, it is necessary to determine the local coordinates of the hanging nodes. Since the hanging nodes are all on the ijmn surface of element A, the local coordinates in the direction of the hanging nodes are equal to either 1 or -1.

Consequently, only the local coordinates in the two other directions of the hanging nodes are needed to be determined. For a given hanging node k on the ijmn surface of element A, the process of determining its location is illustrated in Fig. 3.

Constitutive models

Nonlinear elastic model for rockfill

For a given CFRD, the dam body is described by a nonlinear elastic model proposed by [START_REF] Liu | Practical nonlinear constitutive model for rockfill materials with application to rockfill dam[END_REF]. This constitutive model is capable of describing the dilatancy behavior, the strength nonlinearity of rockfill materials, as well as the influence of the intermediate principal stress. The model has been applied in the 3D FE study of a real CFRD during the construction and impounding phases. The valley where the dam is located is basically symmetrical with slope angles about 4 0 ~4 5  . The riverbed is relatively flat and the dam foundation has 20~30m thick alluvium of sand gravels. The calculated deformation of the dam body and the deflection of the concrete face slabs are in good agreement with the in-situ measurements, and the performance of proposed model is verified. A detailed description of the model can be found in the work by [START_REF] Liu | Practical nonlinear constitutive model for rockfill materials with application to rockfill dam[END_REF]. A rapid review of the above-mentioned constitutive model is given as follows.

The stress-strain relationship of the model in mean stress p-deviatoric stress q space is given by

v v p v q s s p s q d p d q d d d KJ d p d q d d d JG                      (6)
where K is the bulk modulus, G is the shear modulus and J is the coupling modulus. The bulk modulus K, representing the volumetric stiffness with respect to dp, can be expressed as

1 n ba a p K K P P     (7)
where a P is the atmospheric pressure; The shear modulus G, which controls shear strain with respect to dq, and the coupling modulus J, describing the coupling between the volumetric strain produced by an increment dq and the shear strain produced by an increment dp, are given as

2 2 TC TC KG G K D K G       (8) TC TC KG J D K G       (
/ f f u R q q  in which f q
is the shear failure strength and u q is the deviatoric stress in the curve s q   when s  approaches to the infinite strain. In the proposed model, the SMP criterion is adopted and the stress , pq %% are obtained in the new stress transformed stress space. With the introduction of transformed stress tensor ij % [START_REF] Yao | A unified constitutive model for both clay and sand with hardening parameter independent on stress path[END_REF][START_REF] Yao | UH model: three-dimensional unified hardening model for overconsolidated clays[END_REF], the SMP criterion [START_REF] Matsuoka | On the significance of the "spatial mobilized plane[END_REF]) can be represented as an Extended Mises type criterion in the new stress transformed stress space. The transformed stress tensor ij % is given as

    1 1 2 3 1 2 3 () 2 3 / 9 1 c ij ij ij ij c q pp q I q I I I I I I                   % ( 12
)
where I 1 , I 2 and I 3 are the first, second and third stress invariants, respectively.

Elastoplastic model for concrete

In the literature, for ease of use, the linear elastic model is widely used for concrete slabs [START_REF] Zhang | Time-dependent deformation in high concrete-faced rockfill dam and separation between concrete face slab and cushion layer[END_REF][START_REF] Seo | Behavior of concrete-faced rockfill dams during initial impoundment[END_REF][START_REF] Dakoulas | Longitudinal vibrations of tall concrete faced rockfill dams in narrow canyons[END_REF]Wang et al., 2014). However, when the linear elastic model is used to simulate concrete slabs, the obtained tensile stress often exceeds the tensile strength of concrete. Inspired by the research work of [START_REF] Yang | An experimental and numerical investigation of the mechanical behaviour of a concrete and of its permeability under deviatoric loading[END_REF], an elastoplastic model is used for concrete slabs. The plastic yield function of concrete is given as

  0 s f q A C p      ( 13
)
where A is the slope of linear failure surface in the p-q stress space; s C is the cohesion coefficient of concrete;  describes the plastic hardening of concrete and is a function of the current cumulated equivalent plastic deviatoric strain p  :

  00 1 p p B          (14)
It varies from the initial yield threshold 0  to the ultimate failure state 1   . The parameter B controls the rate of plastic hardening.

To describe the volumetric transition between dilatancy and compaction, a non-associated plastic flow rule is adopted:

 ln 0 s cs pC g q A p C p          (15)
where p corresponds to the intersection point between the axis p and the plastic potential surface, and the parameter c  controls the transition boundary between the compaction and dilatancy domain. The model's parameters can be identified by a uniaxial tension test and a uniaxial/triaxial compression test. A detailed description of constitutive model can be found in the work by [START_REF] Yang | An experimental and numerical investigation of the mechanical behaviour of a concrete and of its permeability under deviatoric loading[END_REF].

Numerical verification

To verify the computational accuracy and rationality of the cross-scale finite element method, the dam section of the Tankeng CFRD is simulated, and the obtained numerical results will be compared

FE analysis model

In Fig. 5, two numerical schemes of Tankeng CFRD are presented. In scheme 1, the mesh dimensions of the face slab and rockfill mass are the same and the mesh of CFRD contains 1290 eight-node hexahedron elements, including 44 slab elements (Fig. 5a). In scheme 2, the mesh of face slab and cushion in scheme 1 are further refined, and a mesh discontinuity appears in the CFRD.

Meanwhile, the cross-scale finite element method is adopted to connect two different-size meshs of the face slab and dam body while the arbitrary-node hexahedron element is used for elements in the transition zone. The final mesh of scheme 2 is composed of 1759 elements (Fig. 5b), in which 176 elements are applied for the face slab. In two schemes, the length along the dam axis is both 12m, which is equal to the width of a single face slab. The dam's axial direction is noted as x direction while the direction along the river denotes y direction and the vertical direction represents z direction. For the boundary conditions, all the displacements at the bottom of the foundation and the axial displacement on both sides of the dam are blocked. 1 and2 [START_REF] Liu | Practical nonlinear constitutive model for rockfill materials with application to rockfill dam[END_REF][START_REF] Yang | An experimental and numerical investigation of the mechanical behaviour of a concrete and of its permeability under deviatoric loading[END_REF]. The interface between the concrete face slab and the cushion is described using Goodman elements [START_REF] Goodman | A model for the mechanics of jointed rock[END_REF].

Table 1 The nonlinear model parameters for different materials of Tankeng CFRD

Material

Dry density

(kN/m 3 ) m 0    0    b K 1 n G K 2 n f R

Numerical results

The main numerical results are presented in this section. The distributions of deformation and stress in the dam body are given for different construction and operation stages (Fig. 6), with the consideration of the effects of self-weight and post-impoundment water pressure. One observes that almost the same numerical results are obtained in two cases. For instance, the maximum settlements obtained with schemes 1 and 2 are 121.3cm and 121.4cm, respectively. The maximum settlement is about 0.75% of the dam height, which is of the same order of magnitude as the in-situ measurements observed in most of CFRDs [START_REF] Wen | A statistical review of the behaviour of concrete face rockfill dams based on case histories[END_REF]. Under the hydraulic pressure, the dam movement towards the upstream is partly restricted and the overall trend of the horizontal displacement varies towards the downstream. The maximum horizontal displacements towards upstream in schemes 1 and 2 are 9.02cm and 9.05cm, respectively, and the maximum horizontal displacements towards downstream in schemes 1 and 2 are 26.36cm and 26.34cm, respectively. The difference of the maximum displacements between two schemes is within 0.5%. In two schemes, the maximum principal stress is the same and located at the bottom of the dam. The maximum major principal stress is 2.49 MPa, and the minor principal stress is 1.17 MPa. Fig. 7 compares the dam deformation increments of two schemes, induced by the last step water load. It can be seen that the numerical results of scheme 2 agree well with those of scheme 1. The maximum settlement increments of schemes 1 and 2 are 17.51cm and 17.53cm, respectively. On the other hand, the maximum horizontal displacement increments of schemes 1 and 2 are 16.71cm and 16.85cm, respectively. The difference of the maximum displacement increments between two schemes is within 1%. For two schemes, the obtained deformation and stress fields of the dam body obtained are almost the same, indicating that the cross-scale finite element method can accurately transfer the upstream water load to the dam body. Based on the numerical analysis, one can conclude that the cross-scale finite element method is capable of providing good precision in the load transfer between discontinuous meshes. In the next section, the proposed method will be used to perform a 3D FEM simulation of the Tankeng CFRD.

CFRD 3D finite element analysis and results

FE analysis model

A 3D FEM analysis of the Tankeng CFRD is carried out by using the cross-scale finite element method in this section. The numerical model of the Tankeng CFRD considering the construction and water storage process is shown in Fig. 8. The same construction and water storage processes are used as those given in section 4. The dam crest is 505 m long and 12 m wide. Each concrete face slab is 12 m wide and waterproof vertical joints are used between adjacent slabs. At the bottom of the geometrical model of the Tankeng CFRD, the displacements are fixed to zero along the normal direction. In the numerical simulation, the eight-node hexahedron element is used for face slab and rockfill materials, while the arbitrary-node hexahedron element is used for elements in the transition zone. The studied slab and in the areas near the face slab such as cushion, waterproof joints, and so on (Fig. 9). The face slab elements size is less than 2 m and there are 82308 elements in the refined meshes, about 90 % of the total elements. The concrete slab and rockfill body are described by an elasto-plastic model and a nonlinear elastic model, respectively. Moreover, goodman elements [START_REF] Goodman | A model for the mechanics of jointed rock[END_REF]) are applied on the interface between the concrete face slab and the cushion layer. The vertical and peripheral joints are simulated by a pair of nodes, which could be combined into a single node under compressive stress and separated into two independent nodes under tensile stress.

Results and discussion

Fig. 10 shows the distribution of displacement components in the maximum cross-section after impounding. It can be seen the vertical settlement near the half-dam height is quite high. The maximum settlement is quale to 116 cm, which is about 0.72% of the dam height. This value is in the same order of magnitude as the previous in-situ measurements [START_REF] Hunter | Rockfill modulus and settlement of concrete face rockfill dams[END_REF][START_REF] Özkuzukiran | Settlement behaviour of a concrete faced rock-fill dam[END_REF][START_REF] Wen | A statistical review of the behaviour of concrete face rockfill dams based on case histories[END_REF]. The horizontal displacement of the dam body is affected by the upstream hydraulic pressure, resulting in the maximum displacement near the downstream in the dam body. Due to the gravity of the dam body and hydraulic pressure, the stress field in the dam body increases with the decrease of the dam height, and the maximum stress should occur at the bottom of the dam body (see Fig. 11). The maximum value of the major principal stress is 2.3 MPa, and the maximum value of the minor principal stress is 1.09 MPa. After impounding, the stress evolution in the face slab elements (where the compressive stress is defined as positive) is given in Fig. 13. Since the slabs move to the middle of the face slab along the dam's axial direction, the slabs squeeze each other and the stresses along the axial direction of the dam are compressive in most areas of the face slab. Meanwhile, as the deformation of the face slab is restricted by the valley, tensile stress appears on the face slab near the dam abutment and bank slope, and the maximum tensile stress attains 2.8MPa. The maximum concrete stress of 13.3MPa along the dam axial direction induced in the 0+252m~0+264m slab (elevation of 68 m) is marked by the contour distribution (Fig. 13a). Along the slope direction, compression is observed in most areas of the face slab, while the tensile stress only appears on the face slab near the bank slope with the maximum magnitude of order 2.4MPa (see Fig. 13b). Along the slope direction in the 0+324m slab located at the elevation of 95m, the maximum stress of 5.2MPa is observed in Fig. 13b. Obvious discontinuous distributions of stress along the slope direction can be observed in the face slab. At the bottom of a single slab near the bank slope, the tensile and compressive stresses appear on two sides of the slab respectively. This observation may be related to the development of the bending deformation in two directions of the slab.

Fig. 14 shows the deformation characteristics of the concrete face slab along the slope. Under hydraulic pressure, the bending deformation of the slab located at the dam sides results in large compressive stress in the middle of the slab, and the tensile stress zone occurs at the bottom of the slab near the bank slope, due to the strong restriction from bedrock on the bottom of the slab. At the same time, since the slabs tend to move towards the riverbed (Fig. 12a), the single slab near the bank slope has a bending deformation to the middle of the face slab (Fig. 14). Therefore, the bottom of a single slab near the bank slope is in tension on the side near the middle part of the face slab and in compression on the other side.

Meanwhile, the bending deformation also results in the compressive stress in the middle of a single slab on the side near the riverbed, which is less than that on the other side.

In order to study the local crack/failure in the concrete slab, the distribution of the plastic hardening parameter is illustrated in Fig. 15. The maximum compressive stresses in both the axial and slop directions appear near the middle part of the face slab. However, they are far less than its compressive strength of 39.7MPa. Therefore, the most important plastic hardening observed in concrete face slab does not occur in the middle of the face slab, but appears in the face slab near the bank slope. This observation is coherent with the distribution of tensile stress in the face slab. The maximum value of plastic hardening reaches up to 0.72 in the 0+444m slab at an elevation of 119m where the face slab is in a bi-directional tension state (Fig. 13). To reduce the plasticity of the concrete slab, the reinforcing bars can be added in this zone. As shown in Fig. 15, due to the use of refined meshes, the magnitude of plastic hardening of different regions in each slab is different. Consequently, it is necessary to prevent the face slab damage in the region, having important plastic hardening. Moreover, some discontinuous distribution zones are also observed in the numerical results (Fig. 13 and 15). Although coarse meshes could provide an average effect on the magnitude of stress and plastic hardening of concrete face slabs, some interesting phenomena may be ignored. To sum up, the obtained numerical results are in good agreement with the previous studies [START_REF] Sigtryggsdóttir | Statistical model for dam-settlement prediction and structural-health assessment[END_REF][START_REF] Wen | A statistical review of the behaviour of concrete face rockfill dams based on case histories[END_REF], which confirms the good performance of the proposed numerical method.

Conclusions

In the present work, a cross-scale finite element method, in which the arbitrary-node hexahedron element is used to link the differently scaled elements, is proposed. The performance of the proposed method is verified by comparing the numerical results of the CFRD section using uniform fine mesh and inhomogeneous mesh. The obtained results confirm that the proposed method is capable of getting a good prediction of load transfer between the fine elements of the face slab and the coarse elements of the rockfill dam body.

The proposed method is also applied in the 3D FE calculation of the Tankeng CFRD, in which the mechanical behaviors of rockfill and concrete are described by a nonlinear elastic model and an elastoplastic model, respectively. Since the concrete face slab in CFRDs plays the most important role as the main anti-seepage structure, the meshes comprising the face slab were refined adequately. The plastic hardening parameter  in the elastoplastic model of concrete is used to characterize the damage of the concrete face slab. The cross-scale finite element method proposed in this paper aims to deal with the mesh discontinuity between the concrete face slab and dam body. The obtained results exhibit that the proposed method is capable of providing a satisfactory prediction of the concrete face slab with fine mesh and dam body with coarse mesh. Although different sizes of mesh are applied on face slab and dam body in the proposed method, the stress gradient and the local damage in the concrete slab are the same as those of the conventional finite element analysis with fine mesh. Moreover, the cross-scale finite element method introduces a new element type into the conventional finite element code and possesses strong versatility, practicality and robustness. When the proposed method is used to perform the cross-scale simulation, the calculation efficiency mainly depends on the degree of freedom. As a result, the cross-scale finite element method has high computational efficiency. In the present work, the cross-scale simulation method is only used to analyze the working behavior of concrete face rockfill dams during construction and operation (i.e. water storage) phases. In the near future, the proposed method will be used to study various complex engineering problems, such as the dynamic analysis of earth-rock dams, the stress analysis of tunnel lining, and the structural analysis of underground powerhouses.

Fig. 1 .

 1 Fig. 1. Arbitrary-node hexahedron element

  have a similar impact on the second type of nodes. Because the differential of the special basis ( 1) ( 1

Fig. 2 .

 2 Fig. 2. Explanatory graph for the determination of local coordinates of hanging nodes

  and obtained by fitting the stress-strain curve of an isotropic compression test.

.,

  The stress-dilatancy equation D[START_REF] Liu | A microstructure-based elastoplastic constitutive model for coarse-grained materials[END_REF] can consider the typical stress-strain characteristic and the microstructure change of granular materials, expressed as and m is an experimentally fitting constant; M is the stress ratio   corresponding to the phase transformation point from contraction to dilatation. parameters and can be determined by using a series of triaxial tests; f M is the peak stress ratio under conventional triaxial condition;

  with those of classical FEM analysis. The Tankeng CFRD, 162m in height, is located on the middle reaches of the Oujiang River, Zhejiang Province, China. The typical cross-section of Tankeng CFRD is shown in Fig.4. The upstream slope of the dam body is 1:1.4, while the average slope of the downstream is 1:1.55. The thickness d of the concrete face slab is 0.3 m at the top elevation of 167 m and varies linearly with a function of 0 .3 0 .0 0 3 5 dH  downwards the slope, where H is the vertical

Fig. 4 .

 4 Fig. 4. Typical cross section of Tankeng CFRD

Fig. 5 .

 5 Fig. 5. Numerical schemes of Tankeng CFRD: (a) Scheme 1; (b) Scheme 2

Fig. 6 .Fig. 7 .

 67 Fig. 6. Calculated results of the dam body after water storage (displacement unit: cm; stress unit: MPa): (a) Vertical

Fig. 8 .

 8 Fig. 8. Construction stages and three-dimensional mesh of Tankeng CFRD

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Deformation distribution of dam body at the maximum section after water storage (unit: cm): (a) Vertical

Fig. 12 .

 12 Fig. 12. Deformation distribution of concrete face slab at the maximum section after water storage (unit: cm): (a) Axial

Fig. 13 .Fig. 14 .Fig. 15 .

 131415 Fig. 13. Stress distribution of concrete face slab at the maximum section after water storage (unit: MPa): (a) Stress

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant Nos. 52009036, 51979091).