# Problem with time loop- variables are not updating and being used in the next loop

168 views (last 30 days)

Show older comments

Hi all, I hope you doing well!

i am trying to simulate the motion of ionic species under electric field. I am using a set of equations such as Poisson equation, drift, and difussion equation (where the coefficient of difussion is related to the electrical mobility through the Eistein relationship), and mass balance equation for continuity. My code seems fine and it give me a reasonble results for the first time step. However, nothing change after that and i get the same results for each time step. The ionic species are supposed to move with the electric field and the concentration profile will change. That new concentration profile will be used to calculate the new potential, new electric field, and new charge density. But, anything i do seems to work to be able to see those changes.

##### 12 Comments

Umar
on 26 Jul 2024 at 2:55

Edited: Umar
on 26 Jul 2024 at 3:15

Shallower slope

To modify the code and achieve a shallower slope in the final concentration profiles plot, I had to adjust the initial concentration profile and the diffusion coefficients.

Adjust the initial concentration profile:

Currently, the code sets the initial concentration profile as a linearly increasing function from 0 to c1_0 (3e27 ions/m^3). To achieve a shallower slope, you can modify the initial concentration profile to start from a higher concentration and decrease linearly towards c1_0. For example, you can update the line c1_profile(:, 1) = linspace(0, c1_0, N)'; to c1_profile(:, 1) = linspace(2*c1_0, c1_0, N)';.

Adjust the diffusion coefficients:

The code currently sets the diffusion coefficients (D1) as a linearly increasing function from k_B * Tp * mu1_profile to 2 * k_B * Tp * mu1_profile. To achieve a shallower slope, you can modify the diffusion coefficients to increase less steeply. For example, you can update the line D1 = ((k_B * Tp) / e) * mu1_profile .* linspace(1, 2, N)'; to D1 = ((k_B * Tp) / e) * mu1_profile .* linspace(1, 1.5, N)';.

So, by making these modifications, the initial concentration profile will start from a higher concentration and decrease linearly towards c1_0, resulting in a shallower slope in the final concentration profiles plot. Additionally, the adjusted diffusion coefficients will increase less steeply, further contributing to the desired effect. Here is the updated code,

% Set format to long for higher precision

format long;

% Constants and Parameters

d_microns = 50; % Thickness of glass in micrometers

d = d_microns * 1e-6; % Thickness in meters

N = 500; % Number of spatial points

L = d; % Length of the glass in meters

x = linspace(0, L, N); % Spatial grid

dx = L / (N - 1); % Grid spacing corrected

V0 = 40; % Voltage applied across the glass (in volts)

e = 1.602e-19; % Elementary charge in coulombs

epsilon_0 = 8.854e-12; % Permittivity of free space (F/m)

er = 6; % Relative permittivity of the glass

k_B = 1.381e-23; % Boltzmann constant

E_a = 0.8 * e; % Activation energy

mu_0 = 1e-5; % Pre-exponential factor for mobility (m^2/V/s)

Tp = 250 + 273.15; % Poling temperature in Kelvin

% Initial concentrations (in ions/m^3)

c1_0 = 3e27;

% Charges of the ionic species (z1, z2, z3)

z1 = 1; % Charge of c1

% Function to calculate mobility as a function of temperature

mu1_profile = mu_0 * exp(-E_a / (k_B * Tp));

% Diffusion coefficients using Einstein relationship

D1 = ((k_B * Tp) / e) * mu1_profile .* linspace(1, 1.5, N)'; % Adjusted diffusion coefficients

% Time parameters

dt = 1; % Time step in seconds

T = 10; % Total simulation time in seconds

Nt = ceil(T / dt); % Number of time steps

% Initialize concentration profiles

c1_profile = zeros(N, Nt);

% Set initial concentrations

c1_profile(:, 1) = linspace(2*c1_0, c1_0, N)';

% Initialize potential phi and electric field E

phi = zeros(N, Nt);

E_profile = zeros(N, Nt);

% Set initial potential profile

phi(:, 1) = linspace(V0, 0, N)'; % Linear gradient from V0 to 0

% Solver parameters

maxIter = 1000; % Maximum iterations for potential solver

tol = 1e-4; % Tolerance for potential convergence

% Time evolution loop

for t = 2:Nt

% Solve Poisson equation in 1D using finite differences

for iter = 1:maxIter

% Update charge density rho based on current potential phi

rho = z1 * e * c1_profile(:, t);

% Save current phi for comparison

phi_old = phi(:, t);

% Update potential phi using finite difference method

phi_new = phi(:, t); % Start with the current phi values

phi_new(2:end-1) = 0.5 * (phi(1:end-2, t) + phi(3:end, t) - rho(2:end-1) * dx^2 / (epsilon_0 * er));

% Apply boundary conditions

phi_new(1) = V0;

phi_new(N) = 0;

% Check for convergence

if max(abs(phi_new - phi_old)) < tol

phi(:, t + 1) = phi_new;

break;

end

end

% Update concentrations using drift-diffusion equations

% Calculate electric field E as negative gradient of phi

E_profile(:, t) = -gradient(phi(:, t), dx);

% Calculate drift velocities (assuming linear drift)

v1 = mu1_profile * E_profile(:, t); % Adjusted drift velocities

% Calculate spatial gradients for concentration

dc1_dx = gradient(c1_profile(:, t), dx);

% Diffusion flux based on Fick's law

J_diffusion1 = -D1 .* dc1_dx;

% Drift flux based on drift velocity

J_drift1 = -z1 * e * v1 .* c1_profile(:, t);

% Mass balance equation: d(c)/dt = -div(J)

dc1_dt = -gradient(J_diffusion1 + J_drift1, dx);

% Update concentrations using forward Euler method

c1_profile(:, t + 1) = c1_profile(:, t) + dt .* dc1_dt;

end

% Plotting

figure;

subplot(3, 1, 1);

plot(x * 1e6, phi(:, :), 'LineWidth', 1.5);

xlabel('Position (\mum)');

ylabel('Electric Potential (V)');

title('Electric Potential Profile');

grid on;

legend(arrayfun(@(n) sprintf('t = %.2f s', n * dt), 0:dt:(T - dt), 'UniformOutput', false));

subplot(3, 1, 2);

plot(x * 1e6, E_profile(:, :), 'LineWidth', 1.5);

xlabel('Position (\mum)');

ylabel('Electric Field (V/m)');

title('Electric Field Profile');

grid on;

subplot(3, 1, 3);

plot(x * 1e6, c1_profile(:, 1), 'b-', 'LineWidth', 2); % Initial concentration profile in blue

hold on;

plot(x * 1e6, c1_profile(:, end), 'r-', 'LineWidth', 2); % Final concentration profile in red

hold off;

xlabel('Position (\mum)');

ylabel('Concentration (m^{-3})');

title('Initial and Final Concentration Profiles');

legend('Initial', 'Final', 'Location', 'best');

grid on;

Please see attached plot

Umar
on 26 Jul 2024 at 2:58

### Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!