Implementing the Weak Form in COMSOL Multiphysics

Chien Liu January 6, 2015
Share this on Facebook Share this on Twitter Share this on LinkedIn

This blog post is part of a series aimed at introducing the weak form with minimal prerequisites. In the first blog post, we learned about the basic concepts of the weak formulation. All equations were left in the analytical form. Today, we will implement and solve the equations numerically using the COMSOL Multiphysics simulation software. You are encouraged to follow the steps with a working copy of the COMSOL software.

Recapping the Basic Ideas

Recall that in the previous entry, we studied a simple example of 1D heat transfer at steady state with no heat source, where the temperature T is a function of the position x in the domain defined by the interval 1\le x\le 5.

The weak formulation turns the differential equation for the heat transfer physics into an integral equation, with a test function \tilde{T}(x) as a localized sampling function within the integrand to clamp down the solution. Integrating the weak form by parts provides the numerical benefit of reduced differentiation order. It also provides a natural way to specify boundary conditions in terms of the heat flux. For fixed boundary conditions, in terms of the temperature, the weak formulation uses the same mechanism of test functions and its natural boundary conditions to construct additional terms in the equation system.

In the end, we arrived at an exemplary equation that looks like this:

(1)

\int_1^5 \partial_x T(x) \partial_x \tilde{T}(x) \,dx = -2 \tilde{T}_1 -\lambda_2 \tilde{T}_2 -\tilde{\lambda}_2 (T_2-9)

Here, the integrand on the left-hand side involves only the first derivative of the temperature, the first term on the right-hand side defines that the outgoing flux should be 2 at the left boundary (x=1), and the other two terms on the right-hand side together specify that the temperature should be 9 at the right boundary (x=5).

The Weak Form PDE Interface

To implement Eq. (1) in COMSOL Multiphysics, we use the Model Wizard to create a new 1D model with a Weak Form PDE (w) interface (under Mathematics > PDE Interfaces) and a Stationary study. The dependent variable can be set to T to match the notation in our equation. For the geometry, we make an Interval between 1 and 5. The weak expression under the default “Weak Form PDE 1” node reads: -test(Tx)*Tx+1[m^-2]*test(T), where the first term corresponds to the integrand in our Eq. (1) and the second term corresponds to a heat source, which is not in our simple example and should be removed from the input field.

The weak expression now reads: -test(Tx)*Tx, where Tx is the COMSOL Multiphysics notation for \partial_x T(x), the first derivative of the temperature, and test(Tx) is the first derivative of the test function \partial_x \tilde{T}(x). The negative sign comes from the convention that the input field assumes that the expression is on the right-hand side of the equal sign (as seen in the “Equation” section of the settings window), while the integral in our equation is on the left-hand side.

The Weak Form PDE Interface in COMSOL Multiphysics.

The Weak Contribution Feature

To implement the weak form terms on the right-hand side of Eq. (1) for the boundary conditions, right-click the Weak Form PDE (w) node. We see that there are built-in boundary features such as the Dirichlet Boundary Condition item, which is available in the pop-up menu for your convenience. However, since here we are interested in entering the equation ourselves, we hover the mouse over the item More in the pop-up menu and click on the item Weak Contribution in the next pop-up menu.

The Weak Contribution feature.

In the settings window for the “Weak Contribution 1” node under Boundary Selection, we select boundary 1 at the left end of the domain (at x=1). We then enter the weak expression as: -2*test(T) under the section Weak Contribution in the same settings window. This takes care of the first term on the right-hand side of Eq. (1), which specifies the outgoing flux to be 2 at the boundary x=1.

Fixed Boundary Condition

For the fixed boundary condition at x=5, where the last two terms on the right-hand side of Eq. (1) together specify that T=9, we create another “Weak Contribution” node at boundary 2 at the right end of the domain and an Auxiliary Dependent Variable subnode under it.

We enter lambda2 for the Field variable name in the subnode and then enter the weak expression as the two terms in Eq. (1): -lambda2*test(T)-test(lambda2)*(T-9)

An image highlighting the fixed boundary condition.

Discretization

The COMSOL software discretizes the domain by creating a mesh. Let’s right-click the Mesh 1 node and select Edge and then right-click Edge 1 and select Distribution. Then, we set the “Number of elements” to 4 and click Build All. We intentionally keep the number of elements small to make it easier when we discuss the discretization in more detail later.

A screenshot depicting how to create a mesh.

Also, under the Discretization section in the settings window for the Weak Form PDE (w) interface node, we set the “Element order” to Linear (click on the Show button under Model Builder and then the item Discretization in the pop-up menu to enable the Discretization section):

Choosing linear elements in the settings window.

Compute the Solution in COMSOL Multiphysics

Now we are ready to click Compute and check whether the solution makes sense.

The solution gives a straight line within the domain, which is consistent with the temperature profile at steady state with no heat source. The slope of the line is 2, which is consistent with the boundary condition that the outgoing flux is 2 at x=1. The temperature is 9 at x=5, as specified by the fixed boundary condition. Since there is no heat source, the total heat flux going out of the domain should sum up to zero in the steady state. Thus, the outgoing flux should be -2 at x=5.

We readily verify this by making a point evaluation of the heat flux variable lambda2, as shown in the screenshot below:

A plot shows the results of computing the solution.

Some readers may wonder whether it is always necessary to solve for the auxiliary variable lambda2, the so-called Lagrange multiplier, especially if it is not needed by the modeler and solving for it inevitably requires more computation. As we will see in the following posts, COMSOL Multiphysics provides alternative features and allows the user to decide whether or not to solve for the Lagrange multiplier.

Summary and Next Up

Today, we refreshed the concept of the weak formulation and implemented an exemplary weak form equation (1) in COMSOL Multiphysics. The resulting numerical solution behaves as expected from simple physical arguments.

In future blog posts, we will take a look “under the hood” to see how the weak form equations, such as Eq. (1), are discretized and solved numerically. We will see how the same problem can be solved in different ways and how different boundary conditions can be set up for different types of problems.

Stay tuned!



Comments

  1. Stefano Maffei January 8, 2015   6:03 am

    Thanks a lot for this clear article. I am looking forward for the next ones.
    Do you have a specific reference concerning finite elements/spectral elements that I can use to get more insights into these concepts?

  2. Chien Liu January 8, 2015   10:19 am

    Dear Stefano,

    Thank you for your interest in this blog post. You may find the list of books below of interest.

    Sincerely,

    Chien

    These books might be of interest for getting an in-depth knowledge of finite element analysis:

    * T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications (2000)

    * O. C. Zienkiewicz & R. L. Taylor, The Finite Element Method Set, Butterworth-Heinemann; 7th edition (2011)

    * S.C. Brenner & L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer; 3rd edition (2009)

    For the application of the finite element method to partial differential equation in particular, these books might be of interest:

    * K. Eriksson, D. Estep, P. Hansbo, C. Johnson, Computational Differential Equations, Cambridge University Press; 2nd edition (1996)

    * C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Dover Publications (2009)

    For fluid flow, this book has a section on comparisons between the finite volume method and the finite element method. It is also a good reference for different finite element types used for CFD:

    * P. M. Gresho & R. L. Sani, Incompressible Flow and the Finite Element Method, Volume 2, Isothermal Laminar Flow, John Wiley & Sons (2000)

    For Electromagnetics, especially high-frequency, this is an excellent reference:

    * Jianming Jin, The Finite Element Method in Electromagnetics, Wiley-IEEE Press; 3rd edition (2014)

  3. Camille Spingarn February 3, 2015   5:12 am
  4. Ren Wenxi February 7, 2015   6:48 am

    Dear Liu,
    Thank you so much for your article about the weak form. I have been interested in the application of the weak form in structure analysis. Recently, i am trying to design a simple structure analysis module by weak form. But there is a problem in my code. I hope your help. The following is my code:
    % Displacement field
    % Dependent variables u, v
    % Variables
    E 1.44*10^4[MPa] % Young modulus
    pr 0.2 % Poisson’s ratio
    % Elasticity matrix
    D11 E/(1-pr^2)
    D22 E/(1-pr^2)
    D33 E/(2*(1+pr))
    D12 E*pr/(1-pr^2)
    D21 D12
    D23 0[Pa]
    D32 0[Pa]
    D13 0[Pa]
    D31 0[Pa]
    % Strain
    ex ux
    ey vy
    exy 0.5*(uy+vx)
    % Stress
    sx D11*ex+D12*ey+2*D13*exy
    sy D21*ex+D22*ey+2*D23*exy
    sxy D31*ex+D32*ey+2*D33*exy
    % Equations
    -sx*test(ux)-sxy*test(uy) % X Direction
    -sy*test(vy)-sxy*test(vx) % Y Direction

  5. Ren Wenxi February 7, 2015   6:49 am

    Look forward to your reply~

  6. Chien Liu February 9, 2015   10:37 am

    Hi Ren,

    Thank you for contacting us on this. In this case we recommend that you please contact our support team for help. They can be reached at:

    http://www.comsol.com/support

    Best regards,

    Chien

  7. Stefano Maffei February 18, 2015   5:38 am

    Hi Chien,
    Let’s say I have a system of equations (in my case i have an eigenvalue problem). Then I multiply both equations by their corresponding test functions (say my variables are y and c), integrate by parts and apply the boundary conditions I have. Imagine after all I am left with something like:

    \int [ weak form for y ] dx = [lambda*test(y)*yx]_x1
    \int [ weak form for c ] dx = [test(c)*cx]_x1

    where lambda is the eigenvalue (sorry if this looks complicated, I just think that with an example is better to talk about things). Imagine that somehow i have no informations about the derivatives in x0 (say i x1 everything vanishes).

    My question is: these boundary terms have t be retained, as they might be important in adjusting the behaviour of the solution near x1. How should I treat them? Should I put in the weak contribution subnode
    lambda*test(y)*yx + test(c)*cx
    applied on x1?

    Thanks

  8. Chien Liu February 18, 2015   8:17 am

    Hi Stefano,

    Thank you for contacting us on this. In this case we recommend that you please contact our support team for help. They can be reached at:

    http://www.comsol.com/support

    Thanks!

    Chien

  9. Pu Zhang April 16, 2015   11:54 pm

    Very nice blog post! So is this possibility of inspecting system matrix a new feature in COMSOL 5?

  10. Chien Liu April 17, 2015   11:01 am

    Hello Pu,

    Thank you for the comment. You were probably referring to this entry:

    http://www.comsol.com/blogs/implementing-the-weak-form-with-a-comsol-app/

    The possibility of inspecting the system matrix is not new. In fact it has been available to users from the beginning.

    Of course the Application Builder is a brand new feature in COMSOL 5.0 (and further enhanced in 5.1). We believe the users’ productivity can be greatly improved by COMSOL Apps, as illustrated in the above blog post (and many others).

    Chien

  11. Chien Liu April 17, 2015   1:11 pm

    In the previous message I forgot to add the link to other blog posts on COMSOL Apps. They can be found here:

    http://www.comsol.com/blogs/category/all/applications/

  12. jackshine ma June 18, 2015   9:12 pm

    Dear Chien Liu,
    Very glad to learn from your blog on Weak form of comsol and it’s a very good blog specifying on weak form that merely introduced detailedly anywhere. However, I wonder that why there is a 2-order differential in integrand but only 1 time integral on the left of equation 1. If like this, the calculating result of left term would be in a differential form, which is not like the right polymial.

  13. Chien Liu June 22, 2015   11:47 am

    Hi jackshine ,
    Thank you for the comments.
    I’m not sure I understand your question since there is no time integral in equation 1.
    The integrand should be interpreted as having parenthesis as follows: (dx T) (dx T_test).
    Hope this help clarify the equation.
    Chien

  14. Vishwas Nesarikar October 22, 2015   9:06 am

    Can you please provide derivation of Eqn (1) ? Thx

  15. Chien Liu October 22, 2015   9:12 am

    Dear Vishwas,
    The derivation is in this previous blog post:
    http://www.comsol.com/blogs/brief-introduction-weak-form/
    Hope this helps.
    Sincerely,
    Chien

  16. Parth Swaroop April 11, 2016   5:59 am

    Dear Chien,

    Implementation of this weak formulation and use T_lm is also presented in example of “Tin Melting Front” module provided by comsol. In this example the movement of the melting phase is implemented as :
    T_lm[W/m^2]/(rho*DelH)

    i guess [W/m^2] is used make the units m/s ( normal velocity of the prescribed mesh ).
    how do i know that this expression actually evaluates the difference in flux at the interface and not the temperature profile.

    from physics, the expression for normal velocity is
    ΔH_fus ρ v_n= ϕ_l- ϕ_s

    Parth

  17. Chien Liu April 12, 2016   11:16 am

    Dear Parth,
    In the Tin Melting Front example, the variable T_lm is the heat flux and it plays exactly the same role as the variable lambda2 in the blog post.
    The formula for the normal mesh velocity v can be understood in the following way:
    v * rho = rate of melting
    DelH = latent heat of fusion
    therefore
    v * rho * DelH = heat flux = T_lm
    Does this make sense?
    Sincerely,
    Chien

  18. Jorge Yannie April 17, 2016   8:36 am

    Hello Chien,

    Great post once again! Looking forward to all your Equation-based post! It is my main goal to learn in Comsol!

    Like it lot!

  19. Puneeth Banisetti April 23, 2016   8:35 am

    Hello Chein,
    Thanks a lot for a very informative blog post. I thoroughly enjoyed it. I want to know how we can implement this technique when we have some normal boundary conditions for 3D problems such a parallel plate rheometer. How do we incorporate the boundary condition (T.n).n=0 and what exactly do we take n when we are using cartesian coordinates since n here will be in the radial direction.

    Could you also please post the link for the continuation of this blog post.

    Regards,
    Puneeth

  20. Chien Liu April 26, 2016   10:58 am

    Dear Puneeth,

    You can find related posts by typing in “weak form” in the “Search Blog” search box above.

    For normal variables please refer to the COMSOL Reference Manual > Global and local Definitions > Predefined and Built-in Variables > Geometric Variables and Mesh Variables > Tangent and Normal Variables.

    Sincerely,

    Chien

  21. Edward Smith June 23, 2016   10:36 am

    Hi
    I’ve been trying to implement this example but cannot select the points at either end of the interval. It only lets me select the whole domain. I suspect I have done something wrong in Geometry.
    Thanks

  22. Chien Liu June 23, 2016   10:51 am

    Dear Edward, when creating a weak contribution on a boundary, make sure to select the item from the group of boundary features, not the domain features. It sounds like you selected the latter. You may find the 2nd screenshot in the blog post helpful.
    Sincerely,
    Chien

  23. Edward Smith June 23, 2016   11:32 am

    Hi
    Thanks for the reply.
    What a fool I am 🙂
    Ted

  24. Edward Smith June 23, 2016   11:39 am

    I am trying to use the Laminar flow subset of the N-S equations, and want to prevent any pressures going negative. Can I do this with a weak contribution?

  25. Iris Zhou November 14, 2016   4:36 pm

    Hi Chien,

    Thanks for your blog, and it is really nourishing!

    Right now I am working with Frequency domain Wave-optics Module (ewfd2). I am using electric field as the input boundary condition, and with numerically calculated data from an eigenfrequency solver (ewfd.E). I want to put both electric fields and their derivatives into the boundary condition, and I thought this blog is quite relevant.

    Should I use electric field with Weak Constraints? If so, how should I put numerical derivative, e.g., d(ewfd.Ex,x), to the Weak Constraints?

    Looking forward for your reply! Thanks in advance!

    Iris

  26. Bridget Cunningham November 15, 2016   10:14 am

    Hello Iris,

    Thank you for your comment.
    For questions related to your modeling, please contact our support team.

    Online support center: https://www.comsol.com/support
    Email: support@comsol.com

    Best,
    Bridget

  27. akshay phadnis January 26, 2017   1:13 pm

    Dear Chien,

    Thanks for this article. It’s a great help. There is no other documentation about weak form implementation in COMSOL. I am trying to model the polymer swelling using solid mechanics equation – div(S)+b=0. When I convert this into a weak form, I get the expression in terms of deformation matrix. Now comsol weak form interface cannot take deformation tensor as input, it has to be specified in terms of displacement. Do you have any suggestions on how I can proceed with implementing solid mechanics problem using weak form ?

  28. Bridget Cunningham January 26, 2017   1:37 pm

    Hello Akshay,

    Thank you for your comment.
    For questions related to your modeling, please contact our support team.

    Online support center: https://www.comsol.com/support
    Email: support@comsol.com

    Best,
    Bridget

  29. Manal Badgaish February 19, 2017   1:20 pm

    Dear All,
    I have the following 1D PDE system:
    u_tt= c^2 *u_xx,
    u(x,0)=u_t(x,0)=0, (initials conditions)
    m* u_tt(0,t)= a*PB(t)+rho*c^2*a*u_x(0,t)-D1*u(0,t)- D2*u^2(0,t)-D3, (boundary condition1)
    u_t(Let)=-c*u_x(L,t). (boundary condition2)
    I want to find the solution at x=0 only, i.e I need to solve for u(0,t).

    Is it possible to implement this pdf system in comsol with this nonlinear complicated boundary condition at (x=0)? Can any one tell me the steps of doing it, and how I write the weak form of this complicated boundary conditions.

    Thanks,
    Manal

  30. Bridget Cunningham February 21, 2017   9:57 am

    Hello Manal,

    Thank you for your comment.
    For questions related to your modeling, please contact our support team.

    Online support center: https://www.comsol.com/support
    Email: support@comsol.com

    Best,
    Bridget

  31. Konark B May 16, 2017   12:15 pm

    Hi,
    What is the significance of negative sign in implementing weak form? Should it work equally if we don’t write negative sign since LHS is zero?

  32. Chien Liu May 16, 2017   12:24 pm

    Dear Konark,
    As long as the sign convention in your design is consistent (remember that there can be more than one weak contribution in general), you are free to choose the most convenient one for you.
    Sincerely,
    Chien

  33. eric malitz August 4, 2017   11:14 pm

    If we have multiple weak equations, i.e. a mixed finite element method, how would we code a boundary integral in each equation individually? The “weak contribution” options don’t seem to allow this specification.

  34. Chien Liu August 7, 2017   1:37 pm

    Dear Eric,
    The Weak Contribution feature’s input field allows any combination of multiple field variables, or if you prefer you can use multiple Weak Contribution nodes in the Model Builder tree structure. For any specific question please feel free to contact us at support@comsol.com or http://www.comsol.com/support
    Sincerely
    Chien

  35. Ngo Dinh Nhan August 17, 2017   4:08 am

    Dear Dr.Chien Liu

    I followed your steps in the instructor, but actually, the 1D plot was not same as your plot.
    I have attached my plot in this link ( google photo )
    https://goo.gl/photos/ueqVhn7CY9USQumR6
    Best regards

  36. Chien Liu August 22, 2017   11:15 am

    Dear Ngo,
    It looks like you may have forgotten to set the B.C. at x=1.
    Sincerely,
    Chien

  37. Vaidehi Oruganti February 1, 2018   2:55 pm

    Dear Dr.Chien Liu,

    I am trying to implement weak form modelling to solve a benchmark problem of flow between two flat plates (parabolic profile).

    I have trouble computing the velocity and pressure fields. I am ending up with the error that system matrix is zero and an under-specified problem.

    Method to solve flow between plates:

    (1) Two weak form PDE interfaces. One PDE solves for velocity field (momentum equation) and the other solves pressure field (continuity equation)

    Weak form PDE (1)

    x-comp: mu*(test(ux)*ux+test(uy)*uy)+px*test(u)
    y-comp: mu*(test(vx)*vx+test(vy)*vy)+py*test(v)

    BC(1) : No Slip -u*test(-u)-v*test(-v)
    Constraints are: [-u, -v]
    Constraint Force: [-test(u),-test(v)]

    For inlet and outlet: n.T_stress= -(p_ext)n
    p_ext at inlet is set to 2 [Pa] gauge and p_ext at outlet is zero

    Weak form of continuity:
    test(p)*(ux+vy)

    2 dirichlet conditions , p=2 [Pa] at inlet and p=0[Pa] at outlet are imposed.

    The simulation does not work. Can you throw some light on what I am doing wrong?

  38. Caty Fairclough February 1, 2018   3:06 pm

    Hi Vaidehi,

    Thanks for your comment.

    For your question, please contact our Support team.

    Online Support Center: https://www.comsol.com/support
    Email: support@comsol.com

  39. Vaidehi Oruganti February 1, 2018   3:08 pm

    Dr. Chien Liu,

    I tried contacting the COMSOL support about the issues I have, but I heard that they cannot provide assistance on weak form modelling to people having academic licenses (which I have) Therefore, I request you to throw some light on the problem. I did read the weak form documentation and all of your blog posts on weak form modelling.

    Thank you Sir.

  40. Vaidehi Oruganti February 1, 2018   3:10 pm

    I tried contacting COMSOL support but they say they cannot offer support on weak form modeling. This is the reason I am trying to reach out to you.

    Thank you Sir.

  41. Vaidehi Oruganti February 1, 2018   4:50 pm

    Dear Caty Fairclough,

    The following is the reply I received from COMSOL support just now, about the above question.

    “Dear Vaidehi,

    Thank you for contacting COMSOL Support.

    The question being asked here pertains to equation based modeling.
    Additionally, the application of this equation based modeling is recreating
    an existing fluid flow interface. These types of questions do fall outside
    the scope of the technical support we provide.

    Best Regards,
    COMSOL Technical Support”

    So, you see COMSOL support does not answer any questions about equation based modelling. Therefore, I am trying to get some assistance here.

  42. shi tan April 15, 2019   4:32 pm

    Dear Dr.Liu
    I am using version 5.4, what are detailed steps to get the value of lambda2?

  43. Chien Liu April 16, 2019   9:45 am

    Dear Shi, The steps should be the same. Please contact support if you need further assistance. Sincerely, Chien

Loading Comments...

Categorias


Tags