How to solve a second-order BVP with multiple boundary conditions?

Hello everyone,

I am solving a linear, second-order, nonhomogeneous differential equation of the form

d^2u/dr^2 + (2*m*r+n)/(m*r^2+n*r) * du/dr + (- nu/r^2 - 1/r^2 + nu/r * (2*m*r+n)/(m*r^2+n*r)) * u = alpha*(1+nu)*(dT/dr + ((2*m*r+n)/(m*r^2+n*r) -1) * T) – (rho*r*omega^2*(1-nu^2))/E

as a boundary value problem. The thermal expansion is disregarded for now. The differential equation describes the displacement of a rotating disc, u, of varying thickness, t, with respect to the radius, r. The thickness is varied with respect to the radius and follows a linear function of the type

t = m*r + n

The boundary conditions for this problem are:

at r=r_inner (or left boundary of the domain): 0 = du/dr + nu*u/r – (1+nu)*alpha*T

and

at r=r_outer (or right boundary of the domain): 0 = du/dr + nu*u/r – sigma_rim*(1-nu^2)/E – (1+nu)*alpha*T

I was able to set this up in matlab using the bvp5c solver.

However, now I would like to divide the differential equation into multiple regions (here 5), where the thickness remains either constant or changes linearly according to the above equation. This means the values of m and n, respectively, change with respect to the radius.

The main code is as follows:

%% Initialisation
% disc geometric key parameters
% number of disc sections, in -
nSections = 5;
% radius at stations 1 through 6, in m
rSections = [0.05 0.1 0.15 0.3 0.35 0.4];
% thickness at stations 1 through 6, in m
t = [0.07 0.07 0.01 0.01 0.04 0.04];
% slope of the linear function for the disc thickness
m = (t(2:end)-t(1:end-1))./(rSections(2:end)-rSections(1:end-1));
% constant of the linear function for the disc thickness
n = t(1:end-1) - m .* rSections(1:end-1);

% operational parameters
% disc rotational speed, in 1/s
omega = 5000 * 2*pi/60;
% temperature at disc bore, in K
T_bore = 125 + 273.15;
% temperature at disc rim, in K
T_rim = 550 + 273.15;
% rim load, in N/m^2 (Pa)
sigma_rim = 0;

% material properties
% Young modulus of elasticity, in N/m^2 (Pa)
E = 209e9;
% material density, in kg/m^3
rho = 7990; 
% Poisson's ratio, in -
nu = 0.3;
% coefficient of thermal expansion, in 1/K
alpha = 0;



%% ODE setup and solver
% specify number of regions
region = 1:nSections;

% set initial guess
rmesh = [];
delta_r = 0.01;
for i = 1:nSections
    rmesh = [rmesh rSections(i):delta_r:rSections(i+1)];
end
uinit = [0.001, 0.001];
init = bvpinit(rmesh, uinit);
% solver of the boundary value problem 
sol = bvp5c(...
    @(r,u,region) bvpfcn(r, u, region, m, n, nu, T_bore, T_rim, rSections, E, rho, omega, alpha),...
    @(uL,uR) bcfcn(uL, uR, rSections, nu, alpha, T_bore, T_rim, sigma_rim, E), ...
    init);

with the definition of the differential equation:

function [dudr] = bvpfcn(r, u, region, m, n, nu, T_bore, T_rim, rSections, E, rho, omega, alpha)
 
%% Initialisation
r_bore = rSections(1);
r_rim = rSections(end);
 
%% definition of constants
switch region
    case 1  % inner rim, u in [r1 r2]
        m_var = m(1);
        n_var = n(1);
    case 2  % inner web, u in [r2 r3]
        m_var = m(2);
        n_var = n(2);
    case 3  % outer web, u in [r3 r4]
        m_var = m(3);
        n_var = n(3);
    case 4  % shoulder, u in [r4 r5]
        m_var = m(4);
        n_var = n(4);
    case 5  % outer rim, u in [r5 r6]
        m_var = m(5);
        n_var = n(5);
end
% constant 1
C1 = ((2*m_var*r + n_var)/(m_var*r^2 + n_var*r));
% constant 2
C2 = (-nu/r^2 - 1/r^2 + nu/r * (2*m_var*r+n_var)/(m_var*r^2 + n_var*r));
% temperature distribution as a function of the radius
T = T_bore + (T_rim - T_bore) / (log(r_rim/r_bore)) * log(r/r_bore);
dTdr = (T_rim - T_bore) / (log(r_rim/r_bore)) * 1/r;
% constant 3
C3 = alpha*(1 + nu)*(dTdr + ((2*m_var*r+n_var)/(m_var*r^2+n_var*r) - 1) * T) - rho*omega^2*(1-nu^2)*r/E;
 
%% differential equation 
dudr = [u(2); -C1*u(2)-C2*u(1)+C3];
end

and the boundary condition function:

function [res] = bcfcn(uL, uR, rSections, nu, alpha, T_bore, T_rim, sigma_rim, E)
 
r_bore = rSections(1);
r_rim = rSections(end);

% bore boundary condition (radial stress = 0)
u_left = uL(2,1) + nu*uL(1,1)/r_bore - (1+nu)*alpha*T_bore;
 
% rim boundary condition (radial stress = rim stress)
u_right = uR(2,end) + nu*uR(1,end)/r_rim - sigma_rim * (1-nu^2)/E - (1+nu)*alpha*T_rim;

res = [u_left   
    uR(1,1) - uL(1,2) % compatibility equations at interface: u(i) = u(i+1)
    uR(1,2) - uL(1,3)
    uR(1,3) - uL(1,4)
    uR(1,4) - uL(1,5)
    u_right];
 
end

I get the following error message using the above code:

Error using bvparguments (line 111)
Error in calling BVP5C(ODEFUN,BCFUN,SOLINIT):
The boundary condition function BCFUN should return a column vector of length 10.

I am using 5 sections, meaning I need 2 boundary conditions for the boundaries of the domain and 4 compatibility equations for the interfaces of adjacent segments (displacement, u, being the same of the two differential equations at the interface). This would total 6 boundary conditions. Unfortunately, I wasn't able to solve the problem with the description of how to solve a bvp with multiple boundary conditions (Solve BVP with Multiple Boundary Conditions).

So my questions are:

How do I change the regions' definition and bvp function setup so that only 6 boundary conditions are required?

How is the indexing of the boundaries in the boundary condition function done so that the right interface k is selected as well as the right variable from the system of first-order differentials (u1 = u or u2 = du/dr)?

Any help is greatly appreciated!

Answer 1

[Torsten]

How do I change the regions' definition and bvp function setup so that only 6 boundary conditions are required?

For every equation, you need one interface condition. Further, you need a total of 2 boundary conditions.

Thus the vector you have to return in "bcfun" is of size 2 + 2 + 2 + 2 (interfaces) + 2 (boundaries) = 10.

How is the indexing of the boundaries in the boundary condition function done so that the right interface k is selected as well as the right variable from the system of first-order differentials (u1 = u or u2 = du/dr)?

From the documentation:

uL(i,k) is the solution variable u(i) at the left boundary point of region k.

uR(i,k) is the solution variable u(i) at the right boundary point of region k.

[Lucas]

Thank you very much for your answer, @Torsten!

I've got the code running with the four additional boundary conditions: uR(2,k) = uL(2,k+1).

I reckon I just needed clarification on the indices and the conditions required. So, cheers for that.

[Torsten]

If the transmission conditions between the different regions are - as in your case - continuity of the function and its first derivative - you wouldn't have needed this complicated approach. Taking the radii of the different sections as grid points and an if statement concerning your radial coordinate in function "bvpfcn" would have been sufficient, I guess.

[Lucas]

That does sound much simpler. I'll give it a shot. Thanks for the advice, @Torsten!