.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "full_field_verification_examples/plot_vfm_methods_verification.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_full_field_verification_examples_plot_vfm_methods_verification.py: Virtual Fields Method Verification ================================== In this example, we verify the virtual power residual calculation for a boundary value problem (BVP) with an analytical solution. To conform to our VFM implementation constraints, we chose a thin plate of length :math:`H` loaded in uniaxial tension in the Y direction with a distributed load with an over all magnitude of :math:`L`. Assuming linear elasticity and small deformation, the stress in the Y direction in the plate is .. math:: \sigma_{yy} =\frac{L}{W T} where :math:`W` is the plate width and :math:`T` is the plate thickness. All other stresses are zero. The in-plane displacements in the plate are given by .. math:: u=-\nu\frac{\sigma_{yy}}{E}x and .. math:: v=\frac{\sigma_{yy}}{E}y where :math:`E` is the material elastic modulus, :math:`\nu` is the material Poisson's ratio, and :math:`x` and :math:`y` are the locations on the plate. The internal virtual power becomes .. math:: P_{int}=\frac{L}{Wt}\int_{z=-T/2}^{T/2}\int_{y=-H/2}^{H/2}\int_{x=-W/2}^{W/2}\dot{\mathbf{F}}^*_{yy}dxdydz where :math:`\dot{\mathbf{F}}^*_{yy}` is the :math:`\frac{\partial\mathbf{v}^*}{\partial y}`. The virtual velocity function used is .. math:: \mathbf{v}^*_y=\frac{2y+H}{2H} With this virtual velocity function, the virtual velocity gradient component of interest is .. math:: \frac{\partial\mathbf{v}^*}{\partial y} = 1/h. As a result, the internal virtual power is given as .. math:: P_{int}=L The external virtual power, as expected, results an equivalent value. .. math:: P_{ext} = L\mathbf{v}^*_y = L\frac{2H/2+H}{2H} = L With an analytical solution for the virtual internal and external powers available for the BVP of interest, we will create MatCal VFM objects with the proper input and verify that they return the expect values. We will begin by creating synthetic data that represents the experimental data for this BVP. To do so, we import MatCal's tools and create the variables needed for the analytical solutions above. We will need: #. The materials elastics constants :math:`E` and :math:`\nu`. We choose values that are similar to steel with :math:`E=200` GPa and :math:`\nu=0.25`. #. The plate dimensions :math:`H=15.2` mm, :math:`W=7.6`, and :math:`T=0.1` mm with a discretization in Y and Y. #. The load magnitude as a function of time :math:`L(t)`. We will choose :math:`L(t)` such that the maximum stress in the plate is 100 MPa. #. The number of time steps and the time values. .. GENERATED FROM PYTHON SOURCE LINES 104-124 .. code-block:: Python from matcal import * import numpy as np E = 200e9 nu = 0.25 H = 15.2e-3 W = 7.6e-3 T = 0.1e-3 measured_num_points = 10 measured_xs = np.linspace(0, W, measured_num_points) measured_ys = np.linspace(0, H, measured_num_points) measured_x_grid, measured_y_grid = np.meshgrid(measured_xs, measured_ys) n_time_steps = 10 L = np.linspace(0, 100e6*W*T,n_time_steps).reshape(n_time_steps,1) time = np.linspace(0,1,n_time_steps) .. GENERATED FROM PYTHON SOURCE LINES 125-137 With all the required inputs initialized, we can create the :class:`~matcal.full_field.data.FieldData` objects that are need to have MatCal evaluate the virtual powers for the problem. The :class:`~matcal.full_field.data.FieldData` object requires several fields: #. time #. load #. The X locations of the field data #. The Y locations of the data #. The X displacements for each point (U) #. The Y displacements for each point (V) .. GENERATED FROM PYTHON SOURCE LINES 137-147 .. code-block:: Python sigma_yy = L/(W*T) def displacements(x,y): return -nu*sigma_yy/E*x, sigma_yy/E*y all_x_1D = measured_x_grid.reshape(measured_num_points**2) all_y_1D = measured_y_grid.reshape(measured_num_points**2) u,v = displacements(all_x_1D, all_y_1D) .. GENERATED FROM PYTHON SOURCE LINES 148-155 After creating the values for the fields that we need, they can be combined into a dictionary and converted into a MatCal :class:`~matcal.full_field.data.FieldData` object using the :func:`~matcal.full_field.data.convert_dictionary_to_field_data` function. .. GENERATED FROM PYTHON SOURCE LINES 155-166 .. code-block:: Python data_dict = {'time':time, 'load':L, "X":all_x_1D, "Y":all_y_1D, "U":u, "V":v} field_data = convert_dictionary_to_field_data(data_dict, coordinate_names=["X", "Y"]) .. GENERATED FROM PYTHON SOURCE LINES 167-178 Next, a discretization for the VFM models must be created. MatCal has a simple tool for creating rectangular volumes that we will use, the :func:`~matcal.full_field.TwoDimensionalFieldGrid.auto_generate_two_dimensional_field_grid` function. It requires the number of X and Y discretization points and the field data that the mesh is intended to encapsulate as inputs. It will return MatCal's two dimensional mesh class that can be used as a mesh input to MatCal's VFM model. It currently can only be used for rectangular shapes without holes. .. GENERATED FROM PYTHON SOURCE LINES 178-183 .. code-block:: Python from matcal.full_field.TwoDimensionalFieldGrid import auto_generate_two_dimensional_field_grid auto_mesh = auto_generate_two_dimensional_field_grid(measured_num_points*5, measured_num_points*5, field_data) .. GENERATED FROM PYTHON SOURCE LINES 184-190 The next item needed for the creation of the MatCal VFM models is a SierraSM material file. We will create a simple elastic material model file using Python file tools and then create a MatCal :class:`~matcal.sierra.material.Material` object so the VFM models can use that material. .. GENERATED FROM PYTHON SOURCE LINES 190-206 .. code-block:: Python material_file_string = \ """ begin material steel begin parameters for model elastic youngs modulus = {E} poissons ratio = {nu} end end """ with open("elastic_material.inc", "w") as mf: mf.write(material_file_string) material = Material("steel", "elastic_material.inc", "elastic") .. GENERATED FROM PYTHON SOURCE LINES 207-217 We can now create the VFM models that we want to evaluate. We will look at two VFM models, The :class:`~matcal.sierra.models.VFMUniaxialTensionHexModel` and the :class:`~matcal.sierra.models.VFMUniaxialTensionConnectedHexModel` will both be evaluated in this verification example. The models require a material, a mesh and thickness to be initialized and field data with a two-dimensional displacement field to be added as the boundary condition data. We optionally set the names so their results can be more easily pulled from the results data object returned from our study. .. GENERATED FROM PYTHON SOURCE LINES 217-229 .. code-block:: Python default_hex_model = VFMUniaxialTensionHexModel(material, auto_mesh, T) default_hex_model.set_name("default_VFM_hex_model") default_hex_model.add_boundary_condition_data(field_data) default_hex_model.set_number_of_time_steps(20) connected_hex_model = VFMUniaxialTensionConnectedHexModel(material, auto_mesh, T) connected_hex_model.set_name("VFM_connected_hex_model") connected_hex_model.add_boundary_condition_data(field_data) connected_hex_model.set_number_of_time_steps(20) .. GENERATED FROM PYTHON SOURCE LINES 230-242 The final step is to create the parameters for the study, the objective to be evaluated and the study itself. For our elastic model, the only parameters we need are the elastic modulus and the Poisson's ratio. The study we will use to evaluate the models is the :class:`~matcal.core.parameter_studies.ParameterStudy` since just need the model results at our predetermined steel-like parameter values. Since this study is evaluating a VFM model, we create a VFM model with no inputs because our global data fields of "time" and "load" match the default ones expected by the objective. .. GENERATED FROM PYTHON SOURCE LINES 242-248 .. code-block:: Python E_param = Parameter('E', 100e9, 300e9) nu_param = Parameter('nu', 0.1, 0.5) vfm_objective = MechanicalVFMObjective() vfm_objective.set_name('vfm_objective') .. GENERATED FROM PYTHON SOURCE LINES 249-254 After creating the parameters, objectives, synthetic data and models, we can create the study and add the evaluation sets and parameter values that we want to evaluate. The study is then launched and the results are stored in an object that we can analyze when the study completes. .. GENERATED FROM PYTHON SOURCE LINES 254-262 .. code-block:: Python study = ParameterStudy(E_param, nu_param) study.add_parameter_evaluation(E=E, nu=nu) study.add_evaluation_set(default_hex_model, vfm_objective, field_data) study.add_evaluation_set(connected_hex_model, vfm_objective, field_data) study.set_core_limit(3) results = study.launch() .. rst-class:: sphx-glr-script-out .. code-block:: none You are using exodus.py v 1.21.6 (seacas-py3), a python wrapper of some of the exodus library. Copyright (c) 2013-2023 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. Opening exodus file: matcal_template/default_VFM_hex_model/matcal_default_state/default_VFM_hex_model.g Closing exodus file: matcal_template/default_VFM_hex_model/matcal_default_state/default_VFM_hex_model.g Opening exodus file: matcal_template/default_VFM_hex_model/matcal_default_state/default_VFM_hex_model.g Closing exodus file: matcal_template/default_VFM_hex_model/matcal_default_state/default_VFM_hex_model.g Opening exodus file: matcal_template/default_VFM_hex_model/matcal_default_state/default_VFM_hex_model.g Closing exodus file: matcal_template/default_VFM_hex_model/matcal_default_state/default_VFM_hex_model.g Opening exodus file: matcal_template/default_VFM_hex_model/matcal_default_state/default_VFM_hex_model.g Opening exodus file: matcal_template/default_VFM_hex_model/matcal_default_state/default_VFM_hex_model_exploded.g Closing exodus file: matcal_template/default_VFM_hex_model/matcal_default_state/default_VFM_hex_model_exploded.g Closing exodus file: matcal_template/default_VFM_hex_model/matcal_default_state/default_VFM_hex_model.g Opening exodus file: matcal_template/VFM_connected_hex_model/matcal_default_state/VFM_connected_hex_model.g Closing exodus file: matcal_template/VFM_connected_hex_model/matcal_default_state/VFM_connected_hex_model.g Opening exodus file: matcal_template/VFM_connected_hex_model/matcal_default_state/VFM_connected_hex_model.g Closing exodus file: matcal_template/VFM_connected_hex_model/matcal_default_state/VFM_connected_hex_model.g Opening exodus file: matcal_template/VFM_connected_hex_model/matcal_default_state/VFM_connected_hex_model.g Closing exodus file: matcal_template/VFM_connected_hex_model/matcal_default_state/VFM_connected_hex_model.g .. GENERATED FROM PYTHON SOURCE LINES 263-280 In this verification problem, our goal is to verify that our MatCal VFM objective is accurately evaluating the internal and external virtual powers. As a result, we extract these virtual powers from the study results and compare them to the analytically determined value for the internal power :math:`L\left(t\right)`. We compare them by evaluating their percent error with the following equation .. math:: Pe = 100\frac{P_{eval} - L}{\max\left(L\right)} where :math:`P_{eval}` is the evaluated virtual power (either the internal or external from either model). Since we will evaluate this percent error several times, we make a function to perform the calculation. .. GENERATED FROM PYTHON SOURCE LINES 280-284 .. code-block:: Python def percent_error(prediction, actual): return (prediction - actual)/np.max(actual)*100 .. GENERATED FROM PYTHON SOURCE LINES 285-287 We then extract the results from each of the models and print or plot their errors. .. GENERATED FROM PYTHON SOURCE LINES 287-309 .. code-block:: Python matcal_internal_power_default = results.best_simulation_qois(default_hex_model, vfm_objective, field_data.state, 0)["virtual_power"] matcal_external_power_default = results.get_experiment_qois(default_hex_model, vfm_objective, field_data.state, 0)["virtual_power"] matcal_internal_power_connected = results.best_simulation_qois(connected_hex_model, vfm_objective, field_data.state, 0)["virtual_power"] matcal_external_power_connected = results.get_experiment_qois(connected_hex_model, vfm_objective, field_data.state, 0)["virtual_power"] print(percent_error(matcal_external_power_default, L)) print(percent_error(matcal_external_power_connected, L)) .. rst-class:: sphx-glr-script-out .. code-block:: none [[0.] [0.] [0.] [0.] [0.] [0.] [0.] [0.] [0.] [0.]] [[0.] [0.] [0.] [0.] [0.] [0.] [0.] [0.] [0.] [0.]] .. GENERATED FROM PYTHON SOURCE LINES 310-322 Both of the external powers have no error. This is expected as our experimental data for this study was exact. This result demonstrates in part that all of the MatCal's VFM objective calculations and virtual fields are correctly implemented. We have unit tests that test each of the components to contribute to this individually. The next evaluation compares the internal virtual power to the expected virtual power. This is significantly more involved as it requires the simulation of the VFM model using SierraSM. .. GENERATED FROM PYTHON SOURCE LINES 322-328 .. code-block:: Python internal_power_error_default = percent_error(matcal_internal_power_default.reshape(10,1), L) internal_power_error_connected = percent_error(matcal_internal_power_connected.reshape(10,1), L) print(internal_power_error_default) print(internal_power_error_connected) .. rst-class:: sphx-glr-script-out .. code-block:: none [[ 0. ] [-0.0006654 ] [-0.00260389] [-0.00581517] [-0.01029905] [-0.01608401] [-0.02317037] [-0.03152879] [-0.04115908] [-0.05206094]] [[ 0. ] [-0.0006654 ] [-0.00260389] [-0.00581517] [-0.01029905] [-0.01608401] [-0.02317037] [-0.03152879] [-0.04115908] [-0.05206094]] .. GENERATED FROM PYTHON SOURCE LINES 329-334 When we print the percent errors, it is clear that the some part of the process has introduced errors. The magnitude of the maximum error is 0.05%. To further investigate, we plot the internal power errors as a function time. .. GENERATED FROM PYTHON SOURCE LINES 334-345 .. code-block:: Python import matplotlib.pyplot as plt plt.plot(time, internal_power_error_default, label="default VFM model", marker='o') plt.plot(time, internal_power_error_connected, label="connected hex VFM model") plt.xlabel("time (s)") plt.ylabel("Internal Power Error (%)") plt.legend() plt.show() .. image-sg:: /full_field_verification_examples/images/sphx_glr_plot_vfm_methods_verification_001.png :alt: plot vfm methods verification :srcset: /full_field_verification_examples/images/sphx_glr_plot_vfm_methods_verification_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 346-356 The error is increasing quadratically as the load is increased and the model becomes more deformed. This error is due to our initial assumption of small deformation. SierraSM is formulated as a large deformation code. As a result, our small deformation assumption results in this small, but noticeable error. In later verification examples, we evaluate our methods with synthetic data using large deformation solutions and are able to obtain results with much smaller errors. .. rst-class:: sphx-glr-timing **Total running time of the script:** (2 minutes 1.271 seconds) .. _sphx_glr_download_full_field_verification_examples_plot_vfm_methods_verification.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_vfm_methods_verification.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_vfm_methods_verification.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_vfm_methods_verification.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_