I have 2 variables as a results from the NDSolve[] which are ua and uw.
This is how I do it, and I got the result for ua and uw just fine and I can plot them.
eq1 = Derivative[2, 0][uw][z, t]*cwv -
Derivative[0, 1][ua][z, t]*cw == Derivative[0, 1][uw][z, t];
eq2 = Derivative[2, 0][ua][z, t]*cav -
Derivative[0, 1][uw][z, t]*ca == Derivative[0, 1][ua][z, t];
ic1 = {uw[zmax, t] == ui + p0, uw[z, tmin] == ui};
ic2 = {ua[z, tmin] == uatm};
sol2 = NDSolve[Flatten[{eq1, eq2, ic1, ic2}], {uw, ua}, {z, zmin, zmax}, {t, tmin,
tmax}, MaxSteps -> Infinity, MaxStepSize -> dt/step, InterpolationOrder -> All];`
However, I want to subtract the two variables such that:
result = (ua-101)-uw and I failed.
I use this as a code to subtract
result = ua[z, t] /. sol2 - 101 - uw[z, t] /. sol2;
Can someone please tell me the reason why it is failed and what should I do to fix this?
The complete code is given as follows:
kw = 1*10^(-10);(*m/s*)
ka = 8*10^(-13);(*m/s*)
gw = 10;(*kN/m^3*)
m2w = 2.41*10^(-4);(*coefficient of compressibility of water phase with respect to matric suction*)
m1w = 0.8*10^(-4);(*coefficient of compressibility of water phase with respect to axial stress*)
m2a = 1.11*10^(-4);(*coefficient of compressibility with air phase with respect to matric suction*)
m1a = 0.37*10^(-4);(*coefficient of compressibility of air phase with respect to axial stress*)
p0 = 200;(*applied pressure -kPa*)
tmin = 0;(*minimum or time when the pressure is applied*)
tmax = 5000;(*latest inspected time*)
zmin = 0;(*top elevation*)
zmax = 0.02;(*bottom elevation*)
uatm = 101;(*kPa-atmospheric air pressure*)
step = 1000;(*The higher the step, the higher the accuracy*)
ui = -300;(*Initial matric suction*)
temp = 300;(*temperature in kelvin*)
wa = 28.97;(*kg/kmol*)
s = 0.7887;
n = 1.0696/(1 + 1.0696);
(*CALCULATION*)
r = 8.31432;(*Universal gas constant in j/(mol K*)
g = 9.8;(*gravity acceleration-m2/s*)
da = ka/g;
cwv = kw/(gw*m2w);
cw = ((1 - m2w/m1w)/(m2w/m1w));
ca = (m2a/m1a)/(1 - (m2a/m1a) - (1 - s)*n/(uatm*m1a));
cav = (da/(wa/(r*temp)*uatm))/(m1a*(1 - m2a/m1a) - (1 - s)*n);
dt = tmax - tmin;
(*Assuming that dua/dt is insignificant, Consolidation equation for unsaturated soil is duw/dt=cwv(d^2uw/dy^2)*)
eq = Derivative[2, 0][uw][z, t]*cwv == Derivative[0, 1][uw][z, t];
ic = {uw[zmax, t] == ui + p0, uw[z, tmin] == ui};
sol = NDSolve[Flatten[{eq, ic}], uw, {z, zmin, zmax}, {t, tmin, tmax},MaxSteps -> Infinity, MaxStepSize -> dt/step,InterpolationOrder -> All]
is = 300;(*ImageSize*)
Plot3D[uw[z, t] /. sol, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {z, t, uw[z, t]}, BoundaryStyle -> Thick, ImageSize -> 500]
Column[Manipulate[Plot[uw[z, t] /. sol, {z, zmin, zmax}, PlotRange -> All, AxesLabel -> {z, uw[z, t]}, ImageSize -> is], {t, tmin, tmax}], Manipulate[Plot[uw[z, t] /. sol, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {t, uw[z, t]}, ImageSize -> is], {z, zmin, zmax}], Manipulate[Plot[{uw[z, tmin] /. sol, uw[z, t1] /. sol, uw[z, t2] /. sol, uw[z, t3] /. sol, uw[z, t4] /. sol, uw[z, tmax] /. sol}, {z, zmin,zmax}, AxesLabel -> {z, uw[z, t]}, ImageSize -> is], {t1, tmin + dt/300, tmin + dt/51}, {t2, tmin + dt/50, tmin + dt/16}, {t3, tmin + dt/15, tmin + dt/3}, {t4, tmin + dt/3, tmin + dt/1}]]
(*Assuming that dua/dt is significant*)
eq1 = Derivative[2, 0][uw][z, t]*cwv -
Derivative[0, 1][ua][z, t]*cw == Derivative[0, 1][uw][z, t];
eq2 = Derivative[2, 0][ua][z, t]*cav -
Derivative[0, 1][uw][z, t]*ca == Derivative[0, 1][ua][z, t];
ic1 = {uw[zmax, t] == ui + p0, uw[z, tmin] == ui};
ic2 = {ua[z, tmin] == uatm};
sol2 = NDSolve[Flatten[{eq1, eq2, ic1, ic2}], {uw, ua}, {z, zmin, zmax}, {t,tmin,tmax}, MaxSteps -> Infinity, MaxStepSize -> dt/step, InterpolationOrder -> All];
(*Coupled equation version*)
Column[Plot3D[uw[z, t] /. sol2, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All,AxesLabel -> {z, t, uw[z, t]}, BoundaryStyle -> Thick, ImageSize -> 500],
Plot3D[ua[z, t] /. sol2, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {z, t, ua[z, t]}, BoundaryStyle -> Thick, ImageSize -> 500]]
Column[Manipulate[Plot[uw[z, t] /. sol2, {z, zmin, zmax}, PlotRange -> Full, AxesLabel -> {z, uw[z, t]}], {t, tmin, tmax}],Manipulate[Plot[uw[z, t] /. sol2, {t, tmin, tmax},PlotRange -> Full, AxesLabel -> {t, uw[z, t]}], {z, zmin,zmax}]]Column[Manipulate[Plot[ua[z, t] /. sol2, {z, zmin, zmax}, PlotRange ->Full,AxesLabel -> {z, ua[z, t]}], {t, tmin, tmax}],Manipulate[Plot[ua[z, t] /. sol2, {t, tmin, tmax}, PlotRange -> Full, AxesLabel -> {t, ua[z, t]}], {z, zmin, zmax}]]
Column[Manipulate[Plot[{uw[z, tmin] /. sol2, uw[z, t1] /. sol2, uw[z, t2] /. sol2,uw[z, t3] /. sol2, uw[z, t4] /. sol2, uw[z, tmax] /. sol2}, {z, zmin, zmax}, AxesLabel -> {z, uw[z, t]}, PlotRange -> Full, AxesLabel -> {z, uw[z, t]}], {t1, tmin + dt/1000, tmin + dt/51}, {t2, tmin + dt/50, tmin + dt/16}, {t3, tmin + dt/15,tmin + dt/3}, {t4, tmin + dt/3, tmin + dt/1}], Manipulate[Plot[{ua[z, tmin] /. sol2, ua[z, t1] /. sol2, ua[z, t2] /. sol2, ua[z, t3] /. sol2, ua[z, t4] /. sol2, ua[z, tmax] /. sol2}, {z, zmin, zmax}, AxesLabel -> {z, ua[z, t]}, PlotRange -> Full], {t1, tmin + dt/500, tmin+ dt/51}, {t2, tmin + dt/50, tmin + dt/16}, {t3, tmin + dt/15, tmin + dt/3}, {t4, tmin + dt/3, tmin + dt/1}]]
(*Comparison*)
Manipulate[Plot[{uw[z, t] /. sol, uw[z, t] /. sol2}, {t, tmin, tmax},PlotRange -> All, AxesLabel -> {t, uw[z, t]}], {z, zmin, zmax}];
(*Matric Suction*)
result = ua[z, t] /. sol2 - 101 - uw[z, t] /. sol2;
Plot3D[result, {z, zmin, zmax}, {t, tmin, tmax}, PlotRange -> All, AxesLabel -> {z, t, \[Psi][z, t]}, BoundaryStyle -> Thick, ImageSize -> 500] `
Regards,
