I want to solve a diffusion analysis using finite elements. According to fick's law, governing equation is

$$\frac{\partial h}{\partial t} = D \frac{\partial^{2} h}{\partial x^{2}}$$

. h is relative humidity and D is moisture diffusion coefficient and it is relative humidity dependence.

$$D = 2.26 \times 10^{-6} \Bigg ( 0.05 + \frac{1-0.05}{1+\Big ( \frac{1-h}{1-0.78} \Big )^{4}} \Bigg )$$

Weak form is

$$C \frac{\partial h}{\partial t} + Kh = F$$

When I assemble the diffusivity matrix K, I don't know how I deal with D.

$$B = \begin{bmatrix} \frac{1}{l} & \frac{1}{l} \end{bmatrix}$$

$$K = A \int_{0}^{l} B^{T}DB dx$$

$$K_{e} = \frac{A}{l} \begin{bmatrix} D(h_{1}) + D(h_{2}) & D(h_{1}) + D(h_{2}) \\ D(h_{1}) + D(h_{2}) & D(h_{1}) + D(h_{2}) \end{bmatrix}$$

Is that right?

  • 4
    $\begingroup$ Are you sure the equation is correct? Usually, nonlinear diffusion looks like $\frac{\partial h}{\partial t} = \frac{\partial}{\partial x} \left( D \frac{\partial h}{\partial x} \right)$ in 1D, or $\frac{\partial h}{\partial t} = \nabla \cdot \left( D \nabla h \right)$ in $n$D $\endgroup$
    – cos_theta
    Feb 22, 2021 at 14:49
  • $\begingroup$ What do you mean by $D(h_1)$? Is that the diffusivity evaluated at the first node? $\endgroup$ Feb 22, 2021 at 15:51
  • $\begingroup$ Welcome tho Scicomp! I think @cos_theta is right. The diffusivity is dependent on the local humidity. Which means that it can and will change in time, and spatially. Since it changes spatially, the spatial derivative of D has to be taken into account. (In the wiki article on Fick's law of diffusion) $\endgroup$
    – MPIchael
    Feb 22, 2021 at 19:26
  • $\begingroup$ @cos_theta Aren't those equations the same? $\endgroup$
    – Jaemin Lee
    Feb 23, 2021 at 6:36
  • $\begingroup$ @BillGreene yes ℎ_1 is first node. $\endgroup$
    – Jaemin Lee
    Feb 23, 2021 at 6:37

2 Answers 2


Your $B$ matrix is incorrect. It should be

\begin{equation} B=\left[\begin{array}{cc} -\frac{1}{l} & \frac{1}{l} \end{array}\right] \end{equation}

For this two-node element, $h$ is assumed to vary linearly as a function of $x$. To perform the integration shown in your equation for $K$ it is typically assumed that $D$ is evaluated at the element center


where $h_1$ and $h_2$ are the values of $h$ at the two element nodes. $D_c$ is assumed to be constant over the length of the element.

The alternative of assuming that the appropriate constant $D_c$ is the average of the values at the two nodes

$$D_c = \frac{D(h_1)+D(h_2)}{2}$$

is certainly defensible.

Either approximation for $D$ yields

\begin{equation} K=\frac{AD_c}{l} \left[\begin{array}{cc} 1 & -1 \\ -1 & 1 \end{array}\right] \end{equation}


Let's consider only the spatial operator and neglect the time dependence in the problem. Suppose you have discretized your function space $V$ by means of a basis $\overline{\mathrm{span}\{ \varphi_i : i = 1, \dots, N \}} = V$.

The procedure in FEM is always the same:

  1. Multiply equation with arbitrary test function $\varphi \in V$ and integrate over domain $\Omega \subset \mathbb{R}^d$: $$ \int_\Omega \varphi \nabla \cdot \left( D( h ) \nabla h \right) \,\mathrm{d}x.$$

  2. Integrate by parts: $$ \int_{\partial \Omega} \varphi D(h) \left( \nabla h \right) \cdot \vec{n} \,\mathrm{d}x - \int_\Omega \left( \nabla \varphi \right) \cdot \left( D( h ) \nabla h \right) \,\mathrm{d}x, $$ where $\partial \Omega$ is the boundary of $\Omega$ and $\vec{n}$ is the outward-facing unit normal.

  3. Apply boundary conditions to the boundary integral (possibly eliminating it).

  4. Decompose solution in terms of basis functions, $h = \sum_{j=1}^N h_j \varphi_j$ with coefficients $\vec{h} = \left(h_1, \dots, h_n \right) \in \mathbb{R}^N$, and substitute into equation: $$ \int_{\partial \Omega} \varphi D\left(\sum_{k=1}^N h_k \varphi_k \right) \left( \sum_{j=1}^N h_j \nabla \varphi_j \right) \cdot \vec{n} \,\mathrm{d}x - \int_\Omega \left( \nabla \varphi \right) \cdot \left( D\left( \sum_{k=1}^N h_k \varphi_k \right) \sum_{j=1}^N h_j \nabla \varphi_j \right) \,\mathrm{d}x. $$

  5. The equation must hold for all test functions, in particular for the basis functions $\varphi_i$, $i = 1, \dots, N$: $$ \sum_{j=1}^N h_j \int_{\partial \Omega} D\left(\sum_{k=1}^N h_k \varphi_k \right) \varphi_i \left( \nabla \varphi_j \right) \cdot \vec{n} \,\mathrm{d}x - \sum_{j=1}^N h_j \int_\Omega D\left( \sum_{k=1}^N h_k \varphi_k \right) \left( \nabla \varphi_i \right) \cdot \left( \nabla \varphi_j \right) \,\mathrm{d}x. $$

  6. Assemble the matrix $K(\vec{h}) = G(\vec{h}) + S(\vec{h})$, where $$ \begin{aligned} G_{ij}(\vec{h}) &= \int_{\partial \Omega} D\left(\sum_{k=1}^N h_k \varphi_k \right) \varphi_i \left( \nabla \varphi_j \right) \cdot \vec{n} \,\mathrm{d}x, \\ S_{ij}(\vec{h}) &= \int_\Omega D\left( \sum_{k=1}^N h_k \varphi_k \right) \left( \nabla \varphi_i \right) \cdot \left( \nabla \varphi_j \right) \,\mathrm{d}x.\end{aligned}$$ Note that the index of the test function determines the row of the matrix.

    The assembly can be done by carefully figuring out the support of the integrand and applying a suitable quadrature formula (e.g., Gaussian quadrature). To this end, you need to evaluate $h = \sum_{k=1}^N h_k \varphi_k$ on the support of $\left( \nabla \varphi_i \right) \cdot \left( \nabla \varphi_j \right)$. It may help to transform the integral to the reference element.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.