terzaghi1D completed

Author: JananiPSrinivasan

examples/matlab/terzaghi1D.m

@@ -54,6 +54,9 @@
 mu = 1e-3;        % Dynamic viscosity [Pa·s]
 rho = 1000;       % Fluid density [kg/m^3]
 g = 9.81;         % Gravity [m/s^2]
+Ks = 1e8;          % Bulk modulus [Pa]
+alpha = 1;         % Biot coefficient
+Ss = 1e-5;         % [1/Pa] Specific storage coefficient
 
 %% Numerical (MOLE) Solution
 L = lap(k, m, dx);       % Mimetic Laplacian operator for diffusion
@@ -62,71 +65,145 @@
 % Boundary modifications to Laplacian
 L(1,:) = 0; L(end,:) = 0;
 
-p_numerical = zeros(length(xgrid), length(times_sec)); % Pressure field [Pa]
-q_numerical = zeros(m+1, length(times_sec));  % Darcy flux [m/s] (size = m+1)
+p_numerical = zeros(length(xgrid), length(times_sec));   % Pressure field [Pa]
+q_numerical = zeros(m+1, length(times_sec));             % Darcy flux [m/s] (size = m+1)
+e_numerical = zeros(length(xgrid), length(times_sec));   % Strain field
+u_numerical = zeros(length(xgrid), length(times_sec));   % Displacement field
 
 % Loop over each specified final time
 for i = 1:length(times_sec)
     t_final = times_sec(i);
-    dt = 0.03;
+    dt = 1.0;                       
     nsteps = round(t_final / dt);
 
     % Uniform Initial Condition : p(x,0) = P0
-    p = P0 * ones(size(xgrid))';
+    p = P0 * ones(size(xgrid))';   
 
     % Time integration using Forward Euler
     for step = 1:nsteps
-        p = p + dt * cf * (L * p);
-        p(1) = 0;             % Dirichlet BC at x = 0
-        p(end) = p(end-1);    % Neumann BC at x = L
+        p = p + dt * cf * (L * p);     % Pressure diffusion
+        p(1) = 0;                      % Dirichlet BC at x = 0
+        p(end) = p(end-1);             % Neumann BC at x = L
     end
 
-    p_numerical(:,i) = p;
+    % Compute strain and displacement from final pressure
+    e = (alpha * p - P0) / Ks;
+    u = cumtrapz(xgrid, e); %cumulative integral using the trapezoidal rule.
 
-    % Compute Darcy flux (numerical)
+    % Store results
+    e_numerical(:,i) = e;
+    u_numerical(:,i) = u;
+    p_numerical(:,i) = p;
+    
+    % Compute darcy flux
     dpdx = G * p;
     q = - (K / mu) * (dpdx - rho * g);
-    q_numerical(:,i) = q(1:m+1);  
+    q_numerical(:,i) = q(1:m+1);
+
+end
+
+%% Mass Conservation Residual Evaluation
+
+% Compute dp/dt and de/dt using backward differences
+dpdt = zeros(size(p_numerical));
+dedt = zeros(size(e_numerical));
+
+for i = 2:length(times_sec)
+    dt_local = times_sec(i) - times_sec(i-1);
+    dpdt(:,i) = (p_numerical(:,i) - p_numerical(:,i-1)) / dt_local;
+    dedt(:,i) = (e_numerical(:,i) - e_numerical(:,i-1)) / dt_local;
+end
 
-    % Print numerical results
+% Compute divergence of q using MOLE div()
+divq = zeros(size(p_numerical));  % same shape as pressure field
+for i = 1:length(times_sec)
+    qx = q_numerical(:,i);                 % q has m+1 values
+    divq(:,i) = div(k, m, dx) * qx;        % returns m+2 values (staggered grid)
+end
+
+% Compute full mass balance residual
+residual = Ss * dpdt + alpha * dedt - divq;
+
+%% Combined Numerical Output
+for i = 1:length(times_sec)
     fprintf('\nNumerical results at t = %.2f hr:\n', times_hr(i));
-    fprintf('|     x (m)     | Numerical p [Pa] | Darcy Flux [m/s] |\n');
-    fprintf('|---------------|------------------|------------------|\n');
+    fprintf('|     x (m)     | Numerical p [Pa] | Darcy Flux [m/s] | Residual [1/s]   |\n');
+    fprintf('|---------------|------------------|------------------|------------------|\n');
     for j = 1:m+1
-        fprintf('| %13.6f | %16.6e | %16.6e |\n', xgrid(j), p_numerical(j,i), q_numerical(j,i));
+        if i == 1
+            fprintf('| %13.6f | %16.6e | %16.6e | %16s |\n', ...
+                xgrid(j), p_numerical(j,i), q_numerical(j,i), 'N/A');
+        else
+            fprintf('| %13.6f | %16.6e | %16.6e | %16.6e |\n', ...
+                xgrid(j), p_numerical(j,i), q_numerical(j,i), residual(j,i));
+        end
     end
 end
 
+
 %% Analytical Solution
 N_max = 100; % Number of Fourier series terms
 p_analytical = zeros(length(xgrid), length(times_sec)); % Pressure field [Pa]
-q_analytical = zeros(m+1, length(times_sec));  % Darcy flux [m/s] (size = m+1)
+q_analytical = zeros(m+1, length(times_sec)); % Darcy flux [m/s] (size = m+1)
+e_analytical = zeros(length(xgrid), length(times_sec)); % Strain field
+u_analytical = zeros(length(xgrid), length(times_sec)); % Displacement field
+residual_analytical = zeros(length(xgrid), length(times_sec));
 
 % Loop over all time snapshots
 for i = 1:length(times_sec)
     t = times_sec(i);
     p = zeros(size(xgrid));
 
-     % Compute analytical solution using truncated Fourier sine series
+    % Compute analytical pressure using Fourier series
     for N = 0:N_max
-        n = 2*N + 1; % Odd terms only (satisfies boundary conditions)
+        n = 2*N + 1; 
         term = (1/n) * sin(n*pi*xgrid/(2*l)) .* exp(-(n^2)*(pi^2)*cf*t/(4*l^2));
         p = p + term;
     end
     p = (4/pi) * P0 * p;
     p_analytical(:,i) = p;
 
-    % Compute Darcy flux (analytical)
-    dpdx = gradient(p, dx);  % basic gradient
+    % Compute Darcy flux
+    dpdx = gradient(p, dx);
     q = - (K / mu) * (dpdx - rho * g);
     q_analytical(:,i) = q(1:m+1);
 
-    % Print analytical results
+    % Compute analytical strain and displacement
+    e = (alpha * p - P0) / Ks;
+    u = cumtrapz(xgrid, e);
+    e_analytical(:,i) = e;
+    u_analytical(:,i) = u;
+    
+
+    % Compute mass conservation residual (analytical)
+    if i == 1
+        residual_analytical(:,i) = NaN(size(p_analytical(:,i)));  % undefined at first time step
+    else
+        dt_local = times_sec(i) - times_sec(i-1);
+        dpdt_ana = (p_analytical(:,i) - p_analytical(:,i-1)) / dt_local;
+        dedt_ana = (e_analytical(:,i) - e_analytical(:,i-1)) / dt_local;
+        divq_ana = div(k, m, dx) * q_analytical(:,i);
+        residual_analytical(:,i) = Ss * dpdt_ana + alpha * dedt_ana - divq_ana;
+    end
+
+
+    % Print combined analytical output
     fprintf('\nAnalytical results at t = %.2f hr:\n', times_hr(i));
-    fprintf('|     x (m)     | Analytical p [Pa] | Darcy Flux [m/s] |\n');
-    fprintf('|---------------|-------------------|------------------|\n');
+    fprintf('|     x (m)     | Analytical p [Pa] | Darcy Flux [m/s] | Residual [1/s]   |\n');
+    fprintf('|---------------|-------------------|------------------|------------------|\n');
     for j = 1:m+1
-        fprintf('| %13.6f | %18.6e | %16.6e |\n', xgrid(j), p_analytical(j,i), q_analytical(j,i));
+        if i == 1
+            fprintf('| %13.6f | %18.6e | %16.6e | %16s |\n', ...
+                xgrid(j), p_analytical(j,i), q_analytical(j,i), 'N/A');
+        else
+            fprintf('| %13.6f | %18.6e | %16.6e | %16.6e |\n', ...
+                xgrid(j), p_analytical(j,i), q_analytical(j,i), residual_analytical(j,i));
+        end
+    end
+
+    if i > 1
+        max_res = max(abs(residual_analytical(2:end-1,i)));
+        fprintf('|     Max Residual (interior) at t = %.2f hr: %e |\n', times_hr(i), max_res);
     end
 end
 
@@ -159,17 +236,97 @@
 legend('Location', 'southeast');
 grid on;
 
+
+
+%% Combined Displacement Plot (Numerical vs Analytical)
+figure('Name','Displacement Field (Numerical vs Analytical)');
+set(gcf,'Color','white');
+
+% Top subplot: Numerical displacement
+subplot(2,1,1);
+hold on;
+for i = 1:length(times_sec)
+    plot(xgrid, u_numerical(:,i), '-o', 'DisplayName', [num2str(times_hr(i)) ' hr']);
+end
+xlabel('x (m)');
+ylabel('Displacement u(x,t) [m]');
+title('Numerical Displacement Field');
+legend('Location','southwest');
+grid on;
+
+% Bottom subplot: Analytical displacement
+subplot(2,1,2);
+hold on;
+for i = 1:length(times_sec)
+    plot(xgrid, u_analytical(:,i), '--s', 'DisplayName', [num2str(times_hr(i)) ' hr']);
+end
+xlabel('x (m)');
+ylabel('Displacement u(x,t) [m]');
+title('Analytical Displacement Field');
+legend('Location','southwest');
+grid on;
+
+
+%% Residual Plot (Numerical vs Analytical)
+figure('Name','Mass Conservation Residuals (Numerical vs Analytical)');
+set(gcf,'Color','white');
+
+% Top subplot: Numerical residual
+subplot(2,1,1);
+hold on;
+for i = 2:length(times_sec)
+    plot(xgrid(2:end-1), residual(2:end-1,i), '-o', 'DisplayName', [num2str(times_hr(i)) ' hr']);
+end
+title('Numerical Mass Conservation Residual');
+xlabel('x (m)');
+ylabel('Residual [1/s]');
+legend('Location','northeast');
+grid on;
+
+% Bottom subplot: Analytical residual
+subplot(2,1,2);
+hold on;
+for i = 2:length(times_sec)
+    plot(xgrid(2:end-1), residual_analytical(2:end-1,i), '--s', 'DisplayName', [num2str(times_hr(i)) ' hr']);
+end
+title('Analytical Mass Conservation Residual');
+xlabel('x (m)');
+ylabel('Residual [1/s]');
+legend('Location','northeast');
+grid on;
+
+
 %% Print relative L2 errors
 fprintf('\nRelative L2 Errors (Numerical vs Analytical):\n');
-fprintf('|     Time [hr]    |   Pressure Error   |   Darcy Flux Error |\n');
-fprintf('|------------------|--------------------|---------------------|\n');
+fprintf('|     Time [hr]    | Pressure Error | Darcy Flux Error | Strain Error | Displacement Error | Residual Error |\n');
+fprintf('|------------------|----------------|------------------|--------------|--------------------|----------------|\n');
+
 for i = 1:length(times_hr)
     % L2 error for pressure
     rel_err_p = norm(p_numerical(:,i) - p_analytical(:,i)) / norm(p_analytical(:,i));
-    
+
     % L2 error for flux
     rel_err_q = norm(q_numerical(:,i) - q_analytical(:,i)) / norm(q_analytical(:,i));
-    
-    % Print in table format
-    fprintf('| %16.2f | %18.6e | %19.6e |\n', times_hr(i), rel_err_p, rel_err_q);
+
+    % L2 error for strain
+    rel_err_e = norm(e_numerical(:,i) - e_analytical(:,i)) / norm(e_analytical(:,i));
+
+    % L2 error for displacement
+    rel_err_u = norm(u_numerical(:,i) - u_analytical(:,i)) / norm(u_analytical(:,i));
+
+    % L2 error for mass conservation residual
+    if i == 1
+        rel_err_r = NaN;  % Not defined for first time step
+    else
+        rel_err_r = norm(residual(:,i) - residual_analytical(:,i)) / norm(residual_analytical(:,i));
+    end
+
+    % Print in formatted row
+    if i == 1
+        fprintf('| %16.2f | %14.6e | %16.6e | %12.6e | %18.6e | %14s |\n', ...
+            times_hr(i), rel_err_p, rel_err_q, rel_err_e, rel_err_u, 'N/A');
+    else
+        fprintf('| %16.2f | %14.6e | %16.6e | %12.6e | %18.6e | %14.6e |\n', ...
+            times_hr(i), rel_err_p, rel_err_q, rel_err_e, rel_err_u, rel_err_r);
+    end
 end