The research presented in this thesis is motivated by many applications of high intensity focused ultrasound, ranging from kidney stone treatments, welding, and heat therapy to sonochemistry, where a better understanding of the physical effects through mathematical analysis and optimization is expected to lead to an improvement in precision and reduction of risks. A significant part of the thesis is dedicated to the question of well-posedness for the Westervelt equation with strong nonlinear damping under practically relevant Neumann as well as absorbing boundary conditions. We prove local in time well-posedness through a fixed point approach under the assumption that the initial and boundary data are sufficiently small. We also obtain short time well-posedness for the problem modeling the interface coupling of acoustic regions, since this is a common occurrence in applications, e.g. lithotripsy. Secondly, higher interior regularity for the model in question, as well as for the coupled system, is achieved by employing the difference quotient approach. In the latter case, we also show that the result can be extended up to the boundary of the subdomains provided that the gradient of the acoustic pressure is essentially bounded in space and time on the whole domain. This result is of importance in future numerical approximations of the present problem, as well as in shape optimization problems governed by this model. The last part of the thesis is dedicated to the shape sensitivity analysis for an optimization problem arising in lithotripsy. The goal is to find the optimal shape of a focusing acoustic lens so that the desired acoustic pressure at a kidney stone is achieved. We follow the variational approach to calculating the shape derivative of the cost functional which does not require computing the shape derivative of the state variable; however, assumptions of certain spatial regularity of the primal and the adjoint state are needed to obtain the derivative, in particular to express it in its strong form in terms of boundary integrals.