Simulaciones en Matlab

[V_Dan]

Hola, estoy haciendo un simulador de la Ley de Stokes, pero no se como hacer que cuando llegue a determinada altura sepa que llego al fluido y ahí haga las ecuaciones de fluido
intente con un if, pero esto me daña el programa :S ¿Alguna sugerencia?
Gracias por la atención.
function varargout = Fluido(varargin)
% FLUIDO MATLAB code for Fluido.fig
% FLUIDO, by itself, creates a new FLUIDO or raises the existing
% singleton*.
%
% H = FLUIDO returns the handle to a new FLUIDO or the handle to
% the existing singleton*.
%
% FLUIDO('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in FLUIDO.M with the given input arguments.
%
% FLUIDO('Property','Value',...) creates a new FLUIDO or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before Fluido_OpeningFcn gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to Fluido_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE's Tools me nu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help Fluido

% Last Modified by GUIDE v2.5 31-Jul-2013 16:01:53

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @Fluido_OpeningFcn, ...
'gui_OutputFcn', @Fluido_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before Fluido is made visible.
function Fluido_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to Fluido (see VARARGIN)

% Choose default command line output for Fluido
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes Fluido wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = Fluido_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on selection change in fluido.
function fluido_Callback(hObject, eventdata, handles)
% hObject handle to fluido (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns fluido contents as cell array
% contents{get(hObject,'Value')} returns selected item from fluido


val = get(hObject, 'Value');
global de;
global vis;

switch val
case 2
% El usuario seleccionó la primera opción AGUA
de=1000;
vis=0.001;
case 3
% El usuario seleccionó la segunda opción ACEITE
de=920;
vis=0.03;
case 4
% El usuario seleccionó la 3 opción GLICERINA
de=1200;
vis=1.5;
end
% --- Executes during object creation, after setting all properties.
function fluido_CreateFcn(hObject, eventdata, handles)
% hObject handle to fluido (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end



function radio_Callback(hObject, eventdata, handles)
global r;
% hObject handle to radio (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of radio as text
% str2double(get(hObject,'String')) returns contents of radio as a double
r=(str2double(get(hObject,'String')));
% --- Executes during object creation, after setting all properties.
function radio_CreateFcn(hObject, eventdata, handles)
% hObject handle to radio (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end


% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

global r ;
r=str2double(get(handles.radio,'String'));
global de ;
global vis;

yo=15;
vyo=0;
g=-9.8;

tcaidali=(yo*2/g)^(1/2);


t=0:1/50:5:(((6)^2/7)+20);

%tMax = 0;
%deltaT = 0.50;
%t = tcaidali;



%%Ecuaciones caida libre
x=0*t;
vy = vyo + g*t;
y=yo+vyo*t+(1/2)*g*t.^2;

%%%Material--->Cobre
densidadcobre=8500;
%Entra en el fluido a una altura de cinco metros
if yo==5
%%Ecuaciones Fluidos

m=(densidadcobre*(4/3)*pi*(r.^2));
b=(6*pi*r*vis);
E=(de*(4/3)*pi*(r.^3)*g);

vl=((2*r.^2*g*(densidadcobre-de))/(9*vis));
v=vl*(1-exp((-b*t/m)));
x=vl*(t-(m/b)*(1-exp((-b/m)*t)));



for i=1:length(t) % para animar las grafica
pause(1/50) % tiempo para que se detenga tiempo igual al de arriba


%subplot(2,2,1)% Matriz 2x2 posición1
plot(handles.axes1,t(i),vy(i))
xlabel(handles.axes1,'tiempo t(s)')
ylabel(handles.axes1,'Velocidad en Y - Vy (m/s)')
axis(handles.axes1,[min(t) max(t) min(vy) max(vy)])
hold (handles.axes1,'on')
yo=yo+1;
plot(handles.axes2,x(i),y(i), 'o r')
xlabel(handles.axes2,'Posición x(m)')
ylabel(handles.axes2,'Posición y(m)')
axis(handles.axes1,axis,[-2 2 min(y) max(y)])
hold (handles.axes2,'on')

plot(handles.axes1,t(i),v(i))
xlabel(handles.axes1,'tiempo t(s)')
ylabel(handles.axes1,'Velocidad en Vl - Vl (m/s)')
axis(handles.axes1,axis,[min(t) max(t) min(v) max(v)])
hold (handles.axes1,'on')


end
end
%tMax = 3;
%deltaT = 0.01;
%t = 0;

%Bola caida libre
%while tiempoTranscurrido < tMax
%caídalibre
%x =0*deltaT;
% y = yo+Vyo*deltaT+(1/2)*g*deltaT.^2;

% vy = vy + g*deltaT;

%plot(handles.axes2, t, vy);
% xlabel(handles.axes2,'tiempo en t(s).')
% ylabel(handles.axes2,'Velocidad en Vy(m/s)')
% axis([0 tMax -5 5])
% hold (handles.axes2,'on')

%end

%v=v2*(1-exp((-b/m)*t));
%x=vt*(t-(m/b)*(1-exp((-b/m)*t)));
%vl=(2*r^2*g*(de-df))*9*vis;

%b=(6*Pi*r*vis);
%E=(df*(4/3)*Pi*(r^3)*g);
%m=(de*(4/3)*Pi*(r^2));

%plotear v vs tiempo


% --- Executes during object creation, after setting all properties.
%function axes1_CreateFcn(hObject, eventdata, handles)
% hObject handle to axes1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes1

[Almapa]

No se si a día de hoy ya lo habrás resuelto pero si no explicas un poco más en detalle dónde está exactamente tu problema no te puedo seguir muy bien. Por una parte, en la línea 177 (en mi archivo), donde calculas "m", si se supone que estás calculando la masa, el exponente de "r" para el cálculo del volumen debería ser 3.

Por otra parte, si quieres que utilice las ecuaciones en el fluido a partir de una altura de 5 unidades no puedes pones un simple "if==5" pues solo las aplicará en el punto de y=5 y no en alturas superiores ni inferiores.

Un saludo!

[Almapa]

Estuve haciendo una pequeña simulación y la fuerza de rozamiento obtenida por la expresión de Stokes sale como 4 órdenes de magnitud menos, por lo que prácticamente no afecta a la aceleración de bajada. Posiblemente el error esté en las unidades de algún dato.

Aquí mi simulación:

Código (matlab) [Seleccionar] Expandir

%Datos líquido
de=1000;
vis=1;
%Datos esfera
r=10;
densidadcobre=8500;
m=(densidadcobre*(4/3)*pi*(r.^3));
%Datos iniciales: altura inicial, velocidad inicial, aceleración y altura
%del líquido
yo=15;
vyo=0;
g=-9.8;
yent=5;

%%Ecuaciones caida libre
tent=(-vyo-sqrt(vyo^2-2*g*(yo-yent)))/g;
t=0:0.001:tent;
vy = vyo + g*t;
y=yo+vyo*t+(1/2)*g*t.^2;
y1=y(length(y));
v1=vy(length(vy));
while y1>0
   b=(6*pi*r*vis*v1);
   al=(g*m+b)/m;
   v1= vyo + al*(t(length(t))+0.01);
   y1=yo+vyo*(t(length(t))+0.01)+(1/2)*al*(t(length(t))+0.01)^2;
   y=[y,y1];
   vy=[vy,v1];
   t=[t,t(length(t))+0.01];
end
       


%%Ecuaciones Fluidos
subplot(2,1,1)
plot(t,y)
xlabel('tiempo t(s)')
ylabel('Velocidad en Y - Vy (m/s)')
subplot(2,1,2)
plot(t,vy)

Eso está muy poco optimizado y es un primer ejemplo solo para el caso en el que se tira desde una altitud superior a la entrada en el líquido.