esa-sceva/satcom-chunk-collection huggingface dataset:
## 2 Materials and Methods
### Preparation of Metal-oxide/Carbon Fibrous Mats  
The metal-oxide/carbon composite fibers were prepared by the centrifugal spinning of precursor solutions. The TiO\({}_{2}\)/C, Fe\({}_{2}\)O\({}_{3}\)/C and Fe\({}_{2}\)O\({}_{3}\)/TiO\({}_{2}\)/C composite-fiber anodes were prepared from precursor solutions as follows: For the preparation of TiO\({}_{2}\)/C composite fibers, a PVP solution (15 wt%) in ethanol and acetic acid (10:1 ratio) was prepared by stirring for 4 h. One g of titanium (IV) butoxide was added to the solution and stirred for another 12 h to obtain a homogeneous solution before centrifugal spinning. The Fe\({}_{2}\)O\({}_{3}\)/C composite fibers were prepared from a PVP solution (15 wt%) in a mixture of ethanol and acetic acid with a 10:1 ratio and then stirred for 4 h. Next, 1.5 g of iron (III) acetylacetonate was added to the solution and then stirred for an additional 12 h. Fe\({}_{2}\)O\({}_{3}\)/TiO\({}_{2}\)/C precursor fibers were prepared from a 15 wt% PVP solution in ethanol and acetic acid (10:1 ratio) and stirred for 4 h. 1 g of titanium (IV) butoxide and 1.5 g of iron (III) acetylacetonate were then added to the PVP solution and then stirred for an additional 12 h. The precursor solutions were prepared at room temperature. The TiO\({}_{2}\)/C, Fe\({}_{2}\)O\({}_{3}\)/C, and Fe\({}_{2}\)O\({}_{3}\)/TiO\({}_{2}\)/C composite-fiber anodes were fabricated by the centrifugal spinning of the corresponding precursor solutions followed by thermal treatment. Figure 1 shows a schematic of the centrifugal spinning process and the collected fibers and the as-prepared composite fibers. The CF method relies on applying centrifugal forces to a polymer solution contained in a spinneret rotating at a high rotational speed of up to 10,000 rpm. The polymer solution is forced to exit the nozzle orifice of the spinneret, and the resulting polymer fibers are deposited on the collector as shown in Figure 1A. A description of the CF spinning method was given in a previous work [22; 47; 48]. In brief, 2 mL of the prepared solution was injected into the needle-based spinneret equipped with 30-gauge half-inch regular bevel needles (Figure 1). The rotational speed of the spinneret was kept at 7000 rpm during the CF of all the solutions. The produced fibers were then collected using a 4 \(\times\) 4 in cardboard square covered with a spun-bond polypropylene substrate. Figure 1 shows the fibers formed before and after collection. Once the composite fibrous mat was collected, it was then removed from the substrate and dried at 120 \({}^{\circ}\)C in a vacuum oven for 12 h prior to carbonization. For the preparation of the composite-fiber electrode, the centrifugally spun fibers were dried and then stabilized in air at 200 \({}^{\circ}\)C for 2 h (with a heating rate of 2 \({}^{\circ}\)C/min), followed by carbonization at 550 \({}^{\circ}\)C for 5 h in an argon atmosphere to obtain the final metal-oxide/carbon composite fibrous mat.
## 3 Experiment Results and Discussion  
### Measuring Cases  
As described in the previous subsection, a pair of wing assemblies are manufactured in this paper. The twisting angle of the flapping-wing rotor is set as from \(10^{\circ}\) to \(50^{\circ}\). The average twisting angles are selected to be \(30^{\circ}\). This is because, when the twisting angle maintains at \(30^{\circ}\) during the flapping process, the FWR share the highest lift efficiency [17; 25]. Thus, in order to study the effect of the stroke angle on the lift generation for the flapping-wing rotor, the number of cases under the combination of different input voltage (\(V_{input}\)) and stroke angle (\(\phi\)) are listed in Table 2. The maximum voltage in these cases was 6.5 V since a higher voltage may lead to the stress concentration in the FWR's structure and wreck the structure. In order to ensure that the effect of the asymmetric stroke angle on the FWR's lift is easier to be observed, it is necessary to set the absolute value of the average stroke angle as large as possible in Cases 1 (\(\phi\in[-20,50]\)) and Cases 3 (\(\phi\in[-50,20]\)). However, the structural interference occurs at the U-shape mechanism of the FWR (as shown in Figure 1) when the stroke angle reaches about \(53^{\circ}\). Thus, the stroke angles are set as \(-20^{\circ}\)-\(50^{\circ}\) and \(-50^{\circ}\)-\(20^{\circ}\) in Case 1 and Case 3.
## 2 Algorithm Principle
### INS+LEO2-Satellite Alternate Switching Ranging Integrated Navigation Algorithm
#### 2.2.2 Algorithm Flow  
We assume that the real-time carrier coordinates measured by the INS in the ECEF coordinate system are \((\hat{x}_{u},\hat{y}_{u},\hat{z}_{u})\), which is a known quantity. \(\hat{\rho}_{1}\) and \(\hat{\rho}_{2}\) represent the virtual ranging values of satellites #1 and #2 measured by INS, respectively, the alternate interval is \(\Delta t\), and other parameters are the same as the standard INS+LEO2-satellite integrated navigation algorithm. The corresponding working principle diagram is shown in Figure 2.  
At time \(t_{1}\), using the real ranging value \(\rho_{1}\) of satellite #1, and the virtual ranging value \(\hat{\rho}_{2}\) of satellite #2, the following expressions are obtained:  
\[\left\{\begin{array}{l}\rho_{1}=\sqrt{\left(x_{1}-x_{u}\right)^{2}+\left(y_{ 1}-y_{u}\right)^{2}+\left(z_{1}-z_{u}\right)^{2}}\\ \hat{\rho}_{2}=\sqrt{\left(x_{2}-\hat{x}_{u}\right)^{2}+\left(y_{2}-\hat{y}_{ u}\right)^{2}+\left(z_{2}-\hat{z}_{u}\right)^{2}}\end{array}\right. \tag{3}\]  
At this time, according to Equations (2) and (3), we use the difference between \(\rho_{1}\) and \(\rho_{11}\) and the difference between \(\hat{\rho}_{2}\) and \(\rho_{12}\) as the filter's observations.  
At time \(t_{2}=t_{1}+\Delta t\), using the virtual ranging value \(\hat{\rho}_{1}\) of satellite #1 and the real ranging value \(\rho_{2}\) of satellite #2, the following expressions are obtained:  
\[\left\{\begin{array}{l}\hat{\rho}_{1}=\sqrt{\left(x_{1}-\hat{x}_{u}\right)^ {2}+\left(y_{1}-\hat{y}_{u}\right)^{2}+\left(z_{1}-\hat{z}_{u}\right)^{2}}\\ \rho_{2}=\sqrt{\left(x_{2}-x_{u}\right)^{2}+\left(y_{2}-y_{u}\right)^{2}+ \left(z_{2}-z_{u}\right)^{2}}\end{array}\right. \tag{4}\]  
At this time, according to Equations (2) and (4), we use the difference between \(\beta_{1}\) and \(\rho_{11}\) and the difference between \(\rho_{2}\) and \(\rho_{12}\) as the filter's observations.  
We then loop through this process, at every \(\Delta t\) time alternating switching the real ranging value and virtual ranging value of the satellite, straight to the end of flight time.  
Analysis of Dual-Satellite Alternate Switching Ranging/INS Integrated Navigation Algorithm for Broadband LEO Constellations
## 3 Modeling of City Environment
### Simulation Method  
The simulation is performed using MATLAB. It is divided into two parts: The first is that LoS probability in the street is studied as a function of elevation angle for four types of build areas. Second, the additional path loss resulting from the shadowing effects of buildings is investigated using the uniform theory of diffraction (UTD) method [35; 36]. Our analysis considers the building's walls constructed of concrete and brick with relative permittivity \(\varepsilon\) of 15 and 4, respectively and conductivity \(\sigma\) (0.015 S/m). Diffraction loss for vertical and horizontal polarization is computed separately. Figure 2 shows the geometry that encompasses both LoS and NLoS communication.  
The simulation involves creating four virtual environments within a 1 km \(\times\) 1 km area to mimic different scenarios in MATLAB. The UAV's position is determined along the street. The azimuthal angle ranges from 0 to 360 degrees in increments of 9 degrees. Calculations are carried out for specific elevation angles using the median value obtained from the azimuthal angle. This approach ensures that the results are not dependent on the azimuthal angle, accounting for the uneven spacing of buildings in the real world. The simulation covers elevation angles from 1 to 89 degrees, assuming a maximum UAV altitude of 500 m. This setup reflects the use of the UAV as a base station or relay to meet the demand of 5G and beyond.  
Figure 1: Selected layout for city areas.
# Traffic Signal Time Optimization Based on Deep Q-Network
## 3 Problem Formulation
### Deep Q-Network  
To solve the problem of TSC, this study uses a DQN, which is an RL method. DQN defines the environment, recognizes the current state (\(s_{t}\)) at time \(t\), selects a possible action (\(a_{t}\)), and receives a reward (\(r_{t}\)) for the result. The value for the selected action is defined as the value (\(y_{t}\)) at time \(t\), and the value predicted by the DQN model is defined as the predicted value (\(Q_{s,a}\)):  
\[L_{t}(\theta_{t})=\mathbb{E}_{s,a}\big{[}(y_{t}-Q(s,a;\theta_{t}))^{2}\big{]}. \tag{1}\]  
The Q-network updates the network weights (\(\theta\)) through neural network learning by setting the minimum value of the loss function (\(L\)) to (1). In addition, a DQN consists of neural networks that are Q-networks because DQN combines Q-learning and deep learning. Therefore, it has the advantage of proceeding with learning using high-dimensional input. However, a DQN has a disadvantage in that it collects data sequentially over time in the environment. These data are then used as the input data. Because sequential data have high correlations, experience replay is applied to solve the data correlation problems. Experience replay is a method of storing the experience obtained at each time in the environment and randomly selecting the experience to learn by organizing a mini-batch. This can reduce the high correlation between the data and increase the data reliability.
## 7 Results and Discussion
CNN models can predict correct results for harbour and coastal waters using Gaussian and Rectangular optical pulse shapes, with different pulse width ranges from 0.1 to 0.95, using DGF, WDGF, CEAPF, and BP channel models. The trend in the curves is similar to the identity function, represented by the red line (\(y\)=\(x\)), which means the actual values are close to the predicted ones in both SNR on the left and BER on the right. This shows that the CNN model can decide correctly in various situations involving various types of water, ISI noise, and water environment variations. The relation between SNR and BER, which represents the performance of the optical receiver, is drawn, as represented in Figure 21, using a 0.45 dropout rate and a \(10^{-5}\) learning rate with 28 filters. From the graphs, we can conclude that the suggested CNN models  
Figure 19: Number of filters vs. number of trainable parameters for all CNN models. Figure 20: Performance of the CNN models, the true versus the predicted SNR (**left**) and BER (**right**) values. Various channel models are employed to simulate the behaviour of water in harbours (represented by the colour blue) and coastal areas (represented by the colour green) for different pulse widths. Figure 18: Number of filters vs. loss and RMSE ratios for all CNN models. perform well in making accurate decisions for various instances involving various types of waters, ISI noise, and underwater environment variations. Regarding the high BER values, it means small SNRs are included in our models. This is due to this study's significant noise and channel fluctuations. Although the SNRs are caused by received pulses that pass through noisy channels, the models can predict SNRs accurately.
# A Communication Framework for Image Transmission through LPWAN Technology
## 3 Proposed Framework
### Stage III: Transmission Bandwidth Increase  
Initially, our research looked into the possibility of discovering a method to increase bandwidth. Techniques such as diversity, Multiple In Multiple Out (MIMO), Single In Multiple Out (SIMO), and others were investigated; however, the results indicate that for our framework they are not required. As a result, we focus on implementing LoRa modulation in Software-Defined Radios (SDR) in the transmitter and receiver, allowing us to have more degrees of freedom than LoRa commercial radios in terms of bandwidth, transmission power, spreading factor, preamble sequence definition, and sync sequence chirps definition. LoRa features can be optimized to achieve higher data rates by taking advantage of the flexibility provided by SDR. Test results show that high bandwidth (>500 kHz) results in higher capacity and lower SF values (SF7), resulting in data rates around 50 kbps with good sensitivity at reception, even near noise floor values. This result is adjusted by having a large payload, a short overhead, and using low coding rates [48; 99; 100; 101; 102; 103; 104; 105; 106].
# Enhanced Intrusion Detection Systems Performance with UNSW-NB15 Data Analysis
## 5 Results and Discussion
### Hyper-Parameter Tuning  
After defining the performance matrix, the initial experiment focuses on Logistic Regression. Table 4 details the step-by-step implementation of Logistic regression, considering hyperparameters such as penalty and alpha. Hyperparameter tuning is conducted using the grid search method. Subsequently, the linear support vector classifier is explored, utilizing the same hyperparameters (penalty and alpha) for tuning. Following the assessment of logistic regression and linear support vector machine models, attention shifts to the decision tree model. Despite obtaining a lower AUC compared to the linear support vector model, efforts are made to enhance the decision tree model's accuracy. Three hyperparameters, namely \"min_sample_leaf\", \"min_sample_split\", and \"max_depth\", are considered for tuning. Optimal values for \"min_sample_split\" and \"max_depth\" are determined to be 6 and 10, respectively. Hyperparameter tuning for \"min_sample_leaf\" reveals its optimal value as 9.  
After evaluating the Decision Tree model, the subsequent experiment involves implementing the Random Forest model. Through a train and cross-validation (CV) approach, optimal hyperparameters are identified as criterion = \"gini\", Max_depth = 22, min_sample_split = 6, n_estimators = 300, and n_jobs = \(-\)1. The Gradient Boosted Decision Tree model requires adjustment of several hyperparameters, including \"learning_rate\",  
\begin{table}
\begin{tabular}{l} \hline \hline
**Logistic Regression** \\ Step 1: Import the sklearn libraries \\ Step 2: Define the hyperparameter tuning function to select the best hyperparameters by \\ implementing the grid search method \\ Step 3: Define a function to evaluate the performance of the model in terms of Accuracy, F1-Score, \\ and False Alarm Rate \\ Step 4: Create an instance of SGDC classifier with logistic regression loss function and L2 \\ regularization \\ Step 5: Define a dictionary of hyperparameters to be tuned during hyperparameter tuning \\ Step 6: Create a final LR model with the best hyperparameters. Step 7: Evaluate the performance of the best logistic model on training and testing data \\ \hline \hline \end{tabular}
\end{table}
Table 4: Logistic regression—steps-by-step implementation.  
\"max_depth\", \"colsample_bylevel\", \"subsample\", and \"n_estimators\". Experimentation with various parameter combinations reveals the significant influence of 'learning_rate', and to a lesser extent,'max_depth' and 'n_estimators', while the impact of the other two parameters is minimal. The best hyperparameters for the model are determined as 'n_estimators': 400,'max_depth': 12, 'learning_rate': 0.1, 'colsample_bylevel': 0.5, and'subsample': 0.1.
# Asymmetric Propagation Delay-Aware TDMA MAC Protocol for Mobile Underwater Acoustic Sensor Networks
## 3 Simulation Results and Analysis
### Performance Versus the Network Range  
The \(CD\) and the \(CU\) of the APD-TDMA and the BTB-TDMA are also evaluated under various \(R\) from 500 m to 3000 m when \(N=25\) and \(V\)max = 5 m/s, as illustrated in Figure 9. It can be seen from Figure 9 that the APD-TDMA outperforms the BTB-TDMA with respect to the \(CD\) and the \(CU\) for all network ranges. The average propagation delay from the SN to each UN increases as R increases, and this increased propagation delay degrades the \(CD\) and the \(CU\) of the APD-TDMA and the BTB-TDMA. However, the slope of the APD-TDMA in Figure 9 is significantly smaller than that of the BTB-TDMA, because data packets of the APD-TDMA arrive at the SN in a packet-train manner without collision, and because the idle time at the SN of the APD-TDMA is minimized, which shows the robustness of the APD-TDMA with respect to the network range. On average, the \(CD_{APD}\) is smaller than the \(CD_{BTB}\) by 26%, and the \(CU_{APD}\) is greater than the \(CU_{BTB}\) by 38%.
## 2 Theoretical Derivation and Energy Band Analysis of Mechanical Properties of the Cell
### Theoretical Derivation of Re-Entrant Hexagonal Cell  
A simple schematic diagram of the honeycomb base is shown in Figure 2a, and its standard re-entrant hexagonal honeycomb cell is shown in Figure 2b, where \(h\) is the vertical side length of cell, \(l\) is the hypotenuse length of cell, and \(\theta\) is the angle of cell (\(\theta<0\)).  
Figure 1: Deformation diagram of positive and negative Poisson’s ratio structure under lateral hand tension is given, (**a**) positive Poisson’s ratio structure, (**b**) negative Poisson’s ratio structure.  
According to reference [38], Poisson's ratio and relative density of the re-entrant element with negative Poisson's ratio can be derived as follows. When the re-entrant hexagonal honeycomb configuration is loaded in \(X_{1}\) or \(X_{2}\) direction and deforms in a linear elastic manner, the dimension relation in Figure 2b can be obtained from the geometric relation:  
\[X_{1}=2l\cos\theta \tag{1}\]  
\[X_{2}=2h+2l\sin\theta \tag{2}\]  
It is assumed that the configuration of re-entrant hexagonal honeycomb cell is rigid structure during the micro-deformation process, it means that, \(h\) and \(l\) are constant in the process of deformation, and its deformation only depends on the change of the cell angle \(\theta\). For this structure, [PERSON]'s ratio can be expressed as:  
\[V_{\rm ij}=\left(V_{\rm ij}\right)^{-1}=-\frac{d\varepsilon_{j}}{d\varepsilon _{i}}\therefore i,j=1,\;2 \tag{3}\]  
Figure 2: The re-entrant hexagonal structure with negative Poisson’s ratio is given, (**a**) a simple schematic diagram of the honeycomb base is shown, (**b**) re-entrant hexagonal honeycomb cell.  
where, the \(d\epsilon\) included angle \(\theta\) can be expressed as:  
\[d\epsilon_{i}=\frac{1}{X_{i}}\frac{dX_{i}}{d\theta}d\theta \tag{4}\]  
Substituting Poisson's Ratio into the calculation formula, we can get:  
\[v_{21}=(v_{12})^{-1}=-\frac{d\epsilon_{1}}{d\epsilon_{2}}=-\frac{dX_{1}/d \theta}{dX_{2}/d\theta}\cdot\frac{X_{2}}{X_{1}} \tag{5}\]  
From Equations (1) and (2), it can be known that:  
\[\frac{dX_{1}}{d\theta}=-2l\sin\theta \tag{6}\]  
\[\frac{dX_{2}}{d\theta}=-2l\cos\theta \tag{7}\]  
Substituting it into Equation (5), we can get:  
\[v_{21}=(v_{12})^{-1}=\frac{(h\,/l+\sin\theta)\sin\theta}{\cos^{2}\theta} \tag{8}\]  
According to the geometric relationship, the relative density of cell can be expressed as:  
\[\frac{\rho^{*}}{\rho_{\rm s}}=\frac{(4l+2h)\tau}{4l\cos\theta(h+l\sin\theta)} =\frac{\tau/l(h/l+\sin\theta)}{\cos^{2}\theta} \tag{9}\]
## 2 Models
### Operating Environment Modeling  
Based on the model assumptions, the UAV only conducts climb and descent operations in its takeoff and landing area, and the remaining flight phases are conducted in a two-dimensional plane. Assume that the two-dimensional plane is a rectangle with length \(L\) and width \(W\). The rectangle is divided into \(u\times v\) grids with a square grid with grain size \(l_{g}\), where \(u=\text{int}(Length/l_{g})\), \(v=\text{int}(Width/l_{g})\), and \(\text{int}()\) is the rounding function. The grid center point is the optional air-route point. When planning the air-route of logistics UAVs, the operational risks of UAVs should be considered comprehensively. UAV operational risks can be expressed as:  
1. Collision risk During actual operation, there is no-fly airspace for terrain obstacles such as tall buildings, which should be strictly avoided when planning logistics UAV air routes. To quantitatively describe airspace characteristics, grids in different special areas can be assigned different grid collision-risk values, as expressed in Table 1.  
In Table 1, \(D_{g}\) is the collision-risk value of the dangerous airspace grid, which can be expressed as:  
\[D_{g}=\frac{N_{obstacle(g)}}{N_{surround(g)}} \tag{1}\]  
where \(N_{obstacle(g)}\) refers to the number of no-fly airspace grids around grid \(g\) and \(N_{surround(g)}\) refers to the total number of grids around grid \(g\).  
2. Crash risk The crash risk can be defined by the \"equivalent safety level\" of the UAV, which is the number of casualties per unit of time:  
\[S_{g}=P_{\text{UAV}}\times E_{g}\times\left(1-G_{g}\right)\times N_{g} \tag{2}\]  
where \(S_{g}\) refers to the crash risk of the grid; \(P_{\text{UAV}}\) refers to the probability of UAV crash accident per hour--the current 5% NMAC standard is adopted in this paper, that is, the  
\begin{table}
\begin{tabular}{c c c} \hline \hline
**Grid Category** & **Flight Criteria** & **Collision Risk Value** \\ \hline No-fly airspace & no entry & \(1\) \\ dangerous airspace & flight allowed & \(D_{g}\) \\ free airspace & flight allowed & \(0\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Grid collision risk.  
maximum collision probability between UAVs is 0.05; \(E_{g}\) refers to the ability of a UAV crash to cause damage, which is the sum of gravitational potential energy and dynamic potential energy of the UAV; \(N_{g}\) refers to the number of people affected in the grid, which is assumed to be one over land and zero over water for simplifying the model; and \(G_{g}\) refers to a safety barrier factor related to the geographical characteristics of the grid [23], with the specific settings shown in Table 2.
## 2 Spectral Wave Equation for Uniform Cross-section Waveguides  
### Mechanical Variables  
The mathematical model is presented here for the case of guided waves propagating in a uniform cross-section waveguide made of linear viscoelastic materials. The waveguide is considered in vacuum by assuming a stress-free condition on the waveguide outer surface. A reference Cartesian system is assumed as shown in Figure 1. The time dependent mechanical variables of the problem at a point \(\mathbf{x}\equiv\left(x,y,z\right)\) of the waveguide are the displacement \(\mathbf{u}\), the stress \(\mathbf{\sigma}\), the strain \(\mathbf{\varepsilon}\)  and the load vector \(\mathbf{t}\) with components:  
\[\mathbf{u}(\mathbf{x},t)=\left[u\ v\ w\right]^{\mathrm{T}} \tag{1}\]  
\[\mathbf{\sigma}(\mathbf{x},t)=\left[\sigma_{xx}\ \ \sigma_{yy}\ \ \sigma_{zz}\ \ \sigma_{yz}\ \ \sigma_{xz}\ \ \sigma_{xy}\right]^{\mathrm{T}} \tag{2}\]  
\[\mathbf{\varepsilon}(\mathbf{x},t)=\left[\varepsilon_{xx}\ \ \varepsilon_{yy}\ \ \varepsilon_{zz}\ \ \gamma_{yz}\ \ \gamma_{xz}\ \ \gamma_{xy}\right]^{\mathrm{T}} \tag{3}\]  
\[\mathbf{t}(\mathbf{x},t)=\left[t_{x}\ \ t_{y}\ \ t_{z}\right]^{\mathrm{T}} \tag{4}\]  
where the superscript \(\mathrm{T}\) means a transpose vector. The field variables are chained by the linear viscoelastic constitutive and linear strain compatibility equations:  
\[\mathbf{\sigma}=\mathbf{\tilde{C}}(t)\mathbf{\varepsilon} \tag{5}\]  
\[\mathbf{\varepsilon}=\left[\mathbf{L}_{x}\ \frac{\partial}{\partial x}+\mathbf{L}_{y} \ \frac{\partial}{\partial y}+\mathbf{L}_{z}\ \frac{\partial}{\partial z}\right]\mathbf{u} \tag{6}\]  
where the \(6\times 3\ \ \mathbf{L}_{i}\) operators are defined as:  
\[\mathbf{L}_{x}=\begin{bmatrix}1&0&0\\ 0&0&0\\ 0&0&0\\ 0&0&0\\ 0&0&1\\ 0&1&0\end{bmatrix},\qquad\mathbf{L}_{y}=\begin{bmatrix}0&0&0\\ 0&1&0\\ 0&0&0\\ 0&0&1\\ 0&0&0\\ 1&0&0\end{bmatrix},\qquad\mathbf{L}_{z}=\begin{bmatrix}0&0&0\\ 0&0&0\\ 0&0&1\\ 0&1&0\\ 1&0&0\\ 0&0&0\end{bmatrix} \tag{7}\]  
In Equation (5) \(\mathbf{\tilde{C}}(t)\) is the time-dependent material linear viscoelastic operator with at most 21 independent coefficients.