diffusivity matrix assembly in nonlinear finite element analysis

Question

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?

Answer 1

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

$$D_c=D((h_1+h_2)/2)$$

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}

Answer 2

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.