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

Hi Cesar,

My next question, tell me what are the issues exactly you are facing with “migrationofions.m”, what seems to be the problem exactly or what are your tasks to accomplish with this script file. I did execute the code and end up getting this attached plot. Since you created this script, you know more about it then I do.

Hi Cesar,

Please see my analysis and observations about your code and solution that you were seeking.

I really loved what you did with your code but at the same time it involves lot of complexities and it is quite interesting and intriguing that it provides a simulation of the drift-diffusion process in a glass material and allows for the visualization of the resulting profiles which can be used to study the behavior of charge carriers and their concentrations under the influence of an applied electric field in a glass material. Now, going back to your code, the variable c1_0 represents the initial concentration of the ionic species c1 in ions/m^3 and after some experiments with your code, it was found out that changing the value of c1_0 should affect the concentration profile of c1 over time and, consequently, the behavior of the plots. But before we delve in to parameters that should reflect changes for E and V, c1_0 , Let's explore firstly the expected changes in the plots when modifying the value of c1_0.

Electric Potential Profile plot

The first plot shows the electric potential profile as a function of position. It represents the voltage distribution across the glass. Changing the value of c1_0 will not directly affect the electric potential profile. Therefore, the plot will remain the same regardless of the value of c1_0.

Electric Field Profile plot

The second plot displays the electric field profile as a function of position. The electric field is calculated as the negative gradient of the electric potential. Since changing c1_0 does not affect the electric potential profile, the electric field profile will also remain unchanged.

Final Concentration Profiles plot

The third plot illustrates the concentration profiles of c1 as a function of position. This plot will be affected by changing the value of c1_0. Let's analyze the expected behavior:

Higher c1_0: Increasing the initial concentration c1_0 will result in higher concentrations of c1 throughout the glass. The concentration profile will exhibit a steeper slope, indicating a higher concentration gradient. The peak concentration will also be higher.

Lower c1_0: Decreasing the initial concentration c1_0 will lead to lower concentrations of c1 in the glass. The concentration profile will have a shallower slope, indicating a lower concentration gradient. The peak concentration will also be lower. By changing the value of c1_0, you can observe how the concentration profile of c1 evolves over time. This information is crucial for understanding the behavior of the ionic species in the glass and its impact on the overall system.

My observation based on increase or decrease value of c1_0

I was observing that when you increase or decrease the values of c1_0, the concentration profile in the plot does not show any change in slope. So, this deals with concept of physics and your code. Your code provided simulates the time evolution of the concentration profile in a glass sample subjected to an applied electric field. The glass sample is assumed to have a fixed thickness (d_microns) and a linearly varying electric potential (phi) across its length (L). The concentration profile of a single ionic species (c1_profile) is calculated using drift-diffusion equations. So, I had to break down the code and understand relevant parts:

In your code snippet below, you setup up the necessary constants, parameters, and initial conditions for the simulation. The concentration of the ionic species (c1_0) is set to a specific value (6e27 ions/m^3) at the beginning of the simulation.

% 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 = 6e27;

% 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;

% 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) = c1_0 * ones(N, 1);

Now, the code below performs a time evolution loop to update the concentration profiles at each time step. The Poisson equation is solved using finite differences to update the electric potential (phi) based on the charge density (rho). Then, the concentration profiles are updated using drift-diffusion equations.

% 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
    ...
end

Now, the final piece of code which calculates the electric field (E_profile) as the negative gradient of the electric potential (phi). The drift velocities (v1) are calculated assuming linear drift. The spatial gradients of the concentration (dc1_dx) are calculated using the gradient function. The diffusion flux (J_diffusion1) is calculated based on Fick's law, and the drift flux (J_drift1) is calculated based on the drift velocity. Finally, the concentration is updated using the forward Euler method.

% 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);

% 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;

Now, let's address my observation that why changing the value of c1_0 does not affect the slope of the concentration profile in the plot. This is really important and plays a huge role if you are physics major. The concentration profile is determined by the interplay between diffusion and drift processes. The diffusion flux is proportional to the concentration gradient, while the drift flux is proportional to the drift velocity and the concentration. In this simulation, the drift velocity is assumed to be linearly dependent on the electric field, and the concentration is updated based on the diffusion and drift fluxes. When the value of c1_0 is changed, it affects the initial concentration profile (c1_profile(:, 1)), but it does not directly affect the subsequent concentration profiles. The subsequent concentration profiles are determined by the diffusion and drift processes, which depend on the spatial gradients of the concentration and the electric field. In your provided code, the initial concentration profile is set to c1_0 for all spatial points. As the simulation progresses, the concentration profile evolves based on the diffusion and drift processes. The slope of the concentration profile in the plot is determined by the interplay between these processes and the applied electric field. Changing the value of c1_0 will affect the initial concentration profile, but it may not have a noticeable effect on the slope of the concentration profile in the plot because the subsequent concentration profiles are determined by the diffusion and drift processes, which depend on the spatial gradients and the electric field. So, to observe a significant change in the slope of the concentration profile, you may need to modify other parameters that affect the diffusion and drift processes, such as the diffusion coefficient (D1), the mobility (mu1_profile), or the electric field (phi). In nutshell, changing the value of c1_0 may not have a direct and noticeable effect on the slope of the concentration profile in the plot because the subsequent concentration profiles are determined by the diffusion and drift processes, which depend on other factors such as the spatial gradients and the electric field. To observe a significant change in the slope, you may need to modify other parameters that affect these processes. In my next post, I will provide plots for shallower and steeper slopes.

Cont’d……

Steeper slope

To modify the code and achieve steeper slope in the final concentration profiles plot, I had to adjust the initial concentration profile and the diffusion coefficients. In the updated code below, I had to make two modifications:

Adjusted the diffusion coefficients (D1): We have multiplied the diffusion coefficients by a factor of 5, resulting in steeper diffusion gradients. This change will lead to steeper slopes in the final concentration profiles.

Adjusted the range of the initial concentration profile: We have modified the range of the initial concentration profile to span from 0 to c1_0 (the initial concentration) using the linspace function. This change ensures that the initial concentration profile covers the entire range of concentrations.

Please note that you can further adjust the parameters and coefficients in the code to achieve the desired slope steepness in the final concentration profiles. Experimenting with different values can help you find the optimal settings for your specific requirements. Here's 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, 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(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

Please see attached plot

Cont’d…