Need help to create this equation

Question

Kevin Isakayoga 2021-1-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/722788-need-help-to-create-this-equation

评论： Bjorn Gustavsson 2021-2-18

采纳的回答： Walter Roberson

Hi everyone, I am a little bit confused about this

Could someone helping me to create this? It seems that to many partial derivatives, since I am not good at it.

*you can assume the parameter by yourself or let them be in coefficient

Hope someone could help me to understand clearly how to perform this equation. Thankyou very much!

2 个评论
显示无隐藏无

Kevin Isakayoga 2021-2-17

can somebody help me?

Thanks

Kevin Isakayoga 2021-2-18

编辑：Kevin Isakayoga 2021-2-18

@darova @Bjorn Gustavsson

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2021-2-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/722788-need-help-to-create-this-equation#answer_625567

在 MATLAB Online 中打开

syms rho_eff_hcp

syms p_w(x,t)

syms w(p_w_)

syms w_r(t)

dw_dp_w = diff(w(p_w_),p_w_)

dw_dp_w =

dp_w_dt = diff(p_w(x,t),t);

LHS = 1/rho_eff_hcp * dw_dp_w * dp_w_dt;

dp_w_dx = diff(p_w(x,t),x)

dp_w_dx =

syms K

RHS_temp1 = K*w(p_w_)*dp_w_dx

RHS_temp1 =

RHS_temp2 = diff(RHS_temp1,x)

RHS_temp2 =

dw_r_dt = diff(w_r(t),t);

RHS = -RHS_temp2 + dw_r_dt;

eqn = LHS == RHS

eqn =

Now you are faced with the problem that w should be a function of p_w instead of p_w_ but p_w is a function. You are trying to take the derivative of a function with respect to a function, which requires Calculus of Variations instead of regular calculas.

You also have three functions but only one equation. There is no hope that you would be able to resolve this symbolically, even if the equations were ODE instead of PDE.

9 个评论
显示 7更早的评论隐藏 7更早的评论

Kevin Isakayoga 2021-2-18

编辑：Kevin Isakayoga 2021-2-18

在 MATLAB Online 中打开

clear; clc;
days=36500;
var(1)=0.001;
timeinit=0.001;
counter=2;
while (timeinit<days)
    timeinit=timeinit*1.1;
    var(counter)=timeinit;
    counter=counter+1;
end
timeh=zeros(1,length(var)+1);
timeh1=zeros(1,length(var)+2);
timeh(1)=0; %should be changed
for m=1:length(var)
timeh(m+1)=var(m)*24;
timeh1(m+1)=var(m)*24;
end
timeday=timeh/24;
RH1=linspace(0.005,1,length(timeday));
vmbet=2.72*10^-4;
Cbet=18;
kbet=0.83;
casi=1.8;                    %atomic ratio
%sh20=202.6;                 %m2/g
wc=0.5;         %water to cement ratio
T=20;           %celcius
Tcal=T+273.15;
rhoeff=2.04;                 %density of HCP at 105C g/cm3
sh20=234-297*log(casi);      %m2/g
Sad=sh20*rhoeff*exp(-0.017*(Tcal-293));
cm=10;                    % Prediction length
h=0.1;                    % the interval of distance (cm)
M=(cm)*(1/h);             % Cut-off line of x-axis (days)
E=0.00001;
for mm=1:length(RH1)
    if (RH1(mm)<=0.98)
        dwrht(mm)=(Cbet*RH1(mm)*Sad*kbet^2*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)^2) - (Cbet*Sad*kbet*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)) + (Cbet*RH1(mm)*Sad*kbet^2*vmbet*(Cbet - 1))/((RH1(mm)*kbet*(Cbet - 1) + 1)^2*(RH1(mm)*kbet - 1));
    else
        dwrht(mm)=(Cbet*RH1(mm)*Sad*kbet^2*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)^2) - (Cbet*Sad*kbet*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)) + (Cbet*RH1(mm)*Sad*kbet^2*vmbet*(Cbet - 1))/((RH1(mm)*kbet*(Cbet - 1) + 1)^2*(RH1(mm)*kbet - 1));%50*w0293 - 50*w98ad;
    end
end
ah=1.05-3.8*wc+3.56*(wc)^2;
bh=-14.4+50.4*wc-41.8*(wc)^2;
gamh=31.3-136*wc+162*(wc)^2;
Dh=ah+bh.*(1-2.^(-10.^(gamh*(RH1-1))));
RH=0.766; %relative humadity
RHH=zeros(M+1,length(var));
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% boundary condition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for J=1:length(var)+1
    RHH(:,1)=1;
    if (J==1)
        RHH(1,J)=RH;
    else
        RHH(1,J)=RH;
    end
    RHH(M+1,J)=RH;                %it's a boundary condition
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Crank¡VNicolson method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 r=1/(2*h^2);
for J=1:length(var)
    for I=2:M
        RHH(I,J+1)=RHH(I,J);
    end
    for L=1:200
        LL=L;                   %number of times of calculation
        for I=2:M
            UF(I,J+1)=RHH(I,J+1);
            RHH(I,J+1)=(((r*(timeh(J+1)-timeh(J))/(24))*(1/dwrht(J))*Dh(J)*((RHH(I+1,J+1)-2*RHH(I,J+1)+RHH(I-1,J+1))+(RHH(I+1,J)-2*RHH(I,J)+RHH(I-1,J))+0))+RHH(I,J));
        end
        S=0;
        SO=0;
        for I=2:M
            SO=SO+abs(RHH(I,J+1));
            S=S+abs(RHH(I,J+1)-UF(I,J+1));
        end
        S=S/SO;
        if(S<=E)
            break;
        end
    end
end

maybe this is the code, but my equation didnt show a good results based on the following equation

and I want to get results like I mentioned below

and this is the boundary condition

Kevin Isakayoga 2021-2-18

I'm pretty sure that all of my coefficient such as dh/dw and Dh are correct.However, I just confuse my calculation and maybe my boundary condition was wrong.

Bjorn Gustavsson 2021-2-18

在 MATLAB Online 中打开

No, it doesn't look right. When you write your C-N function you will get to a stage where you have to multiply the RHS with the inverse of the LHS-matrix. You have neither an inv nor a \ in your C-N section. Take a pen and paper, sit down and write out the expressions for the LSH and RHS matrices and then implement those expressions. You will get a step looking something lke this:

I_next = invMLHS*([Mrhs]*I_prev + Q_t_p_half);

HTH

请先登录，再进行评论。