Fick's 2nd Law of Diffusion using FEM (Direct Stiffness)

Direct stiffness finite element method to solve for c (concentration) at given x (1D space) and t (time) values
155 Downloads
Updated 26 Apr 2019

View License

%% Finite Element Method with D = diffusivity applied to Fick's 2nd Law of Diffusion: dc/dt = D*d^2c/dx^2
% by Prof. Roche C. de Guzman
clear; clc; close('all');
%% Given
xi = 0; xf = 0.6; dx = 0.04; % x range and step size = dx [m]
xL = 0; xU = 0.1; % initial value x lower and upper limits [m]
ti = 0; tf = 0.05; dt = 4e-4; % t range and step size = dt [s]
ci = 2; % initial concentration value [ng/L]
cLU = 8; % initial concentration value within x lower and upper limits [ng/L]
D = 1.5; % diffusivity or diffusion coefficient [m^2/s]
%% Calculations
% Independent variables: x and t
X = xi:dx:xf; nx = numel(X); T = ti:dt:tf; nt = numel(T); % x and t vectors and their number of elements
[x,t] = meshgrid(X,T); x = x'; t = t'; % x and t matrices
% Dependent variable: c
c = ones(nx,nt)*ci; % temporary c(x,t) matrix with rows: c(x) and columns: c(t)
% Initial values and Dirichlet boundary
I = find((X>=xL)&(X<=xU)); % index of lower and upper limits
c(I,1) = cLU; % c at t = 0 for lower and upper limits
% FEM (direct stiffness method): [C]*{Y1} + [K]*{Y} = {F}
where C = damping matrix, Y1 = dc/dt, K = stiffness matrix, Y = c, and F = load vector

Cite As

Roche de Guzman (2024). Fick's 2nd Law of Diffusion using FEM (Direct Stiffness) (https://www.mathworks.com/matlabcentral/fileexchange/71353-fick-s-2nd-law-of-diffusion-using-fem-direct-stiffness), MATLAB Central File Exchange. Retrieved .

MATLAB Release Compatibility
Created with R2019a
Compatible with any release
Platform Compatibility
Windows macOS Linux
Categories
Find more on Physics in Help Center and MATLAB Answers
Communities

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
Version Published Release Notes
1.0.0