Maxwell’s construction on van der Waals isotherms

Code: #111-000

File: apps/van_der_waals/p_v_2D.ipynb

Run it online: Binder


Interface

The main interface (main_block_111_000) is divided in two Box: top_block_111_000 and middle_block_111_000.

top_block_111_000 contains the widgets to control the tenperatute range shown and the toggles of the different areas and isother types.

middle_block_111_000 contains the main figure (fig_111_000) and the different controlers to implement Maxwell’s construction: T_slider changes the selected isotherm and p_slider controls the isobaric line.

Some buttons are only available if maxwell_construction_checkbox is activated.

v_slider controls the position of the mark tracer and it is only available if show_tracer_checkbox is activated.

[1]:
from IPython.display import Image
Image(filename='../../static/images/apps/111-000_1.png')
[1]:
../../_images/apps_van_der_waals_p_v_2D_2_0.png

CSS

A custom css file is used to improve the interface of this application. It can be found here.

[2]:
from IPython.display import HTML
display(HTML("<head><link rel='stylesheet' type='text/css' href='./../../static/custom.css'></head>"))
display(HTML("<style>.container { width:100% !important; } .jupyter-button {white-space: normal !important;}</style>"))

Packages

[3]:
from bqplot import *
import bqplot as bq
import bqplot.marks as bqm
import bqplot.scales as bqs
import bqplot.axes as bqa

import ipywidgets as widgets

from scipy import interpolate
from scipy.signal import argrelextrema

import urllib.parse
import webbrowser

import sys

Physical functions

This are the functions that have a physical meaning:

  • get_relative_isotherms
  • experimetal_isotherms
  • get_roots
  • p_indefinite_integral
  • p_definite_integral
  • find_real_fixed_p
  • find_real_fixed_T
[4]:
def get_relative_isotherms(v_range, T_range):
    """This function calculates the theoretical p(v, T) plane
        (in reduced coordinates) according to van der Waals
        equation of state from a given range of volumes
        and tenperatures.

    Args:
        v_range: An array containing the values of v
        (in reduced coordinates)for which the isotherms must be
        calculated.\n
        T_range: An array containing the values of T
        (in reduced coordinates)for which the isotherms must be
        calculated.\n


    Returns:
        isotherms: A list consisted of numpy arrays containing the
        pressures of each isotherm.
    """

    isotherms = []

    for T in T_range:
        p_R = []
        for v in v_range:
            val = (8.0/3.0*T/(v - 1.0/3.0) - 3.0/v**2)
            p_R = np.append(p_R, val)

        isotherms.append(p_R)

    return isotherms
[5]:
def experimental_isotherms(p_range, v_range, T_range, fixed_p, fixed_T):
    """This function calculates the experimental p(v, T) plane
        (in reduced coordinates) according to van der Waals
        equation of state for a given range of volumes
        and tenperatures or for a given range of volumes
        and pressures.

    Args:
        p_range: An array containing the values of p
        (in reduced coordinates)for which the isotherms must be
        calculated. Only used if fixed_p == True.\n
        v_range: An array containing the values of v
        (in reduced coordinates)for which the isotherms must be
        calculated.\n
        T_range: An array containing the values of v
        (in reduced coordinates)for which the isotherms must be
        calculated. Only used if fixed_T == True.\n
        fixed_p: Boolean variable which represents if the isotherms
        must be calculated for a given pressures.\n
        fixed_T: Boolean variable which represents if the isotherms
        must be calculated for a given pressures.\n


    Returns:
        expe_data: A list consisted of numpy arrays containing the
        pressures of each theoretical isotherm.\n
        theo_data: A list consisted of numpy arrays containing the
        pressures of each theoretical isotherm.\n
        v_limits: A list consisted of arrays of the volume limits of
        the phase-transition of each subcritical isotherm.\n
        p_limits: A list consisted of arrays of the pressure limits of
        the phase-transition of each subcritical isotherm.\n
        tenperatures: A list consisted of the tenperatures of the
        isotherms.\n
    """

    if fixed_T:

        theo_data = get_relative_isotherms(v_range, T_range)
        expe_data = []

        v_limits = []
        p_limits = []

        p_range = np.linspace(0.001, 1.0, num=10000)
        pressures, v_isobaric_limits = find_real_fixed_T(p_range, T_range)

        for i in range(len(theo_data)):

            p_expe = []

            if i < len(v_isobaric_limits):

                    v_lim = v_isobaric_limits[i]

                    if len(v_lim) > 1: #check if there is only one point
                        for j in range(len(v_range)):

                            if v_range[j] > v_lim[0] and v_range[j] < v_lim[1]:
                                p_expe.append(pressures[i])

                            else:
                                p_expe.append(theo_data[i][j])

                        v_limits = np.append(v_limits, [v_lim[0], v_lim[1]])
                        p_limits = np.append(p_limits, [pressures[i], pressures[i]])

                    else:
                        p_expe = theo_data[i]
                        v_limits = np.append(v_limits, [1.0])
                        p_limits = np.append(p_limits, [1.0])

            else:

                p_expe = theo_data[i]

            expe_data.append(p_expe)

        tenperatures = T_range

        return expe_data, theo_data, p_limits, v_limits, tenperatures

    elif fixed_p:

        tenperatures, v_isobaric_limits = find_real_fixed_p(p_range, T_range)

        theo_data = get_relative_isotherms(v_range, tenperatures)
        expe_data = []

        v_limits = []
        p_limits = []

        for i in range(len(theo_data)):

            p_expe = []

            if i < len(v_isobaric_limits):

                    v_lim = v_isobaric_limits[i]

                    if len(v_lim) > 1: #check if there is only one point

                        for j in range(len(v_range)):

                            if v_range[j] > v_lim[0] and v_range[j] < v_lim[1]:
                                p_expe.append(p_range[i])

                            else:
                                p_expe.append(theo_data[i][j])

                        v_limits = np.append(
                            v_limits,
                            [v_lim[0],
                             v_lim[1]]
                        )
                        p_limits = np.append(
                            p_limits,
                            [p_range[i],
                             p_range[i]]
                        )

                    else:
                        p_expe = theo_data[i]
                        v_limits = np.append(v_limits, [1.0])
                        p_limits = np.append(p_limits, [1.0])

            else:

                p_expe = theo_data[i]

            expe_data.append(p_expe)

        return expe_data, theo_data, p_limits, v_limits, tenperatures
[6]:
def get_roots(p, T):
    """This function calculates the roots of a van der Waals
    isotherm of a given T and set of pressures.

    Args:
        p: Numpy array consisted of the pressures of the isotherm.\n
        T: Value of the tenperature.\n

    Returns:
        roots_in_range: A list consisted of the real roots.\n
    """


    roots = np.roots([1.0, - 1.0/3.0*(1.0 + 8.0*T/p), 3.0/p, -1.0/p])
    roots_in_range = []

    for root in roots:
        if np.isreal(root):
            root = np.real(root)
            if root > 0:
                roots_in_range.append(root)
    roots_in_range.sort()

    return roots_in_range
[7]:
def p_indefinite_integral(p_0, v_0, T):
    """This function calculates the indefinite integral between
    a van der Waals isotherm and a isobaric line.

    Args:
        p0: Isobaric line's pressure.\n
        v0: Value of the volume.\n
        T: Value of the tenperature.\n

    Returns:
        integral: Value of the indefinite integral between a
        van der Waals isotherm at T and a isobaric line of p0 at a
        volume v0.\n
    """

    integral = 8.0/3.0 * T *np.log(v_0 - 1.0/3.0) + 3.0/v_0 - p_0*v_0

    return integral
[8]:
def definite_integral(p_0, v_range, T):
    """This function 'p_indefinite_integral' function to calculate
    the definite integral between a van der Waals isotherm and a
    isobaric line.

    Args:
        p0: Isobaric line's pressure.\n
        v_range: Tuple or list consisted of volume limits.\n
        T: Value of the tenperature.\n

    Returns:
        integral: Value of the definite integral between a
        van der Waals isotherm at T and a isobaric line of p0 in a
        volume range v_range.\n
    """

    v_0, v_1 = v_range[0], v_range[1]

    integral = p_indefinite_integral(p_0, v_1, T) - p_indefinite_integral(p_0, v_0, T)

    return integral
[9]:
def find_real_fixed_T(p_values, T_values):
    """This function uses Maxwell's construction to find the
       pressures in which phase transition happens given some
       fixed tenperatures.\n

    Args:
        p_values: List of pressures in which the real isotherm is
        searched.\n
        T_values: List of tenperatures of the isotherms.\n


    Returns:
        pressures: List of pressures in which phase transition
        happens.\n
        v_range: Volume limits of phase transition zones.
    """

    eps = 1e-3

    pressures = []
    v_ranges = []

    for T in T_values:

        if T < 1.0:

            for p in p_values:

                roots = get_roots(p, T)

                if len(roots) == 3:

                    v_range = [roots[0], roots[2]]
                    area = definite_integral(p, v_range, T)

                    if abs(area) < eps:

                        pressures.append(p)
                        v_ranges.append(v_range)

                        break

        elif T == 1.0:

            pressures.append(1.0)
            v_ranges.append([1.0])

    return pressures, v_ranges
[10]:
def find_real_fixed_p(p_values, T_values):
    """This function uses Maxwell's construction to find the
       tenperatures in which phase transition happens given some
       fixed pressures.\n

    Args:
        p_values: List of pressures of the isotherms.\n
        T_values: List of tenperatures in which the real isotherm is
        searched.\n


    Returns:
        tenperatures: List of tenperatures in which phase transition
        happens.\n
        v_range: Volume limits of phase transition zones.
    """

    eps = 1e-3

    tenperatures = []
    v_ranges = []

    for p in p_values:

        if p < 1.0:

            for T in T_values:

                roots = get_roots(p, T)

                if len(roots) == 3:

                    v_range = [roots[0], roots[2]]
                    area = definite_integral(p, v_range, T)

                    if abs(area) < eps:

                        tenperatures.append(T)
                        v_ranges.append(v_range)

                        break

        elif p == 1.0:

            tenperatures.append(1.0)
            v_ranges.append([1.0])

    return tenperatures, v_ranges

Main interface

[ ]:
"""

.. module:: p_v_2D.ipynb
    :sypnopsis: This module creates an interface to interact with the
    van der Waals isotherms in the p(v, T) plane.\n

.. moduleauthor:: Jon Gabirondo López (jgabirondo001@ikasle.ehu.eus)

"""

saved = []

#############################
######TOP BLOCK##############
#############################

top_block_111_000 = widgets.VBox(
    [],
    layout=widgets.Layout(align_items='center')
)

change_view_button = widgets.ToggleButton(
    value=False,
    description='Presentation mode (OFF)',
    disabled=False,
    button_style='',
    tooltip='',
    icon='desktop',
    layout=widgets.Layout(
        width='auto'
    )
)

change_view_button.observe(change_view, 'value')

t_min_input = widgets.BoundedFloatText(
    value=0.85,
    min=0.5,
    max=3.0,
    step=0.005,
    description=r'\( T_{min}: \)',
    disabled=False
)

t_max_input = widgets.BoundedFloatText(
    value=1.05,
    min=0.5,
    max=3.0,
    step=0.005,
    description=r'\( T_{max}: \)',
    disabled=False
)

num_input = widgets.IntSlider(
    value=6,
    min=1,
    max=20,
    description=r'\( \# \)',
    disabled=False,
    continuous_update=True,
    orientation='horizontal',
    readout=True,
)

generate_button = widgets.Button(
    description='Generate',
    disabled=False,
    button_style='',
    tooltip='Generate isotherms in the selected T range.',
)

generate_button.on_click(change_isotherms)

#############################
#isotherm_type_block
#############################

theoretical_toggle = widgets.ToggleButton(
    value=True,
    description='Theoretical',
    disabled=False,
    button_style='',
    tooltip='Show theoretical isotherms',
    style = {'description_width': 'initial'},
    layout = widgets.Layout(
        align_self='center',
        width='100%'
    )
)

experimental_toggle = widgets.ToggleButton(
    value=False,
    description='Experimental',
    disabled=False,
    button_style='',
    tooltip='Show experimental isotherms',
    style = {'description_width': 'initial'},
    layout = widgets.Layout(
        align_self='center',
        width='100%'
    ),
)

theoretical_toggle.observe(change_isotherm_type, 'value')
experimental_toggle.observe(change_isotherm_type, 'value')

show_transition_limits_checkbox = widgets.Checkbox(
    value=False,
    description='Show transition limits',
    disabled=False,
    layout = widgets.Layout(width='148px !important')
)

show_transition_line_checkbox = widgets.Checkbox(
    value=False,
    description='Show transition line',
    disabled=False,
    layout = widgets.Layout(width='148px !important')
)

show_existence_limits_checkbox = widgets.Checkbox(
    value=False,
    description='Show existence limits',
    disabled=False,
    layout = widgets.Layout(width='148px !important')
)

show_existence_line_checkbox = widgets.Checkbox(
    value=False,
    description='Show existence line',
    disabled=False,
    layout = widgets.Layout(width='148px !important')
)

show_transition_limits_checkbox.observe(show_limits, 'value')
show_transition_line_checkbox.observe(show_limits, 'value')
show_existence_limits_checkbox.observe(show_limits, 'value')
show_existence_line_checkbox.observe(show_limits, 'value')

isotherm_type_block = widgets.HBox(
    [],
    layout = widgets.Layout(
        align_self='center',
        margin='10px 0 0 0'
    )
)

top_block_111_000.children = [
    change_view_button,
    widgets.HTML('</br>'),
    widgets.HBox([
        t_min_input,
        t_max_input,
        num_input,
        generate_button
    ]),
    widgets.HBox([
        theoretical_toggle,
        experimental_toggle,
    ],
        layout = widgets.Layout(
            margin='20px 0 20px 0',
            width='33%')
    ),
    widgets.HBox([
        show_transition_limits_checkbox,
        show_transition_line_checkbox,
        show_existence_limits_checkbox,
        show_existence_line_checkbox
    ],
        layout = widgets.Layout(
            margin='20px 0 20px 0',
            width='100%')
    )
]

#############################
############MIDDLE BLOCK#####
#############################

middle_block_111_000 = widgets.HBox(
    [],
    layout=widgets.Layout(width='100%')
)

T_values = get_t_range(t_min_input, t_max_input, num_input)

############################
#FIGURE
############################

scale_x = bqs.LinearScale(min = 0.45, max = 5.0)
scale_y = bqs.LinearScale(min = 0.0, max = 2.0)

fig_111_000 = bq.Figure(
    title='Van der Waals isotherms',
    marks=[],
    axes=[],
    animation_duration=0,
    layout = widgets.Layout(
        align_self='center',
        width='75%'
    ),
    legend_location='top-right',
    background_style= {'fill': 'white','stroke': 'black'},
    fig_margin=dict(top=80,bottom=80,left=60,right=30),
    toolbar = True
    )

axis_x = bqa.Axis(
    scale=scale_x,
    tick_format='0.2f',
    tick_style={'font-size': '15px'},
    tick_values=[1.0, 2.0, 3.0, 4.0, 5.0],
    grid_lines = 'none',
    grid_color = '#8e8e8e',
    label='v',
    label_location='middle',
    label_style={'stroke': 'black', 'default-size': 35},
    label_offset='50px'
)

axis_y = bqa.Axis(
    scale=scale_y,
    tick_format='0.2f',
    tick_style={'font-size': '15px'},
    tick_values=[0.0, 1.0, 2.0],
    grid_lines = 'none',
    grid_color = '#8e8e8e',
    orientation='vertical',
    label='p',
    label_location='middle',
    label_style={'stroke': 'red', 'default_size': 35},
    label_offset='50px'
)

fig_111_000.axes = [axis_x, axis_y]


tt = bq.Tooltip(
    fields = ['y', 'x'],
    formats = ['.3f', '.3f'],
    labels = ['p', 'v']
)

marks = []

tracer = bqm.Scatter(
    name = 'tracer',
    x = [1.0],
    y = [1.0],
    scales = {'x': scale_x, 'y': scale_y},
    visible = False,
    colors = ['black'],
    names = [],
    labels=['tracer'],
    tooltip = tt
)

marks.append(tracer)

labels_points = bqm.Scatter(
    name = 'labels',
    x = [],
    y = [],
    scales = {'x': scale_x, 'y': scale_y},
    #opacities = [1.0],
    visible = True,
    colors = ['black'],
    names = [],
    labels=['labels'],
    tooltip = tt,
)

marks.append(labels_points)


isobaric = bqm.Lines(
    x = [],
    y = [],
    scales = {'x': scale_x, 'y': scale_y},
    opacities = [1.0, 0.0],
    visible = False,
    colors = ['green'],
    fill_opacities = [0.35],
    fill = 'none',
    labels = ['isobaric'],
)

marks.append(isobaric)

################################
#LEFT BLOCK
################################

left_block_111_000 = widgets.VBox(
    [],
    layout = widgets.Layout(
        width='15%',
        margin = '70px 0 0 0'
    )
)

maxwell_construction_checkbox = widgets.Checkbox(
    value=False,
    description='Maxwell Construction',
    disabled=False
)

maxwell_construction_checkbox.observe(show_maxwell, 'value')

p_slider = widgets.FloatSlider(
    value=1.0,
    min=0.001,
    max=scale_y.max,
    step=0.001,
    description=r'\( p_r \)',
    disabled=True,
    continuous_update=True,
    orientation='vertical',
    readout=True,
    readout_format='.3f',
    layout = widgets.Layout(
        height = '80%',
        margin = '45px 0 0 0'
    )
)

p_slider.observe(update_isobaric, 'value')

shade_areas_checkbox = widgets.Checkbox(
    value=False,
    description='shade areas',
    disabled=True
)

shade_areas_checkbox.observe(shade_areas, 'value')

integral_value_text = widgets.HTML(
    value="",
)

integral_value_left_text = widgets.HTML(
    value="",
)

integral_value_right_text = widgets.HTML(
    value="",
)

fixed_T_button = widgets.Button(
    description='Calculate real isotherm (fixed T)',
    disabled=True,
    button_style='',
    tooltip="Use Maxwell's construction to" \
            "calculate real isotherm (fixed T)",
    layout = {
        'width' : 'auto',
        'height' : 'auto',
    }
)

fixed_T_button.on_click(find_real_isotherm)

fixed_p_button = widgets.Button(
    description='Calculate real isotherm (fixed p)',
    disabled=True,
    button_style='',
    tooltip="Use Maxwell's construction to" \
            "calculate real isotherm (fixed p)",
    layout = {
        'width' : 'auto',
        'height' : 'auto',
    }
)

fixed_p_button.on_click(add_isotherm_mark_fixed_p)

left_block_111_000.children = [
    maxwell_construction_checkbox,
    shade_areas_checkbox,
    widgets.HTML(value="</br>"),
    widgets.HBox([
        widgets.Label(value="$Left:$"),
        integral_value_left_text,
    ]),
    widgets.HBox([
        widgets.Label(value="$Right:$"),
        integral_value_right_text,
    ]),
    widgets.HBox([
        widgets.Label(value="$Sum:$"),
        integral_value_text,
    ]),
    widgets.HTML(value="</br>"),
    fixed_T_button,
    fixed_p_button
]

################################
#RIGHT BLOCK
################################

right_block_111_000 = widgets.VBox(
    [],
    layout = widgets.Layout(
        width='15%',
        margin = '70px 0 0 0'
    )
)

T_slider = widgets.SelectionSlider(
    options=T_values,
    value=T_values[int(T_values.size/2)],
    description=r'\( T_r \)',
    disabled=False,
    continuous_update=True,
    orientation='vertical',
    readout=True,
    layout = widgets.Layout(
        height = '80%',
        margin = '45px 0 0 0'
    )
)

T_slider.observe(select_isotherm, 'value')

select_visible = widgets.SelectMultiple(
    options=[str(t) for t in T_values],
    value=[str(t) for t in T_values],
    rows=10,
    description='',
    disabled=False,
    layout = widgets.Layout(width = '90%')
)

select_visible.observe(show_isotherm, 'value')

show_all_button = widgets.Button(
    description='Show all',
    disabled=False,
    button_style='',
    tooltip='Show all isotherms',
)

show_all_button.on_click(show_all)

opacity_theo_slider = widgets.FloatSlider(
    value=0.2,
    min=0.0,
    max=1.0,
    step=0.05,
    description="",
    disabled=False,
    continuous_update=True,
    orientation='horizontal',
    readout=True,
    layout = widgets.Layout(width = '90%')
)

opacity_theo_slider.observe(select_isotherm, 'value')

opacity_expe_slider = widgets.FloatSlider(
    value=0.2,
    min=0.0,
    max=1.0,
    step=0.05,
    description="",
    disabled=False,
    continuous_update=True,
    orientation='horizontal',
    readout=True,
    layout = widgets.Layout(width = '90%')
)

opacity_expe_slider.observe(select_isotherm, 'value')

tracer_p = widgets.Label(value=" - ")
tracer_p_base = widgets.Label(value="$p =$")

tracer_v = widgets.Label(value=" - ")
tracer_v_base = widgets.Label(value="$v =$")

tracer_T = widgets.Label(value=" - ")
tracer_T_base = widgets.Label(value="$T =$")

right_block_111_000.children = [
    widgets.Label(value="Select visible isotherms:"),
    select_visible,
    show_all_button,
    widgets.Label(
        value="Theoretical opacity:",
        layout=widgets.Layout(
            margin='20px 0 0 0'
        )
    ),
    opacity_theo_slider,
    widgets.Label(
        value="Experimental opacity:",
        layout=widgets.Layout(
            margin='20px 0 0 0'
        )
    ),
    opacity_expe_slider,
    widgets.HBox([
        tracer_p_base,
        tracer_p,
    ]),
    widgets.HBox([
        tracer_v_base,
        tracer_v,
    ]),
    widgets.HBox([
        tracer_T_base,
        tracer_T,
    ]),
]

################################
#tracer_block
################################

tracer_block_111_000 = widgets.VBox(
    [],
    layout=widgets.Layout(
        align_items="center",
        width="100%"
    )
)

show_tracer_checkbox = widgets.Checkbox(
    value=False,
    description='Show tracer',
    disabled=False,
    layout=widgets.Layout(width='40%')
)

show_tracer_checkbox.observe(show_tracer, 'value')


v_slider = widgets.FloatSlider(
    value=1.0,
    min=0.5,
    max=5.0,
    step=0.001,
    description=r'\( v \)',
    disabled=True,
    continuous_update=True,
    orientation='horizontal',
    readout=True,
    readout_format='.2f',
    layout = widgets.Layout(width = '90%')
)

v_slider.observe(update_tracer, 'value')

label_input = widgets.Text(
    value='',
    placeholder="Label:",
    disabled = True,
)

add_label_button = widgets.Button(
        description='Add label',
        disabled=True,
        button_style='',
        tooltip="Add label in tracer's position",
    )

add_label_button.on_click(add_label_button_clicked)

undo_label_button = widgets.Button(
        description='Undo',
        disabled=True,
        button_style='',
        tooltip="Remove last added label",
    )

undo_label_button.on_click(undo_label_button_clicked)

tracer_block_111_000.children = [
    widgets.HBox([
        show_tracer_checkbox,
    ],
        layout=widgets.Layout(
            align_items="center",
            width="70%")
    ),
    v_slider,
    widgets.HBox([
        label_input,
        add_label_button,
        undo_label_button
    ],
        layout=widgets.Layout(
            align_items="center"
        )
    )
]

#############################
#FIGURE BLOCK
#############################
change_isotherms(None)

figure_block_111_000 = widgets.VBox(
    [],
    layout = widgets.Layout(width='70%')
)

prepare_export_fig_111_000_button = widgets.Button(
        description='Export',
        disabled=False,
        button_style='',
        tooltip='',
    )

prepare_export_fig_111_000_button.on_click(prepare_export)

go_to_export_button = widgets.HTML()

figure_block_111_000.children = [
    widgets.HBox([
        p_slider,
        fig_111_000,
        T_slider
    ]),
    widgets.VBox([
        prepare_export_fig_111_000_button,
        go_to_export_button
        ],
        layout=widgets.Layout(
            align_items="center",
            width="100%"
        )
    ),
    tracer_block_111_000
]

middle_block_111_000.children = [
    left_block_111_000,
    figure_block_111_000,
    right_block_111_000
]

############################
########MAIN BLOCK##########
############################

main_block_111_000 = widgets.VBox(
    [],
    layout=widgets.Layout(
        width='100%',
        align_items='center'
    )
)


main_block_111_000.children = [
    top_block_111_000,
    middle_block_111_000,
]

figures = [
    fig_111_000
]

main_block_111_000