Multiphysics 0471

The code has been written by Louis Denis, Pierre-Yves Goffin, Kévin Maltez and Valentin Vanraes in the context of the Multiphysics Integrated Computational Project (Spring 2022), as part of the Master of Science in Engineering Physics of the Faculty of Applied Sciences at University of Liège, Belgium.

Building the solver

FOR WINDOWS:

To build the coupled FEM/BEM solver:

• go to the srcs folder (cd srcs)
• create a build folder (mkdir build)
• go to this build folder (cd build)
• cmake ..
• make

To execute the solver:

• in build folder: solver.exe ..\hybrid_geo.geo (or any other .geo file present at the root of the srcs folder)

To build the independent FEM solver:

• similar procedure but to be done in the srcs/FEM folder

To build the independent BEM solver:

• similar procedure but to be done in the srcs/BEM folder

Creation of a .geo file

To create a .geo file compatible with the program, some rules must be respected:

• For the elastic FEM part:

  1. The FEM domain should be defined as a Physical Surface with the precise name FEM_domain. First, a Curve Loop gathering the lines representing the boundary of the FEM domain must be created. Then, a Plane Surface containing the corresponding curve loop must be created. Finally, if the plane surface is the first of the .geo file, the physical surface must be created as: Physical Surface("FEM_domain", x) = {1};, where x is the tag of the physical surface. Multiple sub domains can be included in the FEM domain. For example, if plane surfaces have been created for five sub domains, the physical surface must be created as:

     Physical Surface("FEM_domain", x) = {1, 2, 3, 4, 5};. 
     

  2. The mechanical boundary conditions, material properties and volumic forces must be specified for the FEM domain:

  • BOUNDARY CONDITIONS:
    1. The edges on which the user wants to impose a boundary condition must be defined as Physical Curve. As an example, if the user want to impose a condition on the bottom edge related to the Line(1):

       Physical Curve("bottom_edge", x) = {1};

    2. All boundary conditions must be written as SetNumber("Boundary Conditions/name_of_the_physical_curve/...") with name_of_the_physical_curve the physical curve on which the condition is imposed.
    3. The Dirichlet boundary conditions, i.e the imposed horizontal and vertical displacement u_x and u_y of an edge,
       expressed in [m], must be written as Setnumber("Boundary Conditions/name_of_the_physical_curve/ux", desired_value);, same applies for u_y. Note that the horizontal displacement can be imposed on a physical curve without imposing the vertical displacement, and vice-versa. As an example, if the user wants the left edge to be clamped:

       SetNumber("Boundary Conditions/left_edge/ux", 0.);
       SetNumber("Boundary Conditions/left_edge/uy", 0.);

    4. The Neumann boundary conditions, i.e surface traction in the horizontal (t_x) or vertical (t_y) direction imposed on an edge, expressed in [Pa], must be written as SetNumber("Boundary Conditions/name_of_the_physical_curve/tx", desired_value), same applies for t_y. Note that on a specific edge, both horizontal and vertical surface tractions must be prescribed. If the user only wants to prescribe t_x, t_y must also be set to zero.
    5. Dirichlet boundary conditions can also be imposed on a specific point of the domain. In this case, a Physical Point must be created.
  • MATERIAL PROPERTIES: All the material properties must be written as SetNumber("Materials/FEM_domain/name_of_property", value_of_property);. The different properties that must be specified are:
    • Young's modulus, expressed in [Pa], named Young.
    • Poisson's ratio, without physical units, named Poisson.
    • Mass density, expressed in [kg/m3], named rho.
  • VOLUMIC FORCES: Volumic forces in the horizontal (b_x) or vertical (b_y) direction, expressed in [m/s2], must be written as SetNumber("Volumic Forces/FEM_domain/b_x",0);, same applies for b_y.

• For the electrostatic BEM part:

  1. The lines corresponding to the boundary of the BEM surfaces must be created in such a way that the area of the BEM surface is located on the left of the Curve Loop.
  2. Multiple BEM domains can be defined. One BEM domain should be defined as a Physical Surface with the precise name BEM_domain_X, with X the identifier of the BEM domain (X=1 if it is the first domain, X=2 if it is the second one, ...).
  3. Boundary conditions and material properties must then be specified for the BEM domain:
  • BOUNDARY CONDITIONS:
    1. Dirichlet boundary conditions, corresponding to imposing an electric potential expressed in [V], can be imposed on a physical curve as: SetNumber("Boundary Conditions/name_of_the_physical_curve/BEM_domain_X/dirichlet", value_of_potential);.
    2. Neumann boundary conditions, corresponding to imposing the opposite of the exterior normal component of the electric field expressed in [V/m], can be imposed on a physical curve as: SetNumber("Boundary Conditions/name_of_the_physical_curve/BEM_domain_X/neumann", value_of_field);.
  • MATERIAL PROPERTIES: The dielectric permittivity of the material inside the BEM domain must be specified. It is done writing SetNumber("Materials/BEM_domain_X/Epsilon", value_of_Epsilon);.

  Note: These operations must be done for each BEM domain the user wants to create, replacing X by the given identifier of the BEM domain.

• For the coupling between the FEM and BEM domains:

  1. All the lines belonging to the interface between the FEM domain and the different BEM domains must be specified as a single Physical Curve("BEM_FEM_boundary, x)={List_of_the_interface_lines}. As an example, if Line(1), Line(4) and Line(5) belong to the intersection of the BEM domains and the FEM domain:

      Physical Curve("BEM_FEM_boundary", x) = {1, 4, 5};

• For the mechanical or the coupled solver, the user can choose between the linear or non-linear iterative solver, by setting SetNumber("Non_linear_solver",0); for the linear solver, or SetNumber("Non_linear_solver",1); for the non-linear solver.