Ana Sayfa / MATLAB / Python / Problem 8

Kabuk Kırılması — Sonlu Fark mı, Kesin Çözüm mü?

1 m kalınlığında duvar, kısa ve şiddetli darbe (P_r = 5 MPa, t_d = 0,1 ms). Basınç dalgası serbest arka yüzden çekme olarak yansır. Kuram açıktır: kalın duvar ve kısa darbede yansıyan çekme, gelen tepeye eşittir — σ_t/σ₀ = 1.

Kod çalışıyor ama yanlış. Açık sonlu fark şeması keskin darbe cephesinde sahte salınım üretir ve σ_t/σ₀ = 1,211 verir — doğrusal elastik bir dalgada fiziksel olarak imkânsız bir değer: yansıyan çekme gelen tepeyi aşamaz. Kesin d'Alembert çözümü ise 1,000 verir. Bu, bir patlama analiz paketinin kabuk kırılması rutininde bulunan gerçek bir hatadır; site aracı artık kesin çözümü kullanır.

İkinci ders — denetimin kendisi de yanlış olabilir. Bu hatayı ararken önce "ince duvarda çekme σ₀·2t/λ ile sönümlenir" diye analitik bir kontrol yazıldı. O formül yalnızca ilk yansımayı sayıyordu; sönümsüz bir levhada çoklu yansımalar çekmeyi gelen tepeye yaklaştırır. Kesin çözüm hem koddaki hatayı hem de denetimdeki hatayı birlikte ortaya çıkardı.

Kullanılan fiziksel sınır: "yansıyan çekme ≤ gelen tepe". Sayısal bir sonucu bir fiziksel sınırla karşılaştırmak, en ucuz ve en güçlü doğrulama araçlarından biridir.

# spall_fd_vs_exact.py -- Kabuk kirilmasi: acik sonlu fark semasi vs KESIN cozum
# Ayni problem iki yolla. Dogrusal elastik bir dalgada yansiyan cekme gelen
# tepeyi ASAMAZ; kalin duvar + kisa darbede kesin sonuc sigma_t/sigma_0 = 1.
import numpy as np

E, nu, rho = 30e9, 0.2, 2400.0
t  = 1.0                     # m, duvar kalinligi
Pr, ir = 5e6, 250.0          # Pa, Pa.s  -> td = 0.1 ms (kisa darbe)
td = 2*ir/Pr
M  = E*(1-nu)/((1+nu)*(1-2*nu))      # kisitli modul
c  = np.sqrt(M/rho); T = t/c
load = lambda s: np.where((s >= 0) & (s <= td), Pr*(1 - s/td), 0.0)

# --- 1) Acik sonlu fark (orijinal spallCheck.m semasi) ---
nel = 80; he = t/nel; dt = 0.5*he/c
mN = np.full(nel+1, rho*he); mN[[0, -1]] /= 2
u = np.zeros(nel+1); v = np.zeros(nel+1); smin_fd = 0.0
for tt in np.arange(0, td + 4*T, dt):
    s = M*(u[:-1] - u[1:])/he            # basinc +
    f = np.zeros(nel+1); f[0] = load(tt)
    f[:-1] -= s; f[1:] += s
    v += dt*f/mN; u += dt*v
    smin_fd = min(smin_fd, s.min())

# --- 2) KESIN cozum (d'Alembert):  sigma = f(tau - x/c) - f(tau + x/c - 2T),
#        f(tau) = p(tau) + f(tau - 2T)   ->   f(tau) = sum_k p(tau - 2kT)
def F(tau):
    out = np.zeros_like(tau)
    for kk in range(int(np.max(tau)/(2*T)) + 2):
        out += load(tau - 2*kk*T)
    return out
x = np.linspace(0, t, 81); smin_ex = 0.0
for tau in np.linspace(0, td + 6*T, 3001):
    sig = F(tau - x/c) - F(tau + x/c - 2*T)
    smin_ex = min(smin_ex, sig.min())

print(f"dalga hizi c = {c:.0f} m/s   td = {td*1e3:.3f} ms   T = t/c = {T*1e6:.1f} us")
print(f"Sonlu fark : sigma_t/sigma_0 = {-smin_fd/Pr:.3f}")
print(f"Kesin cozum: sigma_t/sigma_0 = {-smin_ex/Pr:.3f}   (kuram: 1.000)")
print(f"Sonlu fark asimi = %{(smin_fd/smin_ex - 1)*100:.1f}")

% spall_fd_vs_exact.m -- Kabuk kirilmasi: acik sonlu fark semasi vs KESIN cozum
% Dogrusal elastik dalgada yansiyan cekme gelen tepeyi ASAMAZ; kalin duvar +
% kisa darbede kesin sonuc sigma_t/sigma_0 = 1.
function spall_fd_vs_exact()
  E = 30e9; nu = 0.2; rho = 2400;
  t = 1.0;  Pr = 5e6; ir = 250;
  td = 2*ir/Pr;
  M  = E*(1-nu)/((1+nu)*(1-2*nu));
  c  = sqrt(M/rho); T = t/c;
  load = @(s) (s >= 0 & s <= td) .* Pr .* (1 - s/td);

  % --- 1) Acik sonlu fark (orijinal spallCheck.m semasi) ---
  nel = 80; he = t/nel; dt = 0.5*he/c;
  mN = rho*he*ones(nel+1,1); mN([1 end]) = mN([1 end])/2;
  u = zeros(nel+1,1); v = zeros(nel+1,1); smin_fd = 0;
  for tt = 0:dt:(td + 4*T)
    s = M*(u(1:end-1) - u(2:end))/he;
    f = zeros(nel+1,1); f(1) = load(tt);
    f(1:end-1) = f(1:end-1) - s;  f(2:end) = f(2:end) + s;
    v = v + dt*f./mN;  u = u + dt*v;
    smin_fd = min(smin_fd, min(s));
  end

  % --- 2) KESIN cozum (d'Alembert) ---
  x = linspace(0, t, 81); smin_ex = 0;
  for tau = linspace(0, td + 6*T, 3001)
    sig = F(tau - x/c, load, T) - F(tau + x/c - 2*T, load, T);
    smin_ex = min(smin_ex, min(sig));
  end
  fprintf('dalga hizi c = %.0f m/s   td = %.3f ms   T = t/c = %.1f us\n', c, td*1e3, T*1e6);
  fprintf('Sonlu fark : sigma_t/sigma_0 = %.3f\n', -smin_fd/Pr);
  fprintf('Kesin cozum: sigma_t/sigma_0 = %.3f   (kuram: 1.000)\n', -smin_ex/Pr);
  fprintf('Sonlu fark asimi = %%%.1f\n', (smin_fd/smin_ex - 1)*100);
end

function out = F(tau, load, T)
  out = zeros(size(tau));
  for kk = 0:(floor(max(tau)/(2*T)) + 1)
    out = out + load(tau - 2*kk*T);
  end
end