Many mathematical models involve \( m+1 \) unknown functions governed by a system of \( m+1 \) differential equations. In abstract form we may denote the unknowns by \( u^{(0)},\ldots, u^{(m)} \) and write the governing equations as $$ \begin{align*} \mathcal{L}_0(u^{(0)},\ldots,u^{(m)}) &= 0,\\ &\vdots\\ \mathcal{L}_{m}(u^{(0)},\ldots,u^{(m)}) &= 0, \end{align*} $$ where \( \mathcal{L}_i \) is some differential operator defining differential equation number \( i \).
There are basically two ways of formulating a variational form for a system of differential equations. The first method treats each equation independently as a scalar equation, while the other method views the total system as a vector equation with a vector function as unknown.
Let us start with the approach that treats one equation at a time. We multiply equation number \( i \) by some test function \( v^{(i)}\in V^{(i)} \) and integrate over the domain: $$ \begin{align} \int_\Omega \mathcal{L}^{(0)}(u^{(0)},\ldots,u^{(m)}) v^{(0)}\dx &= 0, \tag{1}\\ &\vdots \tag{2}\\ \int_\Omega \mathcal{L}^{(m)}(u^{(0)},\ldots,u^{(m)}) v^{(m)}\dx &= 0 \tag{3} \tp \end{align} $$ Terms with second-order derivatives may be integrated by parts, with Neumann conditions inserted in boundary integrals. Let $$ V^{(i)} = \hbox{span}\{\baspsi_0^{(i)},\ldots,\baspsi_{N_i}^{(i)}\},$$ such that $$ u^{(i)} = B^{(i)}(\x) + \sum_{j=0}^{N_i} c_j^{(i)} \baspsi_j^{(i)}(\x), $$ where \( B^{(i)} \) is a boundary function to handle nonzero Dirichlet conditions. Observe that different unknowns live in different spaces with different basis functions and numbers of degrees of freedom.
From the \( m \) equations in the variational forms we can derive \( m \) coupled systems of algebraic equations for the \( \Pi_{i=0}^{m} N_i \) unknown coefficients \( c_j^{(i)} \), \( j=0,\ldots,N_i \), \( i=0,\ldots,m \).
The alternative method for deriving a variational form for a system of differential equations introduces a vector of unknown functions $$ \u = (u^{(0)},\ldots,u^{(m)}),$$ a vector of test functions $$ \v = (u^{(0)},\ldots,u^{(m)}),$$ with $$ \u, \v \in \V = V^{(0)}\times \cdots \times V^{(m)} \tp $$ With nonzero Dirichlet conditions, we have a vector \( \boldsymbol{B} = (B^{(0)},\ldots,B^{(m)}) \) with boundary functions and then it is \( \u - \boldsymbol{B} \) that lies in \( \V \), not \( \u \) itself.
The governing system of differential equations is written $$ \boldsymbol{\mathcal{L}}(\u ) = 0,$$ where $$ \boldsymbol{\mathcal{L}}(\u ) = (\mathcal{L}^{(0)}(\u),\ldots, \mathcal{L}^{(m)}(\u)) \tp $$ The variational form is derived by taking the inner product of the vector of equations and the test function vector: $$ \begin{equation} \int_\Omega \boldsymbol{\mathcal{L}}(\u )\cdot\v = 0\quad\forall\v\in\V\tp \tag{4} \end{equation} $$
Observe that (4) is one scalar equation. To derive systems of algebraic equations for the unknown coefficients in the expansions of the unknown functions, one chooses \( m \) linearly independent \( \v \) vectors to generate \( m \) independent variational forms from (4). The particular choice \( \v = (v^{(0)},0,\ldots,0) \) recovers (1), \( \v = (0,\ldots,0,v^{(m)} \) recovers (3), and \( \v = (0,\ldots,0,v^{(i)},0,\ldots,0) \) recovers the variational form number \( i \), \( \int_\Omega \mathcal{L}^{(i)} v^{(i)}\dx =0 \), in (1)-(3).
We now consider a specific system of two partial differential equations in two space dimensions: $$ \begin{align} \mu \nabla^2 w &= -\beta, \tag{5}\\ \kappa\nabla^2 T &= - \mu ||\nabla w||^2 \tp \tag{6} \end{align} $$ The unknown functions \( w(x,y) \) and \( T(x,y) \) are defined in a domain \( \Omega \), while \( \mu \), \( \beta \), and \( \kappa \) are given constants. The norm in (6) is the standard Euclidean norm: $$ ||\nabla w||^2 = \nabla w\cdot\nabla w = w_x^2 + w_y^2 \tp $$
The boundary conditions associated with (5)-(6) are \( w=0 \) on \( \partial\Omega \) and \( T=T_0 \) on \( \partial\Omega \). Each of the equations (5) and (6) needs one condition at each point on the boundary.
The system (5)-(6) arises from fluid flow in a straight pipe, with the \( z \) axis in the direction of the pipe. The domain \( \Omega \) is a cross section of the pipe, \( w \) is the velocity in the \( z \) direction, \( \mu \) is the viscosity of the fluid, \( \beta \) is the pressure gradient along the pipe, \( T \) is the temperature, and \( \kappa \) is the heat conduction coefficient of the fluid. The equation (5) comes from the Navier-Stokes equations, and (6) follows from the energy equation. The term \( - \mu ||\nabla w||^2 \) models heating of the fluid due to internal friction.
Observe that the system (5)-(6) has only a one-way coupling: \( T \) depends on \( w \), but \( w \) does not depend on \( T \), because we can solve (5) with respect to \( w \) and then (6) with respect to \( T \). Some may argue that this is not a real system of PDEs, but just two scalar PDEs. Nevertheless, the one-way coupling is convenient when comparing different variational forms and different implementations.
Let us first apply the same function space \( V \) for \( w \) and \( T \) (or more precisely, \( w\in V \) and \( T-T_0 \in V \)). With $$ V = \hbox{span}\{\baspsi_0(x,y),\ldots,\baspsi_N(x,y)\}, $$ we write $$ \begin{equation} w = \sum_{j=0}^N c^{(w)}_j \baspsi_j,\quad T = T_0 + \sum_{j=0}^N c^{(T)}_j \baspsi_j\tp \tag{7} \end{equation} $$ Note that \( w \) and \( T \) in (5)-(6) denote the exact solution of the PDEs, while \( w \) and \( T \) (7) are the discrete functions that approximate the exact solution. It should be clear from the context whether a symbol means the exact or approximate solution, but when we need both at the same time, we use a subscript e to denote the exact solution.
Inserting the expansions (7) in the governing PDEs, results in a residual in each equation, $$ \begin{align} R_w &= \mu \nabla^2 w + \beta, \tag{8}\\ R_T &= \kappa\nabla^2 T + \mu ||\nabla w||^2 \tp \tag{9} \end{align} $$ A Galerkin method demands \( R_w \) and \( R_T \) do be orthogonal to \( V \): $$ \begin{align*} \int_\Omega R_w v \dx &=0\quad\forall v\in V,\\ \int_\Omega R_T v \dx &=0\quad\forall v\in V \tp \end{align*} $$ Because of the Dirichlet conditions, \( v=0 \) on \( \partial\Omega \). We integrate the Laplace terms by parts and note that the boundary terms vanish since \( v=0 \) on \( \partial\Omega \): $$ \begin{align} \int_\Omega \mu \nabla w\cdot\nabla v \dx &= \int_\Omega \beta v\dx \quad\forall v\in V, \tag{10}\\ \int_\Omega \kappa \nabla T\cdot\nabla v \dx &= \int_\Omega \mu \nabla w\cdot\nabla w\, v\dx \quad\forall v\in V \tag{11} \tp \end{align} $$
The alternative way of deriving the variational from is to introduce a test vector function \( \v\in\V = V\times V \) and take the inner product of \( \v \) and the residuals, integrated over the domain: $$ \int_{\Omega} (R_w, R_T)\cdot\v \dx = 0\quad\forall\v\in\V \tp $$ With \( \v = (v_0,v_1) \) we get $$ \int_{\Omega} (R_w v_0 + R_T v_1) \dx = 0\quad\forall\v\in\V \tp $$ Integrating the Laplace terms by parts results in $$ \begin{equation} \int_\Omega (\mu\nabla w\cdot\nabla v_0 + \kappa\nabla T\cdot\nabla v_1)\dx = \int_\Omega (\beta v_0 + \mu\nabla w\cdot\nabla w\, v_1)\dx, \quad\forall \v\in\V \tp \tag{12} \end{equation} $$ Choosing \( v_0=v \) and \( v_1=0 \) gives the variational form (10), while \( v_0=0 \) and \( v_1=v \) gives (11).
With the inner product notation, \( (p,q) = \int_\Omega pq\dx \), we can alternatively write (10) and (11) as $$ \begin{align*} (\mu\nabla w,\nabla v) &= (\beta, v) \quad\forall v\in V,\\ (\kappa \nabla T,\nabla v) &= (\mu\nabla w\cdot\nabla w, v)\quad\forall v\in V, \end{align*} $$ or since \( \mu \) and \( \kappa \) are considered constant, $$ \begin{align} \mu (\nabla w,\nabla v) &= (\beta, v) \quad\forall v\in V, \tag{13}\\ \kappa(\nabla T,\nabla v) &= \mu(\nabla w\cdot\nabla w, v)\quad\forall v\in V \tag{14} \tp \end{align} $$
The linear systems governing the coefficients \( c_j^{(w)} \) and \( c_j^{(T)} \), \( j=0,\ldots,N \), are derived by inserting the expansions (7) in (10) and (11), and choosing \( v=\baspsi_i \) for \( i=0,\ldots,N \). The result becomes $$ \begin{align} \sum_{j=0}^N A^{(w)}_{i,j} c^{(w)}_j &= b_i^{(w)},\quad i=0,\ldots,N, \tag{15}\\ \sum_{j=0}^N A^{(T)}_{i,j} c^{(T)}_j &= b_i^{(T)},\quad i=0,\ldots,N, \tag{16}\\ A^{(w)}_{i,j} &= \mu(\nabla \baspsi_j,\nabla \baspsi_i), \tag{17}\\ b_i^{(w)} &= (\beta, \baspsi_i), \tag{18}\\ A^{(T)}_{i,j} &= \kappa(\nabla \baspsi_j,\nabla \baspsi_i), \tag{19}\\ b_i^{(T)} &= \mu((\sum_j c^{(w)}_j\nabla\baspsi_j)\cdot (\sum_k c^{(w)}_k\nabla\baspsi_k), \baspsi_i) \tp \tag{20} \end{align} $$
It can also be instructive to write the linear systems using matrices and vectors. Define \( K \) as the matrix corresponding to the Laplace operator \( \nabla^2 \). That is, \( K_{i,j} = (\nabla \baspsi_j,\nabla \baspsi_i) \). Let us introduce the vectors $$ \begin{align*} b^{(w)} &= (b_0^{(w)},\ldots,b_{N}^{(w)}),\\ b^{(T)} &= (b_0^{(T)},\ldots,b_{N}^{(T)}),\\ c^{(w)} &= (c_0^{(w)},\ldots,c_{N}^{(w)}),\\ c^{(T)} &= (c_0^{(T)},\ldots,c_{N}^{(T)})\tp \end{align*} $$ The system (15)-(16) can now be expressed in matrix-vector form as $$ \begin{align} \mu K c^{(w)} &= b^{(w)}, \tag{21}\\ \kappa K c^{(T)} &= b^{(T)}\tp \tag{22} \end{align} $$
We can solve the first system for \( c^{(w)} \), and then the right-hand side \( b^{(T)} \) is known such that we can solve the second system for \( c^{(T)} \).
Despite the fact that \( w \) can be computed first, without knowing \( T \), we shall now pretend that \( w \) and \( T \) enter a two-way coupling such that we need to derive the algebraic equations as one system for all the unknowns \( c_j^{(w)} \) and \( c_j^{(T)} \), \( j=0,\ldots,N \). This system is nonlinear in \( c_j^{(w)} \) because of the \( \nabla w\cdot\nabla w \) product. To remove this nonlinearity, imagine that we introduce an iteration method where we replace \( \nabla w\cdot\nabla w \) by \( \nabla w_{-}\cdot\nabla w \), \( w_{-} \) being the \( w \) computed in the previous iteration. Then the term \( \nabla w_{-}\cdot\nabla w \) is linear in \( w \) since \( w_{-} \) is known. The total linear system becomes $$ \begin{align} \sum_{j=0}^N A^{(w,w)}_{i,j} c^{(w)}_j + \sum_{j=0}^N A^{(w,T)}_{i,j} c^{(T)}_j &= b_i^{(w)},\quad i=0,\ldots,N, \tag{23}\\ \sum_{j=0}^N A^{(T,w)}_{i,j} c^{(w)}_j + \sum_{j=0}^N A^{(T,T)}_{i,j} c^{(T)}_j &= b_i^{(T)},\quad i=0,\ldots,N, \tag{24}\\ A^{(w,w)}_{i,j} &= \mu(\nabla \baspsi_j,\baspsi_i), \tag{25}\\ A^{(w,T)}_{i,j} &= 0, \tag{26}\\ b_i^{(w)} &= (\beta, \baspsi_i), \tag{27}\\ A^{(w,T)}_{i,j} &= \mu((\nabla\baspsi w_{-})\cdot\nabla\baspsi_j), \baspsi_i), \tag{28}\\ A^{(T,T)}_{i,j} &= \kappa(\nabla \baspsi_j,\baspsi_i), \tag{29}\\ b_i^{(T)} &= 0 \tp \tag{30} \end{align} $$ This system can alternatively be written in matrix-vector form as $$ \begin{align} \mu K c^{(w)} &= b^{(w)}, \tag{31}\\ L c^{(w)} + \kappa K c^{(T)} & =0, \tag{32} \end{align} $$ with \( L \) as the matrix from the \( \nabla w_{-}\cdot\nabla \) operator: \( L_{i,j} = A^{(w,T)}_{i,j} \).
The matrix-vector equations are often conveniently written in block form: $$ \left(\begin{array}{cc} \mu K & 0\\ L & \kappa K \end{array}\right) \left(\begin{array}{c} c^{(w)}\\ c^{(T)} \end{array}\right) = \left(\begin{array}{c} b^{(w)}\\ 0 \end{array}\right), $$
Note that in the general case where all unknowns enter all equations, we have to solve the compound system (23)-(24) since then we cannot utilize the special property that (15) does not involve \( T \) and can be solved first.
When the viscosity depends on the temperature, the \( \mu\nabla^2w \) term must be replaced by \( \nabla\cdot (\mu(T)\nabla w) \), and then \( T \) enters the equation for \( w \). Now we have a two-way coupling since both equations contain \( w \) and \( T \) and therefore must be solved simultaneously Th equation \( \nabla\cdot (\mu(T)\nabla w)=-\beta \) is nonlinear, and if some iteration procedure is invoked, where we use a previously computed \( T_{-} \) in the viscosity (\( \mu(T_{-}) \)), the coefficient is known, and the equation involves only one unknown, \( w \). In that case we are back to the one-way coupled set of PDEs.
We may also formulate our PDE system as a vector equation. To this end, we introduce the vector of unknowns \( \u = (u^{(0)},u^{(1)}) \), where \( u^{(0)}=w \) and \( u^{(1)}=T \). We then have $$ \nabla^2 \u = \left(\begin{array}{cc} -{\mu}^{-1}{\beta}\\ -{\kappa}^{-1}\mu \nabla u^{(0)}\cdot\nabla u^{(0)} \end{array}\right)\tp $$
It is easy to generalize the previous formulation to the case where \( w\in V^{(w)} \) and \( T\in V^{(T)} \), where \( V^{(w)} \) and \( V^{(T)} \) can be different spaces with different numbers of degrees of freedom. For example, we may use quadratic basis functions for \( w \) and linear for \( T \). Approximation of the unknowns by different finite element spaces is known as mixed finite element methods.
We write $$ \begin{align*} V^{(w)} &= \hbox{span}\{\baspsi_0^{(w)},\ldots,\baspsi_{N_w}^{(w)}\},\\ V^{(T)} &= \hbox{span}\{\baspsi_0^{(T)},\ldots,\baspsi_{N_T}^{(T)}\} \tp \end{align*} $$ The next step is to multiply (5) by a test function \( v^{(w)}\in V^{(w)} \) and (6) by a \( v^{(T)}\in V^{(T)} \), integrate by parts and arrive at $$ \begin{align} \int_\Omega \mu \nabla w\cdot\nabla v^{(w)} \dx &= \int_\Omega \beta v^{(w)}\dx \quad\forall v^{(w)}\in V^{(w)}, \tag{33}\\ \int_\Omega \kappa \nabla T\cdot\nabla v^{(T)} \dx &= \int_\Omega \mu \nabla w\cdot\nabla w\, v^{(T)}\dx \quad\forall v^{(T)}\in V^{(T)} \tag{34} \tp \end{align} $$
The compound scalar variational formulation applies a test vector function \( \v = (v^{(w)}, v^{(T)}) \) and reads $$ \begin{equation} \int_\Omega (\mu\nabla w\cdot\nabla v^{(w)} + \kappa\nabla T\cdot\nabla v^{(T)})\dx = \int_\Omega (\beta v^{(w)} + \mu\nabla w\cdot\nabla w\, v^{(T)})\dx, \tag{35} \end{equation} $$ valid \( \forall \v\in\V = V^{(w)}\times V^{(T)} \).
The associated linear system is similar to (15)-(16) or (23)-(24), except that we need to distinguish between \( \baspsi_i^{(w)} \) and \( \baspsi_i^{(T)} \), and the range in the sums over \( j \) must match the number of degrees of freedom in the spaces \( V^{(w)} \) and \( V^{(T)} \). The formulas become $$ \begin{align} \sum_{j=0}^{N_w} A^{(w)}_{i,j} c^{(w)}_j &= b_i^{(w)},\quad i=0,\ldots,N_w, \tag{36}\\ \sum_{j=0}^{N_T} A^{(T)}_{i,j} c^{(T)}_j &= b_i^{(T)},\quad i=0,\ldots,N_T, \tag{37}\\ A^{(w)}_{i,j} &= \mu(\nabla \baspsi_j^{(w)},\baspsi_i^{(w)}), \tag{38}\\ b_i^{(w)} &= (\beta, \baspsi_i^{(w)}), \tag{39}\\ A^{(T)}_{i,j} &= \kappa(\nabla \baspsi_j^{(T)},\baspsi_i^{(T)}), \tag{40}\\ b_i^{(T)} &= \mu(\nabla w_{-}, \baspsi_i^{(T)}) \tp \tag{41} \end{align} $$
In the case we formulate one compound linear system involving both \( c^{(w)}_j \), \( j=0,\ldots,N_w \), and \( c^{(T)}_j \), \( j=0,\ldots,N_T \), (23)-(24) becomes $$ \begin{align} \sum_{j=0}^{N_w} A^{(w,w)}_{i,j} c^{(w)}_j + \sum_{j=0}^{N_T} A^{(w,T)}_{i,j} c^{(T)}_j &= b_i^{(w)},\quad i=0,\ldots,N_w, \tag{42}\\ \sum_{j=0}^{N_w} A^{(T,w)}_{i,j} c^{(w)}_j + \sum_{j=0}^{N_T} A^{(T,T)}_{i,j} c^{(T)}_j &= b_i^{(T)},\quad i=0,\ldots,N_T, \tag{43}\\ A^{(w,w)}_{i,j} &= \mu(\nabla \baspsi_j^{(w)},\baspsi_i^{(w)}), \tag{44}\\ A^{(w,T)}_{i,j} &= 0, \tag{45}\\ b_i^{(w)} &= (\beta, \baspsi_i^{(w)}), \tag{46}\\ A^{(w,T)}_{i,j} &= \mu (\nabla w_{-}\cdot\nabla\baspsi_j^{(w)}), \baspsi_i^{(T)}), \tag{47}\\ A^{(T,T)}_{i,j} &= \kappa(\nabla \baspsi_j^{(T)},\baspsi_i^{(T)}), \tag{48}\\ b_i^{(T)} &= 0 \tp \tag{49} \end{align} $$ The corresponding block form $$ \left(\begin{array}{cc} \mu K^{(w)} & 0\\ L & \kappa K^{(T)} \end{array}\right) \left(\begin{array}{c} c^{(w)}\\ c^{(T)} \end{array}\right) = \left(\begin{array}{c} b^{(w)}\\ 0 \end{array}\right), $$ has square and rectangular block matrices: \( K^{(w)} \) is \( N_w\times N_w \), \( K^{(T)} \) is \( N_T\times N_T \), while \( L \) is \( N_T\times N_w \),
We can reduce the system (5)-(6) to one space dimension, which corresponds to flow in a channel between two flat plates. Alternatively, one may consider flow in a circular pipe, introduce cylindrical coordinates, and utilize the radial symmetry to reduce the equations to a one-dimensional problem in the radial coordinate. The former model becomes $$ \begin{align} \mu w_{xx} &= -\beta, \tag{50}\\ \kappa T_{xx} &= - \mu w_x^2, \tag{51} \end{align} $$ while the model in the radial coordinate \( r \) reads $$ \begin{align} \mu\frac{1}{r}\frac{d}{dr}\left( r\frac{dw}{dr}\right) &= -\beta, \tag{52}\\ \kappa \frac{1}{r}\frac{d}{dr}\left( r\frac{dT}{dr}\right) &= - \mu \left( \frac{dw}{dr}\right)^2 \tp \tag{53} \end{align} $$
The domain for (50)-(51) is \( \Omega = [0,H] \), with boundary conditions \( w(0)=w(H)=0 \) and \( T(0)=T(H)=T_0 \). For (52)-(53) the domain is \( [0,R] \) (\( R \) being the radius of the pipe) and the boundary conditions are \( du/dr = dT/dr =0 \) for \( r=0 \), \( u(R)=0 \), and \( T(R)=T_0 \).
Calculations to be continued...