gistfile1.cpp

void QThermalConductivity::ev_solid_foil()
{
    double l = thickness[2];
    double kapa_s;
    double H = thickness[3];    
    
    T_temp[N[1] - 1] = T[N[1] - 1] + 1.0/chcapacity[2]/ro[2]/H*(-2*kapa[1]*kapa[2]/(kapa[1]*H + kapa[2]*dx)*(T[N[1] - 1] - T[N[1] - 2]) +2/(H+l)*kapa_s*(T[N[1]] - T[N[1] - 1]))*dt;
    
    bool t1 = (!flag && (f <= 0) && (T[N[1]] > T_evt) );
    bool t2 = (!flag && (f >= 1) && (T[N[1]] < T_evt) );
    
    if (t1)
        solder_mode = 1;
    if (t2)
        solder_mode = 2;
    
    if (t1 || t2)
        flag = 1;
    
    if (flag) {
        
        kapa_s = kapa[2]*(H+ l)/H*((kapa_l*f + kapa[3]*(1 - f))/l/(kapa[2]/H + (kapa_l*f + kapa[3]*(1 - f))/l));
        
        T_temp[N[1] - 1] = T[N[1] - 1] + 1.0/chcapacity[2]/ro[2]/H*(-2*kapa[1]*kapa[2]/(kapa[1]*H + kapa[2]*dx)*(T[N[1] - 1] - T[N[1] - 2]) +2/(H+l)*kapa_s*(T[N[1]] - T[N[1] - 1]))*dt;
        
        T_temp[N[1]] = T_evt;
        
        f += dt/(lambda*ro[3]*l)*(mu[1]*(T[N[2] + 1] - T[N[2]]) - kapa_s*(T[N[2]] - T[N[2] - 1])*2/(H + l));
        
        if ((f < 0) || (f > 1))
            flag = 0;
        
    } else {
        if (T[N[1]] < T_evt) {
            
            T_temp[N[1] - 1] = T[N[1] - 1] + 1.0/chcapacity[2]/ro[2]/H*(-2*kapa[1]*kapa[2]/(kapa[1]*H + kapa[2]*dx)*(T[N[1] - 1] - T[N[1] - 2]) +2*kapa[2]*kapa[3]/(kapa[2]*l + kapa[3]*H)*(T[N[1]] - T[N[1] - 1]))*dt;
            
            T_temp[N[1]] = T[N[1]] + 1.0/chcapacity[3]/ro[3]/l*(-2*kapa[2]*kapa[3]/(kapa[2]*l +kapa [3]*H)*(T[N[1]] - T[N[1] - 1]) +
                    mu[1]*(T[N[1] + 1] - T[N[1]]))*dt;
        } else {
            T_temp[N[1] - 1] = T[N[1] - 1] + 1.0/chcapacity[2]/ro[2]/H*(-2*kapa[1]*kapa[2]/(kapa[1]*H + kapa[2]*dx)*(T[N[1] - 1] - T[N[1] - 2]) +  2*kapa[2]*kapa_l/(kapa[2]*l + kapa_l*H)*(T[N[1]] - T[N[1] - 1]))*dt;
            
            T_temp[N[1]] = T[N[1]] + 1.0/chcapacity[3]/ro[3]/l*(-2*kapa[2]*kapa_l/(kapa[2]*l +kapa_l*H)*(T[N[1]] - T[N[1] - 1]) +
                    mu[1]*(T[N[1] + 1] - T[N[1]]))*dt;            
        }        
    }
    
    T_temp[N[1] + 1] = T[N[1] + 1] + 1.0/chcapacity[4]/ro[4]/dx*(-mu[1]*(T[N[1] + 1] - T[N[1]]) + kapa[4]/dx*(T[N[1] + 2] - T[N[1] + 1]))*dt;
    
    if (flag_foil) {
        ddy = sqrt(D0*dg/(k*T_temp[N[1] - 1]*c*(1-c))*qExp(-(Q/(k*T[N[1] - 1])))*dt);
        
        if ((dy_f + ddy) > delta/2) {
            ddy = delta/2 - dy_f;
            dy_f = delta/2;
            flag_foil = false;
        } else {
            dy_f += ddy;
        }
        
        T_temp[N[1] - 1] = T_temp[N[1] - 1] + 2.0/delta*dg/(3*k)*ddy;
    }
}