Newsletter January 2019

Julia Co-Creators Win Wilkinson Prize: Julia co-creators Jeff
Bezanson, Stefan Karpinski and Viral Shah are the winners of the
prestigious James H. Wilkinson Prize for Numerical
Software.
This prize is awarded every four years to recognize innovative software
in scientific computing. The award will be presented at the Society for
Industrial and Applied Mathematics (SIAM) Conference on Computational
Science and Engineering (CSE19) in Spokane, WA Feb 25-Mar 1, 2019.

Julia Growth Metrics

  Cumulative Total as of Jan 2018 Cumulative Total as of Jan 2019 Growth
Number of Julia Downloads Initiated 1.8 million 3.2 million +78%
Julia Packages Available 1,688 2,462 +46%
Number of News Articles Mentioning Julia 93 253 +172%
Julia Discourse Threads + Stack Overflow Questions 8,620 16,363 +90%
GitHub Stars for Julia Language (Not Including Julia Packages) 9,626 19,472 +102%
Published Citations of Julia: A Fresh Approach to Numerical Computing (2017) and Julia: A Fast Dynamic Language for Technical Computing (2012) 613 1,048 +71%

* Note: Julia can also call C, C++, Fortran, Python, R, Java and MPI
libraries

Julia GitHub Stars: Julia GitHub stars have nearly doubled since the
Aug 2018 release of Julia 1.0 and Julia is one of the top 10 languages
on GitHub
measured by stars and forks.

JuliaTeam:
JuliaTeam is an
enterprise solution that makes it as easy to use and develop Julia
packages inside your company as it is in the open source world. The
current release of JuliaTeam integrates with your corporate
authentication systems and eliminates the headaches of installing and
using public Julia packages behind a corporate firewall. JuliaTeam also
gives IT and management insight and control over what packages
developers are using, helping ensure quality and security.

Upcoming
JuliaTeam
releases will include these features as well:

  • Read and search docs for all internal and external packages in a single place

  • Create and manage private package registries

  • Publish and test private packages as easily as public ones, making sure new versions work seamlessly with all the other versions of packages that your teams are using

  • Benchmark your code to make sure it runs as efficiently as possible and stays fast

  • Download a summary of licenses of all the software you depend on

JuliaTeam makes
Julia development at your company as easy, effective and fun as open
source.

For more information, contact us.

Julia Training: Julia Computing offers online and offline
training including Introduction
to Julia, Machine Learning and Artificial Intelligence, Parallel
Computing in Julia and customized courses at every level. Registration
and course descriptions are available
here. Live instructor-led online
courses include:

New Julia Textbook: Introduction to Linear Applied Algebra –
Vectors, Matrices and Least Squares by
Stephen Boyd and Lieven Vandenberghe uses Julia as the language of
instruction and includes a Julia language
companion.
This is the textbook for EE103 (Stanford University) and EE133A (UCLA).
Other Julia textbooks, tutorials and teaching resources are available at
JuliaLang.org/Learning.

Julia for Space Mission Planning: Julia Language Ephemeris and
Physical Constants Reader for Solar System
Bodies was published by
Professor Kaela M. Martin (Embry-Riddle Aeronautical University),
undergraduates Julia Mihaylov and Renee Spear (Embry-Riddle Aeronautical
University) and Damon Landau (National Aeronautics and Space
Administration (NASA) Jet Propulsion Laboratory (JPL), California
Institute of Technology (Caltech)). This research “will help
interplanetary space mission designers more easily locate the relative
positions of planets, moons and small bodies in the solar system at
specific times in order to calculate optimum spacecraft trajectories.”

Julia for School Bus Optimization: The Institute for Operations
Research and Management Sciences (INFORMS) has
selected
MIT Sloan Professor Dimitris Bertsimas, Arthur Delarue and Sebastien
Martin as Franz Edelman award finalists for their use of Julia and
JuMP.jl to make Boston Public Schools bus routes 20% more efficient and
save $5 million per year which is being reinvested in schools.

JuMP-dev Annual Workshop: The third annual JuMP-dev
workshop takes place
March 12-14 in Santiago, Chile.
Registration
is free and talk proposal
submissions
are encouraged.

Julia and Julia Computing in the News

Julia Blog Posts

Upcoming Julia Events

Recent Julia Events

Julia Jobs, Fellowships and Internships

Do you work at or know of an organization looking to hire Julia
programmers as staff, research fellows or interns? Would your employer
be interested in hiring interns to work on open source packages that are
useful to their business? Help us connect members of our community to
great opportunities by sending us an
email, and we’ll get the word out.

There are more than 300 Julia jobs currently listed on
Indeed.com, including jobs at Accenture,
Airbus, Amazon, AstraZeneca, Barnes & Noble, BlackRock, Capital One,
Charles River Analytics, Citigroup, Comcast, Cooper Tire & Rubber,
Disney, Facebook, Gallup, Genentech, General Electric, Google, Huawei,
Johnson & Johnson, Match, McKinsey, NBCUniversal, Nielsen, OKCupid,
Oracle, Pandora, Peapod, Pfizer, Raytheon, Zillow, Brown, Emory,
Harvard, Johns Hopkins, Massachusetts General Hospital, Penn State, UC
Davis, University of Chicago, University of Virginia, Argonne National
Laboratory, Lawrence Berkeley National Laboratory, Los Alamos National
Laboratory, National Renewable Energy Laboratory, Oak Ridge National
Laboratory, State of Wisconsin and many more.

Contact Us: Please contact us if
you wish to:

  • Purchase or obtain license information for Julia products such as JuliaTeam, JuliaPro or JuliaBox

  • Obtain pricing for Julia consulting projects for your organization

  • Schedule Julia training for your organization

  • Share information about exciting new Julia case studies or use cases

  • Spread the word about an upcoming conference, workshop, training, hackathon, meetup, talk or presentation involving Julia

  • Partner with Julia Computing to organize a Julia meetup, conference, workshop, training, hackathon, talk or presentation involving Julia

  • Submit a Julia internship, fellowship or job posting

About Julia and Julia Computing

Julia is the fastest high performance open
source computing language for data, analytics, algorithmic trading,
machine learning, artificial intelligence, and other scientific and
numeric computing applications. Julia solves the two language problem by
combining the ease of use of Python and R with the speed of C++. Julia
provides parallel computing capabilities out of the box and unlimited
scalability with minimal effort. Julia has been downloaded more than 3
million times and is used at more than 1,500 universities. Julia
co-creators are the winners of the 2019 James H. Wilkinson Prize for
Numerical Software. Julia has run at
petascale on
650,000 cores with 1.3 million threads to analyze over 56 terabytes of
data using Cori, one of the ten largest and most powerful supercomputers
in the world.

Julia Computing was founded in 2015
by all the creators of Julia to develop products and provide
professional services to businesses and researchers using Julia.

Winter warm-up: toy models for heat exchangers

By: Julia on μβ

Re-posted from: https://matbesancon.xyz/post/2018-12-27-heat-exchanger/

Enjoying the calm of the frozen eastern French countryside for the last week of 2018,
I was struck by nostalgia while reading a SIAM news article [1] on a
near-reversible heat exchange between two flows and decided to dust off my
thermodynamics books (especially [2]).

Research in mathematical optimization was not the
obvious path I was on a couple years ago. The joint bachelor-master’s program
I followed in France was in process engineering, a discipline crossing
transfer phenomena (heat exchange, fluid mechanics, thermodynamics), control,
knowledge of the matter transformations at hand
(chemical, biochemical, nuclear reactions) and industrial engineering
(see note at the end of this page).

Hypotheses Throughout the article, we will use a set of flow hypotheses
which build up the core of our model for heat exchange.
These can seem odd but are pretty common in process engineering and
realistic in many applications.

  1. The two flows advance in successive “layers”.
  2. Each layer has a homogeneous temperature; we therefore ignore boundary layer effects.
  3. Successive layers do not exchange matter nor heat. The rationale behind this
    is that the temperature difference between fluids is significantly higher than between layers.
  4. Pressure losses in the exchanger does not release a significant heat compared to
    the fluid heat exchange.
  5. The fluid and wall properties are constant with temperature.

Starting simple: parallel flow heat exchange

In this model, both flows enter the exchanger on the same side, one at a
hot temperature, the other at a cold temperature. Heat is exchanged along the
exchanger wall, proportional at any point to the difference in temperature
between the two fluids. We therefore study the evolution of two variables
$u_1(x)$ and $u_2(x)$ in an interval $x \in [0,L]$ with $L$ the length of
the exchanger.

In any layer $[x, x + \delta x]$, the heat exchange is equal to:
$$\delta \dot{Q} = h \cdot (u_2(x) – u_1(x)) \cdot \delta x$$
with $h$ a coefficient depending on the wall heat exchange properties.

Moreover, the variation in internal energy of the hot flow is equal to
$\delta \dot{Q}$ and is also expressed as:

$$ c_2 \cdot \dot{m}_2 \cdot (u_2(x+\delta x) – u_2(x)) $$
$c_2$ is the calorific capacity of the hot flow and $\dot{m}_2$ its
mass flow rate. The you can check that the given expression is a power.
The same expressions apply to the cold flow.
Let us first assume the following:

$$c_2 \cdot \dot{m}_2 = c_1 \cdot \dot{m}_1$$

import DifferentialEquations
const DiffEq = DifferentialEquations
using Plots

function parallel_exchanger(du,u,p,x)
    h = p[1] # heat exchange coefficient
    Q = h * (u[1]-u[2])
    du[1] = -Q
    du[2] = Q
end

function parallel_solution(L, p)
    problem = DiffEq.ODEProblem(
      parallel_exchanger, # function describing the dynamics of the system
      u₀,                 # initial conditions u0
      (0., L),            # region overwhich the solution is built, x ∈ [0,L]
      p,                  # parameters, here the aggregated transfer constant h
    )
    return DiffEq.solve(problem, DiffEq.Tsit5())
end

plot(parallel_solution([0.0,100.0], 50.0, (0.05)))

$$ u_1(x) = T_{eq} \cdot (1 – e^{-h\cdot x}) $$
$$ u_2(x) = (100 – T_{eq}) \cdot e^{-h\cdot x} + T_{eq} $$

With $T_{eq}$ the limit temperature, trivially 50°C with equal flows.

(Full disclaimer: I’m a bit rusty and had to double-check for errors)

This model is pretty simple, its performance is however low from
a practical perspective. First on the purpose itself, we can compute for two
fluids the equilibrium temperature. This temperature can be adjusted
by the ratio of two mass flow rates but will remain a weighted average.
Suppose the goal of the exchange is to heat the cold fluid, the necessary
mass flow $\dot{m}_2$ tends to $\infty$ as the targeted temperature tends to
$u_2(L)$, and this is independent of the performance of the heat exchanger
itself, represented by the coefficient $h$. Here is the extended model using
the flow rate ratio to adjust the temperature profiles:

import DifferentialEquations
const DiffEq = DifferentialEquations

function ratio_exchanger(du,u,p,x)
    h = p[1] # heat exchange coefficient
    r = p[2] # ratio of mass flow rate 2 / mass flow rate 1
    Q = h * (u[1]-u[2])
    du[1] = -Q
    du[2] = Q / r
end

function ratio_solution(u₀, L, p)
    problem = DiffEq.ODEProblem(
      ratio_exchanger, # function describing the dynamics of the system
      u₀,              # initial conditions u0
      (0., L),         # region overwhich the solution is built, x ∈ [0,L]
      p,               # parameters, here the aggregated transfer constant h
    )
    return DiffEq.solve(problem, DiffEq.Tsit5())
end

for (idx,r) in enumerate((1.0, 5.0, 10.0, 500.0))
    plot(ratio_solution([0.0,100.0], 50.0, (0.05, r)))
    xlabel!("x (m)")
    ylabel!("T °C")
    title!("Parallel flow with ratio $r")
    savefig("parallel_ratio_$(idx).pdf")
end

This model has an analytical closed-form solution given by:
$$ T_{eq} = \frac{100\cdot \dot{m}_2}{\dot{m}_1 + \dot{m}_2} = 100\cdot\frac{r}{1+r} $$
$$ u_1(x) = T_{eq} \cdot (1 – e^{-h\cdot x}) $$
$$ u_2(x) = (100 – T_{eq}) \cdot e^{-h\cdot x \cdot r} + T_{eq} $$

Opposite flow model

This model is trickier because we don’t consider the dynamics of the system
along one dimension anymore. The two fluids flowing in opposite directions
are two interdependent systems. We won’t go through the analytical solution
but use a similar discretization as in article [1].

This model takes $n$ discrete cells, each considered at a given temperature.
Two cells of the cold and hot flows are considered to have exchanged heat
after crossing.

Applying the energy conservation principle, the gain of internal energy
between cell $k$ and $k+1$ for the cold flow is equal to the loss of
internal energy of the hot flow from cell $k+1$ to cell $k$. These differences
come from heat exchanged, expressed as:

$$\dot{Q}_k = h \cdot \Delta x \cdot (u_{2,k+1} – u_{1,k}) $$
$$\dot{Q}_k = \dot{m}_1 \cdot c_1 \cdot (u_{1,k+1} – u_{1,k}) $$
$$\dot{Q}_k = \dot{m}_2 \cdot c_2 \cdot (u_{2,k+1} – u_{2,k}) $$

Watch out the sense of the last equation since the heat exchange is
a loss for the hot flow. Again we use the simplifying assumption of
equality of the quantities:
$$ \dot{m}_i \cdot c_i $$

Our model only depends on the number of discretization steps $n$
and transfer coefficient $h$.

function discrete_crossing(n, h; itermax = 50000)
    u1 = Matrix{Float64}(undef, itermax, n)
    u2 = Matrix{Float64}(undef, itermax, n)
    u1[:,1] .= 0.0
    u1[1,:] .= 0.0
    u2[:,n] .= 100.0
    u2[1,:] .= 100.0
    for iter in 2:itermax
        for k in 1:n-1
            δq = h * (u2[iter-1, k+1] - u1[iter-1, k]) * (50.0/n)
            u2[iter, k]   = u2[iter-1, k+1] - δq
            u1[iter, k+1] = u1[iter-1, k]   + δq
        end
    end
    (u1,u2)
end
const (a1, a2) = discrete_crossing(500, 0.1)
const x0 = range(0.0, length = 500, stop = L)

p = plot(x0, a1[end,:], label = "u1 final", legend = :topleft)
plot!(p, x0, a2[end,:], label = "u2 final")
for iter in (100, 500)
    global p
    plot!(p, x0, a1[iter,:], label = "u1 $(iter)")
    plot!(p, x0, a2[iter,:], label = "u2 $(iter)")
end
xlabel!("x (m)")

We can observe the convergence of the solution at different iterations:

After convergence, we observe a parallel temperature profiles along the
exchanger, the difference between the two flows at any point being reduced
to $\epsilon$ mentioned in article [1]. The two differences between our model
and theirs are:

  • The discretization grid is slightly different since we consider the exchange
    to happen between cell $k$ and cell $k+1$ at the node between them, while they
    consider an exchange between $k-1$ and $k+1$ at cell $k$.
  • They consider two flow unit which just crossed reach the same temperature,
    while we consider a heat exchange limited by the temperature difference
    (the two flows do not reach identical temperatures but tend towards it).

Finally we can change the ratio:
$$\frac{\dot{m}_1\cdot c_1}{\dot{m}_2\cdot c_2}$$ for the counterflow model
as we did in the parallel case.

function discrete_crossing(n, h, ratio; itermax = 50000)
    u1 = Matrix{Float64}(undef, itermax, n)
    u2 = Matrix{Float64}(undef, itermax, n)
    u1[:,1] .= 0.0
    u1[1,:] .= 0.0
    u2[:,n] .= 100.0
    u2[1,:] .= 100.0
    for iter in 2:itermax
        for k in 1:n-1
            δq = h * (u2[iter-1, k+1] - u1[iter-1, k]) * 50.0 / n
            u2[iter, k]   = u2[iter-1, k+1] - δq * ratio
            u1[iter, k+1] = u1[iter-1, k]   + δq
        end
    end
    (u1,u2)
end

Julia tip: note that we do not define a new function for this but
create a method for the function discrete_crossing defined above
with a new signature (n, h, ratio).

We can plot the result:

const x0 = range(0.0, length = 500, stop = L)
p = plot(x0, a1[end,:], label = "u1 ratio 1.0", legend = :bottomright)
plot!(p, x0, a2[end,:], label = "u2 ratio 1.0")
for ratio in (0.1,0.5)
    global p
    (r1, r2) = discrete_crossing(500, 0.1, ratio)
    plot!(p, x0, r1[end,:], label = "u1 ratio $(ratio)")
    plot!(p, x0, r2[end,:], label = "u2 ratio $(ratio)")
end
xlabel!("x (m)")

Conclusion

To keep this post short, we will not show the influence of all parameters.
Some key effects to consider:

  • Increasing $h$ increases the gap between the flow temperatures
  • Increasing the number of steps does not change the result for a step size
    small enough
  • Increasing the exchanger length reduces the gap
  • A ratio of 1 minimizes the temperature difference at every point
    (and thus minimizes the entropy). This very low entropy creation is a positive
    sign for engineers from a thermodynamics point of view: we are not “degrading”
    the “quality” of available energy to perform this heat exchange or in other
    terms, we are not destroying exergy.

Feel free to reach out on Twitter
or via email if you have comments or questions, I’d be glad to take both.

Note on process engineering
The term is gaining more traction in English, and should replace
chemical engineering in higher education to acknowledge the diversity of
application fields, greater than the chemical industry alone.
The German equivalent Verfahrenstechnik has been used for decades and
Génie des Procédés is now considered a norm in most French-speaking
universities and consortia.

Edit: thanks BYP for the sharp-as-ever proofreading

Sources:

[1] Levi M. A Near-perfect Heat Exchange. SIAM news. 2018 Dec;51(10):4.

[2] Borel L, Favrat D. Thermodynamique et énergétique. PPUR presses polytechniques; 2nd edition, 2011.

Image sources:
[3] Geogebra

Winter warm-up: toy models for heat exchangers

By: Julia on μβ

Re-posted from: https://matbesancon.github.io/post/2018-12-27-heat-exchanger/

Enjoying the calm of the frozen eastern French countryside for the last week of 2018,
I was struck by nostalgia while reading a SIAM news article [1] on a
near-reversible heat exchange between two flows and decided to dust off my
thermodynamics books (especially [2]).

Research in mathematical optimization was not the
obvious path I was on a couple years ago. The joint bachelor-master’s program
I followed in France was in process engineering, a discipline crossing
transfer phenomena (heat exchange, fluid mechanics, thermodynamics), control,
knowledge of the matter transformations at hand
(chemical, biochemical, nuclear reactions) and industrial engineering
(see note at the end of this page).

Hypotheses Throughout the article, we will use a set of flow hypotheses
which build up the core of our model for heat exchange.
These can seem odd but are pretty common in process engineering and
realistic in many applications.

  1. The two flows advance in successive “layers”.
  2. Each layer has a homogeneous temperature; we therefore ignore boundary layer effects.
  3. Successive layers do not exchange matter nor heat. The rationale behind this
    is that the temperature difference between fluids is significantly higher than between layers.
  4. Pressure losses in the exchanger does not release a significant heat compared to
    the fluid heat exchange.
  5. The fluid and wall properties are constant with temperature.

Starting simple: parallel flow heat exchange

In this model, both flows enter the exchanger on the same side, one at a
hot temperature, the other at a cold temperature. Heat is exchanged along the
exchanger wall, proportional at any point to the difference in temperature
between the two fluids. We therefore study the evolution of two variables
$u_1(x)$ and $u_2(x)$ in an interval $x \in [0,L]$ with $L$ the length of
the exchanger.

In any layer $[x, x + \delta x]$, the heat exchange is equal to:
$$\delta \dot{Q} = h \cdot (u_2(x) – u_1(x)) \cdot \delta x$$
with $h$ a coefficient depending on the wall heat exchange properties.

Moreover, the variation in internal energy of the hot flow is equal to
$\delta \dot{Q}$ and is also expressed as:

$$ c_2 \cdot \dot{m}_2 \cdot (u_2(x+\delta x) – u_2(x)) $$
$c_2$ is the calorific capacity of the hot flow and $\dot{m}_2$ its
mass flow rate. The you can check that the given expression is a power.
The same expressions apply to the cold flow.
Let us first assume the following:

$$c_2 \cdot \dot{m}_2 = c_1 \cdot \dot{m}_1$$

import DifferentialEquations
const DiffEq = DifferentialEquations
using Plots

function parallel_exchanger(du,u,p,x)
    h = p[1] # heat exchange coefficient
    Q = h * (u[1]-u[2])
    du[1] = -Q
    du[2] = Q
end

function parallel_solution(L, p)
    problem = DiffEq.ODEProblem(
      parallel_exchanger, # function describing the dynamics of the system
      u₀,                 # initial conditions u0
      (0., L),            # region overwhich the solution is built, x ∈ [0,L]
      p,                  # parameters, here the aggregated transfer constant h
    )
    return DiffEq.solve(problem, DiffEq.Tsit5())
end

plot(parallel_solution([0.0,100.0], 50.0, (0.05)))

$$ u_1(x) = T_{eq} \cdot (1 – e^{-h\cdot x}) $$
$$ u_2(x) = (100 – T_{eq}) \cdot e^{-h\cdot x} + T_{eq} $$

With $T_{eq}$ the limit temperature, trivially 50°C with equal flows.

(Full disclaimer: I’m a bit rusty and had to double-check for errors)

This model is pretty simple, its performance is however low from
a practical perspective. First on the purpose itself, we can compute for two
fluids the equilibrium temperature. This temperature can be adjusted
by the ratio of two mass flow rates but will remain a weighted average.
Suppose the goal of the exchange is to heat the cold fluid, the necessary
mass flow $\dot{m}_2$ tends to $\infty$ as the targeted temperature tends to
$u_2(L)$, and this is independent of the performance of the heat exchanger
itself, represented by the coefficient $h$. Here is the extended model using
the flow rate ratio to adjust the temperature profiles:

import DifferentialEquations
const DiffEq = DifferentialEquations

function ratio_exchanger(du,u,p,x)
    h = p[1] # heat exchange coefficient
    r = p[2] # ratio of mass flow rate 2 / mass flow rate 1
    Q = h * (u[1]-u[2])
    du[1] = -Q
    du[2] = Q / r
end

function ratio_solution(u₀, L, p)
    problem = DiffEq.ODEProblem(
      ratio_exchanger, # function describing the dynamics of the system
      u₀,              # initial conditions u0
      (0., L),         # region overwhich the solution is built, x ∈ [0,L]
      p,               # parameters, here the aggregated transfer constant h
    )
    return DiffEq.solve(problem, DiffEq.Tsit5())
end

for (idx,r) in enumerate((1.0, 5.0, 10.0, 500.0))
    plot(ratio_solution([0.0,100.0], 50.0, (0.05, r)))
    xlabel!("x (m)")
    ylabel!("T °C")
    title!("Parallel flow with ratio $r")
    savefig("parallel_ratio_$(idx).pdf")
end




This model has an analytical closed-form solution given by:
$$ T_{eq} = \frac{100\cdot \dot{m}_2}{\dot{m}_1 + \dot{m}_2} = 100\cdot\frac{r}{1+r} $$
$$ u_1(x) = T_{eq} \cdot (1 – e^{-h\cdot x}) $$
$$ u_2(x) = (100 – T_{eq}) \cdot e^{-h\cdot x \cdot r} + T_{eq} $$

Opposite flow model

This model is trickier because we don’t consider the dynamics of the system
along one dimension anymore. The two fluids flowing in opposite directions
are two interdependent systems. We won’t go through the analytical solution
but use a similar discretization as in article [1].

This model takes $n$ discrete cells, each considered at a given temperature.
Two cells of the cold and hot flows are considered to have exchanged heat
after crossing.


[3]

Applying the energy conservation principle, the gain of internal energy
between cell $k$ and $k+1$ for the cold flow is equal to the loss of
internal energy of the hot flow from cell $k+1$ to cell $k$. These differences
come from heat exchanged, expressed as:

$$\dot{Q}_k = h \cdot \Delta x \cdot (u_{2,k+1} – u_{1,k}) $$
$$\dot{Q}_k = \dot{m}_1 \cdot c_1 \cdot (u_{1,k+1} – u_{1,k}) $$
$$\dot{Q}_k = \dot{m}_2 \cdot c_2 \cdot (u_{2,k+1} – u_{2,k}) $$

Watch out the sense of the last equation since the heat exchange is
a loss for the hot flow. Again we use the simplifying assumption of
equality of the quantities:
$$ \dot{m}_i \cdot c_i $$

Our model only depends on the number of discretization steps $n$
and transfer coefficient $h$.

function discrete_crossing(n, h; itermax = 50000)
    u1 = Matrix{Float64}(undef, itermax, n)
    u2 = Matrix{Float64}(undef, itermax, n)
    u1[:,1] .= 0.0
    u1[1,:] .= 0.0
    u2[:,n] .= 100.0
    u2[1,:] .= 100.0
    for iter in 2:itermax
        for k in 1:n-1
            δq = h * (u2[iter-1, k+1] - u1[iter-1, k]) * (50.0/n)
            u2[iter, k]   = u2[iter-1, k+1] - δq
            u1[iter, k+1] = u1[iter-1, k]   + δq
        end
    end
    (u1,u2)
end
const (a1, a2) = discrete_crossing(500, 0.1)
const x0 = range(0.0, length = 500, stop = L)

p = plot(x0, a1[end,:], label = "u1 final", legend = :topleft)
plot!(p, x0, a2[end,:], label = "u2 final")
for iter in (100, 500)
    global p
    plot!(p, x0, a1[iter,:], label = "u1 $(iter)")
    plot!(p, x0, a2[iter,:], label = "u2 $(iter)")
end
xlabel!("x (m)")

We can observe the convergence of the solution at different iterations:

After convergence, we observe a parallel temperature profiles along the
exchanger, the difference between the two flows at any point being reduced
to $\epsilon$ mentioned in article [1]. The two differences between our model
and theirs are:

  • The discretization grid is slightly different since we consider the exchange
    to happen between cell $k$ and cell $k+1$ at the node between them, while they
    consider an exchange between $k-1$ and $k+1$ at cell $k$.
  • They consider two flow unit which just crossed reach the same temperature,
    while we consider a heat exchange limited by the temperature difference
    (the two flows do not reach identical temperatures but tend towards it).

Finally we can change the ratio:
$$\frac{\dot{m}_1\cdot c_1}{\dot{m}_2\cdot c_2}$$ for the counterflow model
as we did in the parallel case.

function discrete_crossing(n, h, ratio; itermax = 50000)
    u1 = Matrix{Float64}(undef, itermax, n)
    u2 = Matrix{Float64}(undef, itermax, n)
    u1[:,1] .= 0.0
    u1[1,:] .= 0.0
    u2[:,n] .= 100.0
    u2[1,:] .= 100.0
    for iter in 2:itermax
        for k in 1:n-1
            δq = h * (u2[iter-1, k+1] - u1[iter-1, k]) * 50.0 / n
            u2[iter, k]   = u2[iter-1, k+1] - δq * ratio
            u1[iter, k+1] = u1[iter-1, k]   + δq
        end
    end
    (u1,u2)
end

Julia tip: note that we do not define a new function for this but
create a method for the function discrete_crossing defined above
with a new signature (n, h, ratio).

We can plot the result:

const x0 = range(0.0, length = 500, stop = L)
p = plot(x0, a1[end,:], label = "u1 ratio 1.0", legend = :bottomright)
plot!(p, x0, a2[end,:], label = "u2 ratio 1.0")
for ratio in (0.1,0.5)
    global p
    (r1, r2) = discrete_crossing(500, 0.1, ratio)
    plot!(p, x0, r1[end,:], label = "u1 ratio $(ratio)")
    plot!(p, x0, r2[end,:], label = "u2 ratio $(ratio)")
end
xlabel!("x (m)")

Conclusion

To keep this post short, we will not show the influence of all parameters.
Some key effects to consider:

  • Increasing $h$ increases the gap between the flow temperatures
  • Increasing the number of steps does not change the result for a step size
    small enough
  • Increasing the exchanger length reduces the gap
  • A ratio of 1 minimizes the temperature difference at every point
    (and thus minimizes the entropy). This very low entropy creation is a positive
    sign for engineers from a thermodynamics point of view: we are not “degrading”
    the “quality” of available energy to perform this heat exchange or in other
    terms, we are not destroying exergy.

Feel free to reach out on Twitter
or via email if you have comments or questions, I’d be glad to take both.

Note on process engineering
The term is gaining more traction in English, and should replace
chemical engineering in higher education to acknowledge the diversity of
application fields, greater than the chemical industry alone.
The German equivalent Verfahrenstechnik has been used for decades and
Génie des Procédés is now considered a norm in most French-speaking
universities and consortia.

Edit: thanks BYP for the sharp-as-ever proofreading

Sources:

[1] Levi M. A Near-perfect Heat Exchange. SIAM news. 2018 Dec;51(10):4.

[2] Borel L, Favrat D. Thermodynamique et énergétique. PPUR presses polytechniques; 2nd edition, 2011.

Image sources:
[3] Geogebra