Question: How can I use fsolve to interatively solve a problem that has units and uses the ThermoPhysicalData library?

I have been trying to use fsolve to iteratively solve a free convection heat transfer problem for a heated vertical plate using units and the ThermoPhysicalData library.  I have received some help and have managed to get rid of the eror messages but I still can't seem to get it to solve.

Here is the code


with(Units:-Standard); with(ThermophysicalData);

conv_VerticalPlate_proc := proc (T__side, T__amb, Area, L, varepsilon, q__i)
  local sigma, g, T__s, `T__∞`, T__f, k, rho, Cp, Pr, mu, nu, alpha, beta, Ra__L, Nus__L, h, q__h, q__r, q__total;
  g := 9.81*Unit('m'/'s'^2); sigma := 5.6703*Unit('W'/('m'^2*'K'^4))/10^8;
  if type(T__side, with_unit) then
    T__s := convert(T__side, temperature, kelvin)
  elif T__side = 0 then
    T__s := 273.15*Unit('K')
    T__s := T__side*Unit('K')
  end if;
  if type(T__amb, with_unit) then
    'T__∞` := convert(T__amb, temperature, kelvin)
  elif T__amb = 0 then
    'T__∞` := 273.15*Unit('K')
    `T__∞` := T__amb*Unit('K')
  end if;
  T__f := (1/2)*T__s+(1/2)*`T__∞`;
  k := Property(thermal_conductivity, temperature = T__f, pressure = Unit('atm'), air);
  rho := Property(density, temperature = T__f, pressure = Unit('atm'), air);
  Cp := Property(Cpmass, temperature = T__f, pressure = Unit('atm'), air);
  Pr := Property(prandtl, temperature = T__f, pressure = Unit('atm'), air);
  mu := Property(viscosity, temperature = T__f, pressure = Unit('atm'), air);
  if type(k, with_unit) and type(rho, with_unit) and type(Cp, with_unit) and type(mu, has_unit) then
    nu := mu/rho;
    alpha := k/(rho*Cp);
    beta := 1.0/T__f;
    Ra__L := g*beta*abs(T__s-`T__∞`)*L^3/(alpha*nu);
    Nus__L := (.825+.387*(Ra__L^(1/6))/(((1+(.492/Pr)^(9/16))^(8/27))))^2;
    h := Nus__L*k/L; q__h := h*(T__s-`T__∞`)*Area;
    q__r := varepsilon*sigma*(T__s^4-`T__∞`^4)*Area;
    q__total := q__h+q__r
  end if
end proc;
conv_VerticalPlate_proc(31.5*Unit('degC'), 20*Unit('degC'), Unit('m'^2), Unit('m'), .9);

                      100.8567010 Unit(W)
conv_VerticalPlate_proc(x*Unit('degC'), 20*Unit('degC'), Unit('m'^2), Unit('m'), .9);
                         0.001 Unit(W)
fsolve(100*Unit('W')-conv_VerticalPlate_proc(x, 20*Unit('degC'), Unit('m'^2), Unit('m'), .9), x = Unit('degC') .. 100*Unit('degC'));


I would appreciate help wih this.



Please Wait...