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

160 Ansichten (letzte 30 Tage)
Cesar Nieves am 23 Jul. 2024 um 16:39
Bearbeitet: Umar am 26 Jul. 2024 um 3:15
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 Kommentare10 ältere Kommentare anzeigen10 ältere Kommentare ausblenden
Umar am 26 Jul. 2024 um 2:55
Bearbeitet: Umar am 26 Jul. 2024 um 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.
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)';.
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(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
% 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;
Umar am 26 Jul. 2024 um 2:58
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.

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Basic finden Sie in Help Center und File Exchange

R2022b

### Community Treasure Hunt

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

Start Hunting!

Translated by