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

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.

[Umar]

Hi Cesar,

Please ignore my previous comments. I was not finished with my thoughts and accidentally posted them. You attached two script files and mentioned “ anything i do seems to work to be able to see those changes”. So, I am trying to understand the reasons behind attaching two script files and what is your aim to accomplish with these scripts or did you attach both of them by accident.

[Cesar]

Thank you Umar for your response. i am trying to figure out what its different from the modifications you sent and my code but I can't find any obvious differences. Can you point out the modifications to the previous code? Either way, I tried to incorporate your modifiucations but i am still getting the same results. Thank you

[Cesar]

Hi Umar, I attached two of them by accident. Sorry for the confusion i caused.

[Umar]

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.

[Cesar]

Hi Umar, see a diagram below of how the Potential (V), electric field (E), and concentration of c1 should change with time. The migration of ions, in this case c1, should reflect changes in V and E as the ions moves in time and space. I hope this helps.

[Cesar]

c2 and c3 are immboile for now so they shouldnt affect V and E, but eventually they will once I define a mobility for them.

[Umar]

Hi Cesar,

So, you are struggling with the migration of ions, in this case c1, not reflecting changes in V and E in your attached script file and due to that ions are not moving in time and space theoretically. This is the only issue you are dealing with c2 and c3 being exception since they are immboile as you mentioned. Did I comprehend it correctly, please correct me if I am wrong.

[Cesar]

Yes, thats correct

[Umar]

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……

[Umar]

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…

[Umar]

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]

As I mentioned in my analysis and observation part of your code, 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. Hope, now you should be able to resolve the problem.