Cohete propulsado por agua

Procedimientos numéricos

Sistema de ecuaciones diferenciales

Resolver el sistema de ecuaciones diferenciales, una de primer orden y otra de segundo orden

$\begin{array}{l} v_{2}^{2} = f (h) = \frac{\frac{p_{0} (H - h_{0})}{H - h} + ρ g h - p_{a t}}{\frac{1}{2} ρ (1 - \frac{S_{2}^{2}}{S_{1}^{2}})} \\ \frac{d h}{d t} = - \frac{S_{2}}{S_{1}} \sqrt{f (h)} \\ \frac{d^{2} x}{d t^{2}} = \frac{ρ S_{2}}{ρ S_{1} h + m_{u}} f (h) - g \end{array}$

    void setNuevo(int valor, double carga, int fracRadio){ 
//tanto por ciento de agua
       h0=hDeposito*valor/100; //altura inicial del agua
       cargaUtil=carga;
       rOrificio=rDeposito/fracRadio;
       h=h0;
       t=0.0;
       x=0.0;
       v=0.0;
       bTermina=false;
 //tiempo hasta que se agota el combustible
       xMax=0.0;    //posición cuando se acaba el combustible
       vMax=0.0;     //velocidad cuando se acaba el combustible
       tMax=0.0;
//masa total inicial
    nVeces=0;
    opcion=0;
    xEmbolo=0.0;
    pInicial=presion=pAtm;
    nRTinicial=pAtm*Math.PI*rDeposito*rDeposito*(hDeposito-h0);
    den=500*(1.0-rOrificio*rOrificio*rOrificio*rOrificio
/(rDeposito*rDeposito*rDeposito*rDeposito));
  }
 void mover(){
    switch(opcion){
        case 0:  //llenar de aire
            xEmbolo+=dxEmbolo;  
//el volumen de la bomba es de 20 cm3 o 0.02 litros
            presion=(nRTinicial+pAtm*volBomba*((nVeces-1)+xEmbolo))
/(Math.PI*rDeposito*rDeposito*(hDeposito-h0));  //presión
            pInicial=presion;
           if(xEmbolo>1.0){
                //se para;
           }
           break;
        case 1: //vuelo
            rungeKutta();
            presion=pInicial*(hDeposito-h0)/(hDeposito-h);
            xMax=x;
            vMax=v;
            tMax=t;
            if((h<0.0)||(bTermina)){
               opcion=2;
            }
            break;
        case 2:  //movimiento bajo la acción de la gravedad
            x=xMax+vMax*(t-tMax)-4.9*(t-tMax)*(t-tMax);
            v=vMax-9.8*(t-tMax);
            if(v<0.0)  //se detiene el cohete
            t+=dt;
            break;
        default:
            break;
    }
  }
  double f(double y){    //aceleración del cohete
    double masa=1000*Math.PI*rDeposito*rDeposito*y+cargaUtil;
    double temp=0.0;
    temp=1000*Math.PI*rOrificio*rOrificio*velocidad2(y)/masa-9.8; 
//aceleración (empuje-peso)/masa
        if(temp<0) temp=0.0;
    return temp;
  }
  double g(double y){   //para la altura de agua
    double temp=0.0;
    temp=-rOrificio*rOrificio*Math.sqrt(velocidad2(y))/(rDeposito*rDeposito);
    return temp;
  }

 void rungeKutta(){
    double k1, k2, k3, k4;
    double l1, l2, l3, l4;
    double m1, m2, m3, m4;
//altura
    m1=dt*g(h);
    m2=dt*g(h+m1/2);
    m3=dt*g(h+m2/2);
    m4=dt*g(h+m3);
//velocidad
    l1=dt*f(h);
    l2=dt*f(h+m1/2);
    l3=dt*f(h+m2/2);
    l4=dt*f(h+m3);
//posición
    k1=dt*v;
    k2=dt*(v+l1/2);
    k3=dt*(v+l2/2);
    k4=dt*(v+l3);

    x+=(k1+2*k2+2*k3+k4)/6;
    v+=(l1+2*l2+2*l3+l4)/6;
    h+=(m1+2*m2+2*m3+m4)/6;   //altura del agua
    t+=dt;
 }
 //devuelve el valor de v2 al cuadrado
 double velocidad2(double y){
    double num=pInicial*(hDeposito-h0)/(hDeposito-y)+9.8*1000*y-pAtm;
    if(num<0.0){   //no sale toda el agua del depósito
        num=0.0;
        bTermina=true; 
    }
    return (num/den);
  }