\documentclass[a4paper,11pt]{article}

\topmargin-1cm
\textheight=23,5cm
\textwidth=16cm
\oddsidemargin=0.5cm
\usepackage[normalem]{ulem}
\usepackage[utf8]{inputenc}
\usepackage[dvips]{graphicx}
%\usepackage[pdflatex]{graphicx}
\usepackage{subfig}   %to label and display sub-figures
\usepackage{subfloat}
\usepackage{bm}
\usepackage{appendix}
%\usepackage[font=small,labelfont=bf]{caption}
%\usepackage{subcaption}
\usepackage{multirow}
\usepackage{morefloats}
\usepackage[labelfont=bf]{caption}
\usepackage{algorithm2e}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{cases}
\usepackage{moreverb}
\usepackage[colorlinks,bookmarksopen,bookmarksnumbered,citecolor=red,urlcolor=red]{hyperref}
\usepackage[stable]{footmisc}
\usepackage{cleveref}
\usepackage{cancel}
\usepackage{bm}
\usepackage{cite}
\usepackage{times}
%\usepackage[authoryear]{natbib}

\ifpdf
  \DeclareGraphicsExtensions{.pdf,.jpg,.png}
  \graphicspath{{./figuras/png/}}
\else
  \DeclareGraphicsExtensions{.eps,.ps,.jpg,.png}
  \graphicspath{{./figuras/png/}}
\fi

\providecommand{\newoperator}[3]{\newcommand*{#1}{\mathop{#2}#3}} % safely create new math operators
\newoperator{\std}{\mathrm{std}}{\nolimits} % standard deviation operator

\newcommand\BibTeX{{\rmfamily B\kern-.05em \textsc{i\kern-.025em b}\kern-.08em
T\kern-.1667em\lower.7ex\hbox{E}\kern-.125emX}}

\usepackage{floatrow}
\def\dobi#1{#1\kern0.5pt\hbox{\llap{\it#1}}}

\def\itneg{\fontsize{10}{10}\usefont{OT1}{phv}{bc}{it}\selectfont}

\def\pmb#1{\setbox0=\hbox{$#1$}%
             \kern-.027em\copy0\kern-\wd0
             \kern+.009em\copy0\kern-\wd0
             \kern+.009em\copy0\kern-\wd0
             \kern+.009em\copy0\kern-\wd0
             \kern+.009em\copy0\kern-\wd0
             \kern+.009em\copy0\kern-\wd0
             \kern+.009em\copy0\kern-\wd0
             \kern-.045em\raise+.012em\copy0\kern-\wd0
             \kern+.009em\raise+.012em\copy0\kern-\wd0
             \kern+.009em\raise+.012em\copy0\kern-\wd0
             \kern+.009em\raise-.012em\copy0\kern-\wd0
             \kern+.009em\raise-.012em\copy0\kern-\wd0
             \kern-.018em\copy0\kern-\wd0\raise-.012em\box0}

\def\pder#1#2{{\partial #1 \over \partial #2}}
\def\2pder#1#2{{\partial^2 #1 \over \partial #2}}

\makeatletter
\def\bs{\expandafter\@gobble\string\\}
\def\lb{\expandafter\@gobble\string\{}
\def\rb{\expandafter\@gobble\string\}}
\def\@pdfauthor{C.V.Radhakrishnan}
\def\@pdftitle{elsarticle.cls -- A documentation}
\def\@pdfsubject{Document formatting with elsarticle.cls}
\def\@pdfkeywords{LaTeX, Elsevier Ltd, document class}
\def\file#1{\textsf{#1}\xspace}

%\newfloatcommand{capbtabbox}{table}[][\FBwidth]
\makeindex

%\graphicspath{{latex/gfx/PNG/}{latex/gfx/PDF/}{latex/gfx/}}

% Demo replaces images with black boxes
%\usepackage[demo]{graphicx}
%\usepackage[TABTOPCAP,FIGBOTCAP]{subfigure}

\bibliographystyle{elsarticle-num}
% math shortcuts


\let\ACMmaketitle=\maketitle
\renewcommand{\maketitle}{\begingroup\let\footnote=\thanks \ACMmaketitle\endgroup}

\setcounter{page}{1}

\date{}

\title{\textbf{FIRST ORDER EXPLICIT TIME INTEGRATION SCHEME WITH LARGE TIME STEPS FOR PARABOLIC PROBLEMS USING IRREGULAR GRIDS}}

\author{\textbf{Eugenio O\~{n}ate}$^{1,2}$, \textbf{Francisco Z\'{a}rate}$^{2}$,  \textbf{Juan M. Gimenez}$^{2}$ and \textbf{Sergio R. Idelsohn}$^{2}$\\[.5cm]
$^{1}$ Emeritus Professor at Universitat Polit\`{e}cnica de Catalunya (UPC),\\ Barcelona, Spain, Email: onate@cimne.upc.edu (Corresponding author)\\
$^{2}$ Centre Internacional de M\`{e}todes Num\`{e}rics en Enginyeria (CIMNE),\\ Barcelona, Spain}

\begin{document}
\maketitle


\begin{abstract}

We present a fast first order explicit time integration scheme for solving parabolic problems in mechanics  via standard numerical methods in space using irregular grids, such as unstructured finite element meshes, or grids containing elements or cells of very different sizes. The new scheme extends one of the explicit FIC-Time (EFT) integration methods derived by the authors in \cite{On24} that allow considerable larger time steps than the forward-Euler (FE) scheme. The new EFT scheme overcomes the limitations in the time step size of explicit time integration schemes for irregular grids containing large and small elements. A variable time step is used for eliminating the oscillations near Dirichlet boundaries when large time steps are used. The advantages of the new EFT scheme versus the FE scheme are shown in one-, two- and three-dimensional transient heat conduction problems using irregular finite element grids.

\end{abstract}

\noindent \textbf{Keywords}: Explicit time integration schemes, parabolic equation, large time steps, irregular grid,\\ finite increment calculus, FIC-Time.

\section{INTRODUCTION}

Parabolic problems are common in mechanics. Examples are heat conduction, seepage in porous media, gas diffusion, magnetostatics, Reynolds film lubrication,  creep and viscoplasticity of solids, Stokes flows and compressible gas flows, just to cite a few cases \cite{Belytschko14,Hughes87,ZiTa05,ZiTaNi05}.

The explicit time integration of parabolic problems is restricted by the critical time step which governs the stability of the numerical solution  \cite{Butcher03,Hughes87}. The value of the critical time step is proportional to the inverse of the maximum eigenvalue of the discretized system. This leads to very small time steps when fine or irregular grids are used in the context of the finite element method (FEM), the finite volume method (FVM) and the finite difference method  (FDM)\cite{Hughes87,ZiTa05,ZiTaNi05}.

In O\~{n}ate et al. \cite{On24} the authors presented a family of explicit time integration schemes of first, second and third order accuracy for  parabolic problems that have a higher computational efficiency versus explicit methods of the same order.
The derivation of the new explicit schemes is based on the Finite Element Calculus (FIC) procedure for modelling the transient heat conduction problem \cite{On24}. The stability and time step size of the new \emph{explicit FIC-Time} schemes (or \emph{EFT} schemes in short)  is governed by time stabilization parameters. The examples of application in one (1D) and two (2D) dimensions presented in O\~{n}ate et al. \cite{On24} showed that the new EFT schemes  allow time steps considerably larger than those for the standard first order forward-Euler (FE) scheme and the second and third order  Adams-Bashforth (AB) schemes \cite{Butcher03} while retaining the accuracy of the numerical solution. They also have advantages versus Runge-Kutta (RK) schemes \cite{Butcher03} in terms of the limit time step size (for second order schemes) and the number of matrix-vector multiplications per time step (for second and third order schemes).

In this work the first order EFT12 scheme presented in \cite{On24} is extended for dealing with irregular grids, such as unstructured finite element meshes or grids containing small and large elements (in the FEM) or cells (in the FVM and the FDM). The critical time step for explicit integration schemes drastically reduces in theses cases, as the difference in size between the smallest and largest element/cell in the grid increases. In the paper we show that the EFT12 scheme with a variable time step size overcomes this limitation and, therefore, it allows considerable large time steps than for the FE scheme, while still providing an accurate transient numerical solutions. The variable time step is also essential for obtaining accurate numerical solutions in the vicinity of Dirichlet boundaries using the EFT12 scheme with large time steps.

The lay-out of the paper is the following. In the next section the solution of the semidiscrete form of the parabolic equation with the FE scheme is presented. Next the FIC-time approach used for deriving  the EFT12 scheme is presented. The extension of this scheme for variable time steps is detailed. The good behavior of the EFT12 scheme with large variable time steps is demonstrated in the solution of transient heat conduction problems in one, two and three dimensions using irregular finite element grids.


\section{NUMERICAL SOLUTION OF THE SEMIDISCRETE\\  PARABOLIC EQUATION IN MECHANICS. EXPLICIT FORWARD-EULER SCHEME}\label{first-order-transient}

The semidiscrete form  of the parabolic  equation in mechanics using either FEM, FVM or FDM has the following standard form \cite{Hughes87,ZiTaZhu05}
\begin{equation}\label{first-order-transient-1}
{\bf C}\dot {\bf a} + {\bf K}{\bf a} - {\bf f} = {\bf 0}\qquad \mbox{in } [\Omega , t].
\end{equation}
where ${\bf a}(t)$ is the vector of unknowns of the discretized problem in space, $\dot {\bf a}=\frac{d}{dt}{\bf a}$, ${\bf f}$ is the vector emanating from external actions to the system, and ${\bf C}$ and ${\bf K}$ are square matrices involving physical and geometrical parameters of the problem, $\Omega$ is the space domain where the parabolic equation is solved and $t$ is the time.

Without loss of generality, we will assume here that Eq.\eqref{first-order-transient-1} refers to transient heat conduction problems. In this context, ${\bf a}$ is the vector of the nodal temperatures,  ${\bf C}$ is the heat capacity matrix which terms depend of $\rho c$ where $\rho$ is the density and $c$ is the specific heat, ${\bf K}$ is the conductivity matrix and ${\bf f}$ is the flux vector contributed by internal and external heat sources. Eq.\eqref{first-order-transient-1} is completed with the Dirichlet and Neuman boundary conditions specifying the  values of the temperature and the normal heat flux at the corresponding boundary, respectively, and with the initial temperature field on $\Omega$.
The  matrices and vectors of Eq.\eqref{first-order-transient-1} and the boundary conditions can be found in many text books \cite{Hughes87,ZiTa05}.

In this work we will assume a diagonal form of matrix ${\bf C}$ and  isotropic conductivity depending on a single conductivity parameter k. Also we will assume that all the material properties are independent of time.

The explicit solution of Eq.\eqref{first-order-transient-1} can be written, using the FE scheme \cite{Belytschko83,Hughes87,ZiTa05}, as
\begin{equation}
 \frac{1}{\Delta t} {\bf C}  ({\bf a}^{n+1} -{\bf a}^n )+  {\bf K} {\bf a}^n  - {\bf f}^n =0
\label{first-order-transient-2}
\end{equation}
where $(\cdot )^n$ denotes values at time $t=t_n$.

From Eq.\eqref{first-order-transient-2} we can find
\begin{equation}
{\bf a}^{n+1} = {\bf a}^n -\Delta t {\bf C}^{-1} \left[ {\bf K} {\bf a}^n  - {\bf f}^n\right]
\label{first-order-transient-3}
\end{equation}

Eq.\eqref{first-order-transient-3} shows that for a diagonal form of ${\bf C}$ the time evolution of the nodal temperatures can be obtained without solving a system of equations.

The stability of the FE scheme is found as follows.
The solution of the discretized transient problem expressed by Eq.\eqref{first-order-transient-1}  is assumed of the form
\begin{equation}
{\bf a}({\bf x},t) =\sum\limits_{i=1}^N {\pmb \phi}_i ({\bf x}) y_i (t)
\label{first-order-transient-4}
\end{equation}
where $N$ is the number of nodal temperatures in the discretized system and $y_i (t) $ and ${\pmb \phi}_i$ are, respectively the $i$-th modal amplitude and the $i$-th modal vector of the eigenvalue system
\begin{equation}
\vert {\bf K}-\lambda_i {\bf C}\vert {\pmb \phi}_i=0 \quad , \quad i=1, N
\label{first-order-transient-5}
\end{equation}

In Eq.\eqref{first-order-transient-5} $\lambda_1,\cdots , \lambda_N$ are the eigenvalues of the semidiscrete system ordered from the smallest one, $\lambda_1$, to the largest one, $\lambda_N$.

Substituting Eq.\eqref{first-order-transient-4} into \eqref{first-order-transient-3} for  ${\bf f}^n =0$, and using the orthogonality properties of the modal vectors ${\pmb \phi}_i$ versus the matrices ${\bf C}$ and ${\bf K}$ \cite{Belytschko14,Hughes87}, we can obtain, after little algebra, the expression for the time evolution of the $i$-th modal amplitude as
\begin{equation}
{y}_i^{n+1} = (1 - p_i) y_i^n
\label{first-order-transient-6}
\end{equation}
with
\begin{equation}
p_i = \frac{k_i}{c_i}\Delta t =\lambda_i \Delta t \quad \hbox{with } \lambda_i = \frac{k_i}{c_i}
\label{first-order-transient-7}
\end{equation}
and
\begin{equation}
k_i = {\pmb\phi}_i^T {\bf K} {\pmb\phi}_i\quad , \quad {c}_i = {\pmb\phi}_i^T {\bf C} {\pmb\phi}_i
\label{first-order-transient-8}
\end{equation}


Assuming a stable solution of the type $y_i^{n+1} =\gamma_i y_i^n$, with $\gamma_i$ being the amplification factor for the $i$-th mode, we can obtain from Eq.\eqref{first-order-transient-6},
\begin{equation}
\gamma_i =(1-p_i)
\label{first-order-transient-9}
\end{equation}

The stability condition, ensuring  $\vert \gamma_i\vert <1$ is
\begin{equation}
-1 < (1-p_i)< 1 \qquad \mbox{for } i=1,N
\label{first-order-transient-10}
\end{equation}

From Eqs.\eqref{first-order-transient-7} and \eqref{first-order-transient-10} we deduce the stability and time step conditions for $i=N$, as
\begin{equation}
p_N < 2 \qquad \mbox{and}\qquad \Delta t < \Delta t_l^{\rm FE} =  \frac{2}{\lambda_N}
\label{first-order-transient-11}
\end{equation}
where $\Delta t_l^{\rm FE}$ is the limit time step that guarantees the stability of the FE scheme.


In order to avoid non-physical oscillatory solutions, the following additional condition must be satisfied by the $i$-th amplification factor
\begin{equation}
0< \gamma_i < 1 \qquad , \quad \mbox{for } i =1,m
\label{first-order-transient-12}
\end{equation}
where $m$ is the number of low order modes that govern the solution of the transient problem. We have found that, in practice, it suffices that condition \eqref{first-order-transient-12} is satisfied for the first mode  only, i.e. $0< \gamma_1 <1$. An exception to this rule are the possible oscillations in the space distribution of the unknown near the Dirichlet boundaries during the initial steps of the transient solution. These oscillations can be eliminated by using  time steps smaller than the limit value during the first instants of the transient solution. The alternative procedure, developed in this paper in the context of the new EFT12 scheme, is to use a variable time step, as described in Section \ref{extension-EFT12}.


\section{FIC FORM OF THE SEMIDISCRETE PARABOLIC EQUATION IN MECHANICS} \label{fic-form-transient-equation}

The differential form of the parabolic equation in mechanics is found by establishing the balance of fluxes (heat, momentum, forces, etc.) in an infinitesimal space-time domain belonging to $[\Omega ,t]$ \cite{Belytschko14,Hughes87,ZiTaZhu05}. From the governing differential equation and the boundary conditions, the discretized form of Eq.\eqref{first-order-transient-1} can be found using the FEM, for instance \cite{Hughes87,ZiTa05}.

A higher order form of the parabolic equation can be found using a Finite Element Calculus (FIC) procedure, as described in \cite{On24}. The general FIC-Time form of the semidiscrete parabolic equation can be written as
\begin{equation}\label{fic-form-transient-equation-3}
{\bf r} \pm \delta_{0i}{\bf r}_i \pm  \delta_{1j}\Delta t \dot{\bf r}_j \pm  \delta_{2k}\Delta t^2 \ddot{\bf r}_k \pm \ldots = {\bf 0}
\end{equation}
where $\delta_{0i},\, \delta_{1j},\, \delta_{2k}$, etc.  are  non-dimensional time stabilization parameters, ${\bf r}_i,\, \dot{\bf r}_j,\ddot{\bf r}_k$, etc. denote values of the residual $\bf r$ and its time derivatives associated with the parameters $\delta_{0i},\, \delta_{1j},\, \delta_{2k}$, etc., respectively, with
\begin{equation}\label{fic-form-transient-equation-2}
{\bf r}:= {\bf C} \dot {\bf a} + {\bf K} {\bf a}  - {\bf f}
\end{equation}

Note that {\bf r} vanishes for the exact solution. In Eq.\eqref{fic-form-transient-equation-3} the  sum criterion for repeated indices applies.

For heat conduction problems Eq.\eqref{fic-form-transient-equation-3} can be found by establishing the balance of heat fluxes in a space-time domain of dimensions $[d\Omega , \tau ]$ where $d\Omega = dxdy$ in 2D and $\tau$ is a finite time interval defined as a proportion of the time step size $\Delta t$ of the time integration scheme chosen. In practice, $\tau$ is expressed as $\tau =\delta \Delta t$, where the time stabilization parameter $\delta$  governs the stability of the time integration scheme. Details of the derivation of Eq.\eqref{fic-form-transient-equation-3} can be found in \cite{On98,OnMa99,On22,On24}.

The sign and the terms within the residuals terms in Eq.\eqref{fic-form-transient-equation-3} can be formulated in a number of ways, each one leading to a different FIC-Time integration scheme.  Note that when the time stabilization parameters in Eq.\eqref{fic-form-transient-equation-3} are all zero, the infinitesimal form of the parabolic equation given by Eq.\eqref{first-order-transient-1} is obtained.

The time stabilization parameters $\delta_{pr}$ of Eq.\eqref{fic-form-transient-equation-3} can be expressed in terms of the time interval, $\tau$ and the time increment $\Delta t$ used for solving Eq.\eqref{fic-form-transient-equation-3} as
\begin{equation}\label{fic-form-transient-equation-4}
\delta_{pr} = \frac{1}{(1+p)!}\left(\frac{\tau}{\Delta t}\right)^r \quad , \quad p,r=0,1,2\ldots
\end{equation}

For the EFT12  integration scheme presented in this work the time stabilization parameters $\delta_{pr}$ are identified simply by a single parameter $\delta$.

The FIC approach  has been successfully used in computational mechanics in conjunction with standard numerical methods, such as the FEM and meshless procedures, for finding accurate and stable solutions for steady state and transient problems, such as convection–diffusion–radiation, fluid dynamics and quasi-incompressible solids, among others \cite{Cotela17,Felippa07,Felippa18,Maso22,On98,OnMa99,OnGa01,OnPeMi01,Onetal03,Onate04,
Onate04b,Boroomand05,OnMiHa06,OnVaGa07,OnMiNa16,OnNaMi17,On22,On24,Puigferrat20}.


\section{FIRST ORDER EXPLICIT TIME INTEGRATION SCHEME WITH LARGE TIME STEP STEPS}\label{first-order-explicit-time-integration}

The authors have presented in \cite{On24} a family of first, second and third order explicit FIC-Time (EFT) integration schemes that have advantages versus standard FE, AB and RK schemes of the same order in terms of a larger time step size and a smaller number of matrix-vector-multiplications.

In this work we will focus in the advantages of the first order EFT12 scheme, presented in \cite{On24}, for solving parabolic problems using irregular grids. We note that the same advantages, in terms of accuracy and speed-up versus the FE scheme, can be found in all the first order EFT schemes presented in \cite{On24}, that have a quasi identical behavior to the EFT12 scheme.

\subsection{EFT12 explicit time integration scheme}\label{eft12explicit-time-integration-scheme}

The EFT12 scheme can be written as \cite{On24}
\begin{equation}
{\bf C}\dot {\bf a}^{n+1} + {\bf K} {\bf a}^n  - {\bf f}^n-\delta \left[ {\bf C}\tilde{\dot {\bf a}}^{n+1} + {\bf K} {\bf a}^n  - {\bf f}^n  \right]=0
\label{eft12scheme-1}
\end{equation}

The first derivatives $\dot {\bf a}^{n+1}$ and $\tilde{\dot {\bf a}}^{n+1}$ in Eq.\eqref{eft12scheme-1} are computed as
\begin{equation}
\begin{split}
\dot {\bf a}^{n+1} &= \frac{1}{\Delta t }\left[\frac{3}{2}{\bf a}^{n+1} - 2{\bf a}^n + \frac{1}{2} {\bf a}^{n-1}    \right]\\
\tilde{\dot {\bf a}}^{n+1} &= \frac{1}{2\Delta t } [{\bf a}^{n+1} - {\bf a}^{n-1} ]
\end{split}
\label{eft12scheme-2}
\end{equation}

We note that $\dot {\bf a}^{n+1}$ is computed by a higher order first derivative scheme, while $\tilde{\dot {\bf a}}^{n+1}$ is computed by a backward derivative.

Substituting Eqs.\eqref{eft12scheme-2} into \eqref{eft12scheme-1} gives

\begin{equation}
\frac{1}{\Delta t }{\bf C} \left( \frac{3}{2}{\bf a}^{n+1} - 2 {\bf a}^n  + \frac{1}{2}{\bf a}^{n-1}   \right) + {\bf K} {\bf a}^n  - {\bf f}^n - \delta   \left[\frac{1}{2\Delta t } {\bf C}   ({\bf a}^{n+1} - {\bf a}^{n-1}) + {\bf K} {\bf a}^n  - {\bf f}^n \right] =0
\label{eft12scheme}
\end{equation}

Eq.\eqref{eft12scheme} is a particular case of Eq.\eqref{fic-form-transient-equation-3}, when the first FIC-Time term is retained only.

% Substituting Eq.(\ref{appB-1}) and (\ref{appB-2}) into Eq.\eqref{first-order-explicit-time-integration-0bis} and making $ \tau_{11} =\delta \Delta t$ gives the explicit form
%\begin{equation}\label{appB-3}
%\frac{1}{\Delta t }{\bf C} \left( \frac{3}{2}{\bf a}^{n-1} - 2 {\bf a}^n  + \frac{1}{2}{\bf a}^{n-1}   \right) + {\bf K} {\bf a}^n  - {\bf f}^n - \delta  \Delta t \left[\frac{1}{2\Delta t } {\bf C}   ({\bf a}^{n+1} - {\bf a}^{n-1}) + {\bf K} {\bf a}^n  - {\bf f}^n \right] =0
%\end{equation}

Grouping terms we can express Eq.\eqref{eft12scheme} as
\begin{equation}
a_1 {\bf C}{\bf a}^{n+1} + [b_1 {\bf C}  + \Delta t c_1 {\bf K}] {\bf a}^n + d_1 {\bf C} {\bf a}^{n-1} - \Delta t c_1 {\bf f}^n = {\bf 0}
\label{first-order-explicit-time-integration-1}
\end{equation}
where $a_1,b_1,c_1$ and $d_1$ are parameters that characterize the EFT12 scheme. From observation of Eqs.\eqref{eft12scheme} and \eqref{first-order-explicit-time-integration-1} we deduce the expression of these parameters as
\begin{equation}
a_1 = \frac{1}{2}(3-\delta )\quad , \quad b_1 = -2\quad , \quad c_1 = 1-\delta \quad , \quad d_1 = \frac{1}{2}(1+\delta )
\label{first-order-explicit-time-integration-2}
\end{equation}


The explicit computation of ${\bf a}^{n+1}$ is deduced from Eq.\eqref{first-order-explicit-time-integration-1} as
\begin{equation}
{\bf a}^{n+1}=-\frac{1}{a_1} [b_1 + \Delta t c_1{\bf C}^{-1} {\bf K}]{\bf a}^{n} - \frac{d_1}{a_1} {\bf a}^{n-1} + \Delta t \frac{c_1}{a_1} {\bf C}^{-1} {\bf f}^n
\label{first-order-explicit-time-integration-3}
\end{equation}


Substituting Eq.\eqref{first-order-transient-4} into \eqref{first-order-explicit-time-integration-1} and using the orthogonality properties of the eigenvectors ${\pmb\phi}_i$ versus the matrices ${\bf C}$ and ${\bf K}$,  we can obtain the decoupled equation for the $i$-th modal amplitude as (for ${\bf f}^n={\bf 0}$)
\begin{equation}
a_1{y}^{n+1}_i + [b_1 + c_1 p_i] {y}^n_i + d_1 {y}^{n-1}_i=0
\label{first-order-explicit-time-integration-4}
\end{equation}
where $p_i$ is given by Eq.\eqref{first-order-transient-7}.

The characteristic equation for the modal amplification factor, $\gamma_i$, is obtained as
\begin{subequations}
\begin{equation}
A \gamma_i^2 + B \gamma_i + C=0
\label{first-order-explicit-time-integration-5a}
\end{equation}
where
\begin{equation}
A = a_1\quad , \quad B = b_1 + c_1 p_i \quad , \quad C = d_1
\label{first-order-explicit-time-integration-5b}
\end{equation}
\label{first-order-explicit-time-integration-5}
\end{subequations}

The stability conditions for the coefficients $A,B$ and $C$, ensuring $-1<\gamma_i <1$, can be found via the  Routh-Hurwitz criterion \cite{Belytschko14,Hughes87,Routh-Hurwitz,ZiTaZhu05}. The three stability conditions are the following:

\begin{enumerate}

\item $A+B+C > 0$
\begin{equation}
(a_1+b_1+c_1 + d_1) + c_1p_i >0
\label{first-order-explicit-time-integration-6}
\end{equation}

This gives
\begin{equation}
\begin{split}
a_1 + b_1+c_1 & \ge 0\\
c_1&>0
\end{split}
\label{first-order-explicit-time-integration-7}
\end{equation}

\item $A-C > 0$
\begin{equation}
a_1 - d_1>0
\label{first-order-explicit-time-integration-8}
\end{equation}

\item $A-B+C> 0$
\begin{equation}
 (a_1-b_1+d_1) - c_1 p_i>0
 \label{first-order-explicit-time-integration-9}
\end{equation}

\end{enumerate}

Eq.\eqref{first-order-explicit-time-integration-9} gives the condition for $p_i$ as
\begin{equation}\label{first-order-explicit-time-integration-10}
p_i < \frac{a_1-b_1+d_1}{c_1}
\end{equation}

Condition \eqref{first-order-explicit-time-integration-10} must be satisfied for all the modes. Recalling that $p_i = \Delta t \lambda_i$ and particularizing Eq.\eqref{first-order-explicit-time-integration-10} for the highest mode ($i=N$) yields the allowable time increment for the EFT12 scheme, as

\begin{equation}\label{first-order-explicit-time-integration-11}
\Delta t < \Delta t_l = \frac{p_N}{\lambda_N}=  \frac{a_1-b_1+d_1}{c_1 \lambda_N}
\end{equation}
where $\Delta t_l$ is the limit time step ensuring a stable transient solution.

%Comparing the expression of $\Delta t_l$ of Eq.\eqref{first-order-explicit-time-integration-11} with that of Eq.\eqref{first-order-transient-11} gives the speed-up (SUP) of the new explicit FIC-Time schemes versus the forward Euler scheme (SUP$^{\rm FE}$) as


\subsection{Critical values of the stabilization parameter and of the time step size for the EFT12 scheme}

%The residual vector {\bf r} in \eqref{first-order-explicit-time-integration-0bis} is computed as
%\begin{equation}
% {\bf r} = \frac{1}{\Delta t} {\bf C}\left( \frac{3}{2} {\bf a}^{n+1}- 2{\bf a}^n + \frac{1}{2} {\bf a}^{n-1} \right) + {\bf K} {\bf a}^n  - {\bf f}^n
% \label{appB-1}
%\end{equation}
%
%On the other hand, the  residual $ {\bf r}_1$  in Eq.\eqref{first-order-explicit-time-integration-0bis} is computed as
%
%\begin{equation}\label{appB-2}
%{\bf r}_1 = \frac{1}{2\Delta t} {\bf C} ({\bf a}^{n+1} - {\bf a}^{n-1}) +  {\bf K} {\bf a}^n  - {\bf f}^n
%\end{equation}

Substituting the expression of $a_1$, $b_1$ and $d_1$ of Eq.\eqref{first-order-explicit-time-integration-2} into Eqs.\eqref{first-order-explicit-time-integration-10} and  \eqref{first-order-explicit-time-integration-11} gives
\begin{equation}
\label{critical-values-1}
p_N < \frac{4}{1-\delta}\quad , \quad \Delta t < \frac{4}{(1-\delta ) \lambda_N}
\end{equation}

Comparing Eqs.\eqref{critical-values-1} and \eqref{first-order-transient-11} the speed-up (SUP) of the EFT12 scheme versus the FE scheme is obtained as
\begin{equation}
\label{first-order-explicit-time-integration-12}
\hbox{SUP}^{\rm FE}= \frac{2}{1-\delta}
\end{equation}
and, hence,

\begin{equation}\label{first-order-explicit-time-integration-13}
\Delta t\le  \frac{2}{\lambda_N} \textrm{SUP}^{\rm FE}
\end{equation}

\begin{description}
  \item[{\bf Remark 1}.] The values of the stabilization parameter $\delta$ range from $-1< \delta < 1$. For  $\delta =0$ the limit time step for the EFT12 scheme is $\Delta t < 4/ \lambda_N$, which is twice the critical step value for the FE scheme.
The computational overhead for the EFT12 scheme is the need for introducing in the computation of $\textbf{a}^{n+1}$ the value of the unknowns at time $n-1$ (see Eq.\eqref{first-order-explicit-time-integration-3}), which is a cheap operation.

\end{description}

We will derive next the expression of the  critical value of the stabilization parameter and the corresponding critical time step that ensures a real value of the first amplification factor, as a necessary condition for obtaining a non oscillatory transient numerical solution (see the last paragraph of Section \ref{first-order-transient}).


The characteristic equation for the first mode is obtained from Eq.\eqref{first-order-explicit-time-integration-4} (for $i=1$) using the expressions of $a_1, b_1, c_1$ and $d_1$ of Eq.\eqref{first-order-explicit-time-integration-2}, as
\begin{subequations}\label{appB-5}
\begin{equation}\label{appB-5a}
(3-\delta )\gamma_1^2 - 4R_1^{1/2} \gamma_1 + (1+\delta ) =0
\end{equation}
with
\begin{equation}\label{appB-5b}
R_1 =(1-2r_1)^2
\end{equation}
where the eigenvalue ratio $r_1$ is defined as
\begin{equation}\label{appB-5c}
R_1 =\lambda_1 /\lambda_N
\end{equation}
\end{subequations}

From Eq.(\ref{appB-5}a) we deduce, after little algebra
\begin{equation}\label{appB-6}
\gamma_1 = \frac{2}{(3-\delta )} \left[ R_1^{1/2} + (4R_1 - (1+\delta )(3-\delta ))^{1/2}  \right]
\end{equation}

The critical time stabilization parameter $\delta_c$ is obtained by making zero the discriminant of Eq.\eqref{appB-6}, i.e.
\begin{equation}\label{appB-7}
4R_1 -(1+\delta_c ) (3-\delta_c )=0
\end{equation}
from Eq.\eqref{appB-7} we can obtain
\begin{subequations}
\begin{equation}\label{appB-8a}
\delta_c = 1 -2 (1-R_1) ^{1/2} = 1-2 G_1
\end{equation}
with
\begin{equation}\label{appB-8b}
 G_1 = (1-R_1) ^{1/2}
\end{equation}
\label{appa1-8}
\end{subequations}

The critical value of $\gamma_1$ is obtained from Eq.\eqref{appB-6} substituting $\delta$ by $\delta_c$ as

\begin{equation}\label{appB-9}
(\gamma_1)_c = \frac{2R_1^{1/2} }{(3-\delta_c )}
\end{equation}

The \emph{critical speed-up} versus the FE scheme is obtained as

\begin{equation}\label{appB-10}
\hbox{SUP}_{\rm c}^{\rm FE}= \frac{2}{1-\delta_c} = \frac{1}{G_1}= G_1^{-1}
\end{equation}

It is interesting that the expression of the error towards the steady-state solution per time step, $\varepsilon$, coincides with Eq.\eqref{appB-5a}, substituting $\gamma_1$ by $\varepsilon_1$ \cite{On24}. The minimum error is given by an expression identical to Eq.\eqref{appB-9}, i.e.

\begin{equation}\label{appB-11}
\varepsilon_1 =  \frac{2R_1^{1/2} }{(3-\delta_c )}
\end{equation}

Substituting $\delta_c$ from of \eqref{appB-8a} into \eqref{appB-11} gives
\begin{equation}\label{appB-12}
\varepsilon_1 =  \frac{1- 2r_1 }{1+G_1}
\end{equation}

From the expressions of $R_1$ and $G_1$ of Eqs.\eqref{appB-5a} and \eqref{appB-8b} we can deduce
\begin{equation}\label{appB-13}
G_1 = (1- (1- 2r_1)^2)^{1/2} = 2 (r_1 -r_1^2)^{1/2}
\end{equation}

For small values of $r_1$, then
\begin{equation}\label{appB-14}
G_1  \simeq 2 r_1^{1/2}
\end{equation}

Substituting \eqref{appB-14} into \eqref{appB-12} gives

\begin{equation}\label{appB-15}
\varepsilon_1 =  \frac{1- 2r_1 }{1+2 r_1^{1/2}} = \frac{(1- 2r_1)(1-2 r_1^{1/2})}{1-4r_1}= \frac{(1- 2r_1^{1/2})- 2r_1 (1-2 r_1)}{1-4r_1}\simeq 1 - 2 r_1^{1/2}
\end{equation}

Eq.\eqref{appB-15} coincides with the error towards the steady-state solution in explicit time relaxation schemes, based on solving the second order dynamic equation with a central difference scheme and optimal Rayleigh damping \cite{Underwood83}. This property is common to all the first order EFT schemes \cite{On24}. The advantages of the first order EFT schemes in the context of time relaxation procedures are discussed in \cite{On24}.

\begin{description}
  \item[{\bf Remark 2}.] The critical stabilization parameter $\delta_c$  and the corresponding critical time step for the EFT12 scheme, yield the maximum speed-up versus the FE scheme, while preserving an accurate transient solution. An exception to this is the solution near Dirichlet boundaries, that can show oscillations in the initial time steps of the transient solution, when large time steps are used. These oscillations can be eliminated using the EFT12 scheme with a variable time step, as explained in Section \ref{extension-EFT12}.

\end{description}


\subsection{Extension of the EFT12 scheme to variable time steps}\label{extension-EFT12}

The use of the EFT12 scheme with large time steps can lead to oscillations in the space distribution of the numerical solution in the vicinity of a Dirichlet boundary where the variable is prescribed to a fixed value. This situation is more evident when irregular meshes are used, such as when small elements or cells are used in the vicinity of these boundaries. We have found that these oscillations  vanish if the EFT12 scheme is used with a variable time step. This allows to accurately capture the evolution of the solution during the first steps of the transient analysis, while still using large time steps during most of the transient solution that lead to considerable gains in speed-up versus the FE scheme. We detail below the procedure to use the EFT12 scheme with variable time steps.

\begin{figure}[t]
  \centering
  \includegraphics[width=10cm]{Figure1}
  \caption{Variable time intervals between the solution at times $n-1$, $n$ and $n+1$.}\label{time-intervals}
\end{figure}


Figure \ref{time-intervals} shows two different time intervals $\Delta t_l$ and $\Delta t_2$, between the solution at times $n-1$, $n$ and  $n+1$. The first derivative of a variable $\phi (t)$ at time $n+1$ can be computed as
\begin{equation}
\begin{split}
\dot\phi^{n+1} & = \dot\phi^n + \Delta t_1 \ddot\phi^n = \frac{\phi^{n+1} - \phi^{n-1}}{(\Delta t_1 + \Delta t_2)} + \Delta t_1 \frac{2(\dot\phi^{n+1/2}- \dot\phi^{n-1/2}) }{ (\Delta t_1 + \Delta t_2)} =\\
 & =  \frac{\phi^{n+1} - \phi^{n-1}}{(\Delta t_1 + \Delta t_2)} +\frac{2\Delta t_1}{(\Delta t_1 + \Delta t_2) } \left[  \frac{\phi^{n+1} - \phi^{n}}{ \Delta t_1} - \frac{\phi^{n} - \phi^{n-1}}{ \Delta t_2}  \right]= \\
 & =   \frac{\phi^{n+1} - \phi^{n-1}}{(\Delta t_1 + \Delta t_2)}   + \frac{2}{\Delta t_2(\Delta t_1 + \Delta t_2)} \left[\Delta t_2 (\phi^{n+1} - \phi^{n}) - \Delta t_1 (\phi^{n} - \phi^{n-1})   \right]\nonumber
\label{extensioneft12-1}
 \end{split}
\end{equation}

After grouping terms we can write
\begin{equation}\label{extensioneft12-2}
 \dot\phi^{n+1} = \frac{3}{(\Delta t_1 + \Delta t_2)}\phi^{n+1}  - \frac{2}{\Delta t_2}\phi^n + \left[ \frac{- \Delta t_2 + 2\Delta t_1 }{\Delta t_2 (\Delta t_1 + \Delta t_2)}  \right] \phi^{n-1}
\end{equation}

Similarly for the first order derivative $\tilde{\dot\phi}^{n+1}$ in Eq.\eqref{eft12scheme-2} we can write
\begin{equation}
\tilde{\dot\phi}^{n+1} = \frac{1}{\Delta t_1 + \Delta t_2} [\phi^{n+1} - \phi^{n-1}]
\label{extensioneft12-3}
\end{equation}


Taking into account Eqs.\eqref{extensioneft12-2} and \eqref{extensioneft12-3} for computing the first order derivatives $\dot {\bf a}^{n+1}$ and $\tilde{\dot {\bf a}}^{n+1}$ in Eq.\eqref{eft12scheme-1}, the expression of the EFT12 scheme, accounting for variable time steps, is found as
\begin{equation}\label{extensioneft12-4}
{\bf C}  \left[\frac{3}{\Delta t_{n+1}} {\bf a}^{n+1} - \frac{2}{\Delta t_2} {\bf a}^n +  \frac{1}{\Delta t_{n-1}} {\bf a}^{n-1}  \right] + {\bf K} {\bf a}^n  - {\bf f}^n - \delta \left[ \frac{1}{\Delta t_{n+1}} {\bf C} ({\bf a}^{n+1} - {\bf a}^{n-1}) + {\bf K} {\bf a}^n  - {\bf f}^n     \right]=0
\end{equation}
where
\begin{equation}\label{extensioneft12-5}
 \Delta t_{n+1}  = \Delta t_1 + \Delta t_2 \quad , \quad    \Delta t_{n-1}  = \frac{(\Delta t_1 + \Delta t_2) \Delta t_2}{2\Delta t_1 - \Delta t_2 }
\end{equation}

Clearly, for $\Delta t_1= \Delta t_2 = \Delta t$ we find
\begin{equation}\label{extensioneft12-6}
\dot\phi^{n+1} = \frac{1}{\Delta t} \left[ \frac{3}{2} \phi^{n+1} - 2  \phi^n + \frac{1}{2} \phi^{n-1} \right]\quad \hbox{and} \quad \tilde{\dot\phi}^{n+1} = \frac{1}{2\Delta t}(\phi^{n+1} - \phi^{n-1})
\end{equation}
which are the standard expression of the first derivatives in the EFT12 scheme (see Eq.\eqref{eft12scheme-2}).

The characteristic equation for the modal amplification factor $\gamma_i$ is given by Eq.\eqref{first-order-explicit-time-integration-5a} with  $A$, $B$ and $C$ given by Eq.\eqref{first-order-explicit-time-integration-5b} with
\begin{equation}\label{extensioneft12-8}
a_i =(3-\delta ), \quad b_i=-2 \frac{\Delta t_{n+1}}{\Delta t_2}, \quad c_i = (1-\delta ), \quad d_i = \frac{\Delta t_{n+1}}{\Delta t_{n+1}}+\delta
\end{equation}

Following the procedure described in Section \ref{eft12explicit-time-integration-scheme} the limit time step for $\Delta t_2$ can be found as
\begin{equation}\label{extensioneft12-9}
\Delta t_2= \frac{4}{(1-\delta ) \lambda_N}
\end{equation}
which coincides with the limit time step of the standard EFT12 scheme (Eq.\eqref{critical-values-1}).

\section{EXAMPLES}

\subsection{One dimensional transient heat conduction solved with 2-noded finite elements}

\subsubsection{1D transient heat conduction. Regular mesh of 10 2-noded elements}\label{subsection5.5.1}

The first example shows the performance of the EFT12 scheme in the solution of the transient heat conduction problem in a one-dimensional domain. Figure \ref{geometry-domain} shows the geometry of the domain, the Dirichlet boundary conditions at both ends, the material properties, the initial temperature and the uniform discretization in ten 2-noded linear finite elements. The figure also shows the smallest and largest eigenvalues of the semidiscrete system and their ratio, the limit time step for the forward Euler (FE) scheme ($\Delta t^{\mathrm{FE}}_l =0.005125$s), the critical stabilization parameter (${\delta }_c = 0.374459$) and the critical time step for the EFT12 scheme ($\Delta t^{\mathrm{EFT}}_c=0.016387$s). The maximum speed-up versus the FE scheme is  ${\mathrm{SUP}}^{\mathrm{FE}}_c = 3.197$.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure2}
  \caption{1D domain of unit length ($L= 1$m) with uniform discretization in ten 2-noded linear finite elements. Smallest and largest eigenvalues and their ratio. Limit time step for the FE scheme. Critical stabilization parameter and critical time step for the EFT12 scheme. Reference solution for the time variation of the temperature at nodes $2$, $5$ and $9$ obtained with the FE scheme using  $\Delta t =0.001$s. Temperature increments at the same nodes computed each $0.005$s.}\label{geometry-domain}
\end{figure}

The lower left graph of Figure \ref{geometry-domain} show the reference solution for the time variation of the temperature at nodes $2$, $5$ and $9$ obtained with the FE scheme using $\Delta t =0.001$s. The right lower curves show the temperature increments at the same nodes computed each $0.005$s. These increments show the expected distribution and will be taken as a reference for comparison with the results of the EFT12 scheme.

Figure \ref{1D-uniform-discretization-3} shows the time variation of the temperature at nodes $2$, $5$ and $9$ obtained with the EFT12  scheme using a uniform time step of $\Delta t^{\mathrm{EFT}}=  0.009$s, a value smaller than the critical one, corresponding to   $\delta = - 0.139$. The speed-up versus the FE solution is  ${\mathrm{SUP}}^{\mathrm{FE}} = 1.8$. The temperature increments at the same nodes computed each $0.005$ s show small oscillations near the left Dirichlet boundary due to the large time step used.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure3}
  \caption{1D domain with uniform discretization. Time variation of the temperature at nodes $2$, $5$ and $9$ obtained with the EFT12  scheme using a uniform time step of $\Delta t^{\mathrm{EFT}}=  0.009$s, corresponding to   $\delta = - 0.139$   and  ${\mathrm{SUP}}^{\mathrm{FE}} = 1.8$. The temperature increments at the same nodes computed each $0.005$s show oscillations near the left Dirichlet boundary due to the large time step used.}\label{1D-uniform-discretization-3}
\end{figure}

Figure \ref{1D-uniform-discretization-4} shows the time variation of the temperature at nodes $2$, $5$ and $9$ obtained with the EFT12 scheme with the variable time step pattern shown in the lower graph of the figure that has been implemented as follows
\begin{equation}
\begin{split}
& 0<t\le 0.10s \quad , \quad \Delta t = 0.005 s\\
& 0.010s < t \le 0.20s \quad , \quad \Delta t = 0.005 + 0.05 (t - 0.10)(s)\\
& 0.015s < t \quad , \quad \Delta t =0.0150 s
\end{split}
\end{equation}

The time step increases in the initial stages of the analysis from $0.005$s (which is slightly smaller than the limit time step for the FE scheme) to reach the value of $0.00150$s and a speed-up versus the FE scheme of ${\mathrm{SUP}}^{\mathrm{FE}} = 2.92$.
The upper right figure shows that the temperature increments at nodes $2$, $5$ and $9$ computed each $0.005$s have no oscillations, despite the large time step used during most of the transient analysis. This simple example demonstrates the advantages of using the EFT12 scheme with a variable time scheme in parabolic problems involving Dirichlet conditions. This conclusion  is evidenced in all the more complex problems presented in the next sections.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure4}
  \caption{1D domain with uniform discretization. Time variation of the temperature at nodes $2$, $5$ and $9$ obtained with the EFT12 scheme with a variable time step. The lower figure shows the time step variation along time to reach a time step of $\Delta t^{\mathrm{EFT}}=  0.150$,  corresponding to  $\delta = 0.317$  and   ${\mathrm{SUP}}^{\mathrm{FE}} = 2.92$. The lower right figure shows that the temperature increments at the same nodes computed each $0.005$s have no oscillations, despite the large time step used.}\label{1D-uniform-discretization-4}
\end{figure}


\subsubsection{1D transient heat conduction. Irregular mesh of twenty 2-noded elements with a refined zone adjacent to the left end}

Figure \ref{1D-transient-heat-conduction-irregular-mesh} shows the one dimensional domain of the previous example discretized with a mesh containing ten linear elements of $l=0.1$m and a refined mesh of ten elements of $l=0.01$m near the left end. The figure shows the smallest and largest eigenvalues and the limit time step for the FE scheme ($\Delta t^{\mathrm{FE}}_l = 0.00005125$). Note that the small elements near the left end lead to a limit time step for the FE scheme that is two orders of magnitude smaller than for the regular mesh case (Section \ref{subsection5.5.1}). Figure \ref{1D-transient-heat-conduction-irregular-mesh} also shows the critical stabilization parameter (${\delta }_c = 0.936565$) and the critical time step for the EFT12 scheme ($\Delta t^{\mathrm{EFT}}_c=0.001616$s). The maximum speed-up versus the FE scheme is ${\mathrm{SUP}}^{\mathrm{FE}}_{\mathrm{c}} = 31.53$. The left graph show the EFT12 results with a variable time step  and the reference solution for the temperature evolution at nodes $2$, $5$ and $9$ using the FE scheme with $\Delta t =0.000001$s. Both results are practically coincident.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure5}
  \caption{1D domain discretized with a mesh containing ten linear elements of $l=0.1$m and a refined mesh of ten elements of $l=0.01$m near the left end. Smallest and largest eigenvalues and limit time step for the FE scheme. The left graphs show the EFT12 results with a variable time step and the reference solution for the temperature evolution at nodes $2$, $5$ and $9$ using the FE scheme with $\Delta t =0.000001$s. Both results are  coincident. The lower figure shows the time step variation to reach a time step of $\Delta t^{\mathrm{EFT}}_c=0.001616$s,  corresponding to  ${\delta }_c = 0.936565$  and   ${\mathrm{SUP}}^{\mathrm{FE}} =  31.53$. The upper right figure shows that the temperature increments at nodes $2$, $5$ and $9$ computed each 0.005 s have no oscillations.}\label{1D-transient-heat-conduction-irregular-mesh}
\end{figure}

The lower figure shows the time step variation used here, that has been implemented as follows:
\begin{equation}
\begin{split}
& 0<t\le 0.008s \quad , \quad \Delta t = 0.00005 s\\
& 0.008s < t \le 0.015s \quad , \quad \Delta t = 0.00005 + 0.223714 (t - 0.08)(s)\\
& 0.015s < t \quad , \quad \Delta t =0.001616 s
\end{split}
\end{equation}

The time step increases in the initial stages of the analysis from $0.00005$s (which is slightly smaller than the limit time step for the FE scheme) to reach the critical value of $0.001616$s and a speed-up versus the FE scheme of   ${\mathrm{SUP}}^{\mathrm{FE}} = 31.53$. The upper right figure shows that the temperature increments at nodes $2$, $5$ and $9$ computed each $0.005$s have no oscillations, despite the very large time step used during most of the transient analysis.

\subsubsection{1D transient heat conduction. Regular mesh of ten 2-noded elements containing a very refined interior zone}

Figure \ref{1D-transient-heat-conduction-regular-mesh} shows the 1D domain discretized with a mesh containing ten linear elements of $l=0.1$m and a the zone where a small element of $l=0.01$m has been included. The figure shows the smallest and largest eigenvalues and their ratio, the limit time step for the forward Euler (FE) scheme ($\Delta t^{\mathrm{FE}}_l= 0.000498$), the critical stabilization parameter (${\delta }_c = 0.8028$) and the critical time step for the EFT12 scheme ($\Delta t^{\mathrm{EFT}}_c=0.005051$s). The maximum speed-up versus the FE scheme is  ${\mathrm{SUP}}^{\mathrm{FE}}_{\mathrm{c}} = 10.14$. The left graphs show the EFT12 results with the variable time step pattern detailed below and the reference solution for the temperature evolution at nodes $2$, $5$ and $9$ using the FE scheme with $\Delta t =0.00001$s. There is a good coincidence between both results.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure6}
  \caption{1D domain discretized with a mesh of ten linear elements of $l=0.1$m and a refined central zone containing one small element of $l=0.01$m. Smallest and largest eigenvalues and limit time step for the FE scheme. The left graphs show the EFT12 results with a variable time step and the reference solution for the temperature evolution at nodes $2$, $5$ and $9$ using the FE scheme with  $\Delta t =0.00001$s. Both results are practically coincident.
   The lower figure shows the time step variation to reach a time step of $\Delta t^{\mathrm{EFT}}_c=0.005051$s,  corresponding to  $\delta  = 0.8028$  and   ${\mathrm{SUP}}^{\mathrm{FE}} =  10.14$. The upper right figure shows that the temperature increments at nodes $2$, $5$ and $9$ computed each 0.005s have no oscillations, despite the large time step used.}\label{1D-transient-heat-conduction-regular-mesh}
\end{figure}

The lower graph in Figure \ref{1D-transient-heat-conduction-regular-mesh} shows the time step variation used, that has been implemented as follows:
\begin{equation}
\begin{split}
& 0<t\le 0.008s \quad , \quad \Delta t = 0.00045 s\\
& 0.008s < t \le 0.015s \quad , \quad \Delta t = 0.00045 + 0.065(t - 0.08) (s)\\
& 0.015s < t \quad , \quad \Delta t =0.05050 s
\end{split}
\end{equation}

 The time step increases from $0.00045$s (a value close to the limit time step of the FE scheme to reach  $0.0015$s and a speed-up versus the FE scheme of ${\mathrm{SUP}}^{\mathrm{FE}} = 10.14$ (basically equal to the critical value). The upper right figure shows that the temperature increments at nodes $2$, $5$ and $9$ computed each $0.005$s have no oscillations, despite the large time step used during most of the transient analysis.


\subsection{Two dimensional (2D) transient heat conduction in square domain with internal circular cavity}\label{two-dimensional}

\subsubsection{Transient heat conduction in square domain with an internal circular hole}\label{two-dimensional-unstructured mesh}

	Figure \ref{fig7} shows the geometry, the problem parameters and the unstructured finite element mesh of $870$ 3-noded triangles with an average element size of $0.00$5m. The figure also shows the smallest and largest eigenvalues and their ratio, the limit time step for the FE scheme ($\Delta t^{\mathrm{FE}} =  0.035$s),  the critical stabilization parameter (${\delta }_c  = 0.890$) and the critical time step for the EFT12 scheme ($\Delta t^{\mathrm{EFT}}_c= 0.642$s). The maximum speed-up versus the FE scheme is ${\mathrm{SUP}}^{\mathrm{FE}}_{\mathrm{c}} =  18.34$. Figure \ref{fig7} also shows the contours of the reference solution obtained in this case with the backward Euler scheme using a conjugated gradient solver and $\Delta t =0.0005$s.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure7}
  \caption{Square domain with internal hole. Geometry (dimensions in meters), problem parameters, boundary conditions and FEM mesh. Smallest and largest eigenvalues and their ratio. Limit time step for the FE scheme. Critical stabilization parameter and critical time step for the EFT12 scheme. Contours of the reference solution obtained with the backward Euler scheme using a conjugated gradient solver and $\Delta t =0.0005$s.}\label{fig7}
\end{figure}


Figure \ref{fig8} shows the temperature distribution along lines AB and CD at different times and the temperature evolution at points E and F obtained with the EFT12 scheme and a constant time step, slightly smaller then the critical one of $\Delta t =  0.60$s, corresponding to a value of $\delta=   0.8824$ and a speed-up versus the FE scheme of ${\mathrm{SUP}}^{\mathrm{FE}} =  17.0$. The EFT12 results are
compared with the reference solution of Figure \ref{fig7}. The incorrect results in the initial phase of the solution are due to the effect of the large time step used in the time evolution of the temperature near the lower left Dirichlet boundary.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure8}
  \caption{Square domain with internal hole. Temperature distribution along lines AB and CD at different times and temperature evolution at points E and F obtained with the EFT12 scheme and a constant time step of $\Delta t =  0.60$s corresponding to ${\mathrm{SUP}}^{\mathrm{FE}}_{\mathrm{c}} =  17.0$. The incorrect initial results are due to the effect of the large time steps used in the time evolution of the solution near the lower left Dirichlet boundary.}\label{fig8}
\end{figure}


The numerical results for the temperature evolution in time can be considerably improved, while maintaining a good speed-up versus the FE scheme, using a variable time step.
Figure \ref{fig9} shows the temperature distribution along lines AB and CD at different times and the temperature evolution at points E and F obtained with the EFT12 scheme and a variable time step defined as follows
\begin{equation}
\begin{split}
& 0<t\le 5s \quad , \quad \Delta t = 0.03 s\\
& 5s < t \le 20s \quad , \quad \Delta t = 0.03 + 0.033(t - 5) (s)\\
& 20s < t \quad , \quad \Delta t =0.5250 s
\end{split}
\end{equation}

The time step is kept constant during the first 5 seconds (which represents a 2.5\% of the total analysis time) to $\Delta t=0.03$s, a value close to the limit time step for the FE scheme. Then, the time step increases linearly from $0.03$s  to reach $0,5250$s and a speed- up versus the FE scheme of ${\mathrm{SUP}}^{\mathrm{FE}}_{\mathrm{c}} =  15.0$, as shown in the figure. Good results are obtained for the temperature evolution at all points considered, while maintaining a clear advantage, in terms of a faster numerical solution  versus the FE scheme, throughout most of the transient solution.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure9}
  \caption{Square domain with internal hole. Temperature distribution along lines AB and CD at different times and temperature evolution at points E and F obtained with the EFT12 scheme and a variable time step ranging from $0.030$s to $0.5250$s with ${\mathrm{SUP}}^{\mathrm{FE}}_{\mathrm{c}} =  15.0$. The variable time step leads to a correct distribution of the temperature at all times.}\label{fig9}
\end{figure}


\subsubsection{2D transient heat conduction in square domain with circular hole. Unstructured triangular mesh with a refined zone adjacent to the lower right boundary}

Figure \ref{geometry-boundary-material} shows the geometry, the problem parameters and the unstructured finite element mesh of $870$ 3-noded triangles with  average element size of $0.005$m, containing a very refined zone adjacent to the lower right boundary with average element size of $0.001$m. The figure also shows the smallest and largest eigenvalues and the limit time step for the forward Euler (FE) scheme ($\Delta t^{\mathrm{FE}}_l= 0.00074$s), the critical stabilization parameter (${\delta }_c = 0.984093$) and the critical time step for the EFT12 scheme ($\Delta t^{\mathrm{EFT}}_c=0.09297$s). The maximum speed-up versus the FE scheme is ${\mathrm{SUP}}^{\mathrm{FE}}_{\mathrm{c}} = 125.63$.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure10}
  \caption{Square domain with internal hole. Geometry, problem parameters and FEM mesh with a refined zone adjacent to the lower left boundary. Smallest and largest eigenvalues and their ratio. Limit time step for the FE scheme. Critical stabilization parameter and critical time step for the EFT12 scheme. Contours of the reference solution obtained with the backward Euler scheme using a conjugated gradient solver and $\Delta t =0.0005$s.}\label{geometry-boundary-material}
\end{figure}


Figure \ref{temperature-distribution-constant-time} shows the temperature distribution along lines AB and CD at different times and the temperature evolution at points E and F obtained with the EFT12 scheme and a constant time step of $\Delta t = 0.0850$s corresponding to a value of $\delta =  0.9826$  and ${\mathrm{SUP}}^{\mathrm{FE}}= 115.0$. The EFT12 results are compared with a reference solution obtained using the implicit Backward Euler  scheme with a conjugate gradient solver and $\Delta t =0.0005$s. The incorrect results in the initial phase of the solution are due to the effect of the large time step used in the time evolution of the temperature near the lower left Dirichlet boundary.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure11}
  \caption{Square domain with internal hole with a refined mesh adjacent to the lower left boundary. Temperature distribution along lines AB and CD at different times and temperature evolution at points E and F obtained with the EFT12 scheme and a constant time step of $\Delta t = 0.0850$s corresponding to ${\mathrm{SUP}}^{\mathrm{FE}}= 115.0$. The incorrect  results in the initial times are due to the effect of the large time steps used in the time evolution of the temperature near the lower left Dirichlet boundary.}\label{temperature-distribution-constant-time}
\end{figure}


Figure \ref{temperature-correct-distribution-variable-time} shows the temperature distribution along lines AB and CD at different times and the temperature evolution at points E and F obtained with the EFT12 scheme and a variable time step defined as follows
\begin{equation}
\begin{split}
& 0<t\le 5s \quad , \quad \Delta t = 0.0007 s\\
& 5s < t \le 20s \quad , \quad \Delta t = 0.0007 + 0.0085(t - 5)\\
& 20s < t \quad , \quad \Delta t =0.085 s
\end{split}
\end{equation}

The time step is kept constant during the first 5 seconds (which represents a 2.5\% of the total analysis time) to $\Delta t =0.0007$s, a value close to the limit time step for the FE scheme, Then it increases linearly from $0.0007$s to reach $0.085$s and a speed- up versus the FE scheme of ${\mathrm{SUP}}^{\mathrm{FE}}= 115.0$, as shown in the figure. Good results are obtained for the temperature evolution at all points considered. The EFT12 results can be still improved by reducing the maximum time step, as described next.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure12}
  \caption{Square domain with internal hole with a refined zone adjacent to the lower left boundary. Temperature distribution along lines AB and CD at different times and temperature evolution at points E and F obtained with the EFT12 scheme and a variable time step ranging from $0.0007$s to the critical value of $0.085$s with ${\mathrm{SUP}}^{\mathrm{FE}}= 115$. The variable time step leads to a correct distribution of the temperature at all times.}\label{temperature-correct-distribution-variable-time}
\end{figure}

Figure \ref{temperature-accurate-distribution-variable-time} shows the temperature distribution along lines AB and CD at different times and the temperature evolution at points E and F obtained with the EFT12 scheme and a variable time step defined as follows
\begin{equation}
\begin{split}
& 0<t\le 5s \quad , \quad \Delta t = 0.0007 s\\
& 5s < t \le 20s \quad , \quad \Delta t = 0.0007 +0.0592(t - 5)(s)\\
& 20s < t \quad , \quad \Delta t =0.00592 s
\end{split}
\end{equation}

The time step increases linearly in this case from $0.0007$s to reach $0.00592$s and a speed- up versus the FE scheme of ${\mathrm{SUP}}^{\mathrm{FE}} = 80.0$ (corresponding to ${\delta }_c = 0.9750$).  Accurate results are obtained for the temperature evolution at all points considered, despite the very large time step used during most of the analysis.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure13}
  \caption{Square domain with internal hole with a refined zone adjacent to the lower left boundary. Temperature distribution along lines AB and CD at different times and temperature evolution at points E and F obtained with the EFT12 scheme and a variable time step ranging from $0.0007$s to a maximum value of $0.0592$s corresponding to $\delta =  0.9750$  and ${\mathrm{SUP}}^{\mathrm{FE}}=80.0 $. The variable time step leads to an accurate distribution of the temperature at all times.}\label{temperature-accurate-distribution-variable-time}
\end{figure}


Figure \ref{evolution-global-error-temperature} displays the evolution of the global error in the temperature, computed as shown in the figure, obtained with the EFT12 scheme, using a constant time step of $\Delta t = 0.0850$s and different variable time steps. The variable time step solutions yield a lower error in the initial phase of the solution, as expected. The error for all solutions has the same tendency as the temperature distribution  evolves in time. As expected, the smallest error correspond to the variable time solution with a smaller maximum time step shown in Figure \ref{temperature-accurate-distribution-variable-time}, that reaches a maximum time step of $0.0592$s and  ${\mathrm{SUP}}^{\mathrm{FE}} =80$.

\begin{figure}[h!]
  \centering
  \includegraphics[width=9cm]{Figure14}
  \caption{Square domain with internal circular hole with a refined zone adjacent to the lower left boundary. Evolution of the global error in temperature for the EFT12 scheme using a constant time step of $\Delta t = 0.0850$s and different variable time steps given different speed-ups versus the FE scheme.}\label{evolution-global-error-temperature}
\end{figure}



\subsubsection{Transient heat conduction in square domain with circular hole. Unstructured triangular mesh with an internal refined zone}

Figure \ref{Finite-element-mesh-of-872-3-noded-triangles} shows the geometry, the  problem parameters and the unstructured finite element mesh of $878$ 3-noded triangles with average element size of $0.005$m, containing a  refined zone in the interior part of the analysis domain with triangles of an average element size of $0.001$m. The figure also shows the smallest and largest eigenvalues and the limit time step for the FE scheme ($\Delta t^{\mathrm{FE}}= 0.000235$s), the critical stabilization parameter (${\delta }_c = 0.97163$) and the critical time step for the EFT12 scheme ($\Delta t^{\mathrm{EFT}}_c=0.16543$s). The maximum speed-up versus the FE scheme is ${\mathrm{SUP}}^{\mathrm{FE}}_{\mathrm{c}} = 70.40$.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure15}
  \caption{Square domain with internal circular hole. Geometry (dimensions in meters), problem parameters and mesh of $872$ 3-noded triangles containing an internal zone with very small elements. Smallest and largest eigenvalues and their ratio. Limit time step for the FE scheme. Critical stabilization parameter and critical time step for the EFT12 scheme. Contours of the reference solution obtained with the backward Euler scheme using a conjugated gradient solver and  $\Delta t = 0.0005$s.}\label{Finite-element-mesh-of-872-3-noded-triangles}
\end{figure}


Figure \ref{variable-time-correct-distribution} shows the temperature distribution along lines AB and CD at different times and the temperature evolution at points E and F obtained with the EFT12 scheme and a variable time step defined as follows
\begin{equation}
\begin{split}
& 0<t\le 5s \quad , \quad \Delta t = 0.0007 s\\
& 5s < t \le 20s \quad , \quad \Delta t = 0.0007 +0.0.00953(t - 5) (s)\\
& 20s < t \quad , \quad \Delta t =0.15 s
\end{split}
\end{equation}

After an initial period of five seconds with a constant time step of $0.0007$s  this value increases linearly to reach $0.15$s (a value slight smaller than the critical one of $0.16543$s) (corresponding to ${\delta } = 0.969$)  and  a speed-up versus the FE scheme of ${\mathrm{SUP}}^{\mathrm{FE}} = 63.8$.  Accurate results are obtained for the temperature evolution at all points considered using a quite large time step during most of the analysis.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure16}
  \caption{Square domain with circular hole and refined internal zone. Temperature distribution along lines AB and CD at different times and temperature evolution at points E and F obtained with the EFT12 scheme and a variable time step r ranging from $0.0020$s to a maximum value of $0.150$s corresponding  to  ${\mathrm{SUP}}^{\mathrm{FE}}=63.8$. The variable time step yields a correct distribution of the temperature at all times.}\label{variable-time-correct-distribution}
\end{figure}


Figure \ref{circular-hole-refined-internal-zone} shows that more accurate results for the temperature distribution at all times is obtained with the EFT12 scheme and a variable time step that follows the pattern shown in the figure. The maximum time step is $\Delta t = 0.705$s, that gives a speed-up versus the FE solution of ${\mathrm{SUP}}^{\mathrm{FE}} = 30$. This gives still a remarkable gain in terms of computing time versus the FE scheme.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure17}
  \caption{Square domain with circular hole and refined internal zone. Temperature distribution along lines AB and CD at different times and temperature evolution at points E and F obtained with the EFT12 scheme and a variable time step r ranging from $0.0020$s to a maximum value of $0.0705$s corresponding  to  ${\mathrm{SUP}}^{\mathrm{FE}}=30$. The variable time step yields a correct distribution of the temperature at all times.}\label{circular-hole-refined-internal-zone}
\end{figure}

Figure \ref{internal-circular-hole} displays the evolution of the global error in temperature for the EFT12 scheme, computed using different variable time steps. The error for all the solutions has the same tendency. As expected, the smallest error corresponds to the EFT12 solution with a smaller maximum time step, that corresponds to the solution displayed in Figure \ref{circular-hole-refined-internal-zone} with  a speed-up versus the FE scheme of ${\mathrm{SUP}}^{\mathrm{FE}} =30$.

\begin{figure}[h!]
  \centering
  \includegraphics[width=9cm]{Figure18}
  \caption{Square domain with internal circular hole and refined internal zone. Evolution of the global error in temperature for the EFT12 scheme using different variable time steps.}\label{internal-circular-hole}
\end{figure}


\subsection{Three dimensional (3D) transient heat conduction in a cubic domain with spherical cavity}

\subsubsection{Transient heat conduction in cubic domain with a spherical cavity. Unstructured tetrahedral mesh}

Figure \ref{cubic-domain-geometry-boundary-material} shows the geometry, the problem parameters and the unstructured finite element mesh of $73205$ 4-noded tetrahedral elements with average element size of $0.005$m. The figure also displays the smallest and largest eigenvalues and their ratio, the limit time step for the FE scheme ($\Delta t^{\mathrm{FE}}= 0.00626$s), the critical stabilization parameter (${\delta }_c = 0.95370$) and the critical time step for the EFT12 scheme ($\Delta t^{\mathrm{EFT}}_c= 0.27028$s). The maximum speed-up versus the FE scheme is ${\mathrm{SUP}}^{\mathrm{FE}}_{\mathrm{c}} = 43.12$. Figure \ref{cubic-domain-geometry-boundary-material} also shows the contours of the reference solution obtained with the backward Euler scheme using a conjugated gradient solver and $\Delta t =0.0010$s.


\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure19}
  \caption{Cubic domain with internal spherical cavity. Geometry, problem parameters and  mesh of $73205$ 4-noded tetrahedra. Smallest and largest eigenvalues and their ratio. Limit time step for the FE scheme. Critical stabilization parameter and critical time step for the EFT12 scheme. Contours of the reference solution obtained with the backward Euler scheme using a conjugated gradient solver and  $\Delta t =0.0010$s.}\label{cubic-domain-geometry-boundary-material}
\end{figure}


Figure \ref{cubic-domain-temperature-distribution} shows the temperature distribution along lines AB and CD at different times obtained with the EFT12 scheme and different variable time steps, following the patterns shown in the figure. The speed-ups versus the FE scheme range from  ${\mathrm{SUP}}^{\mathrm{FE}}= 43.20$ to  ${\mathrm{SUP}}^{\mathrm{FE}}= 20$. The results evidence, once more, that the EFT12 scheme with a variable time step yields a correct distribution of the temperature in a considerably faster time than the FE scheme with an increased accuracy as the final time step is reduced. Once again, considerable gains in speed-up versus the FE scheme throughout most of the solution are obtained in all cases.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure20}
  \caption{Cubic domain with internal spherical cavity. Temperature distribution along lines AB and CD at different times obtained with the EFT12 scheme and different variable time steps giving speed-ups versus the FE scheme ranging from   ${\mathrm{SUP}}^{\mathrm{FE}}= 43.20$ to   ${\mathrm{SUP}}^{\mathrm{FE}}= 20$. The variable time step yields a correct distribution of the temperature in all cases with an increased accuracy for  smaller values of the speed-up.}\label{cubic-domain-temperature-distribution}
\end{figure}


\subsubsection{Transient heat conduction in cubic domain with spherical cavity. Unstructured tetrahedral mesh with a refined interior zone}

Figure \ref{cubic-domain-spherical-refined-zone} shows the geometry, the problem parameters and the unstructured finite element mesh of $88569$ 4-noded tetrahedra with an average element size of $0.005$m, containing a refined zone in the interior the analysis domain with smaller  tetrahedral elements of an average size of $0.0005$m. The figure also shows the smallest and largest eigenvalues and their ratio, the limit time step for the forward Euler (FE) scheme ($\Delta t^{\mathrm{FE}}=  0.000023677$s), the critical stabilization parameter (${\delta }_c = 0.99715$) and the critical time step for the EFT12 scheme ($\Delta t^{\mathrm{EFT}}_c= 0.0166129$s). The maximum speed-up versus the FE scheme reaches in this case the extraordinary figure of ${\mathrm{SUP}}^{\mathrm{FE}}_{\mathrm{c}} = 701.64$.


\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure21}
  \caption{Cubic domain with internal spherical cavity and refined internal zone. Geometry, problem parameters and  mesh of $88659$ 4-noded tetrahedra. Smallest and largest eigenvalues and their ratio. Limit time step for the FE scheme. Critical stabilization parameter and critical time step for the EFT12 scheme. Contours of the reference solution obtained with the backward Euler scheme using a conjugated gradient solver and  $\Delta t =0.0010$s.} \label{cubic-domain-spherical-refined-zone}
\end{figure}


Figure \ref{cubic-domain-spherical-refined-zone-unstructured} shows the temperature distribution along lines AB and CD at different times obtained with the EFT12 scheme and different variable time steps, following the patterns shown in the figure.  The speed-ups versus the FE scheme range from the maximum value of ${\mathrm{SUP}}^{\mathrm{FE}}=  701.5$ to  ${\mathrm{SUP}}^{\mathrm{FE}}=  350.5$. The variable time step yields a correct distribution of the temperature in all cases, with an increased accuracy for the smaller values of the final (still very large) time step chosen.

\begin{figure}[h!]
  \centering
  \includegraphics[width=15cm]{Figure22}
  \caption{Cubic domain with internal spherical cavity and refined internal zone. Temperature distribution along lines AB and CD at different times obtained with the EFT12 scheme and different variable time steps giving speed-ups versus the FE scheme ranging from   ${\mathrm{SUP}}^{\mathrm{FE}}=  701.5$ to  ${\mathrm{SUP}}^{\mathrm{FE}}=  350.5$. The variable time step yields a correct distribution of the temperature in all cases, with an increased accuracy for the smaller values of the speed-up.}\label{cubic-domain-spherical-refined-zone-unstructured}
\end{figure}

We highlight that two order of magnitude gains in speed-up versus the FE scheme throughout the transient solution are obtained in all cases, despite the huge time step used during most of the analysis.

\begin{description}
  \item[{\bf Remark 3}.] In the examples presented in the paper a constant time step equal or slightly smaller than the limit value for the FE scheme has been used during the first instants of the solution (5\% of the total analysis time). The time step is then increased linearly, up to an instant equal to 20\% of the total analysis time, to reach a value close or slightly smaller than the critical time step for the EFT scheme. Following that instant, the time step has been kept constant during the rest of the solution. We recommend this practical rule for the numerical solution of similar parabolic problems using the EFT12 scheme with large time steps.
\end{description}

\section{CONCLUDING REMARKS}

The explicit EFT12 time integration scheme with a variable time step allows to obtain accurate transient solutions for the parabolic equation using irregular grids, such as unstructured finite element meshes, or grids containing elements or cells of very different sizes,  using time steps considerable larger than for the  forward Euler (FE) scheme. This overcomes the limitation of explicit time integration schemes in these situations, in terms of the need of using a very small time steps due to a large eigenvalue ratio. The EFT12 scheme with  a variable time step also yields accurate transient solutions near Dirichlet boundaries using much larger time steps than the FE scheme during most of the transient solution.

The examples presented have demonstrated the accuracy of the explicit EFT12 scheme with large time steps for analysis of one, two and three dimensional transient heat conduction problems with Dirichlet boundaries using irregular FEM meshes, containing elements of very different size. The time step sizes allowed by the EFT12 scheme in the problems presented in the paper are between one and two orders of magnitude larger that those required for the FE scheme.

\section*{ACKNOWLEDGEMENTS}

The authors thank Rainald L\"{o}hner for many useful discussions. The authors acknowledge the support from the MCIN/AEI of the Government of Spain and the FEDER programme of the EU for funding this work through project PID2021\_122676NB-100.


\renewcommand{\refname}{REFERENCES}

\begin{thebibliography}{99}

%\bibitem{Alexiades96}  Alexiades, V., Amiez, G., Gremaud, P.A. , Super-time-stepping acceleration of explicit schemes for parabolic problems. Comm. Num. Meth. Eng., 12(1), 31-42, 1996. 10.1002/(SICI)1099-0887(199601)12:1<31::AID-CNM950>3.0.CO;2-5

\bibitem{Belytschko83} Beltyschko, T.  and Hughes T.J.R. (Eds.), Computational methods for transient analysis. Vol. 1 of the Series on Computational Methods in Mechanics, Elsevier, 1983.

\bibitem{Belytschko14} Beltyschko, T., Liu, W.K., Moran, B.  and Elkhodary, K.I., Non linear finite elements for continua and structures, 2d edition, John Wiley \& Sons, 2014.

%\bibitem{Brew71} Brew J. S. and Brotton D. M., Non-linear structural analysis by dynamic relaxation. Int J Numer Meth Engng, 3, 436–483, 1971. 10.1002/nme.1620030403

\bibitem{Butcher03} Butcher, J.C., Numerical methods for ordinary differential equations. John Wiley, 2003.

\bibitem{Cotela17} Cotela, J., Rossi, R., O\~{n}ate, E.,  A FIC-based stabilized finite element formulation for turbulent flows. Comput. Methods Appl. Mech. Engrg.,  315, 607-631, 2017. 10.1016/j.cma.2016.11.020

\bibitem{Felippa07} Felippa C.A., O\~{n}ate E., Nodally exact Ritz discretizations of the 1D diffusion-absorption and Helmholtz equations by variational FIC and modified equation methods. Comput Mech,  39(2), 91-111, 2007. 10.1007/s00466-005-0011-z

\bibitem{Felippa18} Felippa, C.A., O\~{n}ate, E., Idelsohn, S.R.,  Variational framework for FIC formulations in continuum mechanics: High order tensor-derivative transformations and invariants. Arch Comput Methods Eng,  25 (4), 919–963, 2018. 10.1007/s11831-017-9245-0

%\bibitem{Gentzsch78} Gentzsch W., Schl\"{u}ter, \"{U}ber ein Einschrittverfahren mit zyklischer Schrittweiten\"{a}nderung zur L\"{o}sung parabolischer Differentialgleichungen. ZAMM 58, T415-T416, 1978.
%
%\bibitem{Gentzsch80} Gentzsch W.,   \"{U}ber ein verbessertes explizitesEinschrittverfahren zur L\"{o}sung parabolischer Differentialgleichungen. DFVLR Report, 1980.

\bibitem{Hughes87} Hughes, T. J. R., The Finite Element Method – Linear Static and Dynamic Finite Element Analysis. Prentice-Hall, Englewood Cliffs, N.J., 1987.

%\bibitem{Linear multistep method} Linear multistep methods.  https://en.wikipedia.org/wiki/Linear$\_$multistep$\_$method

%\bibitem{Lohner87} Lohner, R., K. Morgan,  An unstructured multigrid  method for elliptic problems. Int. J. Num. Meth Engn., 24, 101-115, 1987.

\bibitem{Maso22} Mas\'{o} M., de-Pouplana I., O\~{n}ate E.,  A FIC-FEM procedure for the shallow water equations over partially wet domains. Comput. Methods Appl. Mech. Engrg.,   389, 114362, 2022. https://doi.org/10.1016/j.cma.2021.114362

%\bibitem{Oakley95} Oakley, D. R., Knight, Jr., N. F., Adaptive dynamic relaxation algorithm for non-linear hyperelastic structures. Part I. Formulation. Comput. Methods Appl. Mech. Engrg., 126, 67-89, 1995. 10.1016/0045-7825(95)00805-

 \bibitem{On98} O\~{n}ate E., Derivation of stabilized equations for numerical solution of advective-diffusive transport and fluid flow problems. {Comput. Methods Appl. Mech. Engrg.},  151, 233-265, 1998. 10.1016/S0045-7825(97)00119-9

\bibitem{OnMa99}O\~{n}ate E., Manzan M., A general procedure for deriving stabilized space-time finite element methods for advective-diffusive problems. {Int J Numer Methods Fluids},  31, 203--221, 1999. 10.1002/(SICI)1097-0363(19990915)31:1<203::AID-FLD964>3.0.CO;2-Z

\bibitem{OnGa01} O\~{n}ate E., Garc\'{\i}a-Espinosa J., A finite element method for  fluid-structure interaction with surface waves using a finite calculus formulation. Comput. Methods Appl. Mech. Engrg.,  191, 635--660, 2001. 10.1016/S0045-7825(01)00306-1

\bibitem{OnPeMi01} O\~{n}ate E., Perazzo F., Miquel J., A finite point method for elasticity problems. Comput Struct,  79, 2151--2163, 2001. 10.1016/S0045-7949(01)00067-0

\bibitem{Onetal03}  O\~{n}ate E.,  Taylor R.L., Zienkiewicz O.C., Rojek J., A residual correction method based on finite calculus. Engineering Computations,  20(5/6), 629--658, 2003. 10.1108/02644400310488790

\bibitem{Onate04} O\~{n}ate E., Possibilities of finite calculus in computational mechanics.  Int J Numer Methods Eng.,   60(1), 255--281, 2004. 10.1002/nme.961

\bibitem{Onate04b} O\~{n}ate E., Rojek J., Taylor R.L. and Zienkiewicz O.C.,  Finite calculus formulation for incompressible solids using linear triangles and tetrahedra. Int J Numer Methods Eng.,  59 (11), 1473-1500, 2004. 10.1002/nme.922

\bibitem{Boroomand05}  Boroomand B., Tabatabaei A. A. and O\~{n}ate E.,  Simple modifications for stabilization of the finite point method, Int. Journal for Numerical Methods in Engineering,  Vol. 63 (3), pp. 351-379, 2005, 10.1002/nme.1278

\bibitem{OnMiHa06} O\~{n}ate E., Miquel J., Hauke G., Stabilized formulation for the advection-diffusion-absorption equation using finite calculus and linear finite elements. Comput. Methods Appl. Mech. Engrg.,  195(33--36), 3926--3946, 2006. 10.1016/j.cma.2005.07.020

\bibitem{OnVaGa07} O\~{n}ate E., Valls A., Garc\'{\i}a-Espinosa J., Modeling incompressible flows at low and high Reynolds numbers via a finite calculus-finite element approach. J Comput Phys,  {224}, 332--351, 2007. 10.1016/j.jcp.2007.02.026

\bibitem{OnMiNa16} O\~{n}ate E., Miquel J., Nadukandi P.,  An accurate FIC-FEM formulation for the 1D advection–diffusion–reaction equation. Comput. Methods Appl. Mech. Engrg., 298, 373–406, 2016. 10.1016/j.cma.2015.09.022

\bibitem{OnNaMi17}    O\~{n}ate E., Nadukandi P., Miquel J.,  Accurate FIC-FEM formulation for the multidimensional steady-state advection–diffusion–absorption equation, Comput. Meth. Appl. Mech. Eng.,  327, 352-368, 2017. 10.1016/j.cma.2017.08.012

\bibitem{On22} O\~{n}ate E., De-Pouplana I., Z\'{a}rate F., Explicit time integration scheme with large time steps for first order transient problems using finite increment calculus, Computer Methods in Applied Mechanics and Engineering,  Vol. 402, 115332, December 2022, 10.1016/j.cma.2022.115332

\bibitem{On24} O\~{n}ate E., Z\'{a}rate F., Gimenez J.M., Lohner R., Idelsohn S.R., Fast explicit time integration schemes for parabolic problems in mechanics, Int J Numer Methods Eng., 125(13), e7477, July 2024, 10.1002/nme.7477

%\bibitem{Papadrakakis81} Papadrakakis, M., A method for the automatic evaluation of the dynamic relaxation parameters. Comput. Methods Appl. Mech. Engrg., 25, 35-48, 1981. 10.1016/0045-7825(81)90066-9
%
%\bibitem{Parabolic problems}Parabolic problems. https://en.wikipedia.org/wiki/Parabolic$\_$partial$\_$differential$\_$equation

\bibitem{Puigferrat20} Puigferrat A., de-Pouplana I., O\~{n}ate E.,  FIC-FEM formulation for the multidimensional transient advection-diffusion-absorption equation. Comput. Methods Appl. Mech. Engrg.,  365, 112984, 1--22, 2020. 10.1016/j.cma.2020.112984

%\bibitem{Runge-Kutta} Runge–Kutta methods. https://en.wikipedia.org/wiki/Runge\%E2\%80\%93Kutta$\_$methods

\bibitem{Routh-Hurwitz} Routh–Hurwitz stability criterion. https://en.wikipedia.org/wiki/Routh\%E2\%80\%93Hurwitz$\_$stability$\_$criterion

%\bibitem{Sauve95}  Sauve, R. G.,  Metzger D. R., Advances in dynamic relaxation techniques for nonlinear finite element analysis. Trans. ASME J, 117(2), 170–176, 1995. 10.1115/1.2842106

\bibitem{Underwood83} Underwood P., Dynamic relaxation. In Computational Methods for Transient Analysis, T. Belytschko and T.J.R. Hughes (Eds.) 1, 245-265, Elsevier, 1983.

%\bibitem{Zienk85} Zienkiewicz O.C., L\"{o}hner R., Accelerated ‘relaxation’ or direct solution? Future prospects for FEM. Int J Numer Methods Eng.,
%    21(1), 1--11, 1985. 10.1002/nme.1620210103

\bibitem{ZiTa05} Zienkiewicz O.C.,  Taylor R.L.,  \emph{The finite element method. Volume 2: Solid and Structural Mechanics}.  6th Ed., Elsevier, 2005.

\bibitem{ZiTaZhu05} Zienkiewicz O.C.,  Taylor R.L., Zhu J.Z. \emph{The finite element method. Volume 1: The basis}.  6th Ed., Elsevier, 2005.

\bibitem{ZiTaNi05} Zienkiewicz O.C.,  Taylor R.L., Nithiarasu P. \emph{The finite element method. Volume 3: Fluid Dynamics}.  6th Ed., Elsevier, 2005.

\end{thebibliography}


\end{document}

