report.tex

\documentclass[twoside,a4paper,12pt]{article}

\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{lmodern}
\usepackage[left=3cm,right=3cm,top=2cm,bottom=2cm,asymmetric]{geometry}

\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage{booktabs}

% fancy notations https://en.wikibooks.org/wiki/LaTeX/Mathematics
\usepackage{mathrsfs}
\usepackage{amssymb}            % for mathfrak -- canonical font

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% hyperlinks in the document
\usepackage{hyperref}
% https://tex.stackexchange.com/questions/823/remove-ugly-borders-around-clickable-cross-references-and-hyperlinks
\hypersetup{
    colorlinks,
    linkcolor={red!80!black},
    citecolor={blue!50!black},
    urlcolor={blue!80!black}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% https://tex.stackexchange.com/questions/135649/make-citemy-reference-show-name-and-year
\usepackage{natbib}             % citations including name and year
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% code listings
\usepackage{listings}

\definecolor{codegreen}{rgb}{0,0.6,0}
\definecolor{codegray}{rgb}{0.3,0.3,0.3}
\definecolor{codepurple}{rgb}{0.58,0,0.82}
\definecolor{backcolour}{rgb}{0.8,0.8,0.8}
\definecolor{BLUE}{rgb}{1.,1.,0.}
 
\lstdefinestyle{mystyle}{
    backgroundcolor=\color{black!10},   
    % commentstyle=\color{black!80},
    commentstyle=\color{codegreen},
    keywordstyle={\bf\ttfamily\color{black}},
    stringstyle={\it\ttfamily\color{black!70}},
    numberstyle=\scriptsize\color{black!70},
    % basicstyle=\footnotesize,
    % breakatwhitespace=false,         
    % breaklines=true,                 
    captionpos=t,                    
    basicstyle=\scriptsize\ttfamily\color{black!90},
    % keepspaces=true,                 
    numbers=left,                    
    numbersep=5pt,   
    stepnumber=5,
    firstnumber=1,
    numberfirstline=false,
    belowcaptionskip=-1em,
    % showspaces=false,                
    % showstringspaces=false,
    % showtabs=false,                  
    % tabsize=2
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% escape underscore environment
% https://tex.stackexchange.com/questions/20890/define-an-escape-underscore-environment
\usepackage{url}
\DeclareUrlCommand\code{\urlstyle{tt}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% strike out 
% https://tex.stackexchange.com/questions/23711/strikethrough-text
\usepackage{soul}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\figref}[1]{Figure~\ref{#1}}
\newcommand{\tabref}[1]{Table~\ref{#1}}

\newcommand{\prog}[1]{\textsf{#1}}

\newcommand{\python}{\prog{Python}}
\newcommand{\numpy}{\prog{NumPy}}
\newcommand{\lapack}{\prog{LAPACK}}

\newcommand{\ie}{{\it i.e.\ }}

% \newcommand{\note}[1]{{\color{red}(#1)}}
% \newcommand{\QM}{{\color{red}(?)}}

\newcommand{\note}[1]{}
\newcommand{\QM}{}
\newcommand{\X}[1]{{\color{brown}\st#1} } % removed text
\newcommand{\A}[1]{{\color{purple}#1}} % added text

\input{nots.tex}

\begin{document}

%% -----------------------------------------------------------------
%% LIST OF CONTENTS
%% -----------------------------------------------------------------

\tableofcontents

%% -----------------------------------------------------------------
%%  LIST OF TABLES, ILLUSTRATIONS, ETC. (IF ANY)
%% -----------------------------------------------------------------

\listoffigures
\listoftables

%% -----------------------------------------------------------------
%% LIST OF ACCOMPANYING MATERIAL (IF ANY)
%% -----------------------------------------------------------------

%% -----------------------------------------------------------------
%% DEFINITIONS; ABBREVIATIONS
%% -----------------------------------------------------------------

%% -----------------------------------------------------------------
%% TEXT (DIVIDED INTO CHAPTERS, SECTIONS, ETC.)
%% -----------------------------------------------------------------
% \mainmatter % start numbering from 1 here -- ommitted due to regulation

\section{Statement of the system}

This program aims to use finite elements methods to solve the boundary
value problem
%
\begin{align}
  \young(x) \area(x) u_{,xx} + f(x) = 0, \label{eq:bvp}
  &&u(L) = \bcg, u_{,x}(0) = -\bch,
\end{align}
where $L$ is the length, $\young(x)$ is Young's modulus and $\area$ is
the cross section area, $u(x)$ is a deformation of the of the object
in consideration, and shorthand notation for a derivative is written
after a coma in a subscript, e.g. $u_{,x} = \partial u /\partial x$.  The
function $f(x)$ represents the body forces. The $\bcg$ and $\bch$ are
essential and natural boundary conditions respectively.

Let $u(x)$ be a solution of system \eqref{eq:bvp} then we may write
%
\begin{align}
  0 = \int_0^L (\young(x) \area(x) u_{,xx} + f(x)) w(x) \dd x 
\end{align}
and $w(x)$ is a chosen \emph{weighting (test) function} that satisfies a
criteria $w(L)=0$.
%
Using integration by parts we get
%
\begin{align}
  \int_0^L \young(x) \area(x) u_{,xx} w(x) \dd x 
  -\int_0^L f(x) w(x) \dd x  =
   \int_0^L u_{,x} (\young(x) \area(x) w(x))_{,x} \dd x -\nonumber \\ -
  (u_{,x} (\young(x) \area(x)  w(x))\big|_0^L -
  \int_0^L f(x) w(x) \dd x =  
  %
  \int_0^L u_{,x} (\young(x) \area(x) w(x))_{,x} \dd x -\nonumber \\ -
  (u_{,x}(L) (\young(L) \area(L)  w(L) -
  u_{,x}(0) \young(0) \area(0)  w(0)) 
  -\int_0^L f(x) w(x) \dd x = 0
\end{align}
%
and using the boundary values and the property of the weighting
function, the previous simplifies as
%
\begin{align}
  \int_0^L u_{,x} (\young(x) \area(x) w(x))_{,x} \dd x =  
  \int_0^L f(x) w(x) \dd x + \bch \young(0) \area(0)  w(0)
  \label{eq:weak-exact}
\end{align}
%
this method is called a \emph{weak (variational)}.

The solution of this system can be found by strong (classical) method
or weak (variational) method. The equivalence with the \emph{strong
  (classical)} method goes under a name \emph{fundamental lemma} which we will not show in here.

To solve the system numerically we will use Galerkin approximation
method. This method uses \eqref{eq:weak-exact} and provides an
approximation of the boundary value problem \eqref{eq:bvp} as
%
\begin{align}
a(\young\area\wwh,\uuh) = (\wwh,f) + \bch \young(0)\area(0)\wwh(0) \label{eq:weak-approx}
\end{align}
%
where we define the operators $a(y,z) = \int_0^L y_{,x} z_{,x} \dd x$
and $(y,z) = \int_0^L y z \dd x$ for arbitrary functions $y = y(x)$
and $z= z(x)$.

We construct the function
\begin{align}
  \uuh  = \vvh + \ggh \label{eq:galerkin-construc}
\end{align}
%
where $\ggh$ is given by satisfying the essential boundary condition
$\ggh(L) = g$ and $\vvh$ corresponds to the displacement. We
substitute \eqref{eq:galerkin-construc} into \eqref{eq:weak-approx} to
obtain
%
\begin{align}
  a(\young\area\wwh,\vvh) = (\wwh, f) + \young(0)\area(0)\wwh(0)\bch - 
  a(\young\area\wwh,\ggh)\label{eq:galerkin}
\end{align}
%
the weighting function can be defined as
%
\begin{align}
  \wwh = \sum_{i=0}^{\nel-1} c_i \shape_i \label{eq:wwh-sum}
\end{align}
%
where $\shape_i$ is so-called a \emph{basis (shape) function} for
element $i$. Similarly, we define
%
\begin{align}
  \vvh = \sum_{i=0}^{\nel-1} d_i \shape_i \label{eq:vvh-sum}
\end{align}
where $c_i, d_i$ are constants.

After substitution of \eqref{eq:wwh-sum} and \eqref{eq:vvh-sum} into
the \eqref{eq:galerkin} and using the bilinearity property of
operators $a(\cdot,\cdot)$ and $(\cdot,\cdot)$ we get an expression
%
\begin{align}
  \sum_{j=0}^{\nel-1} a(\young\area\shape_i, \shape_j) d_j = (\shape_i, f) + \young(0)\area(0)\shape_i(0) \bch - a(\young\area\shape_i,\shape_\nel)\bcg. \label{eq:galerkin-row}
\end{align}
%
Using a notation
%
\begin{subequations}
  \begin{align}
    \stme{ij} =& a(\young\area\shape_i,\shape_j) = \int_0^L \young(x)\area(x)\shape_{i,x}(x),\shape_{j,x}(x) \dd x \label{eq:stme} \\ 
    F_i =& (\shape_i,f) + \young(0)\area(0)\shape_i(0) \bch - a(\young\area\shape_i,\shape_{\nel}) \bcg .
  \end{align}\label{eq:fem-stme-fi}
\end{subequations}
%
Then equation \eqref{eq:galerkin-row} can be rewritten as
%
\begin{align}
  \sum_{j=0}^{\nel-1} \stme{ij} d_j = F_i
\end{align}
or more concisely as
\begin{align}
\stm \vdspl = \vfrc \label{eq:concise-eq}
\end{align}
where $\stm$ is a \emph{stiffness matrix}, $\vdspl$ is a
\emph{displacement vector} and $\vfrc$ is a \emph{force vector}.

\section{First problem}

Consider boundary value problem \eqref{eq:bvp} where $\bcg=\bch =0$ with a
constant force as $f(x) = q$, where $q=-1$, length of the piece is
$L=1$, cross-section area and young modulus are $A(x) = E(x) = 1$.

\subsection{Shape function}

We define the basis function and find its derivative for each element
$i$ of $\nel$ elements. In the middle elements we get
\begin{align}
\shape_i(x) =& \frac{x-x_{i-1}}{h_{i-1}}, &\shape_{i,x}(x) =& \frac{1}{h_{i-1}},&&\text{for } x_{i-1} \leq x \leq x_i \\
\shape_i(x) =& \frac{x_{i+1}-x}{h_{i}},  &\shape_{i,x}(x) =& -\frac{1}{h_{i}}, &&\text{for }x_i \leq x \leq x_{i+1} \\
\shape_i(0) =& 0,                       &\shape_{i,x}(0) =& 0,                &&\text{elsewhere.}
\end{align}
%
whereas for the boundary node we have
%
\begin{align}
\shape_0(x) =& \frac{x_1-x}{h_{0}},             &\shape_0(x) =& -\frac{1}{h_{0}},                          &&\text{for } x_0 \leq x \leq x_{1} \\
\shape_\nel(x) =& \frac{x-x_{\nel-1}}{h_{\nel-1}},&\shape_\nel(x) =& \frac{1}{h_{\nel-1}},             &&\text{for } x_{\nel-1} \leq x \leq x_\nel .
\end{align}

\subsection{\python\ Implementation}

Then we compute the values for the stiffness matrix using numerical
integration. In \python\ this can be done using \code{scipy.integrate}
function with a \code{integrate.quad} for a quadrature rule. The
output array contains a tupple containing the value of integral and
the maximum error of the result.

At first, we have integrated from the limits of $0$ all the way
thorough the length of the element $L$. This caused a problem when the
number of elements were greater than $16$. This was due to the fact
that the shape function is defined piecewise, and the results of
numerical integration were not sufficiently accurate. The way to
address this issue was to integrate only in the immediate surroundings
of each point, where the shape functions are actually
non-zero. However, this raised another issue, because the shape
functions were defined even out of the interval $[0,1]$. To avoid
integration in this area, we have hard-coded constrains directly to
the shape function. This could be obviously corrected, but because
the stiffness matrix will be computed element-wise in future version
of this code, it was neglected at this stage.

After constructing element stiffness matrix and the force vector, we
have proceed to the solution of the system \eqref{eq:concise-eq}. In
the first attempt, we have found the displacement vector by inverting
the stiffness matrix using \code{np.linalg.inv} function (we imported
\numpy\ library as \code{np}). So in the current version, we
solve the system using \code{np.linalg.solve} function, which uses
\lapack\ routine \code{_gesv}. A performance benchmark has shown
that for a $5000$ elements the solution of the system using the
\code{np.linalg.solve} is about four times than using matrix inversion
(\ie $2.5$ seconds compared to $8.6$ seconds on Intel Core i5-3470 CPU
at a clock frequency of 3.20GHz).
%
\note{it would be interesting to see what order of $N$ is the
  computational cost. It should be $\cO{N^2}$ for gaussian elimination
  and $\cO{N^3}$ for matrix inversion. }

The code for the computation is shown as follows.

\lstinputlisting[language=python, style=mystyle, breaklines=true]{galerkin.py}


\subsection{Results}

The results are depicted in \figref{fig:displacement}. The simulation
was done for eight free elements. The material is fixed at the point
$x=1$ at a position $0$. The solution is exact at the nods.

\begin{figure}
  \centering
  \includegraphics[width=1.0\textwidth]{displacement.eps}
  \caption{Dependence of the displacement on the position (x-axis) for
    a computation of eight degrees of freedom. Numerical results at
    the nodes (red dots) are linearly interpolated (red lines) and
    compared to the exact solution (blue line).}
  \label{fig:displacement}
\end{figure}

The output of the program follows.
%
The first part corresponds to the computational cost (in seconds)
spend in the construction of force vector, stiffness matrix and on the
solution of the system. The values are very small for this

Further, the stiffness matrix and the force vector are printed and
finally the solution for the displacement is shown both as a numerical
result and exact solution obtained analytically. As the numerical
results are exact in this case, also the norm of the difference
between the numerical and analytical results is zero.
%
\note{the result at the nodes is 0, but what about the result on the
  midpoints? Hughes seems to talk about relative error in $u_{,x}$
  what is that useful for?  How is the order of convergence?}

\tabref{tab:comp-cost} shows the computational results for the
simulation using $5000$ elements. The remaining parameters are
identical to the code presented. As seen in this example, the
construction of the stiffness matrix is the most time consuming task
despite the fact, that we do not compute the zero elements of the
matrix. The relatively high cost is caused due to numerical
integration.

In the following section we will focus on the computation of the
stiffness matrix element-wise.

\begin{table}
  \centering
  \label{tab:comp-cost}
\caption{Computational cost of the solving of the system for $N=5000$.}
\begin{tabular}{ll}
  \toprule
  Task & cost (seconds)\\
  \midrule
  Constructing stiffness matrix: & 7.161638498306274   \\
  Constructing force vector:     & 0.33005428314208984 \\
  Solving the system:            & 2.460254430770874   \\
  \bottomrule
\end{tabular}
\end{table}

\lstinputlisting[breaklines=true]{galerkin.out}

\section{Element-wise Construction of Stiffness Matrix and Force Vector}

The time of construction of the stiffness matrix for this case has
reduced from $7.16$ seconds to $0.05$ seconds in a computation of
$N= 5000$ nodes (as shown in the \tabref{tab:comp-cost}). This is
probably because now the stiffness is given explicitely, rather then
integrated from the shape functions.

\figref{fig:element} shows the results of the elementwise construction
of the stiffness matrix. This now also allows setting elements of
different length.

\begin{figure}
  \centering
  \includegraphics[width=1.0\textwidth]{element.eps}
  \caption{Element wise construction of the stiffness matrix.}
  \label{fig:element}
\end{figure}

\lstinputlisting[language=python, style=mystyle, breaklines=true]{element.py}

-\note{it remains to explain the element-wise construction of the stiffness matrix}

\section{Combining Solid and Flow Mechanics}

This section aims to develop a model combining the mechanics for solid
and flow. This aims to simulate the formation of the composite
structure of carbon fibre and resin. During the formation the material
is placed into a pressure chamber, where a pressure is applied to the
material in an elevate temperatures. Due to the high temperature, the
resin becomes liquid and the pressure causes the flow of the resin out
of the material such that the carbon fibres become compressed
together. Consequently the resin undergoes a process of cure, that
causes a solidification of the resin so that the material is ``glued''
together. The following text is motivated by an article by
\citet{Hubert1999}.

Here we consider a toy model for a formation of the material as
described.  In this situation the total stress is composed according
to the Terzaghi principle, which states that the porous material
subjected to pressure is opposed by the fluid pressure in the
pores. Mathematicaly, this is described as
%
\begin{align}
\stress(x,t) = \efstress(x,t) - \press(x,t) \label{eq:flow-stress}
\end{align}
%
where $\press(x,t)$ is the resin pressure and the effective stress
$\efstress(x,t)$ that are dependent both on the reference position $x$ and
the time $t$. The effective stress is connected with the
displacement by a constitutive law
\begin{align}
\efstress = u_{,x},
\end{align}
%
where the parameters describing the properties of the material are
neglected and $u= u(x,t)$ represents the displacement.


From the Newton third law (sum of forces equals to zero) we get
%
\begin{align}
  0 =
  \stress_{,x} + f(x) =
  \efstress_{,x}(x,t) - \press(x,t)_{,x} + f(x) =
  % (u_{,x}(x,t) - \press(x,t))_{,x} + f(x) =\nonumber \\
  u_{,xx}(x,t) - \press_{,x}(x,t) + f(x) 
  % , 
  % &&u(0) = \bcg, u_{,x}(L) = -\bch\label{eq:bvp-t}.
\end{align}
%
where the body forces and the inertia of the material are neglected
$f(x) = 0$, as the gravitaional force is small compared to the applied
pressure, and the deformation, and hence the speeds are small.
%
So we obtain a boundary value problem
%
\begin{align}
  u_{,xx}(x,t) - \press_{,x}(x,t) + f(x) \label{eq:bvp-flow}
  = 0 , &&& u(0,t) = 0,\\\nonumber
        &&&u(L,t) =
            \begin{cases}
              -Rt &\text{ for } t < \duration,\\
              -RT &\text{ for } t > \duration.
            \end{cases} 
\end{align}
%
The boundary at the reference point $L$ is compressed by a press with
a constant speed $R$ for in the first $\duration$ time units. Then the press
is stoped at the possition reached after $\duration$.

%

The flow of the resin through the fibre bed will follow Darcy's law.
This law describes the flow of a fluid through a porous material, and
is analogous to Ohm's law from the theory of electric circuits. The
Darcy's law reads as
%
\begin{align}
  \flow(x,t) = -\frac{\permeab(x)}{\viscos(\temp,\cure)} \press_{,x}. \label{eq:darcy}
\end{align}
where the $\flow(x,t)$ is the velocity of the flow, $\permeab$ is a fibre bed
permeability (due to the size of the pores),
%
%
$\viscos$ is a resin viscosity (dependent on temperature $\temp$ and
degree of cure of the resin $\cure$).  For simplicity, let's consider
a case where $\permeab/\viscos=1$.
%
The velocity of the flow constrains the position of the fibre bed,
hence by linking the definition of velocity $\flow = - u_{,t}$ (recall the
notation $u_{,t} = \partial u/\partial t$) with the Darcy's law, we
can write an inital value problem for the internal points of the
material as
%
\begin{align}
  u_{,t}(x,t) =  \press_{,x} , && u(x,0) = 0 , \text{ for } 0 < x < L. \label{eq:velocity}
\end{align}


\subsection{Finite Differences Algorithm}

We write the boundary \eqref{eq:bvp-flow} and initial
\eqref{eq:velocity} value problem in a simplified form, where the
physical units are ignored, as follows
%
\begin{align}
  \label{eq:heat}
  u_{,xx}(x,t) - D u_{,t}(x,t) + f(x) 
  =& 0 ,        &&u(x,0) = 0, u(L,t) =
            \begin{cases}
              -Rt &\text{ for } t < \duration,\\
              -RT &\text{ for } t > \duration,
            \end{cases}\\\nonumber
 &&& u(0,t) = 0.
\end{align}
%
where in this subsection we consider a simplified problem where $D=1$.

The simplest method for solving this system is a finite difference
method. This uses time and space discretisation as
%
\begin{subequations}
  \begin{align}
    x_i = x_0 + i\dx, \text{ for } i = 1, 2, \dots, N_x \\
    t_j = t_0 + j\dt, \text{ for } j = 1, 2, \dots, N_t
  \end{align} \label{eq:discr}
\end{subequations}
%
where $\dx$ is the space step and $\dt$ is the time step, and $x_{N_x} = L$.  So we can
write the displacement at a certain points in time-space as
$\disp{i}{j} \approx u(x_i,t_j)$. Then we find the discretisation as
%
\begin{subequations}
  \begin{align}
    u_{i,t}^j  =& \frac{\disp{i}{j+1} - \disp{i}{j}}{\dt} + \bigO{\dt}, \\
    u_{i,xx}^j =& \frac{\disp{i+1}{j} - 2\disp{i}{j} + \disp{i-1}{j}}{\dx^2} + \bigO{\dx^2},
  \end{align}\label{eq:fin-diff}
\end{subequations}
%
where $u_{i,t}^j = u_{,t}(x_i,t_j)$, $u_{i,xx}^j = u_{,xx}(x_i,t_j)$,
and $\mathcal{O}$ is so-called big-O notation that specifies the order
of the approximation and corresponds to the part of the solution that
is ignored. This equation is valid for the interior point $x_i$ where
$i = 2,\dots, N_x-1$. The finite difference scheme requires to
introduce so-called numerical boundary conditions which represent the
situation at the boundaries. The numerical boundary conditions
interpolate the results at the interior points to approximate the
boundary points. The simplest method is to use linear interpolation to
find the solution at hypothetical points outside the boundary as
%
\begin{subequations}
\begin{align}
  \disp{-1}{j} =& \disp{1}{j} - 2\disp{0}{j}, \\
  \disp{N_x+1}{j} =& 2\disp{N_x}{j} - \disp{N_x-1}{j},
\end{align}
\end{subequations}
%
so that the finite difference gives
%
\begin{align}
  u_{0,xx}^j =& \frac{\disp{1}{j} - 2\disp{0}{j} + \disp{-1}{j}}{\dx^2} + \bigO{\dx^2} =
  \frac{\disp{1}{j} - 2\disp{0}{j} + \disp{1}{j} - 2\disp{0}{j}}{\dx^2} + \bigO{\dx^2} =\nonumber \\
  =&\frac{2(\disp{1}{j} - 2\disp{0}{j})}{\dx^2} + \bigO{\dx^2},
  \\
  u_{N_x,xx}^j =& \frac{\disp{N_x+1}{j} - 2\disp{N_x}{j} + \disp{N_x-1}{j}}{\dx^2} + \bigO{\dx^2} =\nonumber \\
  =&\frac{2\disp{N_x}{j} - \disp{N_x-1}{j} - 2\disp{N_x}{j} + \disp{N_x-1}{j}}{\dx^2} + \bigO{\dx^2} = \bigO{\dx^2}.
     \label{eq:error-bigO}
\end{align}

Substituting \eqref{eq:fin-diff} into \eqref{eq:heat} we get the
numerical algorithm as
%
\begin{align}
  % \frac{\disp{i+1}{j} - 2\disp{i}{j} + \disp{i-1}{j}}{\dx^2} -
  % \frac{\disp{i}{j+1} - \disp{i}{j}}{\dt}  
  % + f(x_i) 
  % = 0 \\
\disp{i}{j+1}  =   \disp{i}{j} + \dt\left(\frac{\disp{i+1}{j} - 2\disp{i}{j} + \disp{i-1}{j}}{\dx^2}   + f(x_i) \right), \label{eq:findiff-alg}
 % +\bigO{\dx^2, \dt}.
\end{align}
%
which has associated error of
\begin{align}
  \errU =\bigO{\dt} + \bigO{\dx^2}.
\end{align}

This numerical method is stable if
\begin{align}
  \frac{\max(D)\dt}{\dx^2} < \frac{1}{2},
\end{align}
which is described in detail in \cite{Strikwerda2004}.
%
\note{derivation from page 47}
%
Using this relation we can choose the space or time step as
%
\begin{align}
  % \frac{\max(D)\dt}{\dx^2} < \frac{1}{2}, \\
  \sqrt{2\max(D)\dt} < \dx, \\
  \dt < \frac{\dx^2}{2\max(D)}.\label{eq:stable-boundary}
\end{align}

\figref{fig:error-analysis}(a) shows the graph of the stability of the
solution dependent on the space step $\dx$. The stable region
corresponds to the values of $\dt$ under the line.

\subsection{Implementation of Simplified Flow Problem}

\begin{figure}
  \centering
  \includegraphics{heat.pdf}
  \caption{Finite difference solution of equation
    \eqref{eq:heat}. Panel (a) shows the dependence of the deformation
    of the fibre bed on time at the boundaries $x=0$ (violet line) and
    $x=1$ (yellow line) and in the interior points $x=0.2$ (green
    line) and $x=0.5$ (blue line) (b) shows the dependence of the
    displacement on the possition in the frame of reference $x$ at the
    initial conditions (violet line), at $t=0.1$ (green line), at
    $t=0.2$ (blue line) and at the end of the simulation $t=0.4$
    (yellow line).  Panel (c) shows colour coded map (left bar) of the
    deformation dependent on the time (vertical axis) and the frame of
    reference (horizontal axis). The a gray scale map from white
    colour which corresponds to small deformation to black which
    corresponds to maximal deformation of $-0.02$ is complemented by a
    contour lines ad four different values as specified in the legend.}
  \label{fig:findiff-heat}
\end{figure}

\figref{fig:findiff-heat} shows the solution of equation
\eqref{eq:heat} using finite difference method
\eqref{eq:findiff-alg}. The time step $\dt=0.001$ and space step
$\dx = 0.1$ such that the space contains $N_x=10$ nodes. The initial
conditions $\disp{i}{0} = 0$. The bottom boundary is fixed at
$\disp{0}{t} = 0$, while the top boundary move with a speed $R=1$ such
that $\disp{N_x}{j} = j \dt R$.


\begin{figure}
  \centering
  \includegraphics{error.pdf}
  \caption{Analysis of the error of the solution. Panel (a) shows the
    boundary between stable and unstable region in the $\dx$--$\dt$
    plane according to \eqref{eq:stable-boundary} where $\max(D) =
    1$. Panel (b) and (c) shows the norm of the error in the
    computation of the heat equation in (b) dependent on $\dt$ at
    $\dx=0.1$, in (c) as dependent on $\dx$ at $\dt=1e-7$.}
  \label{fig:error-analysis}
\end{figure}

\figref{fig:error-analysis} shows the analysis of the error in the
solution of the heat equation. Panel (a) shows the stability region
from equation \eqref{eq:stable-boundary} as specified above. The
panels (b) and (c) show the dependence of the error on the time step
and on the grid step. The norm was computed using the formula
%
\begin{align}
  \errU = \norm{\hat{u} - u} = \frac{1}{N_t N_x}\sqrt{ \sum_{i,j} \left(\hat{u}_{i}^{j} - \disp{i}{j}\right)^2}
\end{align}
where $N_t$ and $N_x$ is the number of time and grid points of in the
stored data files, and $\hat{u}$ corresponds to accurate solution. The
time step of the accurate solution was $\dt = 1\ee{-7}$. The space step
of accurate solution was $\dx=0.1$ and $\dx=1\ee{-3}$ in panel (b) and
(c) respectively.

The simulation time of the solution with $\dx=1\ee{-3}$, $\dt=1\ee{-7}$
was $1$ hour $43$ minutes.

The results confirm the theoretical results from equation
\eqref{eq:error-bigO} giving the first-order in $\dt$ and the
second-order in $\dx$ as the slope in the log-log plot is $1$ and $2$
for panel (b) and (c) respectively.

The implementation of the finite difference code in \python\ is shown
below.
%
\lstinputlisting[language=python, style=mystyle, breaklines=true]{findiff.py}

\subsection{Fibre Consolidation Model}

The goal of the previous text was to develop a simple toy example of a
type of equations describing the system. In the process we have ignored
or largely simplify the physical units to obtain a model with
mathematical characteristics rather then aiming for physical
precision. In the following text we focus to address the assumptions
in order to obtain more physically detailed model. The constants in the
models are shown in \tabref{tab:phys-const}.

The first assumption we made was that the stress is a linear function
of the strain $\varepsilon = u_{,x}$. This is only correct for very
small deformations when the volume fraction of the fibre does not
change much.  The stress is a non-linear variable dependent on the
volume fraction of the carbon fibre in the material. As the fibres are
packed the stiffness will increase non-linearly. The the following
formula shows a good agreement with experiments \citet{Gutowski1987},
\ie it is a heuristic model
% 
\begin{align}
  \efstress(x,t) = \spring
  \frac{\frac{\vfrac(x,t)}{\vzero} -1}{\left(\frac{1}{\vfrac(x,t)} -\frac{1}{\vmax}\right)^4}, \label{eq:stress-nonlin}
\end{align}
%
where $\vzero = \vfrac(x,t_0)$ \note{which is assumed to be constant
  throughout the space $x$}, and $A_s$ is so-called spring
constant\note{reference for the expression of spring constant? In the
  \citet{Gutowski1987} only $A_s$ is specified.},
%
\begin{align}
  A_s = \frac{3\pi E_f}{\beta^4},
\end{align}
%
fitted to experimental data and their values are specified in
\tabref{tab:phys-const}.
%
The fibre volume fraction $\vzero$, and $\vmax$ represents the ratio
between the total volume and the fibre volume at the initial, maximal
packing respectively. The fibre volume fraction is expressed as
\begin{align}
  \vfrac(x,t) = \frac{v_f(x,t)}{v_c(x,t)}
\end{align}
%
where $v_c$ is the volume of the fibre in corresponding volume of
composite $v_f$. As the fibre volume is assumed to remain constant,
which can be expressed as
%
\note{is this assmuption reasonable? What about the compression of the fibre?}
%
\begin{align}
  v_f(x,t) = \vzero v_{c0}. \label{eq:fibre-fraction}
\end{align}
where $\vzero = \vfrac(x,t_0)$ and $v_{c0} = v_c(x,t_0)$.
%
At a time $t$ volumetric strain in the element is defined as
\begin{align}
  \strain_V(x,t) = \frac{v_c(x,t) - v_{c0}}{v_{c0}}, \label{eq:volumetric-strain}
\end{align}
which gives
%
\begin{align}
{v_{c0}}= \frac{v_c(x,t) }{\strain_V(x,t) + 1}. \label{eq:initial-volume}
\end{align}

Substituting \eqref{eq:initial-volume} into \eqref{eq:fibre-fraction}
and then into \eqref{eq:volumetric-strain} yields
%
\begin{align}
  \vfrac(x,t) = \frac{\vzero}{1+\strain(x,t)}, \label{eq:fibre-volume}
\end{align}
%
where $\strain =  u_{,x}$.

\figref{fig:consolidation} shows the change in the properties of the material during the consolidation. All panels shown for the region of the independent variable for which was the model designed. The panel (a) shows the dependency of the volume fraction on the strain as described by the \eqref{eq:fibre-volume}. As the strain decreases, \ie\ fibre bed is compressed, the fibres consolidate which means that the volume fraction of the fibres increase. The panel (b) shows the effective stress in the fibre bed during the consolidation as described by \eqref{eq:stress-nonlin}. In this panel we can observe that as the fibre volume fraction increases, the fibre bed becomes stiffer, which can be observed as higher stress. The last panel (c) shows the combination of the two previous panels (a,b) and can be interpreted as the increase of stiffness with the consolidation.
%
\begin{figure}
  \centering
  \includegraphics{consolidation.pdf}
  \caption{Mechanical properties of the material during the consolidation of the carbon fibres. The combination of dependence of volume fraction on the strain as shown in panel (a), combined with the dependence of stress on volume fraction shown in panel (b), gives the dependence of the stress as a function of strain as shown in panel (c).}
  \label{fig:consolidation}
\end{figure}

In the Darcy's law in equation \eqref{eq:darcy} we assumed the ratio
between fibre bed permeability and the viscosity of the resin
$\permeab/\viscos=1$. However, the permeability will depend on the volume
fraction of the resin. \citet{Dave1987b} suggested a model as
%
\note{I didn't find the model in their paper, this expression is given in \citet{Hubert1999}}
%
\begin{align}
  \permeab = \frac{r_f^2}{4\kozeny} \frac{(1-\vfrac)^3}{\vfrac^2} \label{eq:permeability}
\end{align}
%
where $r_f$ specifies the fibre radius and $\kozeny$ is so called Kozeny
constant (see \tabref{tab:phys-const}). Also, the viscosity is a
function of other variables, such as temperature and the degree of
cure of the resin.

\figref{fig:permeab} shows the dependence of the permeability of the fibre bed during consolidation. Panel (a) shows the permeability as a function of  fibre volume fraction. When the fibre volume fraction is  high, the size of the pores between the material is higher, so the permeability is low. As the material is compressed and the fibre volume fraction reduces, the size of the pores reduces as well which results in the reduction of the permeability.
%
\begin{figure}
  \centering
  \includegraphics{permeab.pdf}
  \caption{Permeability of the fibre bed during the consolidation.
    The figure (a) The combination of dependence of volume fraction on
    the strain as shown in \figref{fig:consolidation}(a), combined
      with the dependence of the permeability on volume fraction shown
      in this figure panel (a), gives the dependence of the
      permeability as a function of strain as shown in panel (b).}
  \label{fig:permeab}
\end{figure}


%
Other assumption which we did the Darcy's law in equation
\eqref{eq:darcy} regards the viscosity of the resin.  The viscosity
which should be substituted to \eqref{eq:darcy} was suggested by  \citet{Kenny1992}
as
%
\begin{align}
  \viscos(\temp,\cure) =
  \viscosinf \exp{\frac{E_\viscos }{R \temp}} \left(\frac{\curegel}{\curegel - \cure}\right)^{A+B\cure} \label{eq:viscosity}
\end{align}
%
where the constants ($\viscosinf$, $\curegel$, $R$, $E_\viscos$, $A$ and
$B$) are explained in \tabref{tab:phys-const}, % $g(\vfrac)$ accounts for
% low viscosity of the air \note{how to call it? what it means?},
$\temp$ is the temperature in Kelvins, which is controlled in the
experiment, $\cure$ is the degree of cure that will be discussed
later and which takes values in the interval $[0,1]$.

\figref{fig:viscos} shows the dependence of the viscosity of the resin
on the temperature and the degree of cure of the resin. Panel (a)
shows the dependence on the temperature. As the temperature increases
the resin ``melts'' which means that it becomes more viscous. Panel
(b) shows the dependence of the viscosity on the degree of cure
$\cure$. Although, the $\cure$ is allowed to get the interval between
$[0,1]$ there is a singularity at a point $\cure = 0.48$, where the
resin becomes nearly solid.

\begin{figure}
  \centering
  \includegraphics{viscosity.pdf}
  \caption{Viscosity of the resin. Panel (a) shows the dependence of the
    viscosity on the temperature for three different values of
    $\cure$: magenta line $\cure=0.1$, green line $\cure=0.2$, cyan
    line $\cure=0.4$, panel (b) shows the dependence of the viscosity
    on the degree of cure for three different temperatures: magenta
    line $\temp = 100$\degree C, green line $\temp = 130$\degree C, and cyan line
    $\temp = 180$\degree C.}
  \label{fig:viscos}
\end{figure}

For ilustration the visosity we show some examples viscosity of
selected fluids in \tabref{tab:viscos-example} obtained from wikipedia
page on viscosity. The highest value of viscosity in the table is a
material called pitch, which is a common name for a viscoelastic
polymers, such as tar or asphalt. At those values of viscosity the
pitch seems to be solid, however in longer observation, the fluidity
can be observed as demonstrated in a pitch drop experiment developed
in the School of Mathematics and Physics at The University of
Queensland, Australia. In this longest running laboratory experiment
pitch was placed in a funnel and allowed to flow down since the year
1930. Due to high viscosity only nine drops have fallen down through
the funel up to a current date.

\begin{table}
  \centering
  \begin{tabular}{ll}
    \toprule
    fluid         & viscosity (Pa$\cdot$s) \\
    \midrule
    air           & $1.86\ee{-5}$          \\
    water         & $8.94\ee{-4}$          \\
    ethanol       & $1.074\ee{-3}$         \\ 
    olive oil     & $0.081$                \\
    honey         & $2$--$10$              \\
    ketchup       & $50$--$100$            \\
    peanut butter & $\sim 250$             \\
    pitch         & $2.3\ee{8}$            \\
    \bottomrule
  \end{tabular}
  \caption{Viscosity values of selected fluids at the room temperature $\sim 25$\degree C (source: wikipedia)}
  \label{tab:viscos-example}
\end{table}

Other interesting table from the wikipedia page on viscosity (not
shown) is on the temperature dependence of the viscosity of water. The
table shows that the viscosity reduced by about one order of magnitude
when heated from $10$ to $100$\degree C.  


The resin cure is described by a differential equation in a form
%
\begin{align}
  \cure_{,t} =  A_\cure \exp{-\frac{E_\cure }{R\temp}} \cure^m (1-\cure)^n \label{eq:cure}
\end{align}
%
where the constants ($A_\cure$, $E_\cure$, $m$, $n$ and $R$) are
given in the table, $\temp$ is the temperature as mentioned above.

This model has two equilibrium points $\cure_{,t}=0$ the $\cure = 0$ is an unstable equilibrium and $\cure=1$ is a stable equilibrium. \figref{fig:cure} shows rate of cure of the resin as a dependence on the temperature (panel a) and cure (panel b). The cure takes strictly positive values, as is also evident from expression \eqref{eq:cure}. 

The dependency of the rate of cure on temperature is almost exponential. From the room temperature, when almost no curing is happening it increases by about four orders of magnintude at the temperature of $200$\degree C.

The dependency of the rate of cure on the degree of cure has a maximum at $\sim 0.2$. At high values of cure, the rate decreases. However, higher values of cure do not need to be considered as from the analysis of the viscosity it is apparent, that the resin is gelated solid for values of cure $\cure > 0.47$.

\begin{figure}
  \centering
  \includegraphics{dcure.pdf}
  \caption{Rate of cure of the resin. Panel (a) shows the dependence of the
    rate of cure on the temperature for three different values of
    $\cure$: magenta line $\cure=0.1$, green line $\cure=0.2$, cyan
    line $\cure=0.4$, panel (b) shows the dependence of the rate of cure 
    on the degree of cure for three different temperatures: magenta
    line $\temp = 100$\degree C, green line $\temp = 130$\degree C, and cyan line
    $\temp = 180$\degree C.}
  \label{fig:cure}
\end{figure}

\begin{table}[t]
  \centering
  \begin{tabular}{llll}
    \toprule
    physical meaning & notation & value & units \\
    \midrule
    Young's modulus               & $E_f$        & $230$         & GPa            \\ 
    data fitting constant         & $\beta$      & $350$         & (unitless)            \\
    initial fibre volume fraction & $\vzero$        & $0.55$        & (unitless)     \\
    maximum fibre volume fraction & $\vmax$        & $0.68$        & (unitless)     \\
    \midrule
    radius of fibre               & $r_f$        & $4\ee{-3}$    & mm             \\
    Kozeny constant              & $\kozeny$          & $0.2$         & 1/s \QM            \\
    \midrule
    viscosity asymptote     & $\viscosinf$ & $3.45\ee{-10}$ & Pa$\cdot$s \QM \\ 
    activation energy \note{of what?}    & $E_\viscos$ & $7.6536\ee{5}$& J/mol\\
    universal gas constant        & $R$          & $8.617$       & eV/K           \\
    degree of cure at  gelation   & $\curegel$   & $0.47$        & (unitless)            \\ 
    data fitting constant        &$A$ & $3.8$&(unitless)\\
    data fitting constant        &$B$ & $2.5$&(unitless)\\
    \midrule
     data fitting constant             &$A_\cure$&$1.53\ee{5}$& 1/s \\
     activation energy             &$E_\cure$& $6.65\ee{4}$  & J/mol          \\
     data fitting constant             &$m$&$0.813$& (unitless) \\
     data fitting constant             &$n$&$2.74$& (unitless) \\
    \bottomrule
  \end{tabular}
  \caption{Physical constants. The parts separated by line explain variable from different equations, the first part corresponds to stress as given in \eqref{eq:stress-nonlin}, the second part corresponds to permeability as given in \eqref{eq:permeability}, the third part corresponds to viscosity as given in \eqref{eq:viscosity}, the fourth part corresponds to cure as given in \eqref{eq:cure}. If the variable is used in more than one equation, it is explained at a place of its first occurrence.}
  \label{tab:phys-const}
\end{table}


\subsection{Procedure of Fabrication}

The resin is cured in a chamber of high pressure and temperature
called autoclave. During this procedure a specific temperature
protocol is followed.  The initial temperature $\temp_0 = 30$\degree C is
increased during $30$ minutes up to the temperature of
$\temp_f = 107.222$\degree C at which value it stays for another half an
hour to allow the resin being expulsed from the fibre bed. After the
first hour of the procedure ($t=3600$~s) the temperature undergoes a
second increase for a duration of $30$ minutes up to the temperature
of $176.66667$\degree C, at which value it stays for another half an
hour. \note{check the duration} The purpose of the increased
temperature is to elicit the cure of the resin.

The simulation of the autoclave cure \eqref{eq:cure} was done using
forward Euler method
%
\begin{align}
  \cure_{i+1} =    \cure_i + \dt \left[A_\cure \exp{-\frac{E_\cure }{R\temp_i}} \cure_i^m (1-\cure_i)^n\right], \label{eq:cure-fe}
\end{align}
%
where $\cure_i \approx \cure(t_i)$ and $\cure(t_0)=\ee{-3}$. The
time step size was choosen as $\dt = 60$~s.
%
\figref{fig:cure-exp}(a) shows the time evolution of the cure of the
resin during a simulation of autoclave temperature protocol, which is
shown in green lines on right axis of (b).  Panel (a) shows that the
cure is negligible for the first $60$--$80$ minutes of the
procedure. Then the cure undergoes a fast upstroke and at $2$ hours
settles at values close to $0.8$.
%
\note{the equation \eqref{eq:cure} might not be realistic for high values of cure}

Panel (b) in \figref{fig:cure-exp} shows the time evolution of
viscosity in the same process as described above. The simulation
results are used to calculate the viscosity using equation
\eqref{eq:viscosity}. As the cure stays close to zero for the first
$80$~minutes, the viscosity can be related to the temperature
profile. After the first $30$~minutes, the resin is in a fluid state,
which allows its flow out of the fibre bed. In the next increase of
temperature, the viscosity further reduces, but then shoots up, as the
cure increases. The graph ends around the point, where the viscosity
passes the gelation point $\curegel$ at which the resin becomes solid
and beyond which the viscosity model is no longer valid.

\begin{figure}
  \centering
  \includegraphics{cure.pdf}
  \caption{Time evolution of the cure of the resin in the autoclave
    experiment (a), and related evolution of the viscosity of the
    resin (b). The ODE \eqref{eq:cure} is subject to the autoclave
    temperature protocol (shown in green lines, labels on the right
    axis in (b) applies to both panels).}
  \label{fig:cure-exp}
\end{figure}


\subsection{Using Viscosity Model}


Implementing the viscosity model into the model of the heat equation
is rather straight forward, as the viscosity is constant throughout
the sample and hence decoupled spatialy. 
%
We use \eqref{eq:darcy} for a generic permeability and viscosity
values and update \eqref{eq:velocity} as
%
\begin{align}
 \press_{,x} =\frac{\viscos(t)}{\permeab(x)}    u_{,t}(x,t) ,  \label{eq:darcy-link}
\end{align}
%
Expression \eqref{eq:darcy-link} is substituted to boundary value problem \eqref{eq:bvp-flow}
so we obtain
%
\begin{align}
   u_{,xx}(x,t) - \frac{\viscos(t)}{\permeab(x)}    u_{,t}(x,t) + f(x) =0.
\end{align}
%
Again, we solve this model using the finite difference method. This
requires using the formulas \eqref{eq:fin-diff} to get a numerical scheme as
%
\begin{align}
  \disp{i}{j+1}  =   \disp{i}{j} + \dt\frac{\permeab(x)}{\viscos_i}\left(\frac{\disp{i+1}{j} - 2\disp{i}{j} + \disp{i-1}{j}}{\dx^2}   + f(x_i) \right). 
\end{align}
%
where viscosity $\viscos_i = \viscos(\temp(t_i),\cure_i)$ is obtained
using expression \eqref{eq:viscosity} and the cure $\cure_i$ is
approximated using equation \eqref{eq:cure-fe}. The permeability
$K=3.56\ee{-6}$ which corresponds to the fibre volume fraction $\vfrac = 0.6$.

\subsection{Fibre Volume Fraction and Effects on Permeability and Stiffness}

Considering the case of static equilibrium $\stress_{,x}(x,t)=0$, and
using \eqref{eq:flow-stress}, \eqref{eq:darcy-link} and assuming that
the Young's modulus $\young = 1$ we can write
%
\begin{align}
  0 = \efstress_{,x} - \frac{\viscos(t)}{\permeab(x)}    u_{,t}, \label{eq:flow-stres_x}
\end{align}
%
where $u_{,t}$ is given in \eqref{eq:fin-diff} allows to derive the
time integration as
\begin{align}
  \disp{}{j+1} = \disp{}{j} + \dt \frac{\permeab(x)}{\viscos(t_j)}\efstress_{,x} + \cO{\dt^2},
\end{align}
%
where $\disp{}{j+1} \approx u(x,t_j)$.

, to derive the
formula for $ \efstress_{,x}$ we need to do a couple of intermediate
calculations.  First, we write the formula for strain as a finite
difference
%
\begin{align}
  u_{,x} =& \frac{\disp{i+1}{j} - \disp{i}{j}}{\dx} + \cO{\dx},
\end{align}
%
which allows to find the fibre volume fraction using
\eqref{eq:fibre-volume} as
%
\begin{align}
  \vfrac(x_i,t_j) = \frac{\vzero}{1+\frac{\disp{i+1}{j} - \disp{i}{j}}{\dx} + \cO{\dx}}
  = \frac{\vzero{\dx}}{{\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}}
  \label{eq:vfrac-findiff}.
\end{align}
%
We substitute this result into \eqref{eq:stress-nonlin} and obtain
%
\begin{align}
  &\efstr{i}{j} = \spring
  \frac{\frac{\vfrac(x_i,t_j)}{\vzero} -1}{\left(\frac{1}{\vfrac(x_i,t_j)} -\frac{1}{\vmax}\right)^4} = \spring
  \frac{\frac{\frac{\vzero}{1+\frac{\disp{i+1}{j} - \disp{i}{j}}{\dx} + \cO{\dx}}}{\vzero} -1}{\left(\frac{1}{\frac{\vzero}{1+\frac{\disp{i+1}{j} - \disp{i}{j}}{\dx} + \cO{\dx}}} -\frac{1}{\vmax}\right)^4} = \spring
  \frac{{\frac{1}{1+\frac{\disp{i+1}{j} - \disp{i}{j}}{\dx} + \cO{\dx}}} -1}{\left({\frac{1+\frac{\disp{i+1}{j} - \disp{i}{j}}{\dx} + \cO{\dx}}{\vzero}} -\frac{1}{\vmax}\right)^4} = \nonumber \\
  &= \spring\frac{\frac{1 -(1+\frac{\disp{i+1}{j} - \disp{i}{j}}{\dx} + \cO{\dx})}{1+\frac{\disp{i+1}{j} - \disp{i}{j}}{\dx} + \cO{\dx}}}{\left({\frac{1+\frac{\disp{i+1}{j} - \disp{i}{j}}{\dx}}{\vzero}} -\frac{1}{\vmax}  + \cO{\dx}\right)^4}
    %
    = \spring\frac{\frac{\frac{\disp{i+1}{j} - \disp{i}{j}}{\dx}}{\frac{\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}{\dx}} + \cO{\dx}}{\left(\frac{1}{\vzero} + \frac{\disp{i+1}{j} - \disp{i}{j}}{\vzero\dx} -\frac{1}{\vmax}  + \cO{\dx}\right)^4}= \nonumber \\
    %
  &= \spring \frac{\frac{\disp{i+1}{j} - \disp{i}{j}}{\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}}%
    {\left(\frac{1}{\vzero} + \frac{\disp{i+1}{j} - \disp{i}{j}}{\vzero\dx} -\frac{1}{\vmax}  + \cO{\dx}\right)^4} + \cO{\dx} = \nonumber \\
    %%%%%%%%%%
  &= \spring {\frac{\disp{i+1}{j} - \disp{i}{j}}{\left(\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}\right){\left(\frac{1}{\vzero} + \frac{\disp{i+1}{j} - \disp{i}{j}}{\vzero\dx} -\frac{1}{\vmax}  + \cO{\dx}\right)^4}}}%
    + \cO{\dx} =  \nonumber \\
  &= \spring \frac{\disp{i+1}{j} - \disp{i}{j}}{\left(\dx + \disp{i+1}{j} - \disp{i}{j}\right){\left(\frac{1}{\vzero} + \frac{\disp{i+1}{j} - \disp{i}{j}}{\vzero\dx} -\frac{1}{\vmax}  + \cO{\dx}\right)^4} + \cO{\dx^2}}
    + \cO{\dx}
     % =  \nonumber \\
  , \label{eq:stress-nonlin-findif}
\end{align}
where $\efstr{i+1}{j} = \efstress(x_i,t_j)$.

That can be used to find $\efstress_{,x}$ as
%
\begin{align}
  &\efstress_{,x} = \frac{\efstr{i+1}{j} - \efstr{i}{j}}{\dx} + \cO{\dx}= \nonumber \\
  &=
  \spring
   \frac{\disp{i+2}{j} - \disp{i+1}{j}}{{\dx}\left(\dx + \disp{i+2}{j} - \disp{i+1}{j}\right){\left(\frac{1}{\vzero} + \frac{\disp{i+2}{j} - \disp{i+1}{j}}{\vzero\dx} -\frac{1}{\vmax}  + \cO{\dx}\right)^4} + \cO{\dx^3}}- \nonumber\\ 
   %%
  &-\spring
  \frac{\disp{i+1}{j} - \disp{i}{j}}{{\dx}\left(\dx + \disp{i+1}{j} - \disp{i}{j}\right){\left(\frac{1}{\vzero} + \frac{\disp{i+1}{j} - \disp{i}{j}}{\vzero\dx} -\frac{1}{\vmax}  + \cO{\dx}\right)^4} + \cO{\dx^3}}%
   + \cO{\dx},
\end{align}
%
which can be substituted to \eqref{eq:flow-stres_x}. Before that we
also need to find the space dependent permeability
$\permeab(\vfrac(x))$ from \eqref{eq:permeability}  and \eqref{eq:vfrac-findiff} as

\begin{align}
  &\permeab(x_i) = \frac{r_f^2}{4\kozeny} \frac{(1-\vfrac)^3}{\vfrac^2}=
    \frac{r_f^2}{4\kozeny} \frac{\left(1-\frac{\vzero{\dx}}{{\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}}\right)^3}{\left(\frac{\vzero{\dx}}{{\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}}\right)^2}=\nonumber \\
  %
  &=\frac{r_f^2}{4\kozeny} \frac{\left(\frac{{\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}-\vzero{\dx}}{{\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}}\right)^3}{\left(\frac{\vzero{\dx}}{{\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}}\right)^2}
    =\frac{r_f^2}{4\kozeny} \frac{\frac{\left({\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}-\vzero{\dx}\right)^3}{\left({\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}\right)^3}}{\frac{\left(\vzero{\dx}\right)^2}{\left({\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}\right)^2}}=\nonumber \\
    %
  &  =\frac{r_f^2}{4\left(\vzero{\dx}\right)^2\kozeny} {\frac{\left({\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}-\vzero{\dx}\right)^3}{{\dx + \disp{i+1}{j} - \disp{i}{j} + \cO{\dx^2}}}}.
\end{align}



\subsection{TODO}

Here we describe the implementation of the models including the
physical units. To simplify the process, the implementation will be
done in several steps as follows
%
\begin{enumerate}
\item constant fibre volume fraction of $\vfrac = 0.6$, temperature $\temp=\QM$, and degree of cure $\cure$; this gives a linear model with appropriate physical units, \label{en:all-const}
\item fibre volume fraction dependent on the deformation according to formula \eqref{eq:fibre-volume}, the remaining variables as in previous case,
\item ramped temperature $\temp_{,t} = q$ which leads to variable degree of cure $\cure$.
\end{enumerate}

Solving the point \ref{en:all-const} we write down the model using the physical units.



%% -----------------------------------------------------------------
%% GLOSSARY
%% -----------------------------------------------------------------

%% -----------------------------------------------------------------
%% BIBLIOGRAPHY
%% -----------------------------------------------------------------
\bibliographystyle{apalike}
\bibliography{references}

%% -----------------------------------------------------------------
%% INDEX (IF ANY)
%% -----------------------------------------------------------------

\end{document}