Problem
For the sake of simplicity, let's consider stationary 1D heat conduction with Neumann boundaries:
$$ \frac{d}{dx}(\kappa \frac{du}{dx}) = 0 \tag{1} $$ $$ \frac{du}{dx}|_0 = F_0 \tag{2} $$ $$ \frac{du}{dx}|_1 = F_1 \tag{3} $$
- $ x \in [0,1] $
- $ u: [0,1] \rightarrow \mathbb{R} $
- $ \kappa = \kappa(u(x)) $
The goal is to end up with the following formulation: $K\hat{u}=f$ , where $K=K(u)$ is the "stiffness" matrix, $\hat{u}$ is the vector of coefficients of the shape functions, and $f=f(u)$ the load vector.
Linear case
In the linear case, $\kappa$ is constant and the (simplified) derivation looks like this:
Choose an appropriate space of ansatz functions $V : [0,1] \rightarrow \mathbb{R}$, and assume that the solution $u$ and test function $v$ can be represented as a linear combination of shape functions $N_i(x) \in V$: $$ u = \sum_i N_i \hat{u}_i \tag{4}$$ $$ v = \sum_i N_i \tag{5}$$
The weak form follows from multiplying the heat equation by the test function and integrating over the domain: $$ \int_0^1{v\frac{d}{dx}(\kappa \frac{du}{dx})dx} = \kappa \int_0^1{v\frac{d^2u}{dx^2}dx} = 0 \tag{6} $$
Integrate by parts to get rid of the second derivative and sneak in the boundaries: $$ \kappa [v\frac{du}{dx}]_0^1 - \kappa\int_0^1{\frac{dv}{dx}\frac{du}{dx}dx} = 0 \tag{7} $$
Substitute $u$, $v$, and the boundaries (writing it in index notation to avoid the sums): $$ \hat{u}_i \kappa \int_0^1{\frac{dN_i}{dx}\frac{dN_j}{dx}dx} = \kappa(N_i(1)F_1 - N_i(0)F_0) \tag{8} $$
The FE components can be named and we're done:
- $K_{ij}:=\kappa \int_0^1{ \frac{dN_i}{dx}\frac{dN_j}{dx}dx }$
- $ f_i := \kappa(N_i(1)F_1 - N_i(0)F_0) $
Nonlinear case
Now, $\kappa$ is a function of $u$, so it cannot be removed from the integral like in $(6)$. Subsequently, integration by parts and the substitution of $u$ and $v$ become impossible (or more complicated?). This is the step I need help with. As I already mentioned, the goal is to end up with the following form:
$$ K\hat{u} = f \tag{9} $$
The question is how to get $K$ and $f$, both of which will be functions of $u$.
In this particular case you can just try to perform a few picard iterations. You start with some initial approximation $u^{(0)}$ and, for $k \geq 0$ you compute $u^{(k+1)}$ by solving $$ \frac{d}{dx}\left(k(u^{(k)}(x)) \frac{d u^{(k+1)}}{dx} \right) =0 \quad \textrm{ + b.c} $$
In each iteration you solve a linear problem which is not hard to formulate, as the diffusion is a known function of $x$.