[ff3d-users] Could you help me confirm a time-dependent PDE?

[Robert Li]

Dear Stephane Del Pino,

  I tested ff3d with a realistic thermal PDE. However,
I found that the result of ff3d was different from
that of SVHEAT and TEMP/W, which are also FEM solution
software. But the other two have close result, or
almost same.

  I am not sure about how to write PDE in input file.
For example, It seems that -div((f2)*grad(T)) = 0; and
-div((f2)*grad(T)) = 0; are different. But I don't
know why?

 Could you help me confirm the following
time-dependent PDE?

1.Real Heat flow Equation
dx(Kx*dx(T)) + dy(Ky*dy(T)) +
dz(Kz*dz(T))=(Cp+Lf*Theat*Theatu1)*(dt(T))
where, T=Temperature,
dt=time

2.I expressed it in ff3d input file as follows:

T-div((f2)*grad(T)) = Tn;

where,
Tn = previously-computed temperature
dt = time step 1 day
Lf=Latent heat of fusion of water
K=conductivity 
Cp=volumetric heat capacity
Theatu1=slope of unfrozen water content
Theat=volumetric water content  
function        factor=((dt)/((Cp)+(Lf)*(Theat)*(Theatu1)));
function        f2=(K*factor);

Is it correct?

Thanks!

Robert
  


I used some selfdefind functions in the following
inputfile for conductivity, etc.

>>>>>>>>>>>input file>>>>>>>>

//load initial parameters
loadpara(".\\Param\\Param.ini");

mesh M = read(medit, "test.1.mesh");

femfunction     Tn(M) = initemperature();                               //T0; 
when
t=0; T
function                Told = 0;                       //T0; when t=0; T
double                  i = 1;                          //loop times
double                  dt = (1);                       //time step 1 day
double                  Lf=334000.0;                                    
//Latent heat of fusion of
water
function        K=(conductivity(Tn,-1) );                       //conductivity
in z direction
function        Cp=(vheatcapacity(Tn,-1));                      //volumetric
heat capacity
function        Theatu1=(slopeofunfwc(Tn, Told, -1));  
//use new slope method
function        Theat=theat();  


function        factor=((dt)/((Cp)+(Lf)*(Theat)*(Theatu1)));
function        f2=(K*factor);

//Get initial value T0, when time t=0;
solve(T) in  M cg(epsilon=1E-6){
         pde(T)

                 -dx((K)*dx(T))-dy((K)*dy(T))-dz((K)*dz(T))  = 0;

                T = -10 on  2;
                T = 2   on  1;

}
//Save result 
save (medit , ".\\temp\\Model9xxxx\\model.000" , T, M,
dos ) ;
save (medit , ".\\temp\\Model9xxxx\\model.000",  M,
dos ) ;

Tn= T;

//solve
do {
cout<<"I-----------------"<<i<<"\n";
         solve(T) in  M cg(epsilon=1E-6){
                 pde(T)
                        T-div((f2)*grad(T)) = Tn;
 
                        T = boundarytemperature(i) * correctvalue(Tn) on 2;
                        T = 2 on 1;

         }

        //Save result 
        save (medit , ".\\temp\\Model9xxxx\\model.00".i , T,
M, dos ) ;
        save (medit , ".\\temp\\Model9xxxx\\model.00".i,  M,
dos ) ;

        i=i+1;

        // correct value here
        Told = Tn;
    Tn = T;

} while (i<=365);

__________________________________
Do you Yahoo!?
Yahoo! Mail SpamGuard - Read only the mail you want.
http://antispam.yahoo.com/tools

//load initial parameters
loadpara(".\\Param\\Param.ini");

mesh M = read(medit, "test.1.mesh");

femfunction     Tn(M) = initemperature();                               //T0; 
when t=0; T
function                Told = 0;                       //T0; when t=0; T
double                  i = 1;                          //loop times
double                  dt = (1);                       //time step 1 day
double                  Lf=334000.0;                                    
//Latent heat of fusion of water
function        K=(conductivity(Tn,-1) );                       //conductivity 
in z direction
function        Cp=(vheatcapacity(Tn,-1));                      //volumetric 
heat capacity
function        Theatu1=(slopeofunfwc(Tn, Told, -1));   //use new slope method
function        Theat=theat();  


function        factor=((dt)/((Cp)+(Lf)*(Theat)*(Theatu1)));
function        f2=(K*factor);

//Get initial value T0, when time t=0;
solve(T) in  M cg(epsilon=1E-6){
         pde(T)

                 -dx((K)*dx(T))-dy((K)*dy(T))-dz((K)*dz(T))  = 0;

                T = -10 on  2;
                T = 2   on  1;

}
//Save result 
save (medit , ".\\temp\\Model9xxxx\\model.000" , T, M, dos ) ;
save (medit , ".\\temp\\Model9xxxx\\model.000",  M, dos ) ;

Tn= T;

//solve
do {
cout<<"I-----------------"<<i<<"\n";
         solve(T) in  M cg(epsilon=1E-6){
                 pde(T)
                        T-div((f2)*grad(T)) = Tn;
 
                        T = boundarytemperature(i) * correctvalue(Tn) on 2;
                        T = 2 on 1;

         }

        //Save result 
        save (medit , ".\\temp\\Model9xxxx\\model.00".i , T, M, dos ) ;
        save (medit , ".\\temp\\Model9xxxx\\model.00".i,  M, dos ) ;

        i=i+1;

        // correct value here
        Told = Tn;
    Tn = T;

} while (i<=365);