Using bvp4c for coupled second order ODEs. Issue in calling bvp4c.
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Isabelle Atkinson
le 23 Jan 2024
Réponse apportée : Torsten
le 23 Jan 2024
I have a system of two 2nd order ODEs. The set up of the ODE system is in the function bvp_function. I get the following error messages when I run the code:
Error using bvparguments
Error in calling BVP4C(ODEFUN,BCFUN,SOLINIT):
The derivative function ODEFUN should return a column vector of length 100.
Error in bvp4c (line 122)
bvparguments(solver_name,ode,bc,solinit,options,varargin);
Error in NewPoiseuilleBaseFlow (line 38)
sol = bvp4c(@(y,X) bvp_function(y,X,b), @(Xa,Xb) boundary_conditons_fcn(Xa,Xb), guess);
Here is the code. Anyone know how to solve this problem?
% Laminar base flow for Poiseuille flow
% using Deka 2023
clear
clc
%% Setting variables
N = 100; % number of collocation points
Mach = 2.0; % Mach number
rho0 = 1.225; % density [kg/m^3]
h = 1; % height of half channel (full channel height = 2h) [m]
gamma = 1.4;
R = 287;
T0 = 288.5; % temperature (and wall temperature) (15 degC) [K]
mu0 = 1.802e-5; % dynamic viscosity [kg m^-1 s^-1]
nu0 = mu0/rho0; % kinematic viscosity [m^2/s]
kappa0 = 0.02476; % thermal conductivity [W/mK]
u0 = Mach*sqrt(gamma*R*T0); % speed [m/s]
Re = rho0*u0*2*h/mu0; % reynolds number
Pr = gamma*R*mu0/((gamma-1)*kappa0); % Prandtl number
b = Pr*(gamma-1)*Mach^2; % constant for setting up ode system
y = linspace(-1,1,N)'; % y vector
T = linspace(T0,T0,N)'; % T vector
u = linspace(1,1,N)'; % initial u vector
%% Solving the x momentum and energy equations using ode45
% initial guess for the solver
uinit = incompVelProfile(N,Mach)';
guess = bvpinit(y, uinit);
%% Solving the system of ODEs using bvp4c
sol = bvp4c(@(y,X) bvp_fn(y,X,b), @(Xa,Xb) BCs(Xa,Xb), guess);
%% Function setting up matrix system of ODEs
function dxdy= bvp_fn(y,X,b)
% X(1)=u, X(2)=u', X(3)=T, X(4)=T'
dxdy(1) = X(2);
dxdy(3) = X(4);
dxdy(4) = -0.5*X(3)*(X(4))^2-b*(X(2))^2;
dxdy(2) = -2*(X(3))^(-3/2)-0.5*(X(3))^(-1)*X(4)*X(2);
dxdy = dxdy(:);
end
%% Function for boundary conditions
function res = BCs(Xa,Xb)
% u=0, T=1 at y=-1 and u=0, T=1 at y=1.
res = [Xa(1)-0;
Xb(1)-0
Xa(3)-1
Xb(3)-1];
end
%% Function for incompressible velocity profile
function [u] = incompVelProfile(N,Mach)
% N = number of collocation points
% Mach = desired Mach number
% returns the incompressible velocity profile for channel flow
% setting variables
rho = 1.225; % density of air [kg/m-3]
nu = 1.45e-5; % kinematic viscosity [m^2/s]
T = 288; % temperature of air, 15 degC [K]
gamma = 1.4; % ratio of specific heat capacities
R = 287; % universal gas constant for air [J/(kg.K)]
% setting domain
y = linspace(-1,1,N); % y points
% working out the equivalent pressure gradient
c = sqrt(gamma*R*T); % speed of sound [m/s]
umax = Mach*c; % max velocity (centre of channel) [m\s]
dPdx = -2*umax*rho*nu; % pressure gradient required for centreline vel [Pa/m]
% obtaining velocity profile
u = zeros(N); % initialising vector
u = -0.5.*1./(rho.*nu).*dPdx.*(1-y.^2);
% have to normalise it
u = u./umax;
end
0 commentaires
Réponse acceptée
Torsten
le 23 Jan 2024
guess = bvpinit(y, uinit);
y must be the spatial mesh,
uinit must either be a function handle that returns a 4x1 vector of initial values for the 4 solution components of X with an y value as input
or
a constant 4x1 vector that is assigned as initial value for the 4 solution components for each y value.
In your case, uinit seems to be an NxN matrix.
Read here about what uinit can be set to:
0 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Boundary Value Problems dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!